The role of 5-methylcytosine regulator-related genes in diagnostic and immune regulatory functions in atherosclerosis

Background: Atherosclerosis (AS) is a chronic inflammatory disease with a poor prognosis, and 5-methylcytosine (m5C) RNA modification plays a significant role in AS-induced cardio-cerebrovascular diseases (CCVDs). However, the effects of m5C modification and related genes in AS remain unclear.

Methods: We analyzed the correlations between m5C RNA modification and its regulatory genes in AS using microarray databases. Specifically, microarray datasets of AS (GSE90074, GSE27034, GSE59421, and GSE159677) were obtained from the Gene Expression Omnibus (GEO) database. Following differential expression and Spearman correlation analyses using GSE90074, m5C regulator-related genes (MRRGs) were screened to construct a protein–protein interaction (PPI) network. A least absolute shrinkage and selection operator (LASSO) logistic regression model was constructed and validated using receiver operating characteristic (ROC) curves in GSE90074 and GSE27034 to identify feature genes. Consensus clustering analysis was then performed to classify AS samples into distinct clusters. In addition, Spearman correlation analysis was used to explore the associations between differentially expressed m5C regulators (DE-MRs) and immune cells based on the identified clusters. Weighted gene co-expression network analysis (WGCNA) was applied to identify AS-related module genes. Subsequently, intersecting genes common to module genes and differentially expressed genes (DEGs) across AS-related clusters were considered candidate biomarkers and were validated by quantitative real-time polymerase chain reaction (qRT-PCR) in human myeloid leukemia mononuclear cells (THP-1). Single-cell RNA sequencing (scRNA-seq) analysis was performed to characterize the immune microenvironment of AS. In vitro experiments and genetic interventions were conducted to investigate the effects of the m5C regulator NSUN3 on macrophage function. Finally, a competitive endogenous RNA (ceRNA) network targeting the identified biomarkers was predicted using the miRNet database.

Results: Based on two differentially expressed m5C regulators (NSUN3 and NSUN5), 546 of 2247 DEGs between AS and control samples were identified as MRRGs for PPI network construction. Twenty hub MRRGs were further incorporated into the LASSO logistic regression model, yielding nine feature genes. Based on these feature genes, AS samples were classified into two clusters, with five immune cell types showing significant differences between clusters. Both NSUN3 and NSUN5 showed the strongest correlations with M0 macrophages. A total of 643 module genes were identified and overlapped with DEGs from the two AS-related clusters, resulting in five biomarkers—MCL1, F13A1, RGS2, Toll-like receptor 8 (TLR8), and TAGAP. The expression patterns of these five biomarkers in the foam cell model were consistent with those observed in public datasets. Furthermore, NSUN3 regulated the production of proinflammatory cytokines in macrophages. Finally, a ceRNA regulatory network was constructed.

Conclusion: Five potential diagnostic biomarkers for AS—MCL1, F13A1, RGS2, TLR8, and TAGAP—were identified. In addition, the m5C regulator NSUN3 plays a critical role in macrophage function, providing experimental evidence that may support the diagnosis and treatment of AS.

1 Introduction

Cardiac-cerebral vascular diseases (CCVDs) are characterized by high morbidity and mortality and are responsible for nearly 20 million deaths annually worldwide (1). CCVDs, such as myocardial infarction and stroke, are primarily caused by atherosclerosis (AS), a complex chronic inflammatory vascular disease involving vascular endothelial dysfunction, inflammatory cell infiltration, oxidative stress, etc. Despite ongoing research, the precise pathogenesis of AS remains unclear. Thus, identifying reliable diagnostic biomarkers for AS and elucidating their associated modulatory networks are critical for early detection, assessment of disease severity, and ultimately, improved prognosis through stratified intervention.

Recent epigenetic studies have drawn attention to 5-methylcytosine (m5C) as a major focus in RNA methylation research. m5C is a ubiquitous modification in both coding and non-coding RNA. Messenger RNA (mRNA) sequencing in mice and human tissues has revealed a significant increase in m5C expression in heart and brain tissues compared with other organs. Gene ontology (GO) enrichment analysis has shown that m5C methylation is highly enriched in mRNA transcripts of mitochondrial-related genes (2). Due to high energy consumption and mitochondrial dependence, m5C methylation plays a crucial role in CCVDs by regulating mitochondrial functional balance.

The target RNA molecules modified by m5C methylation mainly include messenger RNA (mRNA), transfer RNA (tRNA), and ribosomal RNA (rRNA), resulting in diverse effects, such as influencing protein translation and participating in physiological processes associated with various diseases. Moreover, m5C modification sites are often embedded in CG-rich contexts and coding sequence (CDS) regions, thereby regulating diverse biological effects in AS. m5C modification is governed by m5C regulators, including “writers” (methyltransferases), “erasers” (demethylases), and “readers.” Among these, m5C writers mainly include NSUN1–7 and DNMT1–3; erasers include ALKBH1 and the TET family; and ALYREF and YBX1 serve as m5C readers that bind to m5C motifs. Despite some progress, the pathophysiological mechanisms underlying m5C regulators in CCVDs remain poorly understood. Intercellular adhesion molecule 1 (ICAM-1), a crucial factor in upregulating endothelial inflammation and vascular sclerosis, has been reported to be regulated by NSUN2-mediated m5C mRNA methylation (3). In addition, DNMT2, another critical m5C methyltransferase, participates in neuronal growth and development (4). Furthermore, Wang et al. demonstrated that parental zebrafish exposed to brominated diphenyl ethers (BDP) exhibited altered m5C DNA methylation patterns in key m6A regulatory components (methyltransferases, demethylases, and recognition proteins) within gonadal tissues. Notably, these epigenetic modifications were transgenerationally transmitted to offspring through germline inheritance, potentially disrupting vascular developmental processes. Mechanistically, BDP-induced methylation reprogramming impaired vascular morphogenesis by dysregulating epigenetic control of angiogenic signaling pathways, particularly through aberrant spatiotemporal expression of vascular endothelial growth factor (VEGF)-related genes (5, 6). However, the specific mechanisms of m5C regulators and related pathways in AS remain unclear.

In this study, we aimed to identify potential biomarkers for AS by analyzing m5C regulator-related genes (MRRGs). Based on feature genes obtained using a LASSO logistic model, AS samples were classified into distinct clusters. Through overlap analysis of AS-related module genes and differentially expressed genes (DEGs) across clusters, diagnostic biomarkers for AS were further identified. Moreover, the expression of these biomarkers was validated in human myeloid leukemia mononuclear cells (THP-1) using quantitative real-time polymerase chain reaction (qRT-PCR). The immune microenvironment characteristics of AS were explored using single-cell RNA sequencing (scRNA-seq) analysis. In addition, the role of the m5C regulator NSUN3 was investigated in an in vitro foam cell model. Finally, a competitive endogenous RNA (ceRNA) network of biomarkers was constructed using the miRNet database, which may provide insights into underlying mechanisms and serve as potential therapeutic targets for AS.

2 Materials and methods

2.1 Data collection

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).

Table 1

Table 1. Clinical characteristics of the samples in GSE90074 (** p < 0.01, *** p < 0.001).

2.2 Analysis of differentially expressed m5C regulators and related genes

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.

2.3 Construction of a LASSO logistic model

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.

2.4 Consensus clustering analysis and co-expression analysis

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.

2.5 Screening and immune analysis of biomarkers of AS

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.

2.6 Cell cultures and foam cell model construction

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.

2.7 qRT-PCR

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.

2.8 scRNA-seq analysis

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).

2.9 Virus package and cell transfection

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.

2.10 Flow cytometry

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.

2.11 Cytokine array

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.

2.12 Construction of a ceRNA network

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.

2.13 Statistical analysis

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.

3 Results

3.1 –2247 DEGs and two DE-MRs between AS and healthy control samples

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).

Figure 1

Figure 1. Graphical abstract of the current study.

Figure 2

Figure 2. Identification of DEGs and DE-MRs between AS and healthy control samples in the datasets. (A) Flowchart of bioinformatics analyses in GSE90074 database. (B) Volcano plot showed DEGs between AS and control samples. (C) Expression of DE-MRs between AS and control samples.

3.2 –546 MRRGs were obtained via correlation analysis

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.

Figure 3

Figure 3. dentification of MRRGs via correlation analysis. (A) DEGs between AS and control groups showing strong correlation with DE-MRs (|Cor| > 0.3). GO enrichment(B) analysis of MRRGs. (C) KEGG pathway analysis of MRRGs. (D) PPI network based on 20 hub MRRGs and their correlated genes. (E) Venn diagram of 20 AS-related hub genes. (F) Bar plot showing inference scores between the 20 hub MRRGs and AS incidence risk.

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.

3.3 Nine feature genes in the LASSO logistic model of 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).

Figure 4

Figure 4. Feature genes in the logistic regression model of AS. (A) Feature genes selected by the univariate Cox regression analysis. (B) Cross-validation for tuning parameter selection in the LASSO model. (C) Feature genes selected by the LASSO model. (D) Predictive performance of the model in GSE90074 (ROC curves). (E) Predictive performance of the model in GSE27034 (ROC curves).

3.4 Cluster analysis of AS samples and identification of gene modules with WGCNA

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.

Figure 5

Figure 5. Consensus clustering analysis. (A) Consensus clustering cumulative distribution function (CDF) when data were divided into 2–9 clusters. (B) AS samples divided into two clusters based on the expression of nine feature genes (k = 2). (C) PCA validation of cluster separation. (D) tSNE validation of cluster separation. (E) Module clustering dendrogram derived from the TOM. (F) Module–cluster relationships in AS.

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).

Table 2

Table 2. Clinical characteristics of the samples in different sub-clusters (* 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.

Figure 6

Figure 6. Analysis of DEGs between the different clusters. (A) Heatmap of 541 DEGs between cluster 1 and cluster 2 identified by differential analysis. (B, C) GO and KEGG enrichment analyses of the 541 DEGs.

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.

3.5 Immune analysis of AS clusters

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).

Figure 7

Figure 7. Immune analysis of different AS clusters. (A) Differences in immune cell subtypes between cluster 1 and cluster 2. (B) Correlation between NSUN3 and macrophages M0. (C) Correlation between NSUN5 and macrophages M0 (* p < 0.05, ** p < 0.01).

3.6 Identification and validation of five biomarkers

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).

Figure 8

Figure 8. Identification and validation of five biomarkers in AS. (A) Venn diagram of 81 candidate biomarkers. (B, C) Expression trends of five biomarkers in the training and validation sets. (D) Construction of the foam cell model. (E) Verification of biomarkers expression in cell lines by qRT-PCR (*p < 0.05, ****p < 0.0001).

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).

3.7 scRNA-seq analysis revealed the immune microenvironment of AS

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).

Figure 9

Figure 9. scRNA-seq analysis revealing the immune microenvironment of AS. (A) Distributions of cells across different samples visualized by UMAP dimensionality reduction. (B) Composition of different cell types in AS and control groups. (C) Expression level of biomarkers in different cell subtypes. (D) UMAP plots showing expression of two DE-MRs across cell subtypes in AS samples. (E) Expression levels of NSUN3 and NSUN5 across distinct cellular subpopulations between AS and control groups (* p < 0.05, ** p < 0.01, *** p < 0.001, **** p < 0.0001).

3.8 The expression levels of NSUN3 and NSUN5 in the foam cell model

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).

Figure 10

Figure 10. The expression levels of NSUN3 and NSUN5 in the foam cell model. (A) MFI histograms for NSUN3 and NSUN5 in the foam cell model. (B) The expression levels of NSUN3 and NSUN5 in the foam cell model. n=4/group (** p < 0.01).

3.9 The immune and inflammatory effects of NSUN3 in macrophages during AS progression

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).

Figure 11

Figure 11. NSUN3 transfection inhibits apoptosis and promotes M1 to M2 polarization in the AS model. (A) Schematic illustrating the experimental procedure for THP-1 cell treatment. THP-1 cells were infected with control virus or NSUN3 virus for 3 days, followed by PMA-induced differentiation into M0 macrophages for 48 h and ox-LDL treatment for 24 h to generate foam cells. Cells were then harvested for flow cytometry analysis. (B) NSUN3 expression efficiency in ox-LDL–treated THP-1 cells under the indicated conditions (n = 3/group). (C, D) Annexin V/PI apoptosis analysis and CD86^-/CD206^- macrophage polarization in NSUN3-transfected cells and controls. n=3/group (*p < 0.05, **p < 0.01, ****p < 0.0001).

Figure 12

Figure 12. Cytokine profiling of the NSUN3 transfected cells and controls. (A) Heatmap showing relative pixel intensity of cytokine spots in the NSUN3-transfected model and controls. (B) Bar graphs showing the most significantly dysregulated cytokines (**p < 0.01, ***p < 0.001).

3.10 Construction of the ceRNA network of biomarkers

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).

Figure 13

Figure 13. Construction of the ceRNA network of biomarkers. (A) Volcano plot of DE-miRNAs between AS and control samples in GSE59421 and those identified in the TF–mRNA–miRNA network. (B) Venn diagram showing six mature miRNAs.

3.11 Clinical relevance of the five biomarkers

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.

4 Discussion

To date, over 170 known RNA chemical modifications have been identified, with methylation modifications accounting for more than 60% of these, including 6-methyladenosine (m6A) and 5-methylcytosine (m5C) (20, 21). A growing body of research has highlighted the crucial role of methylation in AS and its associated immune regulation (22). The m6A “writer” methyltransferase METTL3 has been reported as a key hub in hemodynamic forces. Oscillatory stress in AS stimulated the upregulation of METTL3 expression, which led to attenuated RelA/p65 phosphorylation, subsequently augmenting nuclear factor kappa B (NF-κB) activity, exacerbating inflammatory responses, and facilitating plaque formation. Increased METTL3 abundance also regulated the levels of nucleotide-binding domain leucine-rich repeat protein 1 (NLRP1) and Kruppel-like factor 4 (KLF4), inducing endothelial inflammation and arterial expansion (23). The RNA-binding protein Matrin-3 exerted anti-inflammatory effects in macrophages by orchestrating the formation of the METTL3–METTL4 methyltransferase complex to promote m6A-dependent mRNA degradation, thereby attenuating ox-LDL-triggered mitogen-activated protein kinase (MAPK) signaling activation and subsequent inflammatory responses (24). Luo et al. showed that NSUN2 methylated ICAM-1 mRNA and increased leukocyte adhesion to vascular endothelial cells, while ICAM-1 expression was significantly diminished in NSUN2^-/^- rats. Mechanistically, TNF-α or homocysteine increased NSUN2 methyltransferase activity by inhibiting NSUN2 phosphorylation. Moreover, donor NSUN2 deficiency hindered the development of allograft arteriosclerosis in vivo (3). In addition, NSUN2 was reported to regulate ALYREF nuclear–cytoplasmic trafficking efficiency, directly influencing cytoplasmic accumulation of m5C-modified mRNAs involved in leukocyte recruitment and vascular smooth muscle cell apoptosis. Mechanistically, NSUN2-mediated m5C modification stabilized ALYREF–mRNA complexes, amplifying inflammatory signaling cascades through enhanced translation fidelity (25). However, diagnostic biomarkers and immune regulatory functions of m5C regulators in AS remain unclear. In this study, we screened AS-related genes significantly correlated with m5C regulators, identified potential biomarkers, constructed a ceRNA network, and conducted a series of experiments to analyze immune and inflammatory effects of m5C regulators on macrophage function. These findings may provide novel insights into therapeutic strategies for AS.

Previously, Nakano et al. reported that knockout of the RNA methyltransferase NSUN3 resulted in marked reductions in mitochondrial protein synthesis and oxygen consumption, leading to impaired mitochondrial activity (26). They demonstrated that the biogenesis pathway of the f5C34 modification in mitochondrial transfer RNA^met (mt-tRNA ^met) begins with NSUN3-mediated catalytic activity. Quantitative analysis of the human mitochondrial proteome revealed a striking codon adaptation: among the 13 mtDNA-encoded polypeptides (containing a total of 207 methionine incorporation sites), AUA codons (167 loci) were 4.2-fold more prevalent than canonical AUG codons (40 loci). This evolutionary adaptation established f5C-dependent AUA decoding as a critical quality control checkpoint. NSUN3 depletion caused failure of wobble base-pair resolution at AUA codons, resulting in translation stalling and collapse of mitochondrial proteome integrity. Mitochondrial dysfunction is recognized as an essential factor in AS initiation and progression. For example, serine/arginine-rich splicing factor 3 (SRSF3) maintained effective clearance of oxidized lipids and apoptotic debris in AS lesions by regulating mitochondrial proteostasis (27–29). Phenotypic consequences of NSUN5 deficiency in mammalian cells include reduced cell proliferation and size, highlighting the importance of m5C regulators in ribosome function and normal cellular physiology (30). In the present study, NSUN3 and NSUN5 were identified as two DE-MRs between AS and healthy control samples. Hub genes in the constructed PPI network included multiple ribosomal proteins (e.g., RPL12, RPL15, RPL18A). As core sensors of cellular stress, alterations in ribosomal function can integrate multiple pathological processes in AS, including endoplasmic reticulum stress, metabolic dysregulation, and inflammatory signaling (31). These findings extend the role of m5C modifications beyond protein synthesis, suggesting broader involvement in the AS pathological environment through regulation of ribosomal function (32).

Atherosclerosis is a chronic inflammatory disease of the arterial vessel wall, with prominent involvement of immune cells, including monocytes, macrophages, dendritic cells (DCs), T cells, and B cells (18, 33, 34). AS is initiated by endothelial dysfunction, which promotes secretion of chemotactic mediators such as CCL2 and chemokine (C-X-C motif) ligand 1 (CXCL1), recruiting circulating myeloid cells, particularly monocytes and neutrophils. Monocytes infiltrate lesions and differentiate into macrophages, which uptake lipoproteins and form foam cells that constitute the central core of atheromas. Neutrophils act as early responders, releasing reactive oxygen species via NADPH oxidase and myeloperoxidase, modifying leukocyte adhesion and promoting endothelial dysfunction. In addition, neutrophils secrete granule proteins that act as extracellular inflammatory mediators, with growing evidence supporting a role for neutrophil extracellular traps (NETs) (35–37). DCs characterized by the CD11c^- human leukocyte antigen-DR (HLA-DR^-) phenotype function as professional antigen-presenting cells via major histocompatibility complex class II (MHC-II)/T cell receptor (TCR) interactions, driving naive T cell activation. In AS, DC subsets (e.g., CD103^- vs. CD8α^-) differentially regulate disease progression through scavenger receptor–mediated uptake of oxidized LDL, antigen presentation, and interleukin (IL)-23–driven T helper 17 (Th17) polarization (38, 39). Mast cells, predominantly localized in the arterial adventitia, promote extracellular matrix degradation through proteolytic enzymes targeting elastic fibers and laminae, thereby contributing to vascular injury and maladaptive remodeling (40). These observations are consistent with our enrichment and immune analyses of MRRGs and AS clusters.

Our CTD analysis revealed that hub MRRGs, including TLR4 and TP53, exhibited higher inference scores. The TLR family comprises pattern-recognition receptors located on the cell membrane. Preclinical studies have shown that inhibition of the TLR4/MyD88/NF-κB axis attenuates neointimal hyperplasia and plaque destabilization (41). Mechanistically, ox-LDL induces TLR4/NF-κB–dependent transcription of matrix metalloproteinase 9 (MMP-9), promoting extracellular matrix degradation and fibrous cap thinning. TLR4 also regulates macrophage polarization; lipopolysaccharide (LPS)–TLR4 interaction promotes M1 polarization and increases proinflammatory cytokines, including TNF-α, IL-1β, IL-6, and monocyte chemoattractant protein 1 (MCP-1) (42). TP53 is a stress-responsive transcription factor activated following genomic damage. TP53-inducible enzymes regulate cellular redox homeostasis and protect against apoptosis by scavenging reactive oxygen species (ROS). These antioxidant functions mitigate oxidative stress–induced DNA damage and genomic instability, promoting cell survival under stress. Increasing evidence suggests that TP53-mediated prosurvival mechanisms contribute to AS pathophysiology through modulation of oxidative stress responses in vascular cells (43, 44).

We constructed a diagnostic prediction model for AS using the LASSO regression algorithm and analyzed regulatory mechanisms of feature genes through a TF–mRNA–miRNA network. PRDM1, regulated by five miRNAs (hsa-miR-223-3p, hsa-miR-181a-5p, hsa-miR-20a-3p, hsa-miR-99a-5p, and hsa-miR-20a-5p) and three TFs (CEBPB, PAX5, and STAT3), may have therapeutic relevance in AS. These findings are consistent with the result reported by Folgado et al., who demonstrated that PRDM1fl/fl Aicda-Cre⁺/ki Ldlr^−/− chimeras increased AS via germinal center-derived plasma cells and/or antibodies (45).

Based on feature gene analysis, AS samples were divided into clusters, with significant differences observed in five immune cell types: activated DCs, macrophages M0, macrophages M2, resting mast cells, and neutrophils. NSUN3 showed the strongest correlation with macrophages M0, which was further validated by flow cytometry. We found that the relative level of NSUN3 was significantly increased in PBMC-derived foam cells compared with the control cells. The apoptosis cells were significantly reduced in NSUN3-transfected macrophages. And the NSUN3 infection could stimulate the transition of foam cells from M1 to M2. Macrophages M0 are the undifferentiated cell type that can be guided to polarized cell types, according to corresponding signals. The switch between polarized macrophages, for instance, from M2 to M1, determined the transformation of the inflammatory microenvironment in the progression of AS. The classically activated M1 subset orchestrated inflammatory pathogenesis and mediated tissue degradation cascades; in contrast, the alternatively activated M2 compartment executed reparative functions, modulated plaque biomechanical stability and repressed pro-inflammatory signaling networks. The knockdown of the m5C regulator NSUN3 has been shown to increase the infiltration of M1 macrophages and decrease the level of M2 macrophages in xenograft models (46). However, the role of NSUN5 in macrophage function remains poorly understood and warrants further investigation.We identified five biomarkers, namely MCL1, F13A1, RGS2, TLR8, and TAGAP, for AS and investigated their expression levels in AS cell lines. Fontaine et al. previously reported that myeloid deletion of MCL-1 promoted apoptosis and lipid accumulation in AS plaques, implicating MCL-1 in the survival and differentiation of leukocytes, particularly neutrophils (47). The coagulation factor XIII A gene (F13A1) is a recognized M2 anti-inflammatory macrophage marker and can be regulated by neuron-derived orphan receptor 1, which is expressed in human AS lesions (48). Genetic polymorphisms in RGS2 have been associated with intima–media thickening of the carotid artery in both hypertensive and general populations, indicating a potential role in AS pathogenesis (49). In addition, Allen et al. demonstrated that TLR8 activation led to enrichment of microbial small RNA (msRNA) on LDL and promoted pro-inflammatory macrophage polarization (50), highlighting the role of TLR8 in AS. However, the potential effects of TAGAP in AS have not yet been fully investigated.

This study established the clinical relevance of five potential biomarkers in AS through systematic correlation analyses with well-established clinical parameters. Our findings demonstrated that these genes were significantly associated with key aspects of AS pathology, including inflammatory mechanisms, hyperlipidemia interaction, disease severity progression, and gender-specific effects. The biomarkers showed distinct clinical utility patterns: MCL1, F13A1, and TAGAP served as early-stage markers, whereas RGS2 and TLR8 were associated with advanced disease, providing complementary information for disease staging. Furthermore, the gender-specific expression patterns, particularly for TLR8, highlighted the importance of considering sex differences in AS biomarker research. Rather than replacing existing diagnostic methods, these transcriptomic biomarkers provide added value by offering molecular-level insights into disease mechanisms, enabling risk stratification, and facilitating early detection.

Chang et al. reported that genetic silencing of MCL-1 augmented aspirin-induced viability loss and apoptosis in glioma cells (51). Higher expression of the hematopoietic transcription factor RUNX1 target F13A1 was associated with acute events in patients with cardiovascular disease receiving aspirin or ticagrelor (52). These findings are consistent with our results showing that aspirin was correlated with both MCL-1 and F13A1. Moreover, we constructed a ceRNA network by identifying six mature miRNAs corresponding to seven miRNAs in the TF–mRNA–miRNA network. We found that lncRNA NEAT1 could regulate MCL-1 expression by targeting hsa-miR-17-5p. Previous studies have shown that NEAT1 regulates the protein kinase B (AKT)/mammalian target of rapamycin (mTOR) signaling pathway in ox-LDL–induced human aortic endothelial cells, accelerating proliferation and inhibiting apoptosis, suggesting that NEAT1 may be a potential therapeutic target in AS (53). Furthermore, reduced levels of miR-17-5p and enhanced Beclin-1–mediated autophagy increased cholesterol efflux from THP-1 macrophage-derived foam cells (54). Downregulation of NEAT1 expression has also been shown to inhibit proliferation and induce apoptosis in rheumatoid arthritis fibroblast-like synoviocytes by targeting miR-17-5p and inhibiting STAT3 (55). However, the mechanisms underlying the lncRNA NEAT1/miR-17-5p/MCL-1 axis in AS require further investigation.

Several limitations of this study should be acknowledged. First, the number of samples available in public databases for both AS and healthy controls was relatively small, which may introduce bias and limit the generalizability of the findings. Prospective validation in longitudinal cohorts with documented cardiovascular outcomes will be essential to establish prognostic value. Second, additional functional studies are required to substantiate the diagnostic and therapeutic relevance of the identified biomarkers and their associated biological pathways, which may guide the development of targeted therapies. Third, the molecular mechanisms of the identified ceRNA networks, such as the lncRNA NEAT1/miR-17-5p/MCL-1 axis, in AS progression remain incompletely understood. Thus, more comprehensive studies are needed to clarify these limitations and provide a more accurate understanding of the molecular mechanisms underlying AS.

5 Conclusions

In conclusion, our findings revealed significant differences in the expression of m5C regulators between AS and healthy control samples, which were used to identify MRRGs and construct a diagnostic model for AS. We further identified five potential diagnostic biomarkers—MCL1, F13A1, RGS2, TLR8, and TAGAP. In addition, our study demonstrated the immune and inflammatory effects of NSUN3 in AS, providing new evidence that may inform the development of novel immunotherapeutic strategies in clinical practice.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Ethics statement

Ethical approval was not required for the studies on humans in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used.

Author contributions

HZ: Writing – review & editing, Writing – original draft, Funding acquisition, Conceptualization. YK: Formal Analysis, Writing – original draft, Data curation, Methodology, Investigation, Validation. PD: Writing – original draft, Data curation. LY: Formal Analysis, Writing – original draft. NL: Conceptualization, Writing – review & editing, Writing – original draft, Funding acquisition.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This research was funded by Shanxi Province Basic Research Program (202203021222353), Shanxi Bethune Hospital Talent Introduction Research Fund Project (2022RC09).

Conflict of interest

The authors declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2025.1636323/full#supplementary-material

References

1. Roth GA, Mensah GA, Johnson CO, Addolorato G, Ammirati E, Baddour LM, et al. Global burden of cardiovascular diseases and risk factors, 1990-2019: Update From the GBD 2019 Study. J Am Coll Cardiol. (2020) 76:2982–3021. doi: 10.1016/j.jacc.2020.11.010

PubMed Abstract | Crossref Full Text | Google Scholar

2. Huang T, Chen W, Liu J, Gu N, and Zhang R. Genome-wide identification of mRNA 5-methylcytosine in mammals. Nat Struct Mol Biol. (2019) 26:380–8. doi: 10.1038/s41594-019-0218-x

PubMed Abstract | Crossref Full Text | Google Scholar

3. Luo Y, Feng J, Xu Q, Wang W, and Wang X. NSun2 deficiency protects endothelium from iInflammation via mRNA methylation of ICAM-1. Circ Res. (2016) 118:944–56. doi: 10.1161/circresaha.115.307674

PubMed Abstract | Crossref Full Text | Google Scholar

4. Rai K, Chidester S, Zavala CV, Manos EJ, James SR, Karpf AR, et al. Dnmt2 functions in the cytoplasm to promote liver, brain, and retina development in zebrafish. Genes Dev. (2007) 21:261–6. doi: 10.1101/gad.1472907

PubMed Abstract | Crossref Full Text | Google Scholar

5. Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature. (2012) 485:201–6. doi: 10.1038/nature11112

PubMed Abstract | Crossref Full Text | Google Scholar

6. Wang C, Xu J, Luo S, Huang J, Ji D, Qiu X, et al. Parental exposure to environmentally relevant concentrations of bisphenol-A bis(diphenyl phosphate) impairs vascular development in offspring through DNA/RNA methylation-dependent transmission. Environ Sci Technol. (2023) 57:16176–89. doi: 10.1021/acs.est.3c03579

PubMed Abstract | Crossref Full Text | Google Scholar

7. Chen YS, Yang WL, Zhao YL, and Yang YG. Dynamic transcriptomic M(5) C and its regulatory role in RNA processing. Wiley Interdiscip Rev RNA. (2021) 12:e1639. doi: 10.1002/wrna.1639

PubMed Abstract | Crossref Full Text | Google Scholar

8. Aanniz T, Bouyahya A, Balahbib A, El Kadri K, Khalid A, Makeen HA, et al. Natural bioactive compounds targeting DNA methyltransferase enzymes in cancer: Mechanisms insights and efficiencies. Chem Biol Interact. (2024) 392:110907. doi: 10.1016/j.cbi.2024.110907

PubMed Abstract | Crossref Full Text | Google Scholar

9. Nombela P, Miguel-Lopez B, and Blanco S. The role of m(6)A, m(5)C and ΨRNA modififications in cancer: Novel therapeutic opportunities. Mol Cancer. (2021) 20:18. doi: 10.1186/s12943-020-01263-w

PubMed Abstract | Crossref Full Text | Google Scholar

10. Dai Q, Ye C, Irkliyenko I, Wang Y, Sun HL, Gao Y, et al. Ultrafast bisulfite sequencing detection of 5-methylcytosine in DNA and RNA. Nat Biotechnol. (2024) 42:1559–70. doi: 10.1038/s41587-023-02034-w

PubMed Abstract | Crossref Full Text | Google Scholar

11. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | Crossref Full Text | Google Scholar

12. Hänzelmann S, Castelo R, and Guinney J. GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinformatic. (2013) 14:7. doi: 10.1186/1471-2105-14-7

PubMed Abstract | Crossref Full Text | Google Scholar

13. Yu G, Wang LG, Han Y, and He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | Crossref Full Text | Google Scholar

14. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303

PubMed Abstract | Crossref Full Text | Google Scholar

15. Friedman J, Hastie T, and Tibshirani R. Regularization Ppaths for generalized linear models via coordinate descent. J Stat Softw. (2010) 33:1–22. doi: 10.18637/jss.v033.i01

PubMed Abstract | Crossref Full Text | Google Scholar

16. Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez JC, et al. pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinf. (2011) 12:77. doi: 10.1186/1471-2105-12-77

PubMed Abstract | Crossref Full Text | Google Scholar

17. Wilkerson MD and Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics. (2010) 26:1572–3. doi: 10.1093/bioinformatics/btq170

PubMed Abstract | Crossref Full Text | Google Scholar

18. Drechsler M, Megens RT, van Zandvoort M, Weber C, and Soehnlein O. Hyperlipidemia-triggered neutrophilia promotes early atherosclerosis. Circulation. (2010) 122:1837–45. doi: 10.1161/circulationaha.110.961714

PubMed Abstract | Crossref Full Text | Google Scholar

19. Aran D, Hu Z, and Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. (2017) 18:220. doi: 10.1186/s13059-017-1349-1

PubMed Abstract | Crossref Full Text | Google Scholar

20. Xu J, Liu Y, Liu J, Xu T, Cheng G, Shou Y, et al. The identification of critical m(6)A RNA methylation regulators as Malignant prognosis factors in prostate sdenocarcinoma. Front Genet. (2020) 11:602485. doi: 10.3389/fgene.2020.602485

PubMed Abstract | Crossref Full Text | Google Scholar

21. Song C, Qiao Z, Cheng J, Zhang Y, Liu W, Yuan S, et al. Optimal definition and risk stratification in prediabetes with stab le coronary heart disease: A Prospective Cohort Study. J Am Heart Assoc. (2025) 14:e037492. doi: 10.1161/JAHA.124.037492

PubMed Abstract | Crossref Full Text | Google Scholar

22. Liu S, Cao Y, and Zhang Y. Regulatory roles of RNA methylation in vascular lesions in ocular and cardiopulmonary diseases. Crit Rev Clin Lab Sci. (2024) 61:726–40. doi: 10.1080/10408363.2024.2370267

PubMed Abstract | Crossref Full Text | Google Scholar

23. Chien CS, Li JY, Chien Y, Wang ML, Yarmishyn AA, Tsai PH, et al. METTL3-dependent N6-methyladenosine RNA modification mediates the atherogenic inflammatory cascades in vascular endothelium. Proc Natl Acad Sci USA. (2021) 118:e2025070118. doi: 10.1073/pnas.2025070118

PubMed Abstract | Crossref Full Text | Google Scholar

24. Sun Z, Chen W, Wang Z, Wang S, Zan J, Zheng L, et al. Matr3 reshapes m6A modification complex to alleviate macrophage inflammation during atherosclerosis. Clin Immunol. (2022) 245:109176. doi: 10.1016/j.clim.2022.109176

PubMed Abstract | Crossref Full Text | Google Scholar

25. He Y, Zhang H, Yin F, Guo P, Wang S, Wu Y, et al. Novel insights into the role of 5-Methylcytosine RNA methylation in human abdominal aortic aneurysm. Front Biosci (Landmark Ed). (2021) 26:1147–65. doi: 10.52586/5016

PubMed Abstract | Crossref Full Text | Google Scholar

26. Nakano S, Suzuki T, Kawarada L, Iwata H, Asano K, and Suzuki T. NSUN3 methylase initiates 5-formylcytidine biogenesis in human mitochondrial tRNA(Met). Nat Chem Biol. (2016) 12:546–51. doi: 10.1038/nchembio.2099

PubMed Abstract | Crossref Full Text | Google Scholar

27. Zheng L, Chen X, He X, Wei H, Li X, Tan Y, et al. METTL4-mediated mitochondrial DNA N6-methyldeoxyadenosine promoting macrophage inflammation and atherosclerosis. Circulation. (2025) 15:946–65. doi: 10.1161/CIRCULATIONAHA.124.069574

PubMed Abstract | Crossref Full Text | Google Scholar

28. Yang X, Zhang X, Tian Y, Yang J, Jia Y, Xie Y, et al. Srsf3-dependent APA drives macrophage maturation and limits atherosclerosis. Circ Res. (2025) 136:985–1009. doi: 10.1161/CIRCRESAHA.124.326111

PubMed Abstract | Crossref Full Text | Google Scholar

29. Kim K, Thome T, Pass C, Stone L, Vugman N, Palzkill V, et al. Multiomic analysis of calf muscle in peripheral artery disease and chronic kidney disease. Circ Res. (2025) 136:688–703. doi: 10.1161/CIRCRESAHA.124.325642

PubMed Abstract | Crossref Full Text | Google Scholar

30. Heissenberge C, Liendl L, Nagelreiter F, Gonskikh Y, Yang G, EM S, et al. Loss of the ribosomal RNA methyltransferase NSUN5 impairs global protein synthesis and normal growth. Nucleic Acids Res. (2019) 47:11807–25. doi: 10.1093/nar/gkz1043

PubMed Abstract | Crossref Full Text | Google Scholar

31. Song H, Zhang J, Liu B, Xu J, Cai B, Yang H, et al. Biological roles of RNA m5C modification and its implications in Cancer immunotherapy. biomark Res. (2022) 10:15. doi: 10.1186/s40364-022-00362-8

PubMed Abstract | Crossref Full Text | Google Scholar

32. Lin G, Cai H, Hong Y, Yao M, Ye W, Li W, et al. Implications of m5C modifications in ribosomal proteins on oxidative stress, metabolic reprogramming, and immune responses in patients with mid-to-late-stage head and neck squamous cell carcinoma: Insights from nanopore sequencing. Heliyon. (2024) 10:e34529. doi: 10.1016/j.heliyon.2024.e34529

PubMed Abstract | Crossref Full Text | Google Scholar

33. Li C, Zhu L, Dai Y, Zhang Z, Huang L, Wang TJ, et al. Diet-induced high serum levels of Trimethylamine-N-oxide enhance the cellular inflammatory response without exacerbating acute intracerebral hemorrhage injury in mice. Oxid Med Cell Longev. (2022) 16:1599747. doi: 10.1155/2022/1599747

PubMed Abstract | Crossref Full Text | Google Scholar

34. Choi YY, Kim A, and Seong KM. Chronic radiation exposure aggravates atherosclerosis by stimulating neutrophil infiltration. Int J Radiat Biol. (2021) 97:1270–81. doi: 10.1080/09553002.2021.1934750

PubMed Abstract | Crossref Full Text | Google Scholar

35. Ma X, Zhao X, Yang Y, Yan J, Shi X, Wu H, et al. Paeonol inhibits NETs-mediated foam cell inflammation through the CitH3/NLRP3/caspase-1 signaling pathway in atherosclerosis. Int Immunopharmacol. (2025) 151:114340. doi: 10.1016/j.intimp.2025.114340

PubMed Abstract | Crossref Full Text | Google Scholar

36. Wu M, Chen M, Zhao Y, Zhang X, Ding X, Yuan J, et al. Neutrophil hitchhiking-mediated delivery of ros-scavenging biomimetic nanoparticles for enhanced treatment of atherosclerosis. Small Methods. (2025) 4:e2402019. doi: 10.1002/smtd.202402019

PubMed Abstract | Crossref Full Text | Google Scholar

37. Ma Y, Cao H, Chen B, Xu X, Zhang Q, Chen H, et al. Simultaneous in vivo imaging of neutrophil elastase and oxidative stress in atherosclerotic plaques using a unimolecular photoacoustic probe. Angew Chem Int Ed Engl. (2024) 63:e202411840. doi: 10.1002/anie.202411840

PubMed Abstract | Crossref Full Text | Google Scholar

38. Zhang L, Al-Ammari A, Zhu D, Zhang H, Zhou P, Zhi X, et al. A nanovaccine for immune activation and prophylactic protection of atherosclerosis in mouse models. Nat Commun. (2025) 16:2111. doi: 10.1038/s41467-025-57467-5

PubMed Abstract | Crossref Full Text | Google Scholar

39. Clement M, Raffort J, Lareyre F, Tsiantoulas D, Newland S, Lu Y, et al. Impaired autophagy in CD11b⁺ dendritic cells expands CD4⁺ regulatory T cells and limits atherosclerosis in mice. Circ Res. (2019) 25:1019–34. doi: 10.1161/CIRCRESAHA.119.315248

PubMed Abstract | Crossref Full Text | Google Scholar

40. Wang M, McGraw KR, Monticone RE, and Pintus G. Unraveling elastic fiber-derived signaling in arterial aging and related arterial diseases. Biomolecules. (2025) 15:153. doi: 10.3390/biom15020153

PubMed Abstract | Crossref Full Text | Google Scholar

41. Wang X, Chen L, Teng Y, Xie W, Huang L, Wu J, et al. Effect of three oral pathogens on the TMA-TMAO metabolic pathway. Front Cell Infect Microbiol. (2024) 14:1413787. doi: 10.3389/fcimb.2024.1413787

PubMed Abstract | Crossref Full Text | Google Scholar

42. Liu W, Wang J, Yang H, Li C, Lan W, Chen T, et al. The metabolite indole-3-acetic acid of bacteroides ovatus improves atherosclerosis by restoring the polarisation balance of M1/M2 macrophages and inhibiting inflammation. Adv Sci (Weinh). (2025), 12e2413010. doi: 10.1002/advs.202413010

PubMed Abstract | Crossref Full Text | Google Scholar

43. Wu X, Wei D, Zhou Y, Cao Q, Han G, Han E, et al. Pesticide exposures and 10-year atherosclerotic cardiovascular disease risk: Integrated epidemiological and bioinformatics analysis. J Hazard Mater. (2025) 485:136835. doi: 10.1016/j.jhazmat.2024.136835

PubMed Abstract | Crossref Full Text | Google Scholar

44. Nagareddy PR and Murphy AJ. TP53 clonal hematopoiesis promotes atherosclerosis. Nat Cardiovasc Res. (2023) 2:102–3. doi: 10.1038/s44161-022-00212-8

PubMed Abstract | Crossref Full Text | Google Scholar

45. Martos-Folgado I, Del Monte-Monge A, Lorenzo C, Busse CE, Delgado P, Mur SM, et al. MDA-LDL vaccination induces athero-protective germinal-center-derived antibody responses. Cell Rep. (2022) 41:111468. doi: 10.1016/j.celrep.2022.111468

PubMed Abstract | Crossref Full Text | Google Scholar

46. Jin S, Li J, Shen Y, Wu Y, Zhang Z, and Ma H. RNA 5-Methylcytosine regulator NSUN3 promotes tumor progression through regulating immune infiltration in head and neck squamous cell carcinoma. Oral Dis. (2024) 30:313–28. doi: 10.1111/odi.14357

PubMed Abstract | Crossref Full Text | Google Scholar

47. Fontaine MAC, Westra MM, Bot I, Jin H, Franssen A, Bot M, et al. Low human and murine Mcl-1 expression leads to a pro-apoptotic plaque phenotype enriched in giant-cells. Sci Rep. (2019) 9:14547. doi: 10.1038/s41598-019-51020-3

PubMed Abstract | Crossref Full Text | Google Scholar

48. De Paoli F, Eeckhoute J, Copin C, Vanhoutte J, Duhem C, Derudas B, et al. The neuron-derived orphan receptor 1 (NOR1) is induced upon human alternative macrophage polarization and stimulates the expression of markers of the M2 phenotype. Atherosclerosis. (2015) 241:18–26. doi: 10.1016/j.atherosclerosis.2015.04.798

PubMed Abstract | Crossref Full Text | Google Scholar

49. Kamide K, Kokubo Y, Yang J, Takiuchi S, Horio T, Matsumoto S, et al. Association of intima-media thickening of carotid artery with genetic polymorphisms of the regulator of G-protein signaling 2 gene in patients with hypertension and in the general population. Hypertens Res. (2011) 34:740–6. doi: 10.1038/hr.2011.25

PubMed Abstract | Crossref Full Text | Google Scholar

50. Allen RM, Michell DL, Cavnar AB, Zhu W, Makhijani N, Contreras DM, et al. LDL delivery of microbial small RNAs drives atherosclerosis through macrophage TLR8. Nat Cell Biol. (2022) 24:1701–13. doi: 10.1038/s41556-022-01030-7

PubMed Abstract | Crossref Full Text | Google Scholar

51. Chang CY, Pan PH, Li JR, Ou YC, Wang JD, Liao SL, et al. Aspirin induced glioma apoptosis through noxa upregulation. Int J Mol Sci. (2020) 21:4219. doi: 10.3390/ijms21124219

PubMed Abstract | Crossref Full Text | Google Scholar

52. Guan L, Voora D, Myers R, Del Carpio-Cano F, and Rao AK. RUNX1 isoforms regulate RUNX1 and target genes differentially in platelets-megakaryocytes: association with clinical cardiovascular events. J Thromb Haemost. (2024) 22:3581–98. doi: 10.1016/j.jtha.2024.07.032

PubMed Abstract | Crossref Full Text | Google Scholar

53. Zhang X, Guan MX, Jiang QH, Li S, Zhang HY, Wu ZG, et al. NEAT1 knockdown suppresses endothelial cell proliferation and induces apoptosis by regulating miR-638/AKT/mTOR signaling in atherosclerosis. Oncol Rep. (2020) 44:115–25. doi: 10.3892/or.2020.7605

PubMed Abstract | Crossref Full Text | Google Scholar

54. Huang C, Yu XH, Zheng XL, Ou X, and Tang CK. Interferon-stimulated gene 15 promotes cholesterol efflux by activating autophagy via the miR-17-5p/Beclin-1 pathway in THP-1 macrophage-derived foam cells. Eur J Pharmacol. (2018) 827:13–21. doi: 10.1016/j.ejphar.2018.02.042

PubMed Abstract | Crossref Full Text | Google Scholar

55. Duan B, Yu Z, Liu R, Li J, Song Z, Zhou Q, et al. Tetrandrine-induced downregulation of lncRNA NEAT1 inhibits rheumatoid arthritis progression through the STAT3/miR-17-5p pathway. Immunopharmacol Immunotoxicol. (2022) 44:886–93. doi: 10.1080/08923973.2022.2092748

PubMed Abstract | Crossref Full Text | Google Scholar