Integrated Transcriptome Analysis Reveals Molecular Subtypes and ceRNA Networks in Multiple Sclerosis

Introduction

MS is an autoimmune disease characterized by diffuse and focal inflammation, demyelination, gliosis, and neuronal damage in the optic nerve, brain, and spinal cord.^1,2 Globally, MS affects over 2 million people, predominantly young adults aged 20–40, especially women. As the disease progresses, patients often develop cognitive deficits, mobility impairments, and loss of sphincter control, severely impacting their quality of life.^2–4 Clinically, MS is traditionally categorized into relapsing–remitting multiple sclerosis (RRMS), primary progressive multiple sclerosis (PPMS), and secondary progressive multiple sclerosis (SPMS), which guide treatment strategies.^4,5 However, despite these clinical categories, increasing evidence suggests that this simplistic classification does not fully capture the underlying complexity of MS pathology.

MS is a complex immune-mediated disease involving a combination of genetic, environmental, and epigenetic factors. Approximately 25% of MS risk is attributed to genetic heritability, while the remaining risk is influenced by environmental exposures, gene-environment interactions, and epigenetic modifications.^6,7 This intricate interplay among various factors highlights the need for a more refined understanding of the disease, particularly at the molecular level.

To advance this understanding, identifying new molecular subtypes of MS based on transcriptomic data has become crucial. Transcriptomics allows for a more objective, reliable, and reproducible classification system, addressing the limitations of traditional clinical subtyping.⁸ By analyzing gene expression profiles, we can identify distinct molecular signatures that may reflect different disease mechanisms and response to treatment. This approach has been successfully applied to various diseases, including cancer,^9–12 and holds promise for uncovering novel molecular pathways and therapeutic targets in MS.

A growing body of evidence underscores the importance of non-coding RNAs (ncRNAs) in regulating gene expression. NcRNAs, including miRNAs, circRNAs, and long non-coding RNAs (lncRNAs), play crucial roles in both immune and neural systems by modulating mRNA stability and translation.¹³ While lncRNA-miRNA-mRNA networks have been extensively studied, emerging evidence highlights that circRNAs are equally critical in regulating gene expression, particularly through their unique covalently closed loop structure.¹⁴ CircRNAs act as molecular sponges, sequestering miRNAs and thereby regulating downstream gene expression involved in processes such as immune modulation and neuroinflammation.^15–17 For instance, circINPP4B promotes Th17 differentiation in RRMS patients, reflecting MS pathogenesis,¹⁸ while circ_0000518 exacerbates MS progression by modulating macrophage/microglia polarization.¹⁹ These findings underscore the indispensable role of circRNA-miRNA-mRNA regulatory networks in disease progression, suggesting the need for further investigation into their specific contributions to MS molecular subtyping and therapeutic targeting.

MiRNAs, small (~22 nucleotides) single-stranded RNAs, regulate biological processes by binding to the 3’ UTR of target mRNAs to suppress their translation.^20,21 MiR-548a-3p has been identified as a biomarker for the absence of disease activity in MS patients treated with fingolimod.²² Additionally, miR-126a-5p levels in microglia correlate with blood-brain barrier (BBB) integrity, highlighting its role in early MS lesions.²³

Given the involvement of non-coding RNAs in MS pathogenesis, this study integrates transcriptomic data to develop a novel MS classification scheme. By combining circRNA and miRNA datasets, we construct a competing endogenous RNA (ceRNA) network, enabling a more comprehensive investigation of the molecular mechanisms underlying MS. This approach provides a promising framework for uncovering novel biomarkers and therapeutic targets that may lead to more precise and individualized treatment strategies for MS.

Materials and Methods

Data Collection

Gene expression profile datasets for peripheral blood samples from individuals with MS and healthy controls (HC) were retrieved from the Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo). The mRNA datasets included three series: GSE17048 with 99 MS cases and 45 HC; GSE190847 with 93 MS cases and 28 HC; and GSE173789 with 23 MS cases and 14 HC. Together, these datasets contain 215 MS patient samples and 87 HC samples. The circRNA dataset GSE161196 consists of samples from 10 MS cases and 10 HC. The miRNA dataset, GSE17846, contains samples from 21 MS cases and 20 HC. Finally, the GSE235357 dataset for human MS disease and the GSE253318 dataset for the mouse model of experimental autoimmune encephalomyelitis (EAE) were downloaded and used to test the results of the analysis.

Cluster Analysis for Subtyping

Prior to the elimination of potential batch effects, the mRNA dataset underwent normalization via log2(x+1) transformation, as facilitated by the SVA package’s combat function. Subsequently, 215 MS patients were classified into subtypes using the ConsensusClusterPlus package, and the resulting clustering was visualized using principal component analysis (PCA).

Identifying Molecular Markers for MS Subtypes

Differentially expressed genes (DEGs) for each subtype compared to HC were identified using the limma package, with |logFC| > 1 and P < 0.05 as the thresholds. The coefficient of variation (CV) of the DEGs was calculated across samples for each subtype. Among the top ten genes with the smallest CV values, the gene with the highest differential expression was selected as the stable and specific molecular marker for each subtype.

GSEA Enrichment Analysis and Immunological Profiling for MS Subtypes

Utilizing pathway information from the bioinformatics database Kyoto encyclopedia of genes and genomes (KEGG) as a functional gene set file, gene set enrichment analysis (GSEA) was conducted using the clusterProfiler package. Protein-protein interaction (PPI) networks were analyzed using the STRING (https://string-db.org/) database with default parameters. All network visualizations in this study were further optimized using Cytoscape software. Immune scores for each sample in the combined MS dataset were calculated using the estimate package, and one-way analysis of variance (ANOVA) was performed on the scores grouped by subtype to evaluate statistical significance of differences. Single-sample gene set enrichment analysis (ssGSEA) is a computational method utilized to estimate the presence and activity levels of various biological processes or cell populations within complex tissue samples.²⁴ This technique has been particularly valuable in quantifying the relative infiltration of a diverse array of immune cell types into disease.

Functional Analysis of DEmiRNAs and Their Target Genes

Using the GEO2R tool (http://www.ncbi.nlm.nih.gov/geo/geo2r/), differentially expressed microRNAs (DEmiRNAs) were identified based on screening criteria of |logFC| > 1 and P < 0.05. The targeting relationships between DEmiRNAs and mRNAs were explored by the miRTarBase database. Subsequently, a KEGG enrichment analysis was performed on the target genes of DEmiRNAs to uncover their impacts on biological processes.

Analysis of DEcircRNAs and the Functional Characterization of Their Target Genes

Using the edgeR package, differentially expressed circRNAs (DEcircRNAs) were identified with screening criteria of |logFC| > 0.5 and P < 0.05. Annotate DEcircRNAs using the circBase database. Download the spliced sequence lengths of the annotated DEcircRNAs and compare them to the reference genome from the Ensembl database (GRCh37.p13, GCA_000001405.14) to identify DEcircRNA types. We conducted a KEGG pathway enrichment analysis on the source genes of the annotated DEcircRNAs. By integrating predictions from miRanda and TargetScan along with circBank data, we constructed the interaction network between DEcircRNAs and miRNAs.

Machine Learning-Based Screening for Feature Target Genes in MS Subtypes

We identified feature immune genes specific to each subtype using random forest (RF) algorithms and least absolute shrinkage and selection operator (LASSO) regression. LASSO regression was implemented using the glmnet package to reduce data dimensionality, selecting the minimum lambda value as the optimal parameter. Initially, the RF model was constructed with the randomForest package, setting “ntree” at 500 for exploratory analysis. Subsequently, an optimized analysis was performed based on the “ntree” value that resulted in the lowest error rate, to identify feature genes. Ultimately, the feature target genes for each subtype were defined by the overlap of gene lists obtained from both analytical approaches.

Results

mRNA Analysis Reveals MS Subtypes and Biomarkers

Integrating mRNA datasets, we derived co-expression profiles for 9848 genes across 302 samples, comprising 215 MS cases and 87 HC. Subsequently, consensus clustering analysis was conducted, initially dividing the MS samples into k clusters (k = 2–9). The consensus score matrix, in conjunction with the cumulative distribution function (CDF) curve of the proportion of ambiguous clustering (PAC) statistics, suggested an optimal cluster number of k = 3. Accordingly, we classified MS into three subtypes: C1, C2, and C3, consisting of 127, 43, and 45 MS patient samples, respectively (Figure 1A). PCA analysis confirmed the clear differentiation between these subtypes (Figure 1B).

Figure 1 Disease typing of MS and screening for stable molecular markers. (A) Consensus score matrix for all samples when k=3. The higher the consensus score between two samples, the more likely they are to be grouped into the same cluster across different iterations. (B) PCA analysis of samples post-clustering; (C–E) Volcano plots for differential expression of each subtype, sequentially for C1, C2, C3; (F) Venn diagram illustrating the overlap of DEGs across subtypes; (G–I) Boxplots of the expression levels of the selected stable molecular markers. **P < 0.01.

For each subtype, differential expression analysis was performed to identify DEGs. Subtype C1 featured 2 downregulated genes and 141 upregulated genes (Figure 1C); subtype C2 presented 44 downregulated genes and 27 upregulated genes (Figure 1D); subtype C3 displayed 1102 upregulated genes (Figure 1E). A Venn diagram analysis of the DEGs across subtypes disclosed that subtype C1 possessed only 1 unique DEG, subtype C2 held 34 unique DEGs, and subtype C3 contained 954 unique DEGs (Figure 1F).

To ascertain stable molecular markers, we computed the CV for HC and MS patient samples within each subtype’s expression matrices. Among the top ten uniquely DEGs with the highest |logFC| values per subtype, those with the lowest CV values were chosen. Ultimately, FCRL1 (logFC = 1.04, CV = 0.043) was designated as the molecular marker for subtype C1, BTG1 (logFC = −1.03, CV = 0.039) for subtype C2, and RPL38 (logFC = 1.15, CV = 0.009) for subtype C3 (Figure 1F). Hence, the three subtypes were labeled as MS-FCRL1, MS-BTG1, and MS-RPL38 (Figure 1A-I).

GSEA Analysis Explores Pathway Signatures in MS Subtypes

We chose the KEGG database as our functional gene set and performed GSEA for each subtype, referencing the HC. Among the significant outcomes, 27 shared pathways across the three subtypes (Figure 2A) fell into five main categories, primarily linked to organismal systems and human diseases. Further analysis of secondary classifications revealed that these pathways predominantly concerned the immune system and viral infections, aligning with the etiology of MS. Noteworthy, key MS-related pathways in the KEGG database include “antigen processing and presentation”, “Th1 and Th2 cell differentiation”, and “Th17 cell differentiation”.

Figure 2 Pathway analysis of MS molecular subtypes. (A) KEGG pathways commonly engaged in the GSEA analysis across the three subtypes. (B–D) KEGG pathways uniquely involved in the GSEA analysis for each subtype, sequentially for C1, C2, and C3.

Moreover, each MS subtype displayed uniquely enriched pathways significantly divergent from others (Figure 2B–D), indicative of specific underlying influences on each subtype. For subtype C1, noteworthy pathways include immunity-related ones such as “Fc gamma R-mediated phagocytosis” and “N-Glycan biosynthesis”, along with the virus-specific “viral life cycle - HIV-1” (Figure 2B). Subtype C2 featured immunity-related pathways, including “primary immunodeficiency”, and virus-specific pathways like “Virion - Herpesvirus” (Figure 2C).

Subtype C3 encompassed a broad spectrum of pathways, including immunity (eg, toll-like receptor signaling), metabolism (eg, citrate cycle), viruses (eg, human cytomegalovirus infection), as well as nervous system-related pathways (dopaminergic synapse) and those linked to neurodegeneration (Figure 2D). Thus, we hypothesize that subtype C3 patients might experience more severe disease manifestations.

Investigating the Hub Gene Sets of Pathways in MS Subtypes

Our subsequent analysis centered on PPI among proteins encoded by hub genes within annotated pathways for each subtype, while simultaneously contrasting the disparities in hub gene sets among the subtypes (Figure 3A–C). Specifically, in subtype C1, RFXAP played a singular role in antigen processing and presentation, while NFATC2 was uniquely involved in Th17 cell differentiation. Distinctive core genes were notably enriched at the apex of the running enrichment score (RES) curves (Figure 3D–F), suggesting potential increased activity in pathways pertinent to MS, thereby promoting pathway activation.

Figure 3 Analysis of hub gene sets involved in pathways for MS subtypes. (A-C) Interaction networks of proteins encoded by hub genes participating in pathways, with genes uniquely contributing to pathways in each subtype annotated. (D–F) GSEA normalized enrichment score plots for pathways in C1, C2, and C3, respectively, with pathway colors consistently applied across the three figures. (G) GSEA normalized enrichment score plot for TCA cycle in C3. (H) Annotative chart of hub genes involved in the citrate cycle for C3.

The citrate (TCA) cycle is of paramount significance in organismal metabolism. It functions not only as the primary source for energy acquisition but also as a critical nexus connecting carbohydrate, fat, and protein metabolism. Moreover, it provides small-molecule precursors essential for other metabolic pathways. Our findings showed that in subtype C3, 23 genes were actively involved, stimulating the TCA cycle (Figure 3G and H).

Deciphering the Immune Cell Infiltration Characteristics of MS Subtypes

The etiology of multiple sclerosis (MS) is highly complex and closely associated with dysregulated immune system activity. To explore this phenomenon, we generated a heatmap illustrating 28 distinct immune cell clusters, providing a clear visualization of the relative abundance of infiltrating immune cell types (Figure 4A). The analysis revealed that subtype C3 had significantly higher immune cell density compared to the other subtypes, a finding supported by the immunoscore data. Specifically, subtype C3 exhibited significantly elevated immunoscores, reflecting a stronger overall immune infiltration (Figure 4D). Collectively, these findings suggest that subtype C3 is characterized by a unique and highly active immune microenvironment compared to other subtypes.

Figure 4 Immunophenotyping analysis of MS molecular subtypes. (A) ssGSEA was performed to identify the relative infiltration of immunocyte clusters in 215 MS samples using expression data. The relative infiltration of each cell type was normalized as a z-score. (B) Anti-disease immunity (ActCD4, ActCD8, TcmCD4, TcmCD8, TemCD4, TemCD8, Th1, Th17, ActDC, CD56brightNK, NK, NKT) and cell types executing pro-disease, immune-suppressive functions (Treg, Th2, CD56dimNK, immDC, TAM, MDSC, Neutrophil, pDC) were analyzed. The Pearson’s correlation coefficient (R) was calculated. The shaded area indicates the 95% confidence interval. (C) Violin plot depicting the immunological scores. (D) Boxplots were drawn to statistically group the results from C according to the defined subtypes. *P < 0.05, **P < 0.01, ***P < 0.001.

Analysis of immune cell infiltration levels revealed that most cells in subtype C3 differed significantly in score from those in the other subtypes. Exceptions included: activated B cells differed significantly between subtypes C1 and C2; activated dendritic cells and memory B cells, between subtypes C2 and C3; MDSCs and effector memory CD8 T cells, across all three subtypes (Figure 4C).

Pearson correlation analysis of the infiltration densities of activated immune cells (eg, activated CD8+ T cells, type 1 T helper cells, and activated dendritic cells) and those providing immunosuppressive functions (eg, MDSCs, regulatory T cells, immature dendritic cells, and neutrophils) revealed a positive correlation in their abundances within the local environment (Figure 4B). This indicates a potential feedback regulatory mechanism within the immune system of MS patients: activated immune cells might not only directly assault self-antigens in the CNS through cytokine release and other signaling molecules, but also indirectly stimulate the recruitment and differentiation of immunosuppressive cells.²⁵

Analysis of miRNAs in MS Patients

Through analysis of miRNA expression profiles, we identified 54 DEmiRNAs—19 downregulated and 35 upregulated (Figure 5A). Querying miRTarbase with these DEmiRNAs yielded predictions for 13,115 target genes. KEGG pathway enrichment analysis of the predicted targets uncovered 32 pathways related to immunity, neurology, metabolism, and viral infection (Figure 5C), involving 778 genes in total.

Figure 5 Analysis of miRNAs in MS patients. (A) Volcano plot depicting differential expression of miRNAs. Each point represents an individual miRNA, with the x-axis showing the log2 fold change and the y-axis showing the negative log10 of the p-value. (B) Heatmap illustrating the expression profiles of target genes for miRNAs. Rows represent different genes, and columns represent samples. The color scale reflects the normalized expression level of each gene across the samples. (C) KEGG pathway enrichment analysis chart for target genes of DEmiRNAs. Bar charts indicate the significance of enrichment, and line graphs represent the ratio of genes involved in each pathway to the total number of genes in that pathway.

Of the 778 genes, 52 miRNAs regulate them (Figure 5B), with pronounced upregulation noted in subtype C3, then subtype C1, and minimal differences seen in subtype C2 relative to HC. This suggests DEmiRNAs influence critical pathways for disease progression in subtype C3 MS patients by binding to target mRNAs. Thus, priority should be given to validating how downregulated miRNAs affect disease progression in subtype C3 MS patients and assessing their potential as therapeutic targets.

Concurrently, it is worth noting that reduced gene expression in subtype C2 might indicate actions of upregulated miRNAs, suggesting this subtype may have unique pathophysiological mechanisms requiring additional scrutiny.

Analysis of circRNAs in MS Patients

Through circRNA expression analysis, we identified 158 DEcircRNAs, including 71 downregulated and 87 upregulated ones (Figure 6A). The annotation of DEcircRNAs was performed in accordance with the chromosomal position information provided by circBase. In light of the sequencing error, this paper employs a screening method wherein the length of sequenced fragments is aligned with that of known circRNA sequences. The ratio of the overlapping region to the combined length of both sequenced and known fragments is set at a minimum of 99%, resulting in a total of 145 fragments annotated, with an annotation rate of 91.8%. This represents a notable enhancement over the findings of the original article.²⁶

Figure 6 Analysis of cirRNAs in MS patients. (A) Volcano plot displaying differential expression analysis of circRNAs. (B) Scatter plot depicting the expression fold change of highly expressed circRNAs in MS patients versus their parental genes. (C) Scatter plot showing the expression fold change of lowly expressed circRNAs in MS patients compared to their parental genes. (D) Bar chart presenting the KEGG pathway enrichment analysis of DEcircRNAs. (E) Schematic diagram illustrating the splicing pattern of circRNAs acting as miRNA sponges. White sections indicate the splice segments derived from introns, whereas colored sections denote splice fragments originating from exons.

When comparing logFC values of DEcircRNAs and their parent genes across subtypes, subtype C3 showed the highest differential expression, followed by C1, and C2 displayed the smallest deviations (Figure 6B and C). KEGG pathway analysis of parent genes disclosed engagement in numerous pathways, notably the MAPK signaling pathway (Figure 6D) and long-term potentiation, crucial for neuronal signal transduction.

CircRNAs are classified into three categories based on their spliced sequence origins: exonic circRNAs (EciRNAs) from exons alone, intronic circRNAs (IciRNAs) from introns exclusively, and exon-intron circRNAs (EIciRNAs) from both exons and introns combined. Generally, EciRNAs are enriched in the cytoplasm,²⁷ whereas IciRNAs and EIciRNAs mainly accumulate in the nucleus.²⁸

Concentrating on the competitive ceRNA network involving circRNAs and miRNAs, we specifically analyzed pairs displaying inverse expression patterns among DEcircRNAs and DEmiRNAs. Our results suggested that 18 DEcircRNAs (6 EciRNAs and 12 IciRNAs) serve as molecular sponges for 22 DEmiRNAs (Figure 6E).

Regulation of Immune Cell Feature Genes by DEmiRNAs

We found that the 22 previously identified DEmiRNAs regulate 161 feature genes in 28 distinct immune cell types (Figure 7A). Analysis of the mean expression levels of these genes across MS subtypes and HC revealed that most genes were significantly upregulated in subtype C3, with subtype C1 showing a lesser extent of upregulation. Conversely, subtype C2 showed high variability in gene expression, with negligible changes in overall expression levels (Figure 7B). Hence, it is imperative to highlight the effect of downregulated miRNAs on target gene expression regulation and their likely critical role in MS pathogenesis. Our results showed that 9 downregulated miRNAs control the expression of 97 feature genes in 24 immune cell types, underscoring their extensive impact, thus demanding additional research to uncover their mechanistic underpinnings.

Figure 7 Regulation of immune cell feature genes by DEmiRNAs. (A) miRNA regulatory network map. Squares represent miRNAs, with the depth of color indicating the magnitude of fold change. Circles denote immune cell feature genes, with identically colored circles belonging to the same type of immune cell, and the color coding is consistent with that of the heatmap. (B) Heatmap of average expression levels of immune cell feature genes across various subtypes and in the HC.

Identification of Hub Target Genes for MS Subtypes

Within the MS co-expression profile, we identified 3501 target genes associated with the 22 DEGs. Referencing the HC, we utilized LASSO regression and RF methodologies to identify hub target genes unique to each MS subtype. Using the LASSO logistic regression algorithm, we identified 12, 11, and 10 hub genes for subtypes C1, C2, and C3, respectively (Figure 8A and B, E and F, I and J). Employing the RF model, we discerned 92, 102, and 40 hub genes for subtypes C1, C2, and C3, respectively (Figure 8C and D, G and H, K and L).

Figure 8 Identification of hub target genes for MS subtypes through machine learning. (A–D) C1 subtype analysis. (A) Ten-fold cross-validation for tuning parameter selection in the LASSO model for C1, with each curve representing an individual gene. (B) LASSO coefficient profiles for C1, with dashed lines indicating the optimal λ. (C) Relationship between the number of trees in the RF model and the error rate for C1. (D) Gene importance ranking in the RF tree for C1. (E–H) C2 subtype analysis, Details as ABCD. (I–L) C3 subtype analysis, Details as ABCD. (M–O) Venn diagrams showing overlap of marker genes identified by both LASSO and RF algorithms in C1, C2, and C3, respectively. (P) Venn diagram showing the intersection of overlapping genes from (M–O).

Venn diagram analyses revealed the intersection of hub target genes identified by both LASSO and RF methods, resulting in 5, 2, and 2 shared genes for subtypes C1, C2, and C3, respectively (Figure 8M–O). Remarkably, one shared gene, EEF1D, was detected between subtypes C1 and C3 (Figure 8P), hinting at possible shared mechanisms among subtypes.

Screening the Regulatory Network of Hub Target Genes

The findings of this study were integrated to construct a ceRNA regulatory network of pivotal functional genes across the spectrum of MS subtypes. The results demonstrated that 16 circRNAs functioned as molecular sponges for 10 miRNAs, competitively binding to these miRNAs and thereby influencing their expression and regulation of hub target genes (Figure 9A). Our analysis and expression data (Figure 9B) from MS patient samples indicate that circ_0045537, an EciRNA, may influence the regulation of TUBA1A in the cytoplasm by competitively binding to miR-196a-5p. Concurrently, the human MS disease dataset and the mouse EAE model dataset were downloaded once more for further analysis. The disease samples exhibiting high expression of RPL38 were designated as the C3 type, whereas the low-expression healthy samples were classified as the control group. TUBA1A expression was then observed, and the results were found to align closely with those of the previously analyzed data set, indicating the reliability of the findings (Figure 9C). This regulatory pathway may affect the dynamic balance of microtubules in neuronal and immune cells in MS patients with type C3 (Figure 9D), thus contributing to disease progression.

Figure 9 The ceRNA regulatory network of hub target genes. (A) Regulatory network map of hub genes. (B) The following heat map depicts the average expression of hub genes across subtypes and HC. (C) Expression heatmap of the validation dataset. In the human dataset, three samples were selected separately, while only two samples each were selected in the mouse dataset. (D) The following schematic diagram illustrates the potential regulatory mechanism of circ_0045537/miR-196a-5p/TUBA1A in C3-type MS. The red and green lines indicate the possible regulatory pathways that may be involved in this process. The red and green squares illustrate the expression trends of the three in patients.

Discussion

The intricate mechanisms of MS and the constraints of conventional classification have prompted the formulation of novel typing strategies to facilitate precision medicine. The ConsensusCluster classification method has been successfully applied in a variety of diseases,²⁹ including pancreatic ductal adenocarcinoma (PDAC),⁹ CRC,³⁰ and breast cancer,²⁹ demonstrating robust molecular typing capabilities. It is capable of identifying stable molecular subtypes and is therefore well suited to MS research, facilitating the elucidation of disease mechanisms. This study employed MS transcriptome data to delineate subtypes and, for the first time, classified MS into three molecular subtypes: MS-FCRL1 (C1) type, MS-BTG1 (C2) type, and MS-RPL38 (C3) type. This facilitates the accurate diagnosis of MS.

The GSEA analysis demonstrated that each subtype is significantly associated with a multitude of pathways pertaining to the immune system and viral infections, particularly those involving Epstein-Barr virus (EBV) infection. Epidemiological, serological, and virological evidence collectively indicate that EBV plays a pivotal role in the aetiology of MS.³¹ It is noteworthy that subtype C3 demonstrates a robust correlation with pathways linked to diverse neurodegenerative disorders, indicating potential parallels in regulatory mechanisms and suggesting a more pronounced neurological impairment compared to other subtypes. MS is typified by immunological dysregulation, leading to immune cell infiltration into the CNS.^2,6 The elevated immune cell infiltration density observed in subtype C3 not only substantiates the hypothesis of a more severe disease course but also underscores the pivotal role of immune-mediated neuronal injury in this subtype.

Pathways associated with MS, such as Th1 and Th2 cell differentiation (hsa04658), Th17 cell differentiation (hsa04659), and antigen processing and presentation (hsa04612), play a crucial role in MS pathogenesis.^32,33 Our study reveals that these pathways are not only commonly involved in all MS subtypes, but the specific genes driving their activation differ across subtypes, indicating their potential contribution to the distinct characteristics of each subtype. For instance, the involvement of HLA-DOA and IL1B in pathway activation is unique to subtype C3, suggesting that these genes may play a key role in the immune dysregulation and inflammatory responses that are characteristic of this subtype. HLA-DOA, a member of the human major histocompatibility complex (MHC) family, has been identified as a significant genetic risk factor for MS,^34,35 with variations in the MHC locus, including HLA-DRB1, being strongly associated with disease susceptibility.^36,37 Similarly, IL1B, a crucial cytokine involved in inflammatory processes, is likely central to the immune response and exacerbation of inflammation in MS.³⁸ These findings highlight the importance of specific gene-pathway interactions in driving the molecular heterogeneity of MS subtypes and underscore the potential of targeting these pathways for subtype-specific therapeutic strategies.

In this study, we identified 18 DEcircRNAs that act as molecular sponges for 22 DEmiRNAs, with nine downregulated miRNAs controlling signature genes in various immune cell types, possibly contributing to immune cell infiltration in subtype C3. Previous research has highlighted the roles of some of these miRNAs in other diseases, although direct evidence in MS remains lacking. For example, miR-548b-5p, elevated in subacute sclerosing panencephalitis (SSPE), is involved in immune modulation by targeting antiviral responses.^39,40 MiR-20b-5p, known for its role in Alzheimer’s disease (AD), regulates neuroinflammation and neuronal survival.^41,42 Similarly, miR-196a-5p, upregulated in Huntington’s disease (HD), is involved in neuroprotective mechanisms and may influence neuronal stress responses relevant to MS.^43,44 Additionally, miR-574-5p regulates immune responses by activating TLR8, and miR-576-3p plays a role in inflammatory pathways that could intersect with MS-related immune processes.^45,46 MiR-1-3p regulates Th17 differentiation, a key process in MS immunopathogenesis,⁴⁷ and may influence macrophage polarization through the circATP8A1/miR-1-3p/STAT6 axis.⁴⁸ Although these miRNAs have been implicated in various immune and neurological processes, further validation in MS models is needed to confirm their role in disease progression. Future studies should focus on experimental validation of these miRNAs and their corresponding circRNA networks in MS, as well as their potential as biomarkers for disease activity and therapeutic targets.

After identifying DEmiRNA target genes, machine learning refined our selection, highlighting EEF1D and TUBA1A as potential hub targets in subtype C3. EEF1D is involved in protein synthesis by catalyzing GDP-GTP exchange, which is critical for eEF1A activity.^49,50 EEF1D isoforms are expressed in the brain and spinal cord, suggesting a role in neural function and development, and mutations in EEF1D have been linked to neurological conditions.^49,51 TUBA1A codes for α-tubulin, a core component of microtubules that are essential for neuronal cell processes, such as migration and cargo transport, and abnormalities in TUBA1A are associated with neurological disorders.^52–55 Additionally, microtubules play a crucial role in immune cell integrity and immune responses, including cytokine release and synapse formation during immune activation.^56–59 These findings suggest that EEF1D and TUBA1A may influence immune and neural cell function in MS, but further experimental validation is required to clarify their specific role in MS pathogenesis. Future studies should focus on confirming the involvement of EEF1D and TUBA1A in MS, especially in the context of immune cell infiltration and neuroinflammation.

In conclusion, differential activation of immune-related pathways across MS subtypes suggests that specific gene expression drives disease heterogeneity and may provide targets for subtype-specific therapies. We also propose a circ_0045537/miR-196a-5p/TUBA1A regulatory axis in subtype C3, which modulates microtubule dynamics in immune and neuronal cells, potentially exacerbating MS severity.