In this study, we utilized RNA-seq data related to AS obtained from the GEO cohort (https://www.ncbi.nlm.nih.gov/geo/) via the R package “GEOquery” (version 2.62.2). Specifically, the GSE90074 microarray, consisting of 93 AS and 50 healthy control blood samples, was used as the training set (Table 1), while the GSE27034 microarray, containing 19 peripheral artery disease (PAD) and 18 control peripheral blood mononuclear cell (PBMC) samples, was utilized as the validation set. The GSE159677 high-throughput sequencing dataset, containing three atherosclerotic core plaques and patient-matched adjacent portions of carotid artery tissue from patients undergoing carotid endarterectomy, was utilized as the single-cell sequencing dataset. In addition, we downloaded the GSE59421 microarray dataset from the Gene Expression Omnibus (GEO), which included microRNA (miRNA) expression data from 33 AS and 63 control blood samples. To identify m5C regulators and their associated genes, we referenced published articles and selected 17 m5C regulators, including 11 writers (NOP2, NSUN2, NSUN3, NSUN4, NSUN5, NSUN6, NSUN7, DNMT1, DNMT3A, DNMT3B, and TRDMT1), four erasers (TET1, TET2, TET3, and ALKBH1), and two readers (ALYREF and YBX1) (7–10).

To identify differentially expressed genes (DEGs) between AS and control samples, we utilized the R package “limma” (version 3.50.1) with a p-value threshold of < 0.05 (11). We then visualized the expression of the top 10 upregulated and downregulated DEGs in a heatmap and compared the expression of m5C regulators between AS and healthy control samples to obtain DE-MRs.

Single-sample gene set enrichment analysis (ssGSEA) was used to calculate the score of each sample to evaluate the overall activity of the m5C-related regulatory network using the R package “GSVA” (version 1.42.0) (12). Correlations between DEGs and DE-MRs were calculated using Spearman correlation analysis with |Cor| > 0.3 and p < 0.05 to screen for m5C regulator-related genes (MRRGs). Moreover, the R package “clusterProfiler” (version 4.2.2) was used to enrich MRRGs into Gene Ontology (GO) functions and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (13). Then, the STRING database (https://cn.string-db.org/) and the software “Cytoscape” (version 3.8.2) were used to generate a protein–protein interaction (PPI) network with a confidence score of 0.9 (14). The CytoHubba tool was applied to screen the top 20 hub MRRGs with the highest correlation degree from the network. Finally, the Comparative Toxicogenomics Database (CTD) (http://ctdbase.org/) was used to detect inference scores between the 20 hub MRRGs and AS incidence risk.

The R package “glmnet” (version 4.1-2) was applied to construct a logistic regression model using MRRGs to identify genes significantly associated with AS (p < 0.001 and OR≠1) (15). These genes were defined as feature genes, followed by further filtering using the minimum value of lambda (λ) in the LASSO algorithm. Next, the R package “pROC” (version 1.18.0) was used to plot receiver operating characteristic (ROC) curves based on both the training and validation sets to evaluate the prognostic performance of the model (16). Moreover, to identify miRNAs and transcription factors (TFs) regulating feature genes, a TF–mRNA–miRNA network was constructed using the miRNet database (https://www.mirnet.ca/) based on gene expression levels in the model.

Consensus clustering analysis was performed using the R package “ConsensusClusterPlus” (version 1.58.0) based on the expression of nine feature genes and the optimal k value to classify AS samples in the training set into distinct clusters (17). To verify the separability of sample distributions between clusters, principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) analysis were conducted using the R packages “factoextra” (version 1.0.7), “FactoMineR” (version 2.4.0), and “Rtsne” (version 0.15).

Next, the R package “Weighted Gene Co-expression Network Analysis (WGCNA)” (version 1.70-3) was employed to identify modules with the highest correlation with different clusters using the training set. A hierarchical clustering tree was generated to remove outlier samples, and the most suitable soft-threshold power (β) was selected to construct an adjacency matrix. The adjacency matrix was then transformed into a topological overlap matrix (TOM), and genes with similar expression patterns were clustered into the same module. Furthermore, relationships between co-expression modules and clusters were calculated using Pearson’s correlation coefficient and visualized using a heatmap.

Differentially expressed genes between clusters were obtained using the R package “limma” based on p < 0.05 and |log₂FC| > 0.5 (18). Enrichment analysis was performed to identify GO functions and KEGG pathways associated with these DEGs. Moreover, immune cell composition in AS samples was inferred using the R package “CIBERSORT,” and correlations between immune cell levels and DE-MRs were examined (19).

Candidate biomarkers were identified by intersecting DEGs with genes from WGCNA-selected modules. Expression levels of candidate biomarkers were compared between AS and control samples in both the training and validation sets, and correlation analysis was conducted between feature genes in the LASSO logistic model and biomarkers. In addition, correlations between biomarkers and drugs were analyzed using the Drug–Gene Interaction Database (DGIdb) (http://www.dgidb.org/). Finally, a TF–mRNA–miRNA network was generated using the miRNet database.

Human myeloid leukemia mononuclear cells (THP-1) were obtained from America Type Culture Collection (ATCC, USA). Cells were cultured in RPMI 1640 medium containing 10% FBS and 1% penicillin-streptomycin (Beyotime, Shanghai, China) in a humidified atmosphere of 5% CO₂ at 37°C. THP-1 cells were seeded into 6-well plates at a density of 1×10⁶cells per well and designated as the control group. To construct the foam cell model, THP-1 cells were exposed to phorbol myristate acetate (100nM, Sigma-Aldrich, Merck KGaA) for 48h, and then incubated with 100μg/ml ox-LDL for 24h.

Total RNA was extracted from cultured cells using the High Purity Total RNA Extraction Kit (BioTeKe, Beijing, China) according to the manufacturer’s instructions. RNA concentration and purity were measured using a NanoDrop 3000 spectrophotometer (Thermo Fisher Scientific, Scotts Valley, CA, USA). Complementary DNA was synthesized using the High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher Scientific, Waltham, MA, USA). Primers used in this study were purchased from Sangon Biotech Co., Ltd. (Beijing, China). Transcripts were quantified using SYBR Green (Roche, Mannheim, Germany) on an ABI 7900 Fast Real-Time PCR System (Applied Biosystems, Foster City, CA, USA) for 45 cycles. Relative mRNA expression was analyzed using the 2-△△Ct method. The following primer sequences were used (F: forward; R: reverse):

MCL1 F: 5’-GCTTCGGAAACTGGACATCA-3.

MCL1 R: 5’-GGAAGAACTCCACAAACCCATC-3.

F13A1 F: 5’-GTCACGAGCGTTCACCTGTT-3.

F13A1 R: 5’-TTGTTCTCCTGTGGGTAGCG-3.

RGS2 F: 5’-ACCACAGAGCCTCATGCTAC-3.

RGS2 R: 5’-GACACCACGTTCAGACCACC-3.

TLR8 F: 5’-TGACTTTACATCTTCCCTTCGG-3.

TLR8 R: 5’-AAGTGCGGATTTGTTGATTGTT-3.

TAGAP F: 5’-TGTCATTTGAAGCCCAGAAGG-3.;

TAGAP R: 5’-ATCAGGGTCGTTGCTGTCGTA-3.

The GSE159677 dataset was processed using the R package “Seurat” (version 5.1.0). Quality control was performed by retaining cells with the number of detected genes (nFeature) between 200 and 6000, the total number of RNA molecules (nCount) between 500 and 20000, and the percentage of mitochondrial genes (percent.mt) below 10%. Logarithmic normalization was subsequently applied, and 2000 highly variable genes were selected using the FindVariableFeatures function for downstream analysis. Principal component analysis (PCA) was performed using the ElbowPlot function to determine the optimal number of principal components. The FindNeighbors function was applied to construct a cell neighborhood graph, and the FindClusters function (resolution = 0.4) was used to identify cell subpopulations for fine-grained clustering analysis. Based on the clustering results, the uniform manifold approximation and projection (UMAP) algorithm was employed for dimensionality reduction and visualization. The FindAllMarkers function was then used to identify cluster-specific highly expressed genes (min.pct = 0.25, logfc.threshold = 0.1). Cell clusters were annotated into different cell types using the R package “SingleR” (version 2.2.0) combined with the CellMarker database, and the annotated cell types were visualized using UMAP plots. Furthermore, the expression patterns of five biomarkers and two DE-MRs (NSUN3 and NSUN5) across the identified cell types were visualized using the R package “ggplot2” (version 3.3.5).

293T cells (3.8×10⁶ per 100 mm dish) were seeded into culture plates. After culturing for 24 h, a mixture of plasmids, including the target plasmids (PCDH plasmid and PCDH-NSUN3 plasmid) and additional packaging plasmids (Addgene, Cambridge, MA, USA), was transfected into the cells using Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA). Viruses were harvested and concentrated at 72 h post-transfection. THP-1 cells were seeded into 6-well plates at a density of 2×10⁵ cells per well. The virus containing NSUN3 was used to infect THP-1 cells to induce NSUN3 overexpression, while a control virus served as the negative control. After 24 h of incubation, the viral medium was replaced with fresh medium, and THP-1 cells were cultured for an additional 48 h. Phorbol 12-myristate 13-acetate (PMA) was added to induce differentiation of THP-1 cells into macrophages (M0) for 48 h, followed by ox-LDL treatment to stimulate foam cell formation.

Anti-human NSUN3 antibody (Thermo Fisher Scientific, Lenexa, KS, USA), anti-human NSUN5 antibody (Thermo Fisher Scientific, Lenexa, KS, USA), Annexin V-FITC/PI Staining Detection Kit (BD Pharmingen, San Diego, CA, USA), anti-human CD86 antibody (BioLegend, Beijing, China), and anti-human CD206 antibody (BioLegend, Beijing, China) were used according to the manufacturer’s instructions. After harvesting and washing cells with PBS three times, cells were digested with EDTA-free trypsin. Cells were stained with NSUN3 antibody, NSUN5 antibody, Annexin V-FITC, propidium iodide (PI), CD86 antibody, and CD206 antibody for 25 min in the dark. FlowJo software was used for data acquisition and analysis.

Cells were collected and lysed using the Human XL Cytokine Array Kit (R&D Systems, Minneapolis, MN, USA). Protein extracts were adjusted to a concentration of 1 mg/mL for each test according to the manufacturer’s instructions. After incubation with antibodies and streptavidin–HRP conjugate, dot images were captured using a Bio-Rad ChemiDoc XRS system. Spot intensities were assessed using ImageJ software and normalized to the average values of control spots.

Differentially expressed microRNAs (DE-miRNAs) between AS and control samples were identified using the R package “limma,” with p < 0.05 as the significance threshold (11). Long noncoding RNAs (lncRNAs) with reciprocal relationships were screened by intersecting miRNAs obtained from the miRNet database with DE-miRNAs. Finally, a ceRNA network was constructed based on miRNA and lncRNA interaction information.

Statistical analyses were performed using R software (version 4.1.0). The R packages “ggplot2” (version 3.3.5) and “pheatmap” (version 1.0.12) were used to generate volcano plots, box plots, violin plots, and heatmaps. Forest plots were generated using the R package “forestplot” (version 2.0.1), and ROC curves were plotted using the R package “geomROC” (version 1.0). scRNA-seq data were analyzed using the R packages “Seurat” (version 5.1.0) and “SingleR” (version 2.2.0). Data are presented as mean ± standard deviation (SD), and cellular experiments were repeated 3–5 times. Differences between two groups were analyzed using a two-tailed Student’s t-test. Nonparametric tests were applied for non-normally distributed variables. Spearman correlation analyses between candidate gene expression levels and key clinical continuous parameters were conducted using IBM SPSS Statistics (version 26.0). Statistical analyses were also performed using GraphPad Prism 10.0 (CA, USA). A p-value < 0.05 was considered statistically significant.

The graphical abstract is shown in Figure 1. We performed differential expression analysis on microarray data from GSE90074 (AS = 93 and control = 50 blood mononuclear samples) and identified 2247 DEGs, including 1131 upregulated and 1116 downregulated genes (Figure 2A). The volcano plot in Figure 2B and Supplementary Table S1 display the top 10 DEGs ranked by logFC, including TXLNG2P, KDM5D, RPS4Y2, DDX3V, RPS4Y1, EIF1AY, CEACAM21, BTNL8, NCRNA00185, and UTY as upregulated genes, and NBEA, EFHB, CUL3, XIST, TEKT4P2, LOC641518, SCGB3A1, KCNH8, DMRTC1, and STAP1 as downregulated genes. We also compared the expression of 17 m5C regulators between AS and control samples and identified two m5C regulators, NSUN3 and NSUN5, that showed significant differences (Figure 2C).

To evaluate the overall activity of the m5C-related regulatory network, we used the ssGSEA method to calculate the score of each sample based on GSE90074 (Supplementary Table S2). This method calculated enrichment scores for potential target gene sets regulated by NSUN3 and NSUN5 within each sample, providing a quantitative measure of their regulatory activity. Spearman correlation analysis was then performed between DEGs and DE-MRs, and 546 DEGs showing strong correlations (|Cor| > 0.3) with DE-MRs were identified (Figure 3A). These genes were defined as m5C regulator-related genes (MRRGs). The MRRGs were enriched in 305 Gene Ontology biological process (GO-BP) terms, including “leukocyte migration” and “leukocyte activation in immune response”; 30 GO cellular component (GO-CC) terms, such as “secretory granule membrane” and “ficolin-1-rich granule”; and 15 GO molecular function (GO-MF) terms, including “structural constituent of ribosome” and “immune receptor activity” (Figure 3B). In addition, these MRRGs were significantly enriched in 241 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, including “lipid and atherosclerosis,” “cytokine–cytokine receptor interaction,” “neutrophil extracellular trap formation,” and “chemokine signaling pathway” (Figure 3C). A comprehensive summary of the enrichment analysis is provided in Supplementary Tables S3 and Supplementary Tables S4.

To explore the molecular mechanisms underlying AS incidence, we constructed a protein–protein interaction (PPI) network based on the 546 MRRGs (Supplementary Figure S1). From this network, the top 20 MRRGs ranked by connectivity degree (k) were identified (Supplementary Figure S2; Supplementary Table S5) and defined as hub MRRGs. These included RPS27A, LYN, SYK, TP53, FCER1G, RPS12, RPL18A, RPL15, RPL39, RPL36, RPL27, RPL35A, RPL30, RPL34, RPL29, FGR, Toll-like receptor 4 (TLR4), RPL12, ITGAM, and TYROBP. Notably, three hub MRRGs—RPS27A, LYN, and SYK—were each correlated with more than 18 genes (Figures 3D; Supplementary Figure S3). Although ribosomal proteins frequently appear as hubs in PPI analyses due to their high and coordinated expression, the network also identified several non-ribosomal hub genes with established roles in immune and inflammatory processes, including LYN, SYK, and TLR4, which are directly relevant to AS.

To further assess the relevance of these hub MRRGs to AS, inference scores between the 20 hub MRRGs and AS incidence were calculated using the Comparative Toxicogenomics Database (CTD). TP53 and TLR4 exhibited relatively high inference scores, both exceeding 150 (Figures 3E, F).

Overall, PPI network analysis of m5C-related genes confirmed multiple genes with established roles in immune and inflammatory processes in AS, such as LYN, SYK, and TLR4, and further demonstrated strong associations between core genes, including TP53 and TLR4, and AS based on CTD analysis. From a systems biology perspective, these results indicate that the m5C-related gene network is deeply involved in immune and inflammatory mechanisms underlying AS, providing a foundation for subsequent immune analyses and enabling the identification of hub genes closely associated with both m5C regulation and AS.

To construct a diagnostic model for AS, feature selection was performed on the 546 MRRGs using a LASSO logistic regression model. This analysis identified nine feature genes: DOK1, RRAGC, EFHD2, PRDM1, WIPI1, HLA-B, URB1, PRKX, and VPS52 (Figures 4A–C). To evaluate the diagnostic performance of the model, receiver operating characteristic (ROC) curves were generated for both the training set (GSE90074) and the validation set (GSE27034). The area under the curve (AUC) values for both datasets exceeded 0.7, indicating that the model has a reasonable ability to distinguish peripheral arterial disease (PAD) samples from control samples (Figures 4D, E).

We utilized consensus clustering analysis to classify AS samples into distinct subclusters based on the expression of nine feature genes. Using the consensus clustering heatmap and cumulative distribution function (CDF) curve, we determined an optimal k value of 2, which separated the samples into cluster 1 and cluster 2 (Figures 5A, B). Further, the distribution of samples in cluster 1 was significantly different from that in cluster 2, as observed by principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE) analysis (Figures 5C, D). To identify genes strongly correlated with cluster assignment based on GSE90074, we first generated a hierarchical clustering tree to remove 14 outlier samples, resulting in 79 remaining samples. After setting the soft-thresholding power (β) to 11, 10 gene modules were identified (Figures 5E; Supplementary Figure S4). Correlation analysis showed that the MEyellow module, comprising 643 genes, had the strongest correlation with cluster 1 and cluster 2 (Figure 5F). Therefore, the MEyellow module was selected for further analysis. This finding indicates that the biological processes represented by this module are closely associated with the fundamental heterogeneity among AS patients. We reasoned that genes within this module, which are central to AS pathophysiology, may harbor key drivers underlying intercluster differences.

Subsequently, baseline characteristics between the two clusters were compared to assess potential clinical relevance. As shown in Table 2, significant differences were observed in gender distribution (p = 0.018) and hyperlipidemia prevalence (p = 0.047) between clusters. Specifically, cluster 1 exhibited a higher proportion of females (52.5% vs. 28.3%, p = 0.018) and a lower prevalence of hyperlipidemia (70.0% vs. 86.8%, p = 0.047) than cluster 2. In addition, a significant difference was observed in CXCL5_rank levels (74.5 vs. 52.0, p = 0.031). No significant differences were identified for other clinical parameters, including ancestry, CAD_class, diabetes, hypertension, or BMI (p > 0.05).

These results indicate that molecular subtyping based on the nine-gene signature stratifies AS patients into two subgroups with distinct clinical characteristics. The identification of these clusters supports the biological relevance of the selected feature genes and reveals clinically meaningful patient subgroups.

To further explore intercluster heterogeneity, differentially expressed genes between cluster 1 and cluster 2 were identified. A total of 541 DEGs were detected, including 292 upregulated and 249 downregulated genes. The heatmap illustrates the top 10 genes with the highest levels of upregulation or downregulation (Figure 6A; Supplementary Table S6). Functional enrichment analysis showed that these DEGs were enriched in 145 GO biological process terms, such as immune response–regulating signaling pathways and positive regulation of cytokine production; 12 GO cellular component terms, including secretory granule lumen and cytoplasmic vesicle lumen; and 10 GO molecular function terms, such as immune receptor activity and phosphoric ester hydrolase activity (Figure 6B). In addition, five KEGG pathways—transcriptional misregulation in cancer, tumor necrosis factor (TNF) signaling pathway, mitogen-activated protein kinase (MAPK) signaling pathway, lipid and AS, and cytokine–cytokine receptor interaction—were significantly enriched (Figure 6C). Detailed enrichment results are provided in Supplementary Tables S7, S8.

To further clarify biological differences between the two clusters, KEGG pathway enrichment analysis was performed on cluster-specific DEGs (Supplementary Tables S9, S10). Focusing on pathways related to inflammation and hyperlipidemia, seven significantly upregulated pathways and one downregulated pathway were identified.

Genes upregulated in cluster 1 were predominantly enriched in immune and inflammatory pathways, including phagosome, Fc gamma R–mediated phagocytosis, chemokine signaling, neutrophil extracellular trap (NET) formation, and Toll-like receptor signaling pathways This coordinated enrichment underscored a pathogenesis centered on robust innate immune activation, neutrophil-driven inflammation, and sustained inflammatory responses, providing a molecular basis for the elevated CXCL5 levels observed clinically in this subgroup.

In contrast, cluster 2 showed downregulation of ribosomal protein–related pathways, suggesting alterations in fundamental cellular processes such as protein synthesis.

Collectively, these results indicate mechanistic divergence between the two clusters, with cluster 1 representing an immune–inflammatory–driven subtype and cluster 2 characterized by metabolic and ribosomal dysregulation, highlighting heterogeneity in AS pathogenesis.

To investigate the immunological basis underlying molecular heterogeneity, immune cell proportions were inferred computationally. CIBERSORT analysis revealed significant differences in the relative abundances of five immune cell types, including activated dendritic cells, macrophages M0, macrophages M2, resting mast cells, and neutrophils, between cluster 1 and cluster 2 (Figure 7A). Next, correlations between immune cell levels and DE-MRs were evaluated. Both NSUN3 and NSUN5 showed the strongest correlations with macrophages M0 (R = 0.35 and R = −0.15, respectively) (Figures 7B, C).

To pinpoint core biomarkers driving subtype differences, we identified 541 DEGs between cluster 1 and cluster 2. To refine this list, we prioritized genes that not only differed between clusters but also resided within the MEyellow co-expression network—the network most central to AS heterogeneity. We reasoned that intersecting the 541 DEGs with the 643 genes of the MEyellow WGCNA module would yield a high-confidence set of biomarkers that were both differentially regulated and integral to the key AS-related network (Figure 8A). Differential expression analysis between AS and control samples based on the training and validation sets identified five biomarkers—MCL1, F13A1, RGS2, TLR8, and TAGAP—that were significantly different between the two groups. Notably, expression trends of these biomarkers were highly consistent between the training and validation sets (Figures 8B, C).

To verify biomarker expression in AS, THP-1 cells were treated with oxidized low-density lipoprotein (ox-LDL) to generate macrophage-derived foam cells. qRT-PCR analysis revealed that expression levels of MCL1, F13A1, RGS2, TLR8, and TAGAP were significantly increased in the ox-LDL group (p < 0.05) (Figures 8D, E). Moreover, correlation heatmap analysis demonstrated strong correlations between feature genes and biomarkers (Supplementary Figure S5). Drug–gene interaction analysis identified 53 drugs, including aspirin and docetaxel, showing reciprocal associations with biomarkers including MCL1, F13A1, RGS2, and TLR8 (Supplementary Figure S6). Finally, a TF–mRNA–miRNA network comprising 14 transcription factors (TFs), four biomarkers, and 16 miRNAs was constructed. Notably, hsa-miR-20a-5p, hsa-miR-93-5p, hsa-miR-106b-5p, and hsa-miR-155-5p regulated two or more biomarkers, and 12 TFs directly regulated MCL1 (Supplementary Figure S7).

scRNA-seq analysis of tissue samples from the GSE159677 dataset was performed to resolve the cellular composition of the AS immune microenvironment and assess expression patterns of two DE-MRs across cell subsets. Following quality control, 2000 highly variable genes were selected for downstream analysis (Figure 9A). Principal component analysis indicated that 40 principal components captured the main biological variation and were used for clustering. UMAP analysis identified 23 distinct cell subpopulations, which were annotated into nine major cell types: macrophages, CD4^- T cells, CD8^- T cells, B cells, natural killer (NK) cells, mast cells, endothelial cells, vascular smooth muscle cells (VSMCs), and fibroblasts (Figure 9B). These results clearly delineated the immune cell landscape in the carotid atherosclerotic plaques. The five biomarkers showed high expression in macrophages (Figure 9C), suggesting a key role for macrophages in AS pathogenesis. Further analysis showed that NSUN3 expression was significantly increased in AS samples compared with control samples (p < 0.05) (Figures 9D, E).

The relative exprssion level (measured by mean fluorescence intensity, MFI) of NSUN3 was statistically increased in the foam cell group (896.50±79.29) compared with the control cell group (526.50±107.28) (p < 0.05), whereas no significant change was observed in NSUN5 (832.25±98.05) vs. (894.00±52.33) (p > 0.05) (Figures 10A, B).

To investigate the functional role of NSUN3, THP-1 macrophages were transfected with an NSUN3-overexpressing vector or a control vector, followed by ox-LDL treatment to establish a foam cell model (Figures 11A, B). Flow cytometry analysis showed that the apoptosis rate in the NSUN3 vector group (15.79% ± 1.12%) was significantly lower than that in the control vector group (29.17% ± 2.63%) (independent Student’s t-test, n=3, p = 0.001). In addition, the proportion of proinflammatory (M1) macrophages indicated by CD86⁺ was significantly reduced, whereas the proportion of anti-inflammatory (M2) macrophages indicated by CD206⁺ was markedly increased in the NSUN3 vector group compared with the control group (Figures 11C, D). Cytokine profiling using a Proteome Profiler Human XL Cytokine Array demonstrated that NSUN3 overexpression significantly increased anti-inflammatory cytokines, including interleukin (IL)-10 and C–C motif chemokine ligand (CCL)17, while decreasing proinflammatory cytokines, including IL-1β, IL-6, and tumor necrosis factor α (TNF-α) (Figures 12A, B; Supplementary Figure S8).

A total of 105 differentially expressed miRNAs (DE-miRNAs) were identified between AS cases (AS = 33) and controls (control = 63) in GSE59421, including 67 upregulated and 38 downregulated miRNAs (Supplementary Table S11). These DE-miRNAs were intersected with miRNAs identified from the TF–mRNA–miRNA network and visualized using a volcano plot (Figure 13A). After deduplication, 88 DE-miRNAs were obtained and further intersected with 13 miRNAs from the miRNet database, resulting in six mature miRNAs (Figure 13B). These mature miRNAs corresponded to seven miRNAs—hsa-miR-17-5p, hsa-miR-17-3p, hsa-miR-93-5p, hsa-miR-106b-5p, hsa-miR-18a-5p, hsa-miR-25-3p, and hsa-miR-146a-5p—in the TF–mRNA–miRNA network and were used to construct the ceRNA network. Among the identified regulatory axes, lncRNA NEAT1/hsa-miR-17-5p/MCL1 was highlighted (Figure 13C).

To assess the clinical relevance of the five biomarkers (MCL1, F13A1, RGS2, TLR8, and TAGAP), we analyzed associations between their expression levels and key clinical parameters of AS. Based on baseline characteristics of the GSE90074 dataset (Table 1), four parameters significantly associated with disease status were identified: gender, hyperlipidemia, disease severity (CAD_class), and the inflammatory marker CXCL5.

Spearman correlation analysis within AS patients showed that RGS2, TLR8, and TAGAP expression levels were negatively correlated with CXCL5 (p < 0.05) (Supplementary Table S12). Stratification by disease status and hyperlipidemia revealed that MCL1, F13A1, TLR8, and TAGAP were significantly upregulated in patients with both AS and hyperlipidemia compared with healthy individuals without hyperlipidemia (p < 0.05) (Supplementary Table S13). This indicated that the upregulation of these genes was closely coupled with the high-risk clinical phenotype of “disease + key risk factor”, demonstrating their potential for patient risk stratification.

In terms of disease severity, stratification analysis by coronary artery disease severity indicated a direct correspondence between gene expression levels and the degree of anatomical lesions (Supplementary Table S14). Among them, MCL1, F13A1, and TAGAP were significantly upregulated early in the disease (CAD_class 2) (p < 0.05), suggesting their value as early diagnostic markers, while RGS2 and TLR8 showed the highest expression at the most severe disease stage (CAD_class 4) (p < 0.05), indicating their association with disease progression.

From the perspective of demographic characteristics, stratification analysis by gender and disease status revealed specific patterns of gene expression (Supplementary Table S15). Several genes (e.g., F13A1, TLR8, TAGAP) showed significantly lower expression in Control Males compared to other groups (p < 0.05), while TLR8 exhibited significant baseline gender differences even in the healthy population (p < 0.05Overall, these findings indicate that the five biomarkers are associated not only with AS presence but also with inflammation, risk stratification, disease progression, and gender-specific characteristics, supporting their potential clinical utility.