Research Paper Volume 15, Issue 13 pp 6361—6379

Integrating the characteristic genes of macrophage pseudotime analysis in single-cell RNA-seq to construct a prediction model of atherosclerosis

Zemin Tian1, , Shize Yang2, ,

Received: April 12, 2023       Accepted: June 19, 2023       Published: July 8, 2023      

https://doi.org/10.18632/aging.204856

Copyright: © 2023 Tian and Yang. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Abstract

Background: Macrophages play an important role in the occurrence and development of atherosclerosis. However, few existing studies have deliberately analyzed the changes in characteristic genes in the process of macrophage phenotype transformation.

Method: Carotid atherosclerotic plaque single-cell RNA (scRNA) sequencing data were analyzed to define the cells involved and determine their transcriptomic characteristics. KEGG enrichment analysis, CIBERSORT, ESTIMATE, support vector machine (SVM), random forest (RF), and weighted correlation network analysis (WGCNA) were applied to bulk sequencing data. All data were downloaded from Gene Expression Omnibus (GEO).

Result: Nine cell clusters were identified. M1 macrophages, M2 macrophages, and M2/M1 macrophages were identified as three clusters within the macrophages. According to pseudotime analysis, M2/M1 macrophages and M2 macrophages can be transformed into M1 macrophages. The ROC curve values of the six genes in the test group were statistically significant (AUC (IL1RN): 0.899, 95% CI: 0.764-0.990; AUC (NRP1): 0.817, 95% CI: 0.620-0.971; AUC (TAGLN): 0.846, 95% CI: 0.678-0.971; AUC (SPARCL1): 0.825, 95% CI: 0.620-0.988; AUC (EMP2): 0.808, 95% CI: 0.630-0.947; AUC (ACTA2): 0.784, 95% CI: 0.591-0.938). The atherosclerosis prediction model showed significant statistical significance in both the train group (AUC: 0.909, 95% CI: 0.842-0.967) and the test group (AUC: 0.812, 95% CI: 0.630-0.966).

Conclusions: IL1RNHigh M1, NRP1High M2, ACTA2High M2/M1, EMP2High M1/M1, SPACL1High M2/M1 and TAGLNHigh M2/M1 macrophages play key roles in the occurrence and development of arterial atherosclerosis. These marker genes of macrophage phenotypic transformation can also be used to establish a model to predict the occurrence of atherosclerosis.

Introduction

Atherosclerosis (AS) is a local chronic inflammatory disease of the artery wall [1]. In AS, the inflammatory response is maintained for decades by the inflow, proliferation, and activation of immune cells [2, 3]. Among the many types of immune cells, macrophages have attracted great interest in AS because of their complex functions and numerous subtypes [4, 5]. Under the stimulation of the atherosclerotic microenvironment, macrophages can polarize into different types of macrophages, which involves changes in gene expression profile and cell functions [6]. In previous studies, polarized macrophages were mainly divided into M1 macrophages and M2 macrophages according to their phenotype and functions [7, 8]. M1 macrophages, as proinflammatory macrophages, can release a variety of chemokines and proinflammatory cytokines, including CCL2 (MCP-1), CCL3 and IL-1β, IL-6, IL-12, IL-23 and TNF, which can induce inflammatory reactions and can also secrete reactive oxygen intermediates (ROI), nitric oxide (NO) and lysosomal enzymes to kill and remove pathogens. However, ROS can induce tissue damage, leading to irreparable tissue damage, which may promote the formation of atherosclerosis and reduce the stability of atherosclerotic plaques [9]. In contrast, activated macrophages (M2) play a preventive role in the progression of human and mouse atherogenesis. M2 macrophages, as anti-inflammatory macrophages, mainly secrete anti-inflammatory cytokines, typically including IL-10 and TGF-β [10, 11]. These factors counteract the proinflammatory effects caused by M1 macrophages, thus promoting tissue repair and reducing tissue damage. In atherosclerotic plaques, M2 macrophages can inhibit the formation of atherosclerosis and maintain plaque stability [12, 13].

With the improved understanding of atherosclerosis by researchers, traditional lipid-lowering and anti-inflammatory treatments are being replaced by targeted treatments [14]. Research on the prevention and treatment of atherosclerosis has also become refined in this era [14]. These technologies include nanotechnology drug treatment, regulation of the ratio of M1/M2 macrophages, control of macrophage phenotype conversion, and regulation of the plaque microenvironment [8].

In this study, we analyzed the phenotypic transformation of macrophages in atherosclerosis through pseudotime analysis. We identified a class of macrophages that simultaneously express M1 and M2 macrophage marker genes and can transform into M1 or M2 macrophages. Finally, we selected the genes with the most obvious differential expression to build an atherosclerosis prediction model based on the results of macrophage pseudotime analysis. These genes may provide new ideas for the targeted treatment of atherosclerosis.

Materials and Methods

Bulk sequencing data processing

We downloaded GSE43292 (training group) and GSE28829 (test group) from the GEO database (https://www.ncbi.nlm.nih.gov/). We used the “WGCNA” package [15] to screen disease-characteristic genes from GSE43292. We utilized the “limma” package and the “VennDiagram” package to identify 223 intersecting genes (Supplementary Table 5) obtained from macrophage cell trajectory analysis. We then took the intersection of these genes with the 1086 AS-characteristic genes (Supplementary Table 4) obtained from GSE43292, resulting in 50 macrophage-related genes (Supplementary Table 6). Next, we compared the residual value and AUC curve of random forest (RF) and support vector machine (SVM) models. Finally, we used the RF method to select disease-characteristic genes again. Independent heatmaps of intersecting genes were drawn in the training and test groups. We selected the important genes with an RF score>1 in the training group and drew the AUC (area under the curve) curve in the test group for verification. Next, we used these genes to construct a disease prediction model. The “rms” and “rmda” packages [16] were used to draw the nomogram, calibration curve, decision curve analysis (DCA) [17], and clinical impact curve of DEGs predicting atherosclerosis in patients.

scRNA sequencing data processing

We downloaded the GSE159677 single-cell dataset from the GEO database (https://www.ncbi.nlm.nih.gov/). We removed cells with fewer than 200 genes, more than 7,000 genes, or more than 10% mitochondrial genes. Analysis was performed on 49576 filtered cells. Using the “LogNormalize” method, gene expression was normalized and scaled. In each sample, the top 2000 highly variable genes (HVGs) were identified using the “vst” method after data normalization. After identifying significant principal components (PCs), PCA was applied. Batch correction was performed using the “Harmony” R package (version 0.1.0) to avoid batch effects resulting from sample identity that could disrupt downstream analysis. Finally, 50 PCs were selected for t-distributed stochastic neighbor embedding (t-SNE) analysis. We set “FindClusters” with a resolution of 2.0 to divide cells into 45 different clusters, which were divided into 9 cell types, and conducted manual inspection according to the results of “FindAllMarkers”.

We used the same method to set “FindClusters” to divide macrophages into 13 different clusters with a resolution of 0.05, divided these clusters into three cell types with marker genes, and used the results of “FindAllMarkers” for manual inspection. The results of uniform manifold approximation and projection for dimension reduction (UMAP) were visualized, and a heatmap of macrophage subtypes was drawn.

We analyzed the trajectory of M1 macrophages and M2 macrophages with the “Monocle” package [18] and drew a heatmap. We used the same method to analyze M2 macrophages and M2/M1 macrophages, M1 macrophages and M2/M1 macrophages.

Consent for publication

The authors agree for publication.

Availability of data and materials

The data for this study were derived from the GEO database (https://www.ncbi.nlm.nih.gov/geo/).

Results

Analysis the heterogeneity of macrophages in the atherosclerotic single-cell transcriptome

Forty-five clusters could be assigned to known cell lineages according to the marker genes defined in a previous study (Figure 1A and Supplementary Table 1). We used t-SNE analysis to visualize the 9 clusters (Figure 1B) and identify marker genes in the 9 cell-type populations (Supplementary Table 2). The expression of cell type marker genes is shown in the dot plot (Figure 1C). We observed 9 cell clusters (CD4T/NK cells: clusters 2, 3, 6, 14, 16, 28, 30, 32 and 37, expressing CD3D [19], CD7 [20], and IL7R (interleukin 7 receptor) [21]; CTLs: clusters 0, 4 and 20, expressing CD8A [22]; Endothelial cells: clusters 5, 15, 21, 22, 33, 39, 43 and 44, expressing vWF and CLDN5 [23, 24]; Fibroblasts/vSMCs: clusters 7, 8, 13, 17, 18, 19, 23, 29, 31 and 35, expressing LUM (Lumican) [25] and ACTA2 (alpha actin 2, smooth muscle); vSMCs: clusters 2, 5, 14, 15, 21, 28, 30, 34, 41 and 42, expressing ACTA2 (alpha actin 2, smooth muscle) [26, 27]; B/Plasma_B cells: clusters 9, 24 and 40, expressing CD79A [28]; Macrophages: clusters 1, 10, 11, 12, 25, 36 and 38, expressing FCGR3A [29]; Mast cells: cluster 27, expressing CPA3 [30]; Monocytes: clusters 34, expressing FCGR3B (https://panglaodb.se/); and DCs: clusters 26, 41 and 42, expressing CD1C [31], CLEC9A [32], LILRA4 (https://panglaodb.se/).

Single-cell transcriptome data: (A, B) The single cells were divided into 45 groups on the basis of their transcriptome data and then ultimately divided into 9 cell populations. (C) Circle chart: The X-axis represents the marker genes that define the cells, and the Y-axis represents the different cell populations.

Figure 1. Single-cell transcriptome data: (A, B) The single cells were divided into 45 groups on the basis of their transcriptome data and then ultimately divided into 9 cell populations. (C) Circle chart: The X-axis represents the marker genes that define the cells, and the Y-axis represents the different cell populations.

Thirteen clusters could be assigned to known cell lineages on the basis of marker gene expression (Figure 2A). We used t-SNE analysis to visualize the 9 clusters (Figure 2B) and identify marker genes in the 3 cell-type populations (Supplementary Table 3). We observed 3 cell clusters (M1: clusters 0, 3, 4, 5, 8 and 12, expressing MACRO [33] and IL1B [34]; M2: clusters 2, 7, expressing MRC1 [35, 36]; and M2/M1: clusters 1, 6, 9, 10 and 11, expressing both MRC1 and IL1B. Difference expressed genes analysis of the macrophage subtypes indicated that M2/M1 macrophages may be more similar to M2 macrophages than to M1 macrophages (Figure 2C, 2D).

(A, B) The macrophage population was divided into 13 clusters, which were finally categorized into 3 cell populations (M1, M2, M2/M1). (C) Differential gene analysis among different subtypes of macrophages: The X-axis represents the macrophage subtypes, and the Y-axis represents the differentially expressed genes. (D) Violin plot of macrophage marker genes.

Figure 2. (A, B) The macrophage population was divided into 13 clusters, which were finally categorized into 3 cell populations (M1, M2, M2/M1). (C) Differential gene analysis among different subtypes of macrophages: The X-axis represents the macrophage subtypes, and the Y-axis represents the differentially expressed genes. (D) Violin plot of macrophage marker genes.

Macrophage trajectory analysis

According to the trajectory analysis of M1 macrophages and M2 macrophages, we found that M2 macrophages can transform into M1 macrophages over time and that this trajectory can be divided into 5 states (Figure 3A3C). By analyzing the cell trajectory (cell trajectory direction: from left to right), we divided the genes into two clusters. The expression of MRC1, the marker gene of M2 macrophages, decreased with time in cluster 2; MARCO and IL1B, the marker genes of M1 macrophages, increased with time in cluster 1 (Figure 3D3G). KEGG analysis was performed on the genes of cluster 1, and a total of 11 signaling pathways were obtained. Among them, the IL17 signaling pathway, HIF-1 signaling pathway, and PPAR signaling pathway were closely related to the formation of atherosclerosis.

(A–C) Trajectory analysis of M1 macrophages and M2 macrophages. (D) The trajectory analysis of the heatmap of M1 macrophages versus M2 macrophages: The X-axis represents the timeline of trajectory analysis, the left Y-axis represents the KEGG enrichment results, and the right Y-axis represents the differentially expressed genes between the two clusters. (E–G) Pseudotime analysis of genes (MRC1, IL1B, and MARCO): The X-axis represents the cell of trajectory analysis, and the Y-axis represents the relative expression of the gene.

Figure 3. (AC) Trajectory analysis of M1 macrophages and M2 macrophages. (D) The trajectory analysis of the heatmap of M1 macrophages versus M2 macrophages: The X-axis represents the timeline of trajectory analysis, the left Y-axis represents the KEGG enrichment results, and the right Y-axis represents the differentially expressed genes between the two clusters. (EG) Pseudotime analysis of genes (MRC1, IL1B, and MARCO): The X-axis represents the cell of trajectory analysis, and the Y-axis represents the relative expression of the gene.

According to the trajectory analysis of M1 macrophages and M2/M1 macrophages, we found that M2/M1 macrophages can transform into M1 macrophages over time and that the trajectory could be divided into 5 states (Figure 4A4C). By analyzing the cell trajectory (cell trajectory direction: from left to right), we divided the genes into two clusters. The expression of MRC1, the marker gene of M2 macrophages, decreased with time in cluster 2; MARCO and IL1B, the marker genes of M1 macrophages, increased with time in cluster 1 (Figure 4D4G). KEGG analysis was performed on the genes of cluster 1, and a total of 11 signaling pathways were obtained. Among them, the IL17 signaling pathway, HIF-1 signaling pathway, PPAR signaling pathway and fluid shear stress and atherosclerosis were closely related to the formation of atherosclerosis.

(A–C) Trajectory analysis of M1 macrophages and M2/M1 macrophages. (D) The trajectory analysis of the heatmap of M1 macrophages versus M2/M1 macrophages. The X-axis represents the timeline of trajectory analysis, the left Y-axis represents the KEGG enrichment results, and the right Y-axis represents the differentially expressed genes between the two clusters. (E–G) Pseudotime analysis of genes (MRC1, IL1B, and MARCO). The X-axis represents the cell of trajectory analysis, and the Y-axis represents the relative expression of the gene.

Figure 4. (AC) Trajectory analysis of M1 macrophages and M2/M1 macrophages. (D) The trajectory analysis of the heatmap of M1 macrophages versus M2/M1 macrophages. The X-axis represents the timeline of trajectory analysis, the left Y-axis represents the KEGG enrichment results, and the right Y-axis represents the differentially expressed genes between the two clusters. (EG) Pseudotime analysis of genes (MRC1, IL1B, and MARCO). The X-axis represents the cell of trajectory analysis, and the Y-axis represents the relative expression of the gene.

According to the trajectory analysis of M2 macrophages and M2/M1 macrophages, we found that M2/M1 macrophages can transform into M2 macrophages over time and that the trajectory could be divided into 5 states (Figure 5A5C). By analyzing the cell trajectory (cell trajectory direction: from left to right), we divided the genes into two clusters. The expression of MRC1, the marker gene of M2 macrophages, decreased with time in cluster 2; IL1B, the marker gene of M1 macrophages, increased with time in cluster 1 (Figure 5D5F). KEGG analysis was performed on the genes of cluster 1, and a total of 6 signaling pathways were obtained. These signaling pathways are mainly related to chemokines. In the single-cell transcriptome, IL1RN was upregulated during the transformation of M2 macrophages into M1 macrophages; SPARCL1, TAGLN and EMP2 were downregulated during the transformation of M2/M1 macrophages into M1 macrophages; IL-1RN was upregulated during the transformation from M2/M1 macrophages to M1 macrophages; NRP1 was upregulated during the transformation from M2/M1 macrophages to M2 macrophages; and SPARCL1, TAGLN and ACTA2 were downregulated during the transformation of M2/M1 macrophages into M2 macrophages.

(A–C) Trajectory analysis of M2 macrophages and M2/M1 macrophages. (D) The trajectory analysis of the heatmap of M2 macrophages versus M2/M1 macrophages. The X-axis represents the timeline of trajectory analysis, the left Y-axis represents the KEGG enrichment results, and the right Y-axis represents the differentially expressed genes between the two clusters. (E, F) Pseudotime analysis of genes (MRC1, IL1B, and MARCO). The X-axis represents the cell of trajectory analysis, and the Y-axis represents the relative expression of the gene.

Figure 5. (AC) Trajectory analysis of M2 macrophages and M2/M1 macrophages. (D) The trajectory analysis of the heatmap of M2 macrophages versus M2/M1 macrophages. The X-axis represents the timeline of trajectory analysis, the left Y-axis represents the KEGG enrichment results, and the right Y-axis represents the differentially expressed genes between the two clusters. (E, F) Pseudotime analysis of genes (MRC1, IL1B, and MARCO). The X-axis represents the cell of trajectory analysis, and the Y-axis represents the relative expression of the gene.

Disease prediction model results

A total of 1086 AS-characteristic genes (module membership vs. gene significance cor=0.87, p<1e-200) from WGCNA and 223 DEGs in the macrophage trajectory analysis were obtained. Of these, we identified 50 intersecting genes (Figure 6A6C). The comparison of the accuracy of SVM and RF in screening disease-characteristic genes indicated that the accuracy of RF was higher (Figure 6D6F). We used the RF analysis method to obtain a total of 6 disease-characteristic genes (IL1RN, TAGLN, SPARCL1, NRP1, EMP2, and ACTA2) (Figure 6G, 6H).

(A) The leftmost color block represents the module, and the rightmost color bar represents the correlation range. In the heatmap in the middle part, the darker the color is, the higher the correlation. Red indicates a positive correlation, and blue indicates a negative correlation. The numbers in each cell indicate relevance and significance. The X-axis represents the sample type. (B) A scatterplot of gene significance (GS) for treat vs. module membership in the turquoise module. There is a highly significant correlation between GS and MM in the module. (C) The left circle represents the disease-characteristic genes screened using the WGCNA method, and the right circle represents the characteristic genes that change most clearly over time between macrophage subtypes. The intersection of the two circles represents the intersecting genes. (D, E) Boxplot of the residual and reserve cumulative distribution of the residual. (F) The ROC curve shows the difference between SVM and RF. (G, H) RF analysis results and screening for important genes.

Figure 6. (A) The leftmost color block represents the module, and the rightmost color bar represents the correlation range. In the heatmap in the middle part, the darker the color is, the higher the correlation. Red indicates a positive correlation, and blue indicates a negative correlation. The numbers in each cell indicate relevance and significance. The X-axis represents the sample type. (B) A scatterplot of gene significance (GS) for treat vs. module membership in the turquoise module. There is a highly significant correlation between GS and MM in the module. (C) The left circle represents the disease-characteristic genes screened using the WGCNA method, and the right circle represents the characteristic genes that change most clearly over time between macrophage subtypes. The intersection of the two circles represents the intersecting genes. (D, E) Boxplot of the residual and reserve cumulative distribution of the residual. (F) The ROC curve shows the difference between SVM and RF. (G, H) RF analysis results and screening for important genes.

The nomogram indicated that the probability of atherosclerosis increased to 90% when the total score of 6 disease-characteristic genes reached 160 and 10% when it reached 120 (Figure 7A). The calibration curve suggested that the model had reasonable accuracy in predicting the incidence of atherosclerosis (Figure 7B). DCA once again proved that the model has clinical utility (Figure 7C). The clinical impact curve showed that the benefit rate of the model was higher when the number of high-risk factors for atherosclerosis was smaller (Figure 7D).

(A) A nomogram was created to represent the disease model. It uses the X-axis to display the expression of a single gene, as well as the score scale of a single gene, the total score scale of all genes, and the disease incidence scale. Meanwhile, the Y-axis shows individual genes, points, total points, and risk of disease. (B) A graph was used to plot the predicted event rate (Predicted Probability) on the abscissa and the observed actual event rate (Actual Rate) on the ordinate, ranging from 0 to 1. This can be interpreted as the event rate in percentage. The diagonal dashed line serves as the reference line, representing the scenario where the predicted value equals the actual value. (C) The DCA graph employs the threshold probability (ThresholdProbability) on the abscissa and the net profit rate after subtracting the disadvantages on the vertical axis. (D) A graph was used to represent the high-risk threshold and benefit rate on the abscissa, and the number of high risks on the ordinate.

Figure 7. (A) A nomogram was created to represent the disease model. It uses the X-axis to display the expression of a single gene, as well as the score scale of a single gene, the total score scale of all genes, and the disease incidence scale. Meanwhile, the Y-axis shows individual genes, points, total points, and risk of disease. (B) A graph was used to plot the predicted event rate (Predicted Probability) on the abscissa and the observed actual event rate (Actual Rate) on the ordinate, ranging from 0 to 1. This can be interpreted as the event rate in percentage. The diagonal dashed line serves as the reference line, representing the scenario where the predicted value equals the actual value. (C) The DCA graph employs the threshold probability (ThresholdProbability) on the abscissa and the net profit rate after subtracting the disadvantages on the vertical axis. (D) A graph was used to represent the high-risk threshold and benefit rate on the abscissa, and the number of high risks on the ordinate.

The NRP1 and IL1RN genes were significantly upregulated in both the training group and test group, while TAGLN, SPARCL1, EMP2 and ACTA2 were significantly downregulated in both the training group and test group (Figure 8A, 8B). The ROC curve values of the six genes in the test group were statistically significant (AUC (IL1RN): 0.899, 95% CI: 0.764-0.990; AUC (NRP1): 0.817, 95% CI: 0.620-0.971; AUC (TAGLN): 0.846, 95% CI: 0.678-0.971; AUC (SPARCL1): 0.825, 95% CI: 0.620-0.988; AUC (EMP2): 0.808, 95% CI: 0.630-0.947; AUC (ACTA2): 0.784, 95% CI: 0.591-0.938) (Figure 8C8H). We constructed an AS prediction model using these six genes, which showed significant statistical significance in both the train group (AUC: 0.909, 95% CI: 0.842-0.967) and test group (AUC: 0.812, 95% CI: 0.630-0.966).

(A) Heatmap of intersecting genes in the training group. (B) Heatmap of intersecting genes in the test group. (C–H) ROC curves of disease signature genes in the test group. (I, J) The area under the curve of the AS prediction model in the train group and test group.

Figure 8. (A) Heatmap of intersecting genes in the training group. (B) Heatmap of intersecting genes in the test group. (CH) ROC curves of disease signature genes in the test group. (I, J) The area under the curve of the AS prediction model in the train group and test group.

Discussion

Atherosclerosis plays an important role in cardiovascular disease. Its pathological mechanisms include dysfunction of the endothelial cell barrier, lipid accumulation, abnormal angiogenesis, and local chronic inflammatory reaction caused by immune cell infiltration (including T cells, B cells, myeloids), among others [3739]. Among numerous immune cells, macrophages are undoubtedly the most important in the evolution of atherosclerosis, and their polarization acts as a double-edged sword in the occurrence and development of atherosclerosis [40]. The polarization of macrophages can also regulate the stability of plaques [41]. In our study, according to the results of pseudotime analysis, we found that M2 macrophages can transform into M1 macrophages, and M2/M1 (MRC1(+) IL1B (+)) macrophages can transform into M1 or M2 macrophages. According to the RF method, a total of six disease-characteristic genes of atherosclerosis were identified: in the single-cell transcriptome, IL1RN was upregulated during the transformation of M2 macrophages into M1 macrophages; SPARCL1, TAGLN and EMP2 were downregulated during the transformation of M2/M1 macrophages into M1 macrophages; IL-1RN was upregulated during the transformation from M2/M1 macrophages to M1 macrophages; NRP1 was upregulated during the transformation from M2/M1 macrophages to M2 macrophages; and SPARCL1, TAGLN and ACTA2 were downregulated during the transformation of M2/M1 macrophages into M2 macrophages. The changes in these genes in the bulk transcriptome were consistent.

The IL1RN gene has been found to have variable numbers of an 86 base pair (bp) tandem repeat in intron 2 [42], and this polymorphism is associated with various inflammatory diseases (systemic lupus erythematosus [43], type 2 diabetes mellitus (T2DM) [44]) and ischemic stroke [45]. A recent study showed that in thoracic aortic dissection, IL1RNHigh macrophages, as a subtype of proinflammatory macrophages, are the main source of MMPs and inflammatory cytokines, which strengthens the formation of thoracic aortic aneurysm [46]. It was also predicted that these macrophages would be considered target cells for the treatment of patients with thoracic aortic aneurysm in the future. This is in line with our research results; in our study, M2 and M2/M1 macrophages were transformed into IL1RNHigh M1 macrophages.

Neuron-1 (NRP1) belongs to the neurotransmitter family [47]. NRP1 is known to be overexpressed by macrophages in acute and chronic inflammation-related diseases (sepsis [48], type II diabetes [49, 50], and metabolic syndrome [51], etc.). Its function is mainly to act as a common receptor of Semaphorin 3A (Sema3A) and vascular endothelial growth factor a (VEGF-A165), thus regulating development and pathological angiogenesis, arteriogenesis and vascular permeability [52]. In atherosclerotic plaques, microvascular growth, or angiogenesis, destroys the stability of the plaque, thus increasing the risk of rupture [53]. In addition, some studies have proven (mainly in the context of cancer) that M2 macrophages can promote the formation of pathological microvessels [54, 55]. In our study, NRP1 was identified as a disease-characteristic gene of atherosclerosis, which indicates that NRP1High M2 macrophages may be a major risk factor for the formation and stability of atherosclerosis.

According to previous studies, ACTA2 is a marker gene of vascular smooth muscle cells [56]. Some studies have shown that VSMCs can take up lipids and become foam cells [57], and undergo phenotypic transfer to macrophage like status in atherosclerotic plaques [58, 59]. In this study, we obtained ACTA2High M2/M1 macrophages and found that they can transform into M1 and M2 macrophages and that the expression of ACTA2 decreased significantly during this period. This indicates that ACTA2High M2/M1 macrophages may be involved in the development of atherosclerosis.

EMP2 is expressed at a high level in a large number of human tissues, including adult ovarian, heart, lung and intestinal tissues and fetal lung tissues [60]. The protein produced by EMP2 is a four-transmembrane protein that has been found to mediate a variety of vascular reactions [61]. In cancer-related studies, we observed that the higher the expression level of EMP2 in various tumor models in vitro, the more obvious the expression of pathological angiogenesis-related pathways was [6265]. Some studies have also shown that the proportion of M1 macrophages in the placenta is significantly increased in EMP2 knockout rats [61]. Based on the above results, we boldly speculate that EMP2High M2/M1 macrophages can not only promote the formation of pathological blood vessels, resulting in the destruction of plaque stability, but also become proinflammatory macrophages (EMP2Low M1 macrophages) through phenotypic transformation.

Secreted acidic cysteine-rich protein-1 (SPARCL1) belongs to the matrix cell protein SPARC family and is an extracellular matrix (ECM) glycoprotein [66]. At present, research on SPARCL1 has mainly been performed in the context of cancers, including lung cancer [67, 68], prostate cancer [69], colon cancer [70], etc. Little is known about its function in the context of cardiovascular disease. One study showed that SPARCL1 can inhibit angiogenesis and support vascular morphogenesis and integrity [71]. We predict that SPACL1High M2/M1 macrophages may play a role in inhibiting the formation of carotid atherosclerosis and maintaining the stability of carotid plaques by inhibiting angiogenesis.

It is well known that the smooth muscle cell marker TAGLN (SM22α) is downregulated in the pathogenesis of atherosclerosis, restenosis, abdominal aortic aneurysm and other arterial diseases [72]. Here, we speculate that TAGLNHigh macrophages (in the M2/M1 macrophage subgroup) may be similar to ACTA2High macrophages and originate from vascular smooth muscle cells. TAGLNHigh M2/M1 macrophages can transform into TAGLNLow M1 and TAGLNLow M2 macrophages, and the expression of TAGLN decreases significantly during this period.

We identified 6 kinds of macrophages that may be related to the occurrence and development of atherosclerosis: IL1RNHigh M1, IL1RNHigh M1, NRP1High M2, ACTA2High M2/M1, EMP2High M1/M1, SPACL1High M2/M1 and TAGLNHigh M2/M1 macrophages. This study is the first to use these macrophage marker genes to construct a model for predicting atherosclerosis. The expression trend of these genes in bulk transcriptome data is consistent with the expression trend of macrophage pseudotime analysis, which again proves the accuracy of these genes as disease-characteristic genes in atherosclerosis.

Conclusions

IL1RNHigh M1 macrophages, NRP1High M2 macrophages, ACTA2High M2/M1 macrophages, EMP2High M1/M1 macrophages, SPACL1High M2/M1 macrophages and TAGLNHigh M2/M1 macrophages play a key role in the occurrence and development of arterial atherosclerosis. These macrophages may become target cells for treating atherosclerosis in the future. These macrophage marker genes can also be used to build a model to predict atherosclerosis.

Author Contributions

The manuscript has been read and approved by all authors. Shize Yang: Conceptualization, Funding acquisition, Project administration, Supervision, Writing – review and editing, Linguistic editing and proofreading. Zemin Tian: Data curation, Formal analysis, Writing – original draft, Visualization, Software, Methodology, Validation, Writing – review and editing. All authors read and approved the final manuscript.

Acknowledgments

We thank the GEO database for allowing us free access.

Conflicts of Interest

In addition, the authors confirmed that their work was not influenced by competing financial interests or personal relationships. All other authors have no conflicts of interest to declare.

Ethical Statement and Consent

Not applicable. The data for this study were obtained from public databases and no additional ethical approval was required.

Funding

No funding.

References

View Full Text Download PDF


Copyright © 2025 Rapamycin Press LLC dba Impact Journals
Impact Journals is a registered trademark of Impact Journals, LLC