A flow diagram summarizing the entire study design is provided in Figure 1.

Based on differentially expressed genes (log₂FC > 2, adj. p < 0.05 by Wilcoxon Rank Sum test) of all cell types identified above, we conducted biological progress and pathway enrichment analysis with Metascape (Zhou et al., 2019) (https://metascape.org/) (Supplementary Figure S1B). Furthermore, we performed GO and KEGG pathway analysis of marker genes for macrophages, fibroblasts, fibromyocytes, endothelial cells, and smooth muscle cells, using clusterProfiler package (v3.14.3) in R (Supplementary Figures S1C, D). We selected GO-enriched annotations that are closely related to atherosclerotic disease, and visualized significant results in R. The results of the significant enrichment of KEGG pathway were also visualized in R.

We gathered 3,000 highly variable genes from macrophages, fibroblasts, fibromyocytes, endothelial cells and smooth muscle cells based on gene expression matrix of human single-cell sequencing dataset, including 2,066, 2,048, 1,178, 2,504, and 1,518 cells, respectively. These five types of cells were applied to predict transcriptional regulators independently, using the SCENIC (version 1.1.2-2) pipeline workflow (Aibar et al., 2017) in R (Figure 3A). After filtering genes and cells with the conditions of default parameters, we collected potential transcription factors (TFs) using GENIE3 (version 1.8.0) and selected potential regulons (direct binding targets) based on the RcisTarget database 1.6.0 (hg19 motif databases of 500 bp upstream and TSS centered 10 kb) (Figure 3A). Next, we used AUCell (version 1.8.0) to evaluate and rank the activities of all predicted transcriptional regulators and regulons (Figure 3B). Thus, highly reliable regulatory factors and corresponding optimal motifs were obtained. We used Cytoscape (version 3.2.1) to visualize the enriched transcription factor regulatory network between different cells (Figure 3A).

Cell-cell communication analysis was performed with 3,000 highly variable genes of all cells from human single-cell sequencing dataset using CellPhoneDB software (version 2.1.2) (Efremova et al., 2020) in python 3.7. We adopted statistical methods to analyze the expression of obtained ligand-receptor pairs among different cell types, and only pairs with significant difference (p < 0.05, one-sided Wilcoxon rank-sum test) in average expression levels were selected for subsequent analysis and visualization (Supplementary Table S2).

Before starting the pseudotime analysis, we extracted fibroblasts and fibromyocytes from the mouse single-cell sequencing datasets. After merging the data of three stages including 0 weeks, 8 weeks, and 16 weeks, we reconstructed a sub-dataset for single-cell trajectories analysis. In downstream analysis, Seurat (version 3.1.2) and Monocle2 package (version 2.14.0) (Qiu et al., 2017) in R were used for finding marker genes (log₂FC > 1.5, adj. p < 0.05, Wilcoxon Rank Sum test) of cell clusters and ordering the cells by the expression of these genes. Finally, cells were divided into different states to imitate the differentiation of cells, and the expression dynamics of specific gene were visualized along the differentiation trajectories (Figure 5).

To learn more information, we performed bulk RNA-seq analysis of human atherosclerotic plaques (GSE120521) and human aortic dissection (AD) tissues (GSE147026). Plaques were obtained at carotid endarterectomy in symptomatic patients, and the visible area where the plaques rupture was dissociated and defined as unstable plaques. Genes with FDR (False Discovery Rate) values less than 0.05 were considered to have significant difference in dissection and normal arteries, and the results were visualized by volcano and heat maps (Figures 6A–C, Supplementary Figures S4A, B). In addition, we imported the FPKM (Fragments Per Kilobase per Million) expression values into R, which were obtained from RNA sequencing of human stable plaques and unstable plaques. Low-expression genes were dropped, and differential expression analysis was performed using the limma package (version 3.42.2) (Figure 6C). Similarly, we defined adj. p < 0.05 (Benjamini-Hochberg adjusted p-value) as the threshold for the significant differential expression of genes in stable and unstable plaques (Figure 6C).

We selected the first 3,000 highly variable genes from human atherosclerotic artery single cell sequencing data, and chose fibroblasts and fibromyocytes for gene set enrichment analysis using GSEA software (Subramanian et al., 2007). Cells of each type were sorted according to the specific gene expression level from high to low, and then divided into specific gene high expression group and low expression group. Then we compared the different expression patterns of these sorted data sets, with the guidance of the defined gene set from KEGG pathway analysis (Figures 6D, 6E; Supplementary Figure S4C).

The development of atherosclerosis and the process of plaque formation involve the participation of various cells. To understand the mechanism differences between all types of cells, we partitioned cellular clusters and constructed a single-cell map based on single-cell RNA sequencing data from human atherosclerotic artery tissues with severely pathological plaques. After critical quality filtering, we obtained 11 756 cells for the following graph-based clustering. We performed shared nearest neighbor (SNN) and t-distributed stochastic neighbor embedding (t-SNE), and divided these cells into 11 cell types (Figure 2A) by comparing differential expressed genes with the previous reported specific markers. Fibroblasts were identified on the basis of abundant lumican (LUM), decorin (DCN), type I collagen (COL1A1 and COL1A2), and fibulin 1 (FBLN1). In particular, clusters that specifically express fibronectin 1 (FN1) were defined as “fibromyocytes,” which were reported to be fibroblast-like cells transformed from smooth muscle cells (Wirka et al., 2019). The clusters labeled “smooth muscle cell” were characterized with the high expressions of smooth muscle actin (ACTA2), myosin heavy chain 11 (MYH11), myosin light chain 9 (MYL9), and tropomyosin 2 (TPM2). We defined endothelial cells based on the expressions of platelet and endothelial cell adhesion molecule 1 (PECAM1), von Willebrand factor (VWF), fatty acid binding protein 4 (FABP4), claudin 5 (CLDN5), and interferon alpha inducible protein 27 (IFI27). With the guidance of immune-related factors, we distinguished five kinds of immune cells, including macrophages, T cells, B cells, NK cells, and mast cells. Macrophages specifically express C1q (C1QA, C1QB and C1QC), CD74, and CXCL8. The high expressions of C-C motif chemokine ligand 4 (CCL4) and C-C motif chemokine ligand 5 (CCL5) indicated the existence of NK cells. Due to the similar marker genes in T cells, B cells, and NK cells, we identified NK cells through the annotation of the original published report (Wirka et al., 2019) and the latest literature (Hudspeth et al., 2012). Tryptase alpha/beta 1 (TPSAB1) and tryptase beta 2 (TPSB2) were distinctively expressed in mast cells. T cells and B cells were distinguished by CD3D and CD52, respectively. Genes encoding immunoglobulin (IGKC, IGHM, IGHA1, IGHG3, IGLC2, and IGLC3) were detected in plasma cells. Neurons were classified in terms of the predominant expression of proteolipid protein 1 (PLP1), which encodes a transmembrane proteolipid protein (the predominant component of myelin). The expressions of the hallmarks used to annotate cell types are shown in a bubble chart (Supplementary Figure S1A).

To gain overall insights into the cellular composition of atherosclerotic arteries, we calculated the proportion of each cell type, which is shown in a pie chart (Figure 2A) along with t-SNE clustering plot. We found that endothelial cells accounted for the largest proportion, up to 21%. Followed by macrophages, fibroblasts, smooth muscle cells and fibromyocytes accounted for 18%, 17%, 13% and 10%, respectively. Endothelial cells are involved in the adhesion of leukocytes (Libby et al., 2011), which is the initial step in the formation of atherosclerosis. Moreover, when the atherosclerosis happens, it contributes to reduce vascular inflammation and inhibits disease development by autophagy (Zhang et al., 2020). Both regulation of immune metabolism by macrophages (Koelwyn et al., 2018) and accumulation of diseased macrophages (Robbins et al., 2013) play the key roles in the formation, growth, and eventual rupture of plaques. Fibroblasts and fibromyocytes accounted for 27% in total, indicating that they may play an important role in maintaining the physiological structure of arteries and promoting the pathological development of arterial disease. This is consistent with the previous reports (Sartore et al., 2001). Vascular smooth muscle cells are essential source of foam cells in atherosclerosis (Pi et al., 2021). The previous research has reported that smooth muscle cells could be transformed into unique “fibromyocytes” in vivo in mouse and human atherosclerotic lesions (Wirka et al., 2019), and the transformation is related to the pathogenesis of atherosclerosis. This resembles with a recent report that both the loss of smooth muscle cells inside the blood vessels and the significant increase in the amount of non-smooth muscle cells are related to calcification and inflammation of the aortic wall (Chen et al., 2020).

We obtained the differentially expressed genes for each cell population by cell grouping with log2FC > 2 and adj. p < 0.05 (Wilcoxon Rank and test), and we selected several differentially expressed genes that are most interesting to display in the bubble (Figure 2B). Furthermore, the biological significance of the above differential expressed genes were explored through gene ontology (GO) and pathway analysis using Metascape (Zhou et al., 2019) (Supplementary Figure S1B). We found that several biological processes related to atherosclerosis including humoral immune response, leukocyte migration, and lymphocyte proliferation were enriched, suggesting the critical role of immune system. In order to understand the biological function of various cells in atherosclerosis, GO gene ontology and KEGG pathway analyses were performed on cell-type-specific genes in smooth muscle cells, endothelial cells, fibroblasts, fibromyocytes, and macrophages (Supplementary Figures S1C, D). The data showed that C1QA, C1QB, and C1QC were involved in the complement coagulation cascade (hsa04610), HLA-DPA1 and CD74 were related to the antigen processing and presentation pathway (hsa04612), and SERPINF1 was connected to the biological process of inflammatory response and cell motility.

We selected 12 more significant or newly discovered genes of interest from the differential expressed genes as the following research objects. LUM was mainly expressed in fibroblasts and fibromyocytes. SERPINF1 and MYOC were enriched in fibroblasts. CD74, SRGN, C1QC, HLA-DPA1, and HLA-DRB5 were abundantly expressed in macrophages. IGKC, IGLC2, IGHG3, and IGHA1 were specifically expressed in plasma cells. At the same time, we also found that SERPINF1 is enriched in fibroblasts in another human atherosclerotic single-cell RNA-seq data analysis (Supplementary Figure S2; Supplementary Table S1). The t-SNE plots (Figure 2C) exhibited the distribution of these genes, and the violin plots (Figure 2D) showed their expression levels in various cells. Previous reports indicated that LUM was expressed in smooth muscle cells in human coronary atherosclerosis and promoted the formation of collagen fibers (Onda et al., 2002). CD74 had also been reported as a new therapeutic target to reduce atherosclerotic inflammation (Martín-Ventura et al., 2009). C1q was involved in regulating macrophage molecular signals and inflammatory responses (Ho et al., 2016). Although the role of SERPINF1 gene in atherosclerosis had not been reported yet, the pigment epithelium-derived factor (PEDF) encoded by SERPINF1 gene had been widely reported in preventing atherosclerosis (Wen et al., 2017b; Wang et al., 2019b). Due to anti-inflammatory, anti-oxidation, anti-angiogenesis, anti-thrombosis, and anti-tumorigenic properties of PEDF, a large number of studies have been conducted in recent years, but the specific mechanism of the anti-atherosclerosis is still unclear. The current literature mainly focused on the expression of SERPINF1 in smooth muscle cells and endothelial cells, but our data showed that SERPINF1 was specifically expressed in fibroblasts. Our another data analysis also found that SERPINF1 was certainly expressed within human fibroblasts (Supplementary Figure S2). With more and more literature suggesting that fibroblasts play the important role in atherosclerosis, our finding will reveal a potential mechanism of SERPINF1 in the atherosclerosis.

To investigate the transcriptional regulatory network underlying the main cell types, the transcription factor (TF) analysis was conducted based on single-cell expression matrix from human atherosclerotic coronary arteries using SCENIC(Aibar et al., 2017). 43 TFs including 19 in macrophages, 18 in endothelial cells, and 6 in fibroblasts were predicted. Subsequently, we focused on genes that were differential expressed between different cell types to construct transcription factor-target gene network (Figure 3A) and identify the best binding motif of highly confident TFs (Figure 3B). The regulatory network showed that some genes dominantly detected in the same cell type were regulated by the same transcription factor. Such as DCN, SERPINF1, and CFD highly expressed in fibroblasts were subject to CEBPD. MAF was in charge of C1QA, C1QB, and C1QC, which were mainly distributed in macrophages. It had been verified in ApoE (−/−) mice that the over expression of DCN could reduce inflammation in atherosclerotic plaques, thereby decelerating the progression of atherosclerosis (Al Haj Zen et al., 2006). Moreover, both genes specifically expressed in fibroblasts and genes predominantly expressed in macrophages are supervised by the same upstream TFs. For example, both LUM and CD74 are driven by JUN. In addition, IRF7 is responsible for ID1 distributed in endothelial cells and SERPINF1 in fibroblasts. Studies had also shown that ID1 might affect the development of atherosclerosis through downregulation of low-density lipoprotein receptor (LDLR), uptaking lipids in endothelial cells (Zhang et al., 2018).

To study the intercellular communication in arterial lesion tissues, we constructed cell-cell communication networks based on the single-cell RNA sequencing dataset of human atherosclerotic coronary arteries. We totally obtained 219 potential ligand-receptor pairs across 11 cell types using CellPhoneDB (Efremova et al., 2020) software, and then the significant ligand-receptor pairs (p < 0.05, one-sided Wilcoxon rank-sum test) were visualized the interaction networks according to the cell types (Supplementary Table S2). As the network shows, there are more ligand-receptor pairs among fibroblasts, fibromyocytes, macrophages, endothelial cells, and smooth muscle cells, compared with other cell types. We also found that the interactions of macrophages with other cell types were significant, and 30 ligand-receptor pairs were detected between fibroblasts and macrophages, indicating their key role in intercellular communication. Combined with the TF analysis, we speculated that fibroblasts played the vital role in atherosclerosis, and it may affect the biological changes of endothelial cells, macrophages, and smooth muscle cells.

Considering that a simple advanced disease model did not reflect the dynamic changes of the proportion of various cells types and the expressions of differential genes during disease development, we conducted similar analysis of single-cell sequencing data from 0, 8, and 16 weeks ApoE^−/− atherosclerotic mice. According to published research, 0-week mice were defined as the baseline time (Wirka et al., 2019). ApoE^−/− mice having a high-cholesterol diet formed foam cells at 8-week, and formed advanced fibrous plaques at 15-weeks (Nakashima et al., 1994). Therefore, the three periods of mice can be a further explanation for the single-cell sequencing analysis of more severely human atherosclerotic lesions tissues. There were 6,836 cells retained after removing the low quality cells and genes in 0 weeks 7,326 and 2,740 cells were gathered in 8 weeks and 16 weeks, respectively. Globally, 9 main cell types were identified in three stages (Figures 4A–C) through supervising the expression levels of canonical cell-type-specific markers. Fibroblasts were identified by the high expression of Pil6, Comp, Smoc2, Clec3b, and Chad. Extracellular matrix (ECM) genes, such as Col1a2, Lum, and Fn1 were abundantly expressed in fibromyocytes. Genes encoding actin, such as Myh11, Myl9, Tpm2, and Acta2 were specifically expressed in smooth muscle cells. Endothelial cells were established on the basis of striking expressions of Pecam1, Vwf, Fabp4, and Cldn5. As expected, C1qa, C1qb, C1qc, and Cd74 were predominantly expressed in macrophages. Cd3g and Cd3d were specifically expressed in T cells. Neuron significantly expressed Plp1 and Egf18. In addition, pericytes were defined by the presence of Notch3. Epithelial cells were confirmed by the specific expression of Upk3b and Cadm4. Both pericytes (Juchem et al., 2010) and epithelial cells (Wilhelmson et al., 2018) had been reported in the pathological process of atherosclerosis.

We showed the expressions of several important genes in the violin diagram (Figures 4A–C) to explain their dynamic changes as the arterial disease develops. The specific distribution of these genes in various disease periods was also described in the scatter plots (Figure 4D). Since we had not defined plasma cells in the single-cell data of mice, the dynamics of related genes remained uncertain. We found that Serpinf1 was expressed in fibroblasts and fibromyocytes at 8 and 16 weeks, and the expression was still dominant in fibroblasts. Lum was highly expressed in fibroblasts and fibromyocytes. Myoc was also mainly expressed in fibroblasts. C1qc, Cd74, and Srgn were highly expressed in macrophages. Moreover, a large amount of Cd74 and Srgn were expressed in T cells, and abundant Srgn was also detected in endothelial cells. Meanwhile, Serpinf1 were also found to be expressed in fibroblasts of additional mice atherosclerotic single-cell RNA-seq data analysis (Supplementary Figures S3A–C; Supplementary Table S3), and was highly expressed at 8 and 22-weeks of atherosclerosis (Supplementary Figure S3D).

Thus, we severally analyzed the cell composition active of each cell type in three stages (Figure 4E). By calculating the ratio of each cell type in different periods (0-week, 8-weeks, and 16-weeks), we found smooth muscle cells (53%, 37%, and 29%), fibroblasts (37%, 34%, and 34%), endothelial cells (5%, 7%, and 6%), and macrophages (3%, 9%, and 9%) were still the dominant cell types, which was completely consistent with the above human results in our study. With the disease progressing, the proportion of fibromyocytes elevated significantly, but the ratio of smooth muscle cells declined obviously. We supposed that this disorder was strongly related to the pathological process of the disease, which might be caused by the phenotype modulation of smooth muscle cells into fibroblast-like cells. Our view was also supported by the report of Robert CW et al. (Wirka et al., 2019).

We sought to gain in-depth insights of the dynamic shifts of the above discovered potential biomarker through conducting further pseudotime analysis of single-cell sequencing dada. According to our study, Serpinf1 was highly expressed in fibroblasts and fibromyocytes, and previous research also depicted that these cells had the important role in atherosclerosis. To obtain independent evidence of the expression dynamics of Serpinf1 in fibroblasts and fibromyocytes, we constructed the predict single-cell differentiation trajectory according to the expression of marker genes, which was composed of pseudo-temporally ordered fibroblasts and fibromyocytes across the three stages (0-week, 8-weeks, and 16-weeks) (Figure 5A). Finally, 3 differentiation branch points were identified on the trajectory tree, and the cells were divided into 7 states (Figure 5C). Combined with the pseudotime of trajectory (Figure 5B), we concluded that state 1 was the initial state, and state 4 was the final state in the differentiation trajectory tree.

Notably, the expression of Serpinf1 decreased in the pseudo-time trajectory, which was most noticeable in state 4 and state 5 (Figure 5E). We further analyzed the cellular components of state 1, state 4, and state 5. We found that cells at 0-week and 8-weeks mainly concentrated in state 1 and state 5. State 4 was enriched with cells at 16-weeks at the end of differentiation trajectory (Figure 5D). In addition, we analyzed the expression of Serpinf1 in the three stages (Figure 5F). The results revealed the fact that Serpinf1 experienced a significantly diminishing expression (p = 1.688e-191, likelihood ratio tests) with the progression of atherosclerosis. We speculated that the downregulation of Serpinf1 might be a possible reason for the intensification of atherosclerosis. Previous studies had reported that PEDF encoded by SERPINF1 gene was a biomarker of atherosclerosis (Tahara et al., 2011). It had been confirmed in both humans and mice model that PEDF played a protective role against atherosclerosis. PEDF could improve the stability of atherosclerotic plaques through inhibiting macrophage inflammatory response (Wen et al., 2017a), and the deficiency of PEDF accelerated the development of atherosclerosis by promoting the absorption of endothelial fatty acid (Wang et al., 2019a). These reports partly corroborate our suggestion, thus we also conclude that the downregulation of Serpinf1 gene expression in fibroblasts and fibromyocytes results in a shortage of PEDF, which aggravate the development of atherosclerosis and accelerate the formation of plaques.

Regarding gene expression differences between human atherosclerotic tissue and normal arterial tissue, the available sequencing data was scarce so far. Studies had shown that patients with severe atherosclerosis had a higher risk of Stanford type B aortic dissection (Barbetseas et al., 2008). Since atherosclerosis is a high risk factor for aortic dissection, we hope to further explore the possible relationship between atherosclerosis and arterial dissection. In order to verify whether the above potential biomarkers were specifically dysregulated in diseased tissues, we performed the analysis of RNA-seq data of human arterial dissection lesion tissues. Moreover, we also monitored the expression characteristics of potential biomarkers in stable plaques and unstable plaques (ruptured plaques), based on the RNA-seq data of human atherosclerotic carotid plaques. The volcano plot showed the differential expressions of genes (Figures 6A, B), and the heat map specifically describes the expressions of differential expressed genes in each sample (Supplementay Figures S4A, B).

We found that both C1QC and SRGN in the arterial dissection were obviously upregulated compared to the normal group, while the change of SERPINF1 was not significant (Figure 6A). We speculated that the inflammatory response mediated by macrophages played an important role in aortic dissection. Our view was supported by previous study that increased macrophage recruitment and inflammatory response could promote blood vessel remodeling and pathological progress of dissection (Wang et al., 2018b). Compared with stable arterial plaques, SERPINF1, SRGN, and IGKC were evidently upregulated in unstable plaques (Figure 6B). Interestingly, this result was contrary to the downregulation of Serpinf1 as the disease progresses that we previously obtained from the analysis of single-cell data in mice. We speculated that SERPINF1 played the complicated role in preventing atherosclerosis. In the early stages of the disease, the downregulation of SERPINF1 mediated the formation and development of the disease. In the advanced stage, the expression of SERPINF1 was upregulated, which reduced the stability of plaques and promoted the rupture. However, there was no report about the role of SRGN and IGKC in unstable arterial plaques. We propose that the upregulation of SRGN may be a novel biomarker for unstable arterial plaques and arterial dissection, and the increase of IGKC in serum could also be considered as a new potential marker for the formation of severe unstable arterial plaques.

Further, we performed the related gene ontology (GO) annotations and KEGG pathway analyses of SERPINF1, SRGN, C1QC, and IGKC using DAVID (Huang et al., 2009) (https://david.ncifcrf.gov/) (Figure 6C). SRGN was involved in granzyme-mediated apoplotic signaling pathway (GO: 0008626). C1QC was related to complement coagulation cascades (hsa04610). IGKC was associated with Fc-epsilon receptor signaling pathway (GO: 0038095). SERPINF1 participated in the negative regulation of angiogenesis (GO: 0016525).

Finally, we explored the possible signaling pathways involved in atherosclerosis. We performed gene set enrichment analysis of genes with specific expression pattern based on human single-cell data. The significant enrichment results (FDR<0.05) were shown with normalized enrichment score (NES) (Figure 6E; Supplementary Figure S4C). We detected that several signaling pathways were enriched in fibroblasts and fibromyocytes with high expression of SERPINF1, including Jak-STAT signaling pathway, Wnt signaling pathway, NOD-like receptor signaling pathway, TGF-beta signaling pathway, Toll-like receptor signaling pathway, and MAPK signaling pathway. Among them, the Jak-STAT signaling pathway is significantly enriched in fibroblasts expressing SERPINF1. Moreover, some common cardiovascular diseases were enriched in cells with low expression of SERPINF1, such as arrhythmogenic right ventricular cardiomyopathy (ARVC), dilated cardiomyopathy (DCM) and hypertrophic cardiomyopathy (HCM) (Figure 6D).

Our findings were supported by many literatures. In the model of atherosclerotic rabbits, the phosphorylation levels of JAK2 and STAT3 were increased, and the expression of SOCS3 was also upregulated, indicating that inhibition of JAK2/STAT3/SOCS3 signaling can alleviate atherosclerosis (Yang et al., 2020). The latest research indicated that lipopolysaccharide (LPS)-activated scavenger receptors in mouse macrophages were regulated by JAK-STAT-dependent pathways (Hashimoto et al., 2020). Studies had confirmed that Wnt/β-catenin signaling caused VSMC proliferation during intimal thickening in atherosclerosis (Tsaousi et al., 2011). Rui Wang et al. (Wang et al., 2018a) reported that the activation of NLRP3 (Nod-like receptor family) inflammatory bodies promoted the formation of foam cells in vascular smooth muscle cells. There are many reports about TGF-beta signaling pathway in atherosclerosis. TGF-beta1 is an anti-atherosclerotic cytokine with vascular protection. The over expression of TGF-beta1 impeded the progress of plaque growth, stabilized plaque structure, and prevented aortic dilation (Frutkin et al., 2009). In addition, in vivo experiments in mice have also verified the important role of Toll-like receptor signaling in atherosclerosis (Karadimou et al., 2020).