Characteristics of participants and study design

This study enrolled 46 HC (13 males, 33 females) and 52 patients diagnosed with NMOSD (9 males, 43 females). The mean age was 42 ± 2.42 years in the HC group and 47.8 ± 2.11 years in the NMOSD group, with no significant age difference (P = 0.07) or sex difference (P = 0.23) between the groups. In the NMOSD cohort, the EDSS was 3.16 ± 0.34 and the MSSS was 5.38 ± 0.44. The mean duration of disease was 42.88 ± 9.74 months. Out of the NMOSD group, 31 patients exhibited an AQP4 antibody titer exceeding 100:1, whereas 21 patients had a titer below 100:1, including 13 who tested negative for the antibody. All NMOSD patients were negative for MOG antibodies. Detailed participant details information can be found in Table S1. For conducting a multi-omics analysis, blood samples were collected from both groups to investigate alterations in immune cell phenotypes, plasma cytokines, and metabolites in the NMOSD, analyzed using bioinformatics tools as depicted in Fig. 1A.

Fig. 1

Study design and changes of blood immune cell subsets in NMOSD. A Study design. B-M Box plots representing expression values of differential expressed blood immune cell subsets in HC (n = 46) and NMOSD (n = 48), with the central line representing the median, the box denoting the interquartile range (IQR), and the whiskers extending to the furthest data point within 1.5 times the IQR from the box. FDR corrected q-value: * P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.001. mDC, myeloid dendritic cell; pDC, plasmacytoid dendritic cell

Alterations of blood immune cell subtypes in NMOSD

To explore the influence of peripheral immune cells on the progression of NMOSD, we assessed thirty immune cell subtypes using mass cytometry (Fig. S1A). We identified many significantly changed immune cell subtypes in NMOSD patients as compared to HC group. Specifically, the NMOSD group exhibited an increased proportion of monocytes (Fig. 1B), along with decreased proportion of regulatory T cells (Treg) and its subtypes—activated Treg and secreting Treg, when compared to the HC group (Fig. 1C–E). Additionally, we observed reductions in dendritic cells (DC), encompassing both plasmacytoid DC (pDC) and myeloid DC (mDC), as well as natural killer (NK) cells, T cells, CD4+T cells, CD28+T cells, CD28+CD4T cells, CD28+CD8T cells, and Th17 cells (Fig. 1F–M, Fig. S1B-C). Additionally, we investigated potential correlations between blood immune cell subtypes and the MSSS; however, our findings indicated no significant correlations.

Alterations of plasma cytokines in NMOSD

Plasma cytokines play a pivotal role in modulating peripheral inflammation and may interact with peripheral immune cells, influencing NMOSD pathogenesis. To gain a comprehensive understanding of the role of plasma cytokines in NMOSD pathogenesis, we examined the expression profiles of 440 cytokines using a cytokine array. We found significant elevations in several cytokines in the NMOSD group compared to the HC group, including hepatocyte growth factor activator inhibitor type 2 (HAI-2), S100 calcium binding protein A8 (S100A8), growth differentiation factor 15 (GDF15), and Prostasin. In contrast, several cytokines were significantly decreased in NMOSD, such as stem cell factor (SCF), neurotrophin-3 (NT-3), bone morphogenetic protein 2 (BMP-2), galectin-2, and beta-nerve growth factor (b-NGF) (Fig. 2A).

Fig. 2

Changes of plasma cytokines in NMOSD and correlation with MSSS. A Bar plot of significantly differential expressed cytokines between HC (n = 46) and NMOSD (n = 46), determined by linear regression analysis controlling for age and sex. red: cytokines increased in NMOSD, blue: cytokines decreased in NMOSD. B Scatter plot showing the top 6 cytokines that have a negative correlation with the MSSS, as determined by spearman analysis, n = 45. C Scatter plot showing the top 6 cytokines that have a positive correlation with the MSSS, as determined by spearman analysis

Furthermore, we investigated correlations between plasma cytokine levels and the MSSS in NMOSD patients. Interestingly, our results indicated a higher number of cytokines negatively correlated with MSSS than those positively correlated (Table S2). Specifically, cytokines such as osteoactivin, alpha-fetoprotein (AFP), CD99, FAP, surfactant protein D (SP-D), and brain-derived neurotrophic factor (BDNF) displayed strong negative correlations with MSSS (Fig. 2B). On the contrary, angiogenin, insulin-like growth factor-binding protein 6 (IGFBP-6), C–C motif chemokine 5 (RANTES), thrombospondin-1 (TSP-1), chitinase-3-like protein 1 (CHI3L1), and dickkopf-3 (Dkk-3) demonstrated the most prominent positive correlations with MSSS (Fig. 2C).

Alterations of plasma metabolites in NMOSD

To gain insights into the alteration of plasma metabolites in NMOSD, we analyzed 630 plasma metabolites using targeted metabolomics. We employed liquid chromatography-mass spectrometry (LC–MS) for small molecular metabolites and flow injection analysis (FIA) for larger molecules such as cholesteryl esters, glycerophospholipids, glycerol esters, sphingolipids, and hexoses. The distinct profiles observed in NMOSD versus the HC group were notably separated in the LC–MS mode and in the PC2 dimension of FIA mode (Fig. S2A), indicating significant alterations in plasma metabolites in NMOSD.

Further analysis of differential metabolites revealed increased levels of small molecule metabolites, including bile acids (e.g., tauroursodeoxycholic acid [TDCA], taurocholic acid [TCA], glycodeoxycholic acid [GDCA], glycocholic acid [GCA]), hypoxanthine, lactate (Lac), sarcosine, p-Cresol sulfate (p-Cresol-SO4), α-aminobutyric acid (AABA), γ-aminobutyric acid (GABA), cysteine (Cys), and homocysteine (HCys) in NMOSD patients. Conversely, decreases were observed in dehydroepiandrosterone (DHEAS), homoarginine (HArg), octadecadienoic acid (FA[18:2]), serotonin, and docosahexaenoate (DHA) (Fig. 3A).

Fig. 3

Changes of plasma metabolites in NMOSD and correlation with disease severity. A Bar plot illustrates the log2 fold change of LC–MS metabolites with significant differential expression between HC (n = 46) and NMOSD patients (n = 47), identified through statistical analysis involving linear regression models and FDR correction for multiple comparisons. B Box plots representing expression values of differential expressed sphingolipids in HC and NMOSD, with the median and interquartile range shown. FDR corrected q-value: **P < 0.01, ****P < 0.001. C Bar plot displayed enriched pathways of significantly differential expressed LC–MS and FIA metabolites using MetaboAnalyst enrichment analysis. D Scatter plots showing small molecular metabolites that have a negative correlation with the MSSS, as determined by spearman analysis

Lipid-related alterations included an increase in triglycerides (TGs), whereas levels of phosphatidylcholines (PCs), cholesterol esters (CEs), ceramides (Cer), and hexosylceramides (HexCers) decreased (Fig. S2B). Notably, various sphingomyelin (SM) subclasses consistently decreased in NMOSD, including SM.C16:0, SM.C16:1, SM.C18:0, SM.C18:1, SM.C20:2, SM.C24:0, SM.C26:0, SM(OH)C14:1, SM(OH)C16:1, SM(OH)C22:1, SM(OH)C22:2, and SM(OH)C24:1 (Fig. 3B), highlighting these metabolites' reduction in plasma as a significant NMOSD feature, potentially useful as diagnostic markers.

Pathway enrichment analysis of differential metabolites pointed to significant enrichment in numerous pathways, including aminoacyl-tRNA biosynthesis, arginine biosynthesis, taurine and hypotaurine metabolism, cysteine and methionine metabolism, primary bile acid biosynthesis, starch and sucrose metabolism, and sphingolipid metabolism (Fig. 3C). These results implicated the contribution of these altered metabolic pathways to NMOSD pathogenesis. Furthermore, correlation analysis between differential metabolites and the MSSS showed taurine, hippuric acid (HipAcid), creatinine, glutamate (Glu), ornithine (Orn), and Lac levels negatively correlated with MSSS (Fig. 3D, Table S3). Additionally, negative correlations with MSSS were also noted for lipid metabolites, including PCs, TGs, and various sphingolipid subclasses (Table S4), suggesting their potential as biomarkers for monitoring NMOSD progression.

Associations among immune cell subsets, cytokines, and metabolites

To elucidate the potential interactions among differential immune cell subsets, cytokines, and metabolites in relationship to the development of NMOSD, we conducted Spearman correlation analysis. This analysis identified significant correlations: cytokines such as HAI-2, S100A8, and GDF15 were positively correlated with monocytes, yet exhibited negative correlations with Tregs, DCs, NK cells, and various T cell subtypes. In contrast, cytokines including IL-3, DNAX accessory molecule 1 (DNAM-1), and Cadherin-11 displayed reverse correlation patterns (Fig. S3A). Furthermore, the analysis extended to correlations between immune cell subtypes and metabolites, uncovering that DHEAS, FA(18:2), and HArg negatively correlated with monocytes and positively with Tregs, DCs, NKs, and various T cell subtypes. Conversely, metabolites such as hypoxantine, xanthine, and methionine sulfoxide (Met-SO) showed positive correlations with monocytes, but negative correlations with Tregs, DCs, NKs, and various T cell subtypes (Fig. S3B). Sphingolipids (e.g., SM, Cer, HexCer) and PCs displayed negative correlations with monocytes and positive correlations with DCs, NKs, and T cell subsets (Th2, Th17, CD28+T, CD28+CD4T, CD28+CD8T), whereas TGs exhibited the opposite correlation pattern (Fig. S3C). These findings suggest intricate functional interactions between varying immune cell subsets and specific cytokines or metabolites.

Further, we used the xMWAS integrative analysis to contrast the multi-omics interactions between NMOSD and HC groups. In the HC group, the analysis selected 9 immune cell subtypes, 11 cytokines, and 14 metabolites, which were organized into seven communities with robust inter-correlations within the integrative network (γ ≥ 0.7, P < 0.05) (Fig. 4A and Table S5). In contrast, the NMOSD group analysis identified 11 immune cell subtypes, 54 cytokines, and 14 metabolites, forming ten communities with similarly high correlations (Fig. 4B and Table S6). Notably, this analysis underscored a significant increase in the number of plasma cytokines associated with NMOSD, highlighting the critical role of plasma cytokines in NMOSD pathogenesis. Moreover, using eigenvector delta centrality measure (DCM), we found that secreting Treg, Treg, TDCA, resting Treg, glycochenodeoxycholic acid (GCDCA), taurochenodeoxycholic acid (TCDCA), hypoxanthine, and GABA showed the largest change in eigenvector centrality between NMOSD and HC (Fig. 4C). This analysis underscores the impactful roles of Treg subtypes and secondary bile acids in NMOSD pathogenesis, demonstrating that immune cell subtypes, cytokines, and metabolite communities may cooperatively contribute to disease development.

Fig. 4

Integrative analysis of blood multi-omics data. A Integrative network plot of cytokines, immune cells, LC–MS and FIA metabolites in HC group (n = 46) by using xMWAS, correlation higher than 0.7, one color representing one community. B Integrative network plot of cytokines, immune cells, LC–MS and FIA metabolites in NMOSD group (n = 34) by using xMWAS, correlation higher than 0.7, one color representing one community. C The identification of nodes that undergo network changes, which is determined based on the Delta centrality (importance) measure (ECMnmo—ECMcontrol) in which EMC represents eigenvector centrality measure

Integrative features discriminating NMOSD from HC

The DIABLO analysis aimed to discern a multi-omics signature that could effectively distinguish NMOSD from HC. By integrating data on immune cells, metabolites, and cytokines, we identified distinct signatures within each category that distinguish NMOSD from HC (Fig. 5A–C). Notably, immune cell subtype analysis revealed higher proportions of monocytes and effector CD8T cells (TE CD8) in NMOSD, whereas HC exhibited higher proportions of CD28+T cells, CD28+CD4T cells, and activated Treg among others (Fig. 5D). In terms of metabolites, NMOSD was characterized by higher levels of GABA, sarcosine, and AABA in NMOSD, contrasting with higher levels of FA (18:2), HArg, DHEAS, SM(OH)C22:1, SM(OH)C24:1, SM(OH)C22:2, and SM.C24:0 in HC (Fig. 5E). For cytokine expression, NMOSD was marked by higher levels of HAI-2, GDF15, p53-associated parkin-like cytoplasmic protein (PARC), S100A8, and T-cell immunoglobulin mucin receptor 3 (TIM-3), while HC had higher levels of SCF, Syndecan-1, and retinol-binding protein 4 (RBP4), among others (Fig. 5F).

Fig. 5

Integrative analysis revealed discriminative features between NMOSD and HC. Omics integration analysis among immune cells, cytokines and metabolites. The highly correlated multi-omics signatures that discriminate HC (n = 46) and NMOSD (n = 34) are identified from supervised model using DIABLO. A–C scatterplot of samples on first 2 components for each block (Omics), colored by group (blue: HC, orange: NMOSD). D–F Loading plot of component 1 from supervised model for each block (omics), important signatures resulted from DIABLO are ordered by absolute importance (x-axis), color indicate the class for which the median expression value is the highest for each feature. G Omics integration analysis among immune cells, cytokines and metabolites. Heatmap of picked important signatures (of component 1 and 2 from each block) with dendrogram computed by Euclidean distance and complete linkage method, high correlations between features (γ > 0.5)

Further analysis involved heat-map clustering to visualize the correlations within the multi-omics signatures differentiating NMOSD from HC. This revealed two primary clusters of highly correlated communities (γ > 0.5) encompassing immune cell subtypes, metabolites, and cytokines (Fig. 5G). Notably, an elevated subcommunity in NMOSD, featuring elements like GABA, choline, Lac, AABA, S100A8, sarcosine, HAI-2, monocytes, TIM-3, PARC, and GDF15, demonstrated strong inter-correlations. Conversely, a subcommunity marked by decreased levels in NMOSD, essential for its differentiation from HC, included RBP4, Th17 cells, Th2 cells, CD28+T cells, CD28+CD4T cells, CD4T cells, Treg, and activated Treg, displaying close connections among its members. This intricate multi-omics DIABLO analysis provides a deep insight into NMOSD's distinguishing biomarker interactions, illuminating the intricate and dynamic interplay essential for understanding this disease.

Differences in peripheral characteristics among clinical subtypes of NMOSD

To elucidate the clinical heterogeneity of NMOSD, we stratified patients based on clinical parameters, including the EDSS, AQP4 antibody titer, disease duration, number of relapses in the past year, and episode history (Table S1). In our analysis, NMOSD patients were dichotomized into high and low EDSS score groups (EDSS-H, EDSS ≥ 4; EDSS-L, EDSS < 4), high and low AQP4 antibody titers (AQP4-H, antibody titer ≥ 100:1; AQP4-L, antibody titer < 100:1), by disease duration into T1 (< 3 years) and T2 (> 3 years), by relapse frequency into low relapse (LR, no relapses in the past year) and high relapse (HR, ≥ 1 relapse in the past year), and finally into first episode versus recurrent groups.

Firstly, the comparison between the EDSS-H and EDSS-L groups revealed that the former exhibited elevated levels of macrophage inflammatory protein 3a (MIP-3a), GRO (CXCL1), and matrix metalloproteinases (MMP-3, MMP-2), and alongside lower levels of IL-10, IFN-γ, and IL-7, etc. (Fig. 6A). Notably, there were no discernible differences in blood immune cell subtypes and metabolites between these subgroups.

Fig. 6

Difference in peripheral features between different clinical subgroups. Bar plot of significantly differential expressed cytokines (P < 0.05) between EDSS-H (n = 16) and EDSS-L (n = 24) (A), between AQP4-L (n = 17) and AQP4-H (n = 23) (B, between T2 (n = 16) and T1 (n = 24) (C), between HR (n = 24) and LR (n = 16) (D), between Recurrence (n = 18) and First episode (n = 22) (E). F The serum level of FAP in first episode (n = 20) and recurrence (n = 20) of NMOSD, Student’s t-test, *P < 0.05. G The serum level of FAP in T1 (< 3 years) (n = 27) and T2 (≥ 3 years) (n = 13) of NMOSD, Student’s t-test, *P < 0.05

Next, the analysis of the AQP4-L versus AQP4-H groups showed lower levels of secreting Treg in AQP4-L (Table S7), alongside significant lower levels in bone morphogenetic protein 4 (BMP-4), BDNF, and neuron-specific enolase (NSE), etc. (Fig. 6B). Metabolomic analysis showed lower levels of ceramides and higher levels of asparagine (Asn) and serine (Ser) in the AQP4-L group (Table S7).

Comparing T2 and T1 subgroups, the former displayed an increased presence of CD28+CD8T cells, osteopontin (OPN), and vascular cell adhesion molecule 1 (VCAM1), alongside reduced CA15-3, MIP-3b, neutrophil-activating protein 78 (ENA-78), FAP, and CCN family member 3 (NOV) (Fig. 6C, Table S8). Metabolic differences included elevated glutamine (Gln) with decreased phenylalanine (Phe) and Glu in T2 (Table S8).

In the HR versus LR comparison, the HR group had increased levels of resting Treg, endoglin, carcinoembryonic antigen-related cell adhesion molecule 1 (CEACAM-1), adiponectin, renin, and matrix extracellular phosphoglycoprotein (MEPE), but lower levels of NOV, FAP, Lac, and GLCAS (Fig. 6D, Table S9).

In the comparison between recurrent and first episodes, recurrent episode was characterized by higher levels of Th17 cells (Table S10), vascular endothelial growth factor receptor 2 (VEGF-R2), VCAM1. In contrast, markers such as CA15-3, kallikrein-5, ENA-78, and FAP were lower (Fig. 6E). Besides, Phe was notably lower in recurrent episodes (Table S10).

We examined the distinctions between the relapse and remission phases. Our findings indicate that the relapse phase is associated with higher levels of memory B cells and monocytes, and lower levels of CD28+CD8T cells and NK cells compared to the remission phase (Fig. S4A). Furthermore, during the relapse, there are elevated levels of SP-D and reduced levels of OPN (Fig. S4B), as well as increased AABA and decreased citrulline (Cit) levels (Fig. S4C). Besides, considering that medication might influence peripheral factors, we also assessed the effects of different treatments, including glucocorticoids, combinations of glucocorticoids with intravenous immunoglobulin (IVIG), and glucocorticoids with immunosuppressants. The analysis indicated no significant changes in the proportions of immune cell subtypes across the groups. Similarly, most cytokines and metabolites did not differ significantly among the various treatment groups. However, TSP-1 and B7 homolog 3 (B7-H3) exhibited changes with glucocorticoids combined with intravenous immunoglobulin (IVIG) or with immunosuppressants, while p-Cresol-SO4 levels varied between the use of glucocorticoids alone and glucocorticoids paired with immunosuppressants (Table S11). This suggests that these treatments have a minimal impact on peripheral factors.

Notably, FAP was differentially expressed in multiple clinical subgroups of NMOSD, such as lower expression in EDSS-H compared to EDSS-L (Fig. 6A), in T2 compared to T1 (Fig. 6C), in HR compared to LR (Fig. 6D), and in the recurrence group compared to the first episode (Fig. 6E), suggesting its potential as a key marker for NMOSD progression. To investigate the link between FAP and disease progression, we collected serum samples from an independent cohort of 40 NMOSD patients, equally divided between 20 with first episode and 20 with recurrent cases (Table S12). We used ELISA method to measure the serum FAP concentrations and compared the levels between the first episode and recurrent cases. Our findings showed that patients with recurrent cases had significantly lower serum FAP concentrations than those experiencing their first episode (Fig. 6F). Furthermore, we found that patients with a long disease duration (T2) had lower serum FAP levels compared to those with a short duration (T1) (Fig. 6G). These results collectively underscore the relevance of FAP as a crucial marker for the progression of NMOSD.