A total of 32 female ApoE-deficient (ApoE^−/–) C57BL/6 mice (10 weeks old) were obtained from the Jackson Laboratory (Bar Harbor, ME). According to method in our previous study (12), all mice underwent combined partial ligation of left renal artery and left common carotid artery (LCCA) to establish the mice model of atherosclerotic plaque rupture with thrombosis. All mice were anesthetized and euthanized 8 weeks after surgery using 4% isoflurane inhalation, in accordance with the AVMA Guidelines for the Euthanasia of Animals (2020 edition). Following euthanasia, the animals were subjected to a whole-body perfusion with phosphate-buffered saline (PBS). Upon examination of the left carotid artery under a microscope, all mice exhibited the presence of atherosclerotic plaques. Notably, thrombosis was observed in 16 of these mice, while the remaining 16 exhibited significant plaque formation without thrombosis. Consequently, the mice were classified into two experimental groups based on the presence of thrombus: the thrombus group (n = 16) and the vulnerable plaque group (n = 16). For further analyses, both groups were subdivided equally into the single-cell sequencing group and the pathological section group, respectively. Carotid artery specimens from both sides (right and left) were collected from all 32 mice for downstream processing. Moreover, the specimens were then distributed into three groups for single-cell RNA sequencing. Right carotid arteries from the 16 mice were combined and served as the control group (A1 group). Left carotid artery samples from eight mice without thrombus formation but with plaques were pooled as A2 group, and left carotid artery samples from eight mice with thrombus formation were served as A3 group. In addition to sequencing, tissue sections from the left carotid arteries of eight mice with plaques but without thrombosis and other eight mice with thrombus formation were embedded for cryosectioning and further histological analysis. All surgeries were performed under the dissecting microscope. Mice were provided with a standard rodent diet and tap water ad libitum during the experiment.

All animal work was performed in accordance with the guidelines on animal care of Shanghai Jiao Tong University School of Medicine. The experimental protocol was approved by the Institutional Animal Care and Use Committee (IACUC) of Shanghai Jiao Tong University School of Medicine (Approved number: JUMC2023-059-B).

Using a dissecting microscope, the isolated vascular tissue was carefully chopped and incubated in a solution containing dissociation enzyme (1 mg/ml collagenase type II, 0.02 mg/ml Deoxyribonuclease I) for 1 h at 37°C. Thereafter, three rounds of centrifugation-washing (500×g for 5 min each) with ice-cold PBS containing 0.04% BSA were performed to remove cell-free RNA released from disrupted cells or residual dissociation reagents. Then, the cell suspension was filtered through a 40 µm sterile cell strainer (JETBIOFIL, Cat#css010040), which not only removed cell clumps but also reduced free-floating RNA fragments. To prepare a single-cell suspension, 0.12% trypsin solution was used to resuspend cells. Based on AOPI Dual-fluoresces counting, cell activity above 90% (primary cell activity above 70%) and concentration between 300 and 600/μl cells were used for sequencing. The single-cell RNA sequencing process was provided by Renji Hospital Affiliated to Shanghai Jiao Tong University School of Medicine. Briefly, the prepared single cell suspension is combined with a mixture of gel beads containing barcode information and enzymes, so as to form GEMs (Gel Bead in EMulsions) containing glue beads (with prefabricated 10× primers), single cells and Master Mix. Single cells were resuspended in PBS with 0.04% BSA and added to each channel. The captured cells were lysed, and the released RNA was barcoded through reverse transcription in individual GEMs. Barcoded cDNA was amplified, and the quality was controlled using Agilent 4200TapeStation System, followed by the scRNA-seq libraries preparation. Then, sequencing was performed on an Illumina Novaseq 6000 sequencer with a pair-end 150 bp (PE150) reading strategy. Finally, the Sequenced Reads were obtained after the sequencing was completed.

The data preprocessing was performed based on sequenced reads obtained above. Briefly, the raw data of single-cell sequencing was obtained in fastq format, followed by quality assessment and alignment to a reference genome. Subsequently, data quality control was conducted to filter out high-quality cells with gene detection counts between 200 and 7,000, as well as a mitochondrial gene ratio below 15%. DecontX algorithm via the celda R package was applied to distinguish between true intracellular gene expression and ambient RNA contamination by leveraging the expression patterns of marker genes and the distribution of gene counts across cells with the contamination probability threshold to 0.05 (13). Then, the data was integrated using the seurat package in R (version: 3.1.2) with Canonical Correlation Analysis (CCA) to correct the batch effects. Next, the data was standardized using global scaling normalization with the LogNormalize method, followed by the linear dimension reduction via Principal Component Analysis (PCA). Finally, the K-Nearest Neighbors (KNN) algorithm from the Seurat package in R was used to cluster cells, followed by Uniform Manifold Approximation and Projection (UMAP) non-linear dimension reduction.

According to genes associated with mouse cell clustering reported in previous literature, the scMCA and CellMarker2.0 software were used to acquire markers associated with mouse cell clustering, followed by annotation. Subsequently, Pearson correlation coefficients were computed between cell subgroups, and a heatmap depicting the correlation was generated to conduct subgroup correlation analysis. To analyze fibroblast cell development, we employed Monocle2 software to perform pseudotime trajectory analysis. Based on prior clustering analyses, three distinct fibroblast subtypes were identified: Fibroblasts, Pi16+ Fibroblasts, and lsg + Fibroblasts. Using Monocle's differential analysis, we selected pseudotime-related genes to capture key expression patterns critical to cell differentiation. By learning the sequence of gene expression changes that each cell undergoes, individual cells were ordered according to their pseudotime values. This allowed us to simulate the dynamic developmental process, aligning each cell along its respective trajectory. Cells were subsequently categorized into multiple differentiation states based on gene expression profiles, and a visual lineage tree was constructed to predict cellular differentiation and developmental pathways. Finally, based on human homologous genes of mouse origin, the iTALK package in R software was employed to analyze receptor-ligand interaction networks among various states and within each state.

Based on different cell cluster, the cluster biomarkers (DEGs) between groups were revealed by using FindMarkers of seurat package in R software. The min.pct = 0.25 and logfc.threshold = 0.25 were used as the cut-off values for cluster biomarkers revealing. Genes with P < 0.05 and | log Fold Change (FC) | > 1 in A1 vs. A2 and A1 vs. A3 were investigated as DEGs in current study. The results of DEGs were visualized by volcano plot using ggplot software (Version: 3.0.0) (14). Finally, the VENN plot analysis was further performed on these DEGs to explore thrombus-related DEGs (TDEGs) using jveen software (15).

GO function and KEGG pathway analyses were performed based on the TDEG using clusterProfiler package (version: 3.16.0) (16) of R. The GO functions including biological process (BP), cellular components (CC), and molecular function (MF). P value < 0.05 was considered as the thresholds for current enrichment analysis. Moreover, according to STING database (version: 11.0) (17), the protein interaction information was extracted, and PPI pairs (median confidence = 0.7) among TDEGs in this study were predicted. Finally, the PPI network was constructed by Cytoscape (version: 3.6.1) software (18).

The common TDEGs in at least 2 cell clusters were considered as hub genes in current study. Then, the miRNAs that target with hub genes were predicted using TargetScan and miRDB in miRWalk software (version: 3.0) (19), followed by the miRNA-hub gene interaction network construction. Finally, the results were visualized by Cytoscape software.

The upstream TFs of hub genes were predicted using iRegulon in Cytoscape software. The P < 0.05 was selected as the thresholds for TFs investigation. Finally, the TFs-hub gene network was visualized by Cytoscape software.

The drugs targeted by homologous human genes of hub genes were screened using Drug-Gene Interaction database (DGIdb, version: 4.0) (20). Based on the drug-target gene relations, the drug-target gene interaction network was constructed using Cytoscape software (version: 3.9.2).

Human carotid artery tissue samples were collected from patients who underwent carotid endarterectomy or autopsy at the Division of Cardiology, Renji Hospital. The study protocol was approved by the Ethics Committee of Renji Hospital, Shanghai Jiao Tong University School of Medicine. Carotid artery tissues from patients without atherosclerotic lesions (n = 8), with stable atherosclerotic plaques (n = 10), and with ruptured atherosclerotic plaques and intraluminal thrombosis (n = 12) were collected.

qPCR was performed using the SYBR Premix Ex Taq II Kit (TaKaRa Bio, Otsu, Japan) on a StepOnePlus Real-Time PCR System (Thermo Fisher Scientific, Waltham, MA, USA) to detect the relative expression levels of the target genes (COL5A1, VCAN, PTGS2, ITGAV, and ITGA8) and the internal reference gene (GAPDH). The relative expression level of each target gene was calculated using the 2^−ΔΔCt method.

The tissue sections of A1, A2 and A3 group used for histological analysis. The result showed that when compared with A1 group (Figure 1A), the plaque and thrombus could be observed in A2 group (Figure 1B) and A3 group (Figure 1C), respectively. The comparative analysis of intima area showed that when compared with A1 group, there was a significant increase of intima area in both A2 group and A3 group (all P < 0.05) (Figure 1D). All results indicated that the mice model established by using LCCA in current study could successfully induce atherosclerotic lesions with plaque disruption associated with lumen thrombosis.

After QC on A1 group, A2 group and A3 group, the high-quality single cells were obtained and kept for downstream analysis. Additionally, we applied Scrublet to identify potential doublets. Overall, the QC process effectively removed low-quality cells and genes, improving the overall quality of current data.

Top 10 principal components (PCs) were revealed by using Elbow-plot analysis (Supplementary Figure 1A). Then, these 10 PCs were used for cell clustering based on UMAP nonlinear dimensionality reduction (Supplementary Figure 1B). The result showed that there were totally 25 clusters revealed by UMAP analysis. The detail information of clustering in each group was showed in Supplementary Figure 1C. Then, the obtained cells were classified based on the marker gene of mice. The results showed that a totally 17 cell subsets were annotated in current study (Figure 2A). In the carotid artery tissues of control (A1), plaque-only (A2), and thrombus (A3) groups. These subsets include not only well-characterized cell types in atherosclerosis, such as endothelial cells, SMCs, Trem2 + macrophages, T cells, but also three fibroblast subtypes, such as Fibroblasts, Pi16+ Fibroblasts, Isg + Fibroblasts, and specialized populations like “Cytokine-simulated Endothelial cells” and “Modulated SMCs.” Furthermore, the detail information for cell types in all three groups (A1, A2 and A3) revealed that pro-inflammatory populations, such as inflammatory macrophages and Isg + Fibroblasts, are nearly absent in A1 but significantly expanded in A2 and A3, while fibroblasts consistently represented the largest cell fraction across all groups, suggesting their central role in both plaque formation and thrombus progression (Figure 2B). In addition, the result of correlation analysis among cell types were showed in Figure 2C.

Fibroblasts are critical for extracellular matrix (ECM) synthesis and tissue repair, but their developmental trajectory and subtype-specific roles in atherosclerotic thrombosis remain unclear. The pseudotime analysis revealed distinct developmental trajectories among the fibroblast cell types. Cells were arranged according to their pseudotime along the inferred trajectories, with early-stage cells situated near the root (Figure 2D). This analysis enabled the identification of several bifurcation points, indicating critical decision points in cell fate determination. The fibroblast cell types displayed differential expression patterns that corresponded to unique developmental states (Figure 2E). The results suggested a clear developmental hierarchy, with distinct fibroblast subtypes, Fibroblasts, Pi16+ Fibroblasts, and lsg + Fibroblasts, distributed along the pseudotime trajectory (Figure 2F). In A2 (plaque-only), most fibroblasts are in the early/intermediate states (Fibroblasts, Pi16+ Fibroblasts), which express high levels of ECM-related genes consistent with plaque structural maintenance. In A3 (thrombus), fibroblasts shift toward the late Isg + Fibroblast state, which is enriched in interferon-stimulated genes and pro-inflammatory cytokines, suggesting a switch from “repair-focused” to “inflammation-amplifying” functions during thrombus formation.

Finally, the receptor-ligand communication relationship among states and in each state were showed in Figure 2G. For example, Fibroblasts/Pi16+ Fibroblasts in A2 secrete ligands like TGF-β1, which signals to SMCs via the TGFBR1/2 receptor to promote ECM synthesis. In A3, Isg + Fibroblasts and Inflammatory Macrophages secrete CXCL10 and CCL2, which bind to CXCR3/CCR2 on T cells and endothelial cells, recruiting immune cells and disrupting endothelial barrier function to facilitate thrombus formation.

The DEGs in each cell subset, compared to other categories, were identified as biological marker genes for corresponding cell subset. The results showed that 812 DEGs were revealed in the A2 vs. A1 comparison, while 1,004 DEGs were identified in the A3 vs. A1 comparison. Meanwhile, a total of 1,004 DEGs were explored between A3 group and A1 group. Due to the large number of genes involved, we presented the results for five representative genes, including COL5A1, VCAN, PTGS2, ITGAV, and ITGA8 (Figure 3A). Based on these DEGs, the VENN plot analysis further revealed 376 TDEGs (Supplementary Figure 2). These TDEGs were mainly assembled in GO functions including cell adhesion (BP, GO: 0007155, Figure 3B), membrane (CC, GO: 0016020, Figure 3C) and protein binding (MF, GO: 0005515, Figure 3D). In addition, these TDEGs were mainly enriched in KEGG pathways like Akt signaling pathway (mmu04151, Figure 3E).

A PPI network was constructed based on TDEGs. The result showed that there were 190 nodes and 319 interactions in current PPI network (Figure 4). According to the degree of these nodes, the Top 20 nodes including IL6, CD4 and ITGAV were selected as hub nodes (Supplementary Table 1 in Supplementary Information).

The common TDEGs in at least 2 cell clusters were enrolled as hub genes, followed by totally 65 miRNAs that target with hub genes predicted. Then, a miRNA-mRNA interaction network was constructed with 4 down-regulated mRNAs, 27 up-regulated mRNAs and 65 miRNAs (Figure 5).

With P < 0.05, a total of 36 TFs of hub genes were predicted. Then, a TFs-hub gene network was visualized by Cytoscape software (Figure 6). The result showed that there were 36 TFs, 23 up-regulated genes and 3 down-regulated genes in current network.

A total 143 drugs targeted by homologous human genes of hub genes were screened using DGIdb. Then, a drug-gene interaction network was constructed with several drug-gene relations such as Aspirin-PTGS2 (Figure 7). The result showed that there were 143 drugs, 13 up-regulated genes and 2 down-regulated genes in current network.

As shown in Figure 8, the levels of five hub genes in clinical subjects were investigated. Compared with the normal group, the expression of all five genes was significantly upregulated in the plaque group, and the upregulation was more pronounced in the plaque rupture with thrombosis group. Among them, COL5A1, VCAN, PTGS2, and ITGAV showed extremely significant increases (P < 0.01); ITGA8 was significantly decreased in the plaque group and thrombosis group (P < 0.01). These findings further confirm that these genes are potential clinical biomarkers for this disease.