Integrated transcriptomics of multiple sclerosis peripheral blood mononuclear cells explored potential biomarkers for the disease.

Introduction

Multiple sclerosis (MS) is characterized by neuroinflammation and demyelination of the central nervous system (CNS) [ 1 ]. Although the exact etiology of MS remains undetermined, autoimmunity is widely acknowledged as the underlying mechanism [ 2 ]. Immune cells in the periphery are implicated in the autoimmune response during progressive and relapsing-remitting courses of MS [ 3 ], particularly autoreactive T lymphocytes, which are substantially involved in the destruction of the myelin sheath [ 4 ].

After the activation in the periphery, T-cells pass through the blood-brain barrier (BBB), reaching the CNS where they undergo further reactivation through interactions with local antigen-presenting cells (APCs) and autoantigen recognition [ 5 ]. Among T-cells, CD4 + lineages, including Th1, Th17, and Th22, significantly contribute to the MS pathogenesis by releasing IFN-γ, IL-17, and IL-22, respectively. These cytokines eventually lead to the recruitment of other immune cells, inflammation, and increased permeability of the BBB [ 6 ]. CD8 + cells, comprising the largest proportion of T-cells in the CNS lesions. They express TNF upon activation, which can directly harm neurons and disrupt the restoration of the myelin sheath [ 7 ]. B lymphocytes not only produce autoantibodies but also, more significantly, act as APCs and express cytokines such as IL-6 and TNF-α [ 8 ]. Furthermore, emerging research indicates the involvement of innate immunity cells like NK cells [ 9 ], monocytes [ 10 ], and neutrophils [ 11 ] in the immunopathology of MS.

Magnetic resonance imaging (MRI) stands as the most common diagnostic procedure for MS [ 12 ]. However, it was recently reported that overdependence on the technique can lead to an increased number of misdiagnoses [ 13 ]. Exploration of peripheral blood mononuclear cells (PBMCs) transcriptome has the potential to identify potential biomarkers in neurodegenerative disorders such as Alzheimer's and Parkinson's disease [ 14 , 15 ].

This study aimed to identify key genes and modules in MS through the integration of PBMCs gene expression data. After integrating GSE21942 and GSE203241 , a total of 106 DEGs were identified. WGCNA showed that 59 genes had a correlation consistent with the direction of their expression levels. Moreover, the analysis singled out a module highly enriched in neurodegenerative disorders. Random forest (RF) and least absolute shrinkage and selection operator (LASSO) machine learning methods identified COPG1 as the most prominent potential biomarker. Although ROC analysis of the integrated data and two validation sets ( GSE141804 and GSE146383 ) demonstrated acceptable diagnostic efficacy for COPG1 , RPN1 , and KDM3B , the inconsistent directions of their alterations suggest that they cannot be considered reliable biomarkers for the disease. Furthermore, assessment of their diagnostic efficacy in the other external datasets ( GSE247181 , GSE59085 , and GSE17393 ) did not indicate diagnostic potential.

Methods

Data collection and preparation

The Gene Expression Omnibus (GEO) database [ 16 ] was reviewed to find datasets incorporating both MS and healthy control samples of PBMCs. The details of seven retrieved datasets, published between 2009 and 2023, are presented in Table 1 .

Table 1 Microarray datasets obtained from PBMCs of MS and healthy controls. Table 1 GEO accession Year Control MS Platform Reference GSE203241 2022 8 15 GPL571 [ 20 ] GSE21942 2010 15 13 GPL570 [ 21 ] GSE247181 2023 10 28 GPL24539 [ 22 ] GSE146383 2020 46 15 GPL571 [ 23 ] GSE141804 2019 7 15 GPL571 Unpublished GSE59085 2014 7 15 GPL570 [ 24 ] GSE17393 2009 7 8 GPL571 [ 25 ]

Microarray datasets obtained from PBMCs of MS and healthy controls.

The quality control stage of data processing included the detection and exclusion of outliers and unwanted samples. Considering that pediatric MS typically presents with a more aggressive onset and disabling symptoms compared to adults [ 17 ] and accounts for a small 3–10 % of MS patients [ 18 ], samples from pediatric patients were excluded from the GSE203241 and GSE146383 datasets. Using principal component analysis (PCA) and hierarchical clustering, GSM545838 from GSE21942 and GSM6165202 , GSM6165212 , and GSM6165227 from GSE203241 were detected as outlier samples and excluded from the subsequent analyses. Similarly, GSM7884845 was detected as an outlier in GSE247181 ( Supplementary Fig. S1 ). In GSE21942 , the samples GSM545843 and GSM545844 , along with GSM545845 and GSM545846 , are technical replicates from identical MS patients that were not merged due to an oversight during data processing. These four technical replicates were not detected as outliers in the PCA and hierarchical clustering dendrogram and were included among samples for further analysis.

After the integration of GSE21942 and GSE203241 , batch-effect removal was conducted using the SVA [ 19 ] package. Integrated data comprised 51 samples, including 23 healthy controls and 28 MS patients, and expression data on 13,237 genes. Although the integration aimed to increase the sample pool, power analysis revealed that the resulting composition yielded a power of 41.38 %, which does not ensure detecting existing differences among the groups.

Demographic information on the collected datasets

After retrieving sample information using the GEOquery [ 26 ] package, the gender and age of the studied population after filtrations were singled out. The gender count, as well as the mean and standard deviation (SD) of the case and control groups, were computed, which are presented in Table 2 . In GSE203241 , adult MS patients were diagnosed with RRMS based on the 2017 McDonald diagnostic criteria. Patients had not been treated with disease-modifying drugs (DMD) or steroids in a month prior to the transcriptomic analysis [ 20 ]. In GSE21942 , patients were diagnosed using the Poser's criteria. Meanwhile, it has not been clearly documented whether the patients were diagnosed with RRMS or SPMS in the study. Of the twelve female patients in the study, four patients had received immunomodulatory treatment, such as cortisone or interferon-beta (β-interferon), when samples were collected [ 21 ].

Table 2 Demographic data of the retrieved datasets. F: Female, M: Male, SD: Standard deviation, NA: Not available. Table 2 GEO accession Number of cases (F/M) Number of controls (F/M) Mean age in MS (SD) Mean age in controls (SD) GSE203241 15 (10/5) 8 (3/5) 31.2 (6.34) 33.125 (5.67) GSE21942 12 (12/0) 15 (15/0) 54.2 (NA) 71.6 (NA) GSE247181 28 (25/3) 10 (5/5) NA NA GSE146383 15 (8/7) 46 (24/22) 31.87 (6.05) 34.44 (7.06) GSE141804 15 (9/6) 7 (3/4) 36.05 (11.11) 33.04 (8.33) GSE59085 15 (NA) 7 (NA) NA NA GSE17393 8 (8/0) 7 (7/0) NA NA

Demographic data of the retrieved datasets. F: Female, M: Male, SD: Standard deviation, NA: Not available.

Differential expression and functional enrichment analyses

DEGs were selected based on an adjusted p-value of less than 0.05 and a |logFC| exceeding 0.5 using the limma [ 27 ] package. Subsequently, the ClusterProfiler [ 28 ] package was utilized to find the enriched Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) annotations. After excluding terms with an adjusted p-value above the significance threshold of 0.05, the top terms were ordered and visualized according to the proportion of involved genes.

Construction of miRNA-mRNA regulatory networks

The STRING database [ 29 ] was employed to retrieve the protein-protein interactions (PPI) of the DEGs using the confidence threshold set at 0.4. Then, nodes in the PPI network were used as an entry in the miRDB [ 30 ] to find miRNAs interacting with the provided list. A target prediction score of above 80, a maximum of 2000 targets in the human genome, and targeting 4 or more nodes in the network were used as criteria for identification of key regulatory miRNAs of the network. Subsequently, the PPI and miRNA-mRNA networks were constructed using Cytoscape v3.10.3 [ 31 ].

WGCNA

WGCNA was carried out to identify MS-specific modules by the WGCNA [ 32 ] package. After the determination of the optimal soft thresholding power, the analysis was performed through a blockwise method. A signed WGCNA analysis was performed, with modules required to have a minimum of 30 genes. Modules were merged if the correlation of their module eigengenes was above the 75 percent threshold. Lastly, modules were labeled numerically, and those with a p-value below 0.01 were deemed important. Then, pathway analysis was conducted for each of the modules using the Enrichr online tool [ [33] , [34] , [35] ].

Moreover, for module 7, which was revealed to be highly enriched in neurodegenerative disorders, the PPI network was retrieved from the STRING database with a high confidence score (0.7) to highlight the genes that not only are co-expressed at the transcript level but also are known to interact at the protein level. After excluding genes with no interaction, a network was constructed using Cytoscape. Highly connected nodes were selected by observational evaluation of the degree of connectivity in the network. Among these, genes with module membership above 0.8 and gene significance greater than 0.2 were discerned as hub genes in the module. KEGG pathways enriched with the hub genes were subsequently visualized by the circlize package [ 36 ].

Similarly, hub genes that met the criteria of having membership and gene significance above 0.8 and 0.2, respectively, were identified separately in each of the singled-out modules. Then, to provide insights on the biological functions of these hub genes, functional enrichment analyses were carried out.

Biomarker identification using machine-based learning methods

The intersection of upregulated genes with modules significantly correlated positively, and downregulated genes with modules significantly correlated negatively, was subjected to machine-based learning methods. After evaluating the optimum number of trees in each split (m try parameter), the RF algorithm was employed to assess the classification potential of genes, using the randomForest [ 37 ] package. Similarly, the LASSO regression was applied to identify the most prominent genes, using the glmnet [ 38 ] package.

At last, using the pROC package [ 39 ], the diagnostic efficacy of the top genes selected by the machine learning methods was assessed by conducting ROC analysis in the integrated data, as well as each of the datasets independently. Furthermore, ROC analysis was also carried out in external datasets, including GSE247181 , GSE146383 , GSE141804 , GSE59085 , and GSE17393 , to examine the robustness of the results.

Characterization of immune cell proportions in case and control groups

To determine whether differences in immune cell composition could confound the identification of DEGs, immune cell deconvolution was performed on the integrated PBMC data using the deconvolution of the absolute immune signal method [ 40 ]. After determining the proportion of specific immune cell types in each sample, the mean and standard deviation of these proportions were computed separately for MS patients and healthy controls. To assess the statistical significance of differences among groups in terms of immune cell proportions, the Wilcoxon rank-sum test was performed. The resulting p-values were adjusted for multiple comparisons using the false discovery rate (FDR) method.

Cell type-specific enrichment of DEGs

The Enrichr [ 41 ] online tool was employed to identify in which cell types of the blood the DEGs are preferentially expressed. The results from both CellMarker 2024 [ 42 ] and PanglaoDB Augmented 2021 [ 43 ] were considered. The results from the former database reflecting cell types for mice were excluded. Lastly, the top 10 enriched cell types, which were underscored by these two databases, were illustrated using the circlize [ 36 ] package.

Results

DEGs and functional enrichment

As Fig. 1 a–c reveal, samples clustered regarding the disease state rather than the study following the ComBat batch-effect removal. Among 13,237 analyzed genes, a total of 106 DEGs, including 71 upregulated and 35 downregulated genes, were identified ( Fig. 1 d). The complete list of log fold changes (logFC) and adjusted p-values for the 71 upregulated genes and 35 downregulated genes are presented in Supplementary Table S1 . By conducting the Gene Ontology (GO) test, the top 10 terms in each GO classification, including biological processes (BP), cellular components (CC), and molecular functions (MF), were visualized. The DEGs were significantly enriched in the immune response-regulating cell surface receptor signaling pathway (GO:0002768), the immune response-activating signaling pathway (GO:0002757), and the immune response-activating cell surface receptor signaling pathway (GO:0002429) biological processes. Similarly, other biological processes reflected the involvement of the DEGs in orchestrating autoimmune processes such as differentiation of lymphocyte and mononuclear cell (GO:0030098 and GO:1903131, respectively), inflammation (GO:0050727), and reactive oxygen species metabolic process (GO:0072593). In the cellular components category, the external side of the plasma membrane (GO:0009897) and blood microparticle (GO:0072562) were the most prominent terms. Tetrapyrrole binding (GO:0046906) and cytokine activity (GO:0005125) emerged as the key terms in the molecular function classification. In line with gene ontology findings, the B cell receptor (BCR) signaling pathway was the most enriched KEGG pathway, followed by the hematopoietic cell lineage, MAPK signaling pathway, and Th1, Th2, and Th17 cell differentiation ( Fig. 1 e).

Differential expression and functional enrichment analyses. ( a ) Datasets before batch-effect correction. ( b ) Datasets after batch-effect correction. ( c ) Samples are clustered according to the disease state after batch-effect correction. ( d ) A total of 106 DEGs were identified (|logFC| > 0.5 and adjusted p-value <0.05). ( e ) Gene ontology classifications and KEGG pathway enrichment analyses.

PPI and miRNA-mRNA regulatory networks

Of the 106 DEGs in MS, a PPI network incorporating 69 nodes and 278 edges was constructed using a 0.4 confidence threshold ( Fig. 2 a). The hub genes in this network were IFNG , MS4A1 , CD22 , CXCL8 , CD79A , BLK , CD79B , JUN , and CD69 , with interactions ranging from 22 for IFNG to 16 for JUN and CD69 . While IFNG , the central hub gene in the network, was downregulated, remaining hub genes were upregulated in MS patients compared to healthy controls. Furthermore, COPG1 , RPN1 , and STT3A appeared as discrete networks. Likewise, FBXO16 and SNRNP40 interacted with HNRNPL and created a separate network.

Networks of the DEGs in MS. ( a ) PPI network. The intensity of the border color reflects the logFC values. ( b ) The miRNA-mRNA network illustrates the key regulators of the PPI network. Blue label reflects downregulation, while red label represents upregulation. The size of the nodes positively correlates with their number of PPIs in the networks.

Construction of the miRNA-mRNA regulatory network for the PPI network demonstrated that members of the miR-92, miR-181, miR-15, miR-200, and miR-450 families are prominent regulators of the network. Furthermore, miR-124-3p, miR-493-5p, and miR-506-3p were also key miRNAs targeting several nodes ( Fig. 2 b). Notably, miR-124-3p and miR-506-3p also regulated STT3A , thereby influencing the discrete network of COPG1 , RPN1 , and STT3A ( Supplementary Table S2 ).

In addition, exploring key regulators for the hub genes further underscored the regulatory role of members of the miR-181 family, including miR-181a-5p, miR-181b-5p, miR-181c-5p, and miR-181d-5p, as they targeted more than one hub gene ( MS4A1 and CD69 ). In addition, miR-493-5p was another multi-target miRNA that was revealed as a regulatory miRNA for CXCL8 and JUN ( Supplementary Table S2 ).

Identification of MS-specific modules by WGCNA

After selecting number 11 as the optimal soft-thresholding power ( Fig. 3 a and b), blockwise module construction assigned 13,237 genes into 15 modules ( Fig. 3 c). Modules with an absolute correlation value above 0.4 and a p-value less than 0.01 were considered significant. These filtrations led to the identification of 3 positively correlated and 5 negatively correlated modules.

WGCNA detection of gene modules correlating with the disease state. ( a ) Clustering of samples after removing outliers. ( b ) A soft-thresholding power of 11 was selected as the optimal power. ( c ) Correlation of module eigengenes with disease state and their corresponding p-values (∗P < 0.05, ∗∗P < 0.01, ∗∗∗P < 0.001). ( d ) The pie chart illustrates the proportion of each module. The most prominent KEGG pathways are presented below each module name based on the combined score in the Enrichr database. ( e ) The overlap of genes between DEGs and WGCNA analyses. Purple indicates positively correlated modules, whereas green indicates negatively correlated modules.

Among the positively correlated modules, module 9 was enriched in the FoxO signaling pathway. Meanwhile, the B cell receptor signaling pathway, hematopoietic cell lineage, and primary immunodeficiency were the most significant KEGG pathways associated with module 12. Pathway enrichment analysis failed to identify any pathway with an adjusted p-value less than 0.05 for module 14.

Of the negatively correlated modules, module 7 was highly enriched in neurodegenerative disorders; therefore, it was singled out for further analysis. Module 1 was functionally associated with lysosome and phagosome pathways. The ubiquitin-mediated proteolysis and the phosphatidylinositol signaling system were the most prominent biological pathways in which genes of module 2 were involved. Modules 3 and 8 were significantly enriched in the MAPK signaling pathway and protein processing in the endoplasmic reticulum, respectively ( Fig. 3 d).

To identify the potential biomarkers, 59 genes, including 31 upregulated and 28 downregulated genes observed in the positively and negatively correlated modules, respectively, underwent machine-based learning methods ( Fig. 3 e). The entire WGCNA statistics for the 59 shared genes between the analyses can be found in Supplementary Table S3 .

Identification of hub genes in the module involved in neurodegeneration

Considering the neurodegenerative aspect of MS, an in-depth analysis was carried out for module 7. Of the 250 genes in module seven, 23 genes were selected as highly connected genes because they stand in the top 10 % of genes in module seven based on their number of interactions in this module [ 44 ]. The cutoff was set because the number of interactions dropped significantly for the genes located after NDUFB2 from 24 for this gene to 14 for ATP5PF ( Fig. 4 a).

In-depth analysis aiming at the identification of hub genes in module seven. ( a ) High-confidence (0.7) PPIs of hub genes for the co-expressed genes in module 7. ( b ) Gene significance against module membership for this module. The 23 central genes in the PPI network are colored. ( c ) Gene ontology analysis of the 23 central genes in the PPI network. ( d ) KEGG pathway analysis of the 23 central genes in the PPI network. Pathways with the subcategory of neurodegenerative disorders are represented by the green color.

In addition to high connectivity, hub genes also meet the criteria of having a module membership above 0.8 and a gene significance for the trait of interest exceeding 0.2 [ 45 ]. By applying these additional cutoff thresholds, 11 genes, namely COX5B , UQCRQ , NDUFS3 , UQCR10 , NDUFA4 , COX7B , NDUFA7 , NDUFS7 , NDUFA2 , NDUFA1 , and NDUFB2 , were determined as hub genes ( Fig. 4 b).

Gene ontology analysis ( Fig. 4 c) showed that the 23 nodes in the PPI network of module 7 were involved in the structure and function of mitochondria. Moreover, the subcategory of neurodegenerative disorders, including Parkinson's, prion, Huntington's, Alzheimer's, and amyotrophic lateral sclerosis, contributed to the largest proportion of KEGG pathways enriched with these genes ( Fig. 4 d).

WGCNA-based identification of hub genes and related pathways in the other modules

The hub genes in the remaining prominent modules were discerned based on module membership and gene significance, and the pathways they are involved in were determined ( Supplementary Fig. S2 ).

Modules 1, 2, 3, 7, 8, and 9 were selected from modules that exhibited a negative correlation with MS. Of the 728 genes in module 1, 102 were identified as hub genes in the module. Considering both GO terms and KEGG pathways, the genes were predominantly involved in innate immunity, inflammation, and myeloid cell-mediated immune defense processes. Among the 354 genes in module 2, 40 genes were identified as hub genes and were enriched in metabolic adaptation and proliferation, including the Ras signaling pathway. Of the 289 genes in module 3, 48 were identified as hub genes, and they contributed to RNA splicing and transport, particularly nucleocytoplasmic transport, which was the only significant pathway (p = 3.10e-2). Of the 200 genes in module 8, 42 were identified as hub genes. Functional enrichment analysis revealed terms reflecting various aspects of homeostasis, including protein processing in the endoplasmic reticulum [ 46 ], ATP-dependent chromatin remodeling [ 47 ], and cellular senescence [ 48 ].

Modules 9 and 12 were positively correlated with MS and identified as modules of interest based on their statistical significance (p < 0.01) and initial enrichment analysis ( Fig. 3 d). Among the 170 genes in module 9, 41 were identified as hub genes. The T cell receptor signaling pathway was the most prominent pathway identified by functional enrichment analysis. Of the 98 genes in module 12, 37 were identified as hub genes, and functional enrichment analysis highlighted the enrichment of the B cell receptor signaling pathway.

Identification of potential biomarkers using machine-based learning models

The 59 genes identified by both differential expression and WGCNA analyses were subjected to RF and LASSO methods. Interestingly, COPG1 was evaluated as the most prominent gene by both machine-based learning algorithms, reflecting its significance in the integrated data ( Fig. 5 ). According to the RF results, RPN1 and KDM3B stood in second and third places, respectively, based on the MeanDecreaseGini and MeanDecreaseAccuracy values. However, in LASSO, ADAM28 was the second most prominent gene, followed by SMG8 , which had a coefficient markedly lower than that of ADAM28 ( Fig. 5 e). Therefore, with a holistic view, COPG1 , RPN1 , ADAM28 , and KDM3B were selected for further analysis, as they were underscored by the machine learning methods.

Employing machine learning methods for the identification of potential biomarkers. ( a ) Error against the number of trees used in the RF classifier. ( b ) Importance of genes based on the mean decrease in accuracy. ( c ) Importance of genes based on the mean decrease Gini. ( d and e ) Coefficient paths and binomial deviance plotted against the logarithm of lambda used in LASSO regression, respectively. ( f ) Genes with a non-zero coefficient were ordered according to the absolute value of the coefficient.

Biomarker validation

By conducting ROC analysis, the diagnostic efficacy of COPG1 , RPN1 , KDM3B , and ADAM28 was independently assessed in GSE21942 , GSE203241 , and the integrated data, as well as in GSE247181 , GSE146383 , GSE141804 , GSE59085 , and GSE17393 as external datasets. As Fig. 6 reveals, COPG1 , RPN1 , and KDM3B had an area under the ROC curve above 0.7 in the integrated data, GSE21942 , GSE203241 , and two of the validation sets ( GSE146383 and GSE141804 ). Meanwhile, these three genes demonstrated low reproducibility of diagnostic efficacy in GSE247181 , GSE59085 , and GSE17393 ( Fig. 6 ).

ROC analysis evaluating the diagnostic efficacy of ( a ) COPG1 , ( b ) RPN1 , ( c ) KDM3B , and ( d ) ADAM28 across the datasets. Colors represent the diagnostic performance of the gene of interest across different datasets.

Expression pattern of selected biomarkers

Although the area under the ROC curve of COPG1 , RPN1 , and KDM3B passed the 0.7 threshold in GSE141804 and GSE146383 , the direction of their gene expression alterations in these two datasets was not consistent with that of in GSE21942 and GSE203241 ( Fig. 7 ). While COPG1 , RPN1 , and KDM3B were downregulated in GSE21942 and GSE203241 ( Fig. 7 a and b), they were upregulated in GSE141804 and GSE146383 ( Fig. 7 c and d). Given this inconsistency, along with their insufficient performance as potential biomarkers in GSE247181 , GSE59085 , and GSE17393 , they were not potential biomarkers for MS.

Expression level of COPG1 , RPN1 , and KDM3B in ( a ) GSE21942 , ( b ) GSE203241 , ( c ) GSE141804 , and ( d ) GSE146383 . Blue and red boxes demonstrate statistically significant downregulation and upregulation, respectively, whereas the alterations of genes with black boxes were not statistically significant.

Immune cell profiles and cell type-specific enrichment of DEGs in cases and controls

Deconvolution of the absolute immune signal from PBMC transcriptomes in the integrated data revealed a greater abundance of naïve B cells and basophils in MS patients compared to healthy individuals. The mean proportion of naïve B cells was 2.76 in the control group, whereas it was substantially higher in MS patients at 7.95 than that of controls. Similarly, the mean proportion of basophils was 1.86 in healthy controls but markedly higher in MS patients at 2.87. These statistically significant differences (adjusted p-value <0.01) in the proportions of these two immune cell types between cases and controls may have partially confounded the findings ( Table 3 ). Considering the enrichment of DEGs in B cell subtypes, this assumption was further strengthened ( Fig. 8 ). In the CellMarker database, the genes that reflected B cell enrichment were upregulated genes. Similarly, enriched cell types underscored by the PanglaoDB Augmented 2021 database mainly involved upregulated genes. Collectively, it is plausible to convey that identification of a proportion of upregulated genes stems from a higher number of B cells in MS patients compared to healthy individuals.

Table 3 Immune cell proportions in patients with MS and healthy controls. NK: Natural Killer cells, pDCs: Plasmacytoid dendritic cells, mDCs: Myeloid dendritic cells, LD: low‐density. Table 3 Cell type Mean in control (SD) Mean in MS (SD) p-value Adjusted p-value T Naïve 18.83 (11.95) 23.73 (18.58) 3.8e-01 5.3e-01 T Memory 24.17 (10.36) 23.35 (9.08) 9.4e-01 9.4e-01 B Naïve 2.76 (2.13) 7.95 (7.44) 1.5e-03 8.4e-03 B Memory 1.73 (0.75) 1.92 (0.93) 3.9e-01 5.3e-01 Plasmablasts 0.36 (0.11) 0.43 (0.15) 7.8e-02 2.2e-01 NK 18 (6.09) 16.39 (4.97) 4.4e-01 5.3e-01 pDCs 2.32 (0.46) 2.27 (0.83) 3.8e-01 5.3e-01 Neutrophils LD 2.08 (0.54) 2.34 (0.62) 1.7e-01 3.8e-01 Basophils LD 1.86 (0.48) 2.87 (1.33) 1.5e-03 8.4e-03 mDCs 0.57 (0.23) 0.53 (0.14) 9.2e-01 9.4e-01 Monocytes 25.92 (5.81) 22.69 (5.26) 5.6e-02 2.0e-01

Immune cell proportions in patients with MS and healthy controls. NK: Natural Killer cells, pDCs: Plasmacytoid dendritic cells, mDCs: Myeloid dendritic cells, LD: low‐density.

Cell-type enrichment of DEGs in MS. Plots showing enriched cell types based on ( a ) the CellMarker 2024 and ( b ) the PanglaoDB Augmented 2021 databases.

Discussion

Despite recent advancements in the diagnosis of MS, the exploration of non-invasive biomarkers could significantly improve diagnostic accuracy. Herein, an integrated microarray analysis of peripheral blood mononuclear cells was conducted to identify potential biomarkers for the disease. However, the biomarkers identified through the analysis of the integrated data of GSE203241 and GSE21942 did not result in the identification of acceptable diagnostic markers. Factors such as small sample size, specific demographic characteristics, MS subtype specificity, undetermined disease stage, pretreatment conditions in some GSE21942 samples, and varying diagnostic criteria across studies may explain the lack of generalizability of the identified biomarkers. Meanwhile, WGCNA identified a module highly enriched in neurodegenerative-based disorders. The hub genes in this module, considering the number of interactions, module membership, and gene significance for MS, were genes encoding multiple components of the respiratory chain complex. While the mitochondrial impairment in brain lesions has long been documented. Its impairment in the immune cells has recently been reported, and the findings of this study further support the idea. While this study did not unravel promising biomarkers, further transcriptomic studies are needed to identify biomarkers from blood samples as valuable non-invasive diagnostic sources.

A total of 106 genes were found to be differentially expressed, which were functionally enriched in well-established pathways involved in the pathogenesis of MS. Gene ontology analysis reflected the contribution of DEGs to various aspects of immune system functions. Highly enriched biological processes were predominantly engaged in lymphocyte differentiation, production of immune mediators, and regulatory signaling pathways. Most significantly, immune response-regulating cell surface receptor signaling pathways, such as TLR [ 49 ], JAK/STAT [ 50 ], PD-1/PD-L1 [ 51 ], and TGFβ [ 52 ], are implicated in the immune responses in MS.

Similarly, the DEGs were associated with numerous immune-related components and complexes. The most significant cellular component term, the external side of the plasma membrane, has been previously reported to be enriched with common DEGs in neurodegenerative disorders [ 53 ]. Inflammation, affecting the MS brain in almost all clinical stages of the disease, has been identified to be linked to neurodegeneration through several plausible mechanisms of action [ 54 ]. Meanwhile, it remains controversial whether neurodegeneration or neuroinflammation is the primary causal factor [ 54 , 55 ]. In line with the ranking of blood microparticles as the second most significant term in this category, elevated levels of blood microparticles originating from endothelial and platelet cells were observed in different types of MS, disrupting the barrier function of endothelial cells in vitro [ 56 ]. In the molecular function category, tetrapyrrole binding and cytokine activity were the most prominent terms, followed by terms mainly associated with immune response and molecular signaling. Tetrapyrrole, heme, and haptoglobin binding-related terms were similarly enriched with genes of the hemoglobin complex, including HBA1 , HBD , HBG1 , and HBB . Although the exact role of these genes in MS pathogenesis has yet to be elucidated, recent advancements in the field have hypothesized that these genes might be involved in functions beyond oxygen transport, such as the role of HBD in counteracting oxidative stress [ 57 ].

KEGG pathways reflected the enrichment of pathways involved in autoimmune disorders, such as the BCR signaling pathway, whose dysregulations have been widely reported in autoimmune disorders, including MS, RA, and SLE [ 58 ]. Previous investigations have demonstrated that B cell subsets contribute to the pathogenesis of MS, and the BCR signaling pathway holds potential as a therapeutic target [ 59 ]. While their exact role remains unclear, B cells contribute to MS through both antibody-dependent and antibody-independent mechanisms. The former is supported by the therapeutic outcome of anti-CD20 (aCD20) therapy. Plasma cells and plasmablasts, the primary antibody-producing cells, do not express CD20 or express it only in small amounts and thus are unaffected by the therapy. Since released antibodies are not impacted by the treatment, the significant reduction in disease activity cannot be attributed to the antibody production role of B cells, suggesting that B cells also play an antibody-independent role in the pathogenesis of MS [ 60 ]. On the other hand, the observed efficacy of the B cell depletion therapy for autoimmune demyelinating disorders, including MS, underscores the role of B cells in MS [ 61 ].

In accordance with the identification of hematopoietic cell lineage as the second most enriched pathway, a recent study demonstrated that bone marrow hematopoietic stem and progenitor cells (HSPCs) exhibit dysregulated activity in MS patients. The dysregulation results in elevated production of myeloid cells, particularly monocytes and neutrophils, which subsequently are able to infiltrate the CNS and contribute to the pathogenesis of MS [ 62 ]. The flow cytometry analysis of peripheral blood mononuclear cells from MS patients identified MAPK signaling, the third most enriched pathway in this study, as a key pathway responsible for the multiplication and survival of immune cells, particularly B cells [ 63 ]. Cellular differentiation of the Th1, Th2, and Th17 were pathways that referred to the well-established role of T helper cells [ 64 ]. Pathogenic differentiation of Th1 and Th17 cells and their cytokine profiles plays a central role in the pathogenesis of MS [ 65 ].

Of the top 10 KEGG pathways, three pathways were classified under the subcategory of infectious diseases, including two parasitic infections, namely Chagas disease and leishmaniasis, as well as a viral infection that was Epstein-Barr virus (EBV) infection. Consistently, it has previously been reported that Chagas disease affects the CNS in both acute and chronic phases, with acute infection often causing meningoencephalitis and chronic infection leading to cognitive issues [ 66 ]. While the co-occurrence of leishmaniasis and MS in natalizumab-treated MS patients has been previously observed, and it was attributed to the drug's suppressive effect on the immune system [ 67 ], the enrichment of DEGs in leishmaniasis suggests potential shared pathogenic signatures. EBV infection is a well-determined risk factor for MS [ 68 ]. While the mechanisms through which EBV favors the onset of MS have yet to be identified [ 69 ], the cellular and viral factors induced by EBV were suggested to potentially stimulate inflammatory and autoimmune responses via interaction with T cells and NK cells [ 70 ].

Based on the number of interactions in the PPI network, the hub genes were IFNG , MS4A1 , CD22 , CXCL8 , CD79A , BLK , CD79B , JUN , and CD69 , which were in the top 10 % of connectivity. IFNG , the gene encoding interferon gamma (IFN-γ), was identified as the central hub gene of the network with 22 interactions. The study of experimental autoimmune encephalomyelitis (EAE) models has revealed that IFN-γ has a stage-specific functionality with contrasting outcomes dependent on the specific stage of the disease. While IFN-γ induces inflammation in the early stages, it suppresses the immune system in the chronic phase [ 71 ]. Therefore, considering IFNG as a therapeutic target depends on the disease stage. CXCL8 was another immune regulator that was marked among the hub genes with a significant upregulation in MS patients compared to healthy controls. By conducting qPCR, a 3-fold upregulation has been observed in the peripheral blood of MS patients when compared with controls [ 72 ]. CXCL8 is expressed by several immune cells, such as T cells, monocytes, and endothelial cells. This chemokine plays a key role in the recruitment of neutrophils, and its upregulation has been observed in astrocytes located in MS lesions. The activation of CXCR2 by CXCL8 contributes to the BBB abnormalities in MS [ 73 ]. On the other hand, multiple markers of B cells, including MS4A1 (CD20), CD22 , CD79A , CD79B , and BLK , were discerned as hub genes of the network. Previous studies have documented that class-switched memory B cells are increased in the blood of MS patients during the relapse stages compared to remission stages [ 74 ]. Given the increased number of CS memory B cells during the active phase of the disease and the higher proportion of naïve B cells in PBMC specimens from MS patients compared to controls in this study, it is plausible that these B cell markers reflect the greater presence of B cells in patients. Furthermore, in the interaction network, CD69 , a conventionally acknowledged marker for T-cell activation expressed shortly after TCR and cytokine stimulation [ 75 ], was detected as a hub gene. This transmembrane protein is expressed in various peripheral blood mononuclear cells, such as B and T lymphocytes, NK cells, and monocytes [ 76 ].

Neurodegeneration in the CNS is influenced by miRNAs, which affect the development of lymphocytes and, in turn, play a key role in MS pathogenesis [ 77 ]. Therefore, key miRNAs were explored by the construction of the miRNA-mRNA regulatory network. Members of the miR-92, miR-181, miR-15, miR-200, and miR-450 families were key regulators of the PPI network by targeting multiple genes.

The four members of the miR-92 family included miR-92a-3p, miR-92b-3p, miR-363-3p, and miR-25-3p [ 78 ]. In a study evaluating the correlation between circulating miRNA levels and neuroimaging results, it was revealed that increased serum levels of miR-92a-3p are associated with the elevated volumes of cervical spine lesions [ 79 ]. It has also been documented that miR-363-3p negatively correlates with T2 lesions in patients with MS, as revealed by MRI [ 80 ]. In a cohort study, miR-25-3p was suggested as an early biomarker associated with the severity of MS [ 81 ].

Four members of the miR-181 family, including miR-181a-5p, miR-181b-5p, miR-181c-5p, and miR-181d-5p, were key regulators of the DEGs. Moreover, these miRNAs targeted CD69 and MS4A1 , which were two of the hub genes in the PPI network ( Supplementary Table S2 ). The current body of literature demonstrates the altered expression of the miR-181 family in neurodegenerative disorders, including MS, Parkinson's disease, and Alzheimer's disease [ 82 ]. Quantitative real-time PCR (qPCR) of PBMCs identified miR-181a-5p and miR-181b-5p as differentially downregulated and upregulated miRNAs, respectively, in MS patients [ 83 , 84 ]. Moreover, targeting genes involved in neuroinflammation, such as MAP2K1 , CREB1 , ATXN1 , and ATXN3 , by miR-181a-5p reflects its therapeutic potential [ 84 ]. The expression of miR-181c-5p is reported to be dysregulated in PBMCs, cerebrospinal fluid (CSF), and serum of relapsing-remitting multiple sclerosis (RRMS) patients [ 85 ]. A study on sporadic and familial forms of MS found miR-181d-5p to be downregulated in PBMCs, as revealed by qPCR analysis [ 86 ].

Members of the miR-15 family identified as regulators were miR-15a-5p, miR-15b-5p, miR-16-5p, and miR-195-5p [ 87 ]. It has been previously demonstrated that miR-15b-5p is a significant upregulated miRNA in the blood of RRMS patients in comparison with healthy individuals [ 88 ]. In a deep learning study, miR-16-5p was suggested to considerably contribute to the pathogenesis of MS [ 89 ]. Similarly, dysregulation of miR-195-5p has been observed in MS, and it was marked as a pathogenic miRNA considering its positive correlation with the severity of MRI clinical data [ 90 ].

Three members of the miR-200 family were discerned as key regulators, including miR-200b-3p, miR-200c-3p, and miR-429 [ 91 ]. The associations of the miR-200 family with multiple neurodegenerative disorders and MS have been previously suggested [ 92 ]. A study on EAE mice demonstrated upregulation of miR-200b-3p in MS brain samples compared to controls. On the contrary, the study reported downregulation for miR-200a-3p and miR-200c-3p [ 93 ].

Assessment of miR-450b-5p serum levels in RRMS and SPMS patients revealed its reduction and suggested it as a potential biomarker for disease progression [ 94 ]. Other key miRNAs, each interacting with multiple mRNAs, included miR-493-5p, miR-124-3p, miR-32-5p, miR-367-3p, and miR-506-3p. Compared to healthy controls, miR-493-5p elevation was observed in the peripheral blood leukocytes of untreated RRMS patients [ 95 ]. Dysregulation of miR-124-3p has been observed in a wide spectrum of neurodegenerative disorders, including MS and Parkinson's disease. It has also been found to have therapeutic potential, as upregulation of miR-124-3p can inhibit the activation of macrophages and microglia, which are involved in the neuroinflammation in MS [ 96 , 97 ]. Evaluation of miR-506-3p expression levels in the peripheral blood of RRMS individuals using qPCR following fingolimod treatment, an immune-regulating therapy approved for these patients, revealed no significant alterations in its levels in MS patients and responders compared to healthy controls and non-responders, respectively [ 98 ].

Assigning 13,237 genes into 15 modules using WGCNA resulted in the identification of a set of co-expressed genes, Module 7, which was negatively correlated with MS (correlation = −0.53, p-value <0.01) and highly enriched in neurological disorders. Functional enrichment analysis of 23 highly connected genes in the PPI network of this module demonstrated that mitochondrial functioning and pathways of neurodegenerative disorders were fully enriched with these genes. Out of 23 highly connected genes in this module, 11 were selected as hub genes, including COX5B , UQCRQ , NDUFS3 , UQCR10 , NDUFA4 , COX7B , NDUFA7 , NDUFS7 , NDUFA2 , NDUFA1 , and NDUFB2 , based on the criteria of module membership >0.8 and gene significance >0.2. These hub genes encode different subunits of electron transport chain complexes, including NADH oxidoreductase (Complex I), ubiquinol-cytochrome c reductase (Complex III), and cytochrome c oxidase (Complex IV), supporting the evidence of defect in mitochondrial functioning in PBMCs of MS patients. Mitochondrial dysfunctions have been widely established as a contributor to MS [ 99 ]. However, at the peripheral level, the reduction of mitochondrial proteins and impairment of this organelle have been reported relatively recently [ 100 ]. According to the current literature on MS, the downregulation of COX5B [ 101 ], NDUFS3 , and NDUFA4 [ 102 ] has been documented in blood samples, while a downward trend in the expression of COX7B [ 103 ], NDUFA4 , NDUFB2 , and UQCRQ [ 104 ] has been reported only in brain lesions. A genome-wide association study (GWAS) reported NDUFA7 as a differentially methylated gene in both immune cells and neurons [ 105 ]. Another association study found NDUFA7 , NDUFS5 , and NDUFS7 genes to be associated with MS [ 106 ].

The WGCNA-based identification of hub genes in the other modules and their functional enrichment also highlighted pathogenic signatures in MS ( Supplementary Fig. S2 ). The hub genes in module 1 were enriched in innate immunity-related terms. While the role of adaptive immunity in the pathogenesis of MS is well-established, emerging research has underscored the contribution of innate immunity in both the onset and progression of the disease [ 107 ]. The hub genes in module 2 were enriched in pathways, a considerable proportion of which were proliferation-related, such as Ras signaling and mitotic-related terms. With a holistic perspective, RAS signaling is involved in autoimmunity, and it has been considered as a therapeutic target for MS [ 108 ]. The hub genes in module 3 were enriched in nucleocytoplasmic transport. Recent findings have revealed dysregulations in this pathway may have a crucial role in the neurodegeneration component of MS [ 109 ]. The hub genes in module 8 reflected dysregulations in homeostasis-related pathways. Notably, protein processing in the endoplasmic reticulum, ATP-dependent chromatin remodeling, and cellular senescence were highlighted by KEGG pathway analysis. On the other hand, the hub genes in modules 9 and 12, which were positively correlated with MS, were enriched in T cell and B cell receptor signaling pathways, respectively.

Among the 59 genes selected from the coordinated results of differential expression analysis and WGCNA ( Fig. 3 ), COPG1 , RPN1 , and KDM3B emerged as prominent candidates based on the converging results of machine-learning methods and ROC analysis ( Fig. 5 , Fig. 6 ). However, the contradictory direction of alterations between the integrated data and external datasets, where their AUC exceeded the 0.7 threshold, refutes their potential as biomarkers for MS ( Fig. 7 ). While their expression demonstrated a downward pattern in the integrated data, they were upregulated in GSE141804 and GSE146383 .

COPG1 is the gamma subunit-encoding gene of the COPI complex, which is involved in vesicle transport from the Golgi to the endoplasmic reticulum (ER) and intra-Golgi traffic [ 110 ]. In a bioinformatics study on MS, COPG1 was ascertained as one of the five hub mRNAs in the circRNA-miRNA-hub mRNA network that met the criterion of having a non-zero coefficient in LASSO regression analysis. Moreover, the expression level of COPG1 in MS was markedly lower than that of the control group [ 111 ]. Similarly, the differential expression analysis of GSE21942 identified COPG1 as one of the most downregulated genes in PBMCs of MS patients [ 112 ], which has also been documented in a recent independent study that utilized the dataset [ 113 ]. In a study on COPA syndrome, which is caused by mutations in the alpha subunit of the COPI complex, applying CRISPR/Cas9 to induce COPG1 deficiency in Th-1 cells resulted in the transcription of inflammatory genes [ 114 ]. This observation suggests a potential link between COPG1 downregulation and inflammation [ 115 ].

RPN1 encodes a regulatory subunit of the 26S proteasome, which is responsible for the ATP-dependent degradation of ubiquitin-marked proteins [ 116 ]. In neurons, defects in this central catalytic component of the ubiquitin-proteasome system (UPS) lead to the accumulation of misfolded proteins and subsequent neurodegeneration in the CNS [ 117 ]. Meanwhile, the 26S proteasome is also involved in oxidative stress, transcription of genes, and the release of neurotransmitters. Moreover, it can indirectly affect T-cell development and migration by regulating dopamine release [ 118 ]. Despite these key regulatory activities of the 26S proteasome, its consideration as a therapeutic target in MS has remained controversial. On the one hand, free proteasomes were detected in the blood specimens of MS patients, and the 26S proteasome was observed to damage the myelin sheath [ 119 ]. On the other hand, 26S proteasome inhibitors result in the exacerbation of neurodegeneration as a consequence of ubiquitinated protein aggregation [ 118 ]. The substantial downregulation of RPN1 in MS has also been underscored in a recent bioinformatic study on MS [ 113 ].

KDM3B encodes lysine demethylase 3B, which is responsible for demethylating histone 3 lysine 9 (H3K9) and is predominantly considered to create a transcriptionally active site [ 120 ]. KDM3B downregulation found by differential expression analysis in the present study is in accord with various reports of hypermethylation in PBMCs of MS patients [ 121 ]. By conducting exome sequencing, pathogenic variants of KDM3B were reported to be associated with a range of cognitive impairments, including behavior problems, attention deficit hyperactivity disorder (ADHD), autism spectrum disorder (ASD), and epilepsy [ 120 ]. Moreover, KDM3B plays a crucial role in the mediation of autophagy through the regulation of autophagy-related genes [ 122 , 123 ]. Autophagy, in turn, is involved in different aspects of MS pathogenesis, notably oxidative stress, myelination, and the activation of immune cells. Meanwhile, it's controversial whether inhibition or amplification of autophagy holds therapeutic potential, as its activation in different cell types is followed by contrasting consequences [ 124 ].

This study encountered several fundamental limitations. Notably, PBMCs are composed of a mixture of cells whose marked disparities in cell proportions between cases and controls could potentially bias the findings. B cell depletion therapies are commonly used in MS treatment [ 125 ]. With this in mind, using the deconvolution of the absolute immune signal method and cell type-specific enrichment of DEGs, the proportion of immune cells was computed in the cases and controls. On the contrary, the proportion of B cells in MS patients was evaluated to be higher than that in healthy individuals, suggesting no obvious effect of the therapy on the results. Secondly, different diagnostic criteria for MS were applied in the investigations, and it was not determined whether the patients in GSE21942 are diagnosed with RRMS or SPMS. Furthermore, all patients in GSE21942 were female, and considering the importance of gender in autoimmune disorders, this factor may influence the results and partially explain the lack of generalizability of the identified biomarkers. There were also four samples in the study that received treatment at the time of sample collection, which could not be excluded as they had not been marked in the sample information that was retrievable by the GEOquery package. Furthermore, the integration of datasets limits the analysis to the genes whose expressions have been commonly recorded across different datasets, potentially concealing data on genes with diagnostic value.

Conclusion

Herein, an integrated transcriptomic approach was employed to make an effort to identify potential non-invasive biomarkers for MS. COPG1 , RPN1 , and KDM3B were initially underscored by machine-based learning methods and exhibited acceptable diagnostic efficacy in the integrated data. However, examination of their expression patterns in external validation datasets did not confirm their diagnostic potential. The failure to identify biomarkers in this study reflects the importance of considering factors such as type of MS, pretreatment conditions, demographic information, and disease stage, as well as employing a large sample size. Furthermore, considerable pathogenic signatures were identified. Notably, the detection of a gene module negatively correlated with MS and enriched in neurodegenerative disorders led to the identification of core genes within this module that encode respiratory chain complexes. This finding, consistent with recent studies [ 100 ], further underscores the dysregulation of mitochondria in circulating immune cells in MS. Although this study could not identify promising biomarkers for MS, further research on blood transcriptomes, a valuable non-invasive source that has gained increasing interest in recent years, is needed to achieve this objective.

Ethics approval and consent to participate

Not applicable.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

The author appreciates the contributions made by the authors of the GSE203241 , GSE21942 , GSE247181 , GSE146383 , GSE141804 , GSE59085 , and GSE17393 .

Footnotes

Supplementary data to this article can be found online at https://doi.org/10.1016/j.bbrep.2025.102022 .

Supplementary data

The following are the Supplementary data to this article:

Multimedia component 1

Multimedia component 2

Multimedia component 3

Data availability

The datasets analyzed during the current study are available in the Gene Expression Omnibus (GEO) database ( https://www.ncbi.nlm.nih.gov/geo/ ).

References

Associated Data

Supplementary Materials

Multimedia component 1

Multimedia component 2

Multimedia component 3

Data Availability Statement

The datasets analyzed during the current study are available in the Gene Expression Omnibus (GEO) database ( https://www.ncbi.nlm.nih.gov/geo/ ).