PBMCs were collected from 10 patients with NMOSD at Osaka University Hospital who were diagnosed according to the 2015 NMOSD criteria (20). The protocol of this study involving human participants was reviewed and approved by the Ethics Committee of Osaka University (11298–36). All patients were positive for anti-AQP4 antibodies. Immunosuppressive drugs other than steroids were limited to azathioprine or the calcineurin inhibitor tacrolimus. Patients using biologics were excluded. The blood samples from a relapsed patient included in this study were taken at the timepoints of 570 days and 506 days before the relapse in the remission phase, and 50 days and 79 days after the relapse in the acute phase respectively. For the control samples, PBMCs derived from three patients positive for anti-acetylcholine receptor antibody, the hallmark of another autoantibody-mediated disease coined as myasthenia gravis, were collected. None of the control patients were taking steroids at the time of sample collection. Detailed information on the cohort is depicted in Table 1 for patients with NMOSD, and in Table 2 for the control samples. PBMCs were isolated with Ficoll-Paque™ PLUS (Cytiva Sweden AB, Uppsala, Sweden), then stored in a freezer at −80°C.

Total RNA was extracted using the RNeasy Micro Kit (Qiagen, Valencia, CA, USA). Full-length cDNA was generated using a SMART-Seq HT Kit (Takara Bio, Kusatsu, Japan). An Illumina library was prepared using a NexteraXT DNA Library Preparation Kit (Illumina, San Diego, CA, USA) according to SMARTer kit instructions. Sequencing was performed on an Illumina NovaSeq 6000 sequencer (Illumina, San Diego, CA, USA) in the 100-base single-end mode. After adapter trimming using Trimmomatic, sequenced reads were mapped to the human reference genome sequence hg19 using TopHat version 2.1.1. The fragments per kilobase of exons per million mapped fragments (FPKM) value was calculated using Cuffnorm version 2.2.1. The raw data have been deposited in the NCBI’s Gene Expression Omnibus database (GEO GSE271586).

Principal component analysis (PCA) and K-means analysis of bulk RNA-seq data were performed using iDEP (an integrated web application for differential expression and pathway analysis of RNA-seq data) (21). For the analysis, iDEP v2.0 was used, and parameters were set at the defaults. We used the top 2000 genes and set the number of clusters to four.

For pathway analysis, we used the Enricher application, and mapped out the pathways for each cluster. Either “Reactome 2022” or “WikiPathway 2023 Human” was used to calculate the data.

A transformed matrix from normalized counts was used for WGCNA (v1.72.5) analysis using R 4.3.2 to identify significant modules that correlated with steroid dosage. Then, we calculated the adjacency using the adjacency function with power = 16. The adjacency matrix was then converted into a Topological Overlap Matrix (TOM) using TOMsimilarity, and a gene tree was calculated using hierarchical clustering against 1 - TOM with the “average” method. A dynamic tree cut was conducted with the following parameters; deepSplit = 0, pamRespectsDendro = FALSE, minModuleSize = 30. The eigengenes of each module were calculated using the moduleEigengenes function, and modules with strongly correlated eigengenes were merged, with MEDissThres set to 0.2. These eigengenes were used for the correlation with clinical information indicating the steroid dosage associated with each sample.

To detect DEGs between relapse and remission samples, we used the edgeR package (v4.0.15) in R. We set the thresholds for significant differential expression at log₂ fold change > 1 and a false discovery rate < 0.01.

Using the four clusters identified by K-means in Patient 10, we created a chronological landscape of gene expression levels for each patient other than Patient 10, on the basis of cluster 2 (downregulated during relapse) and cluster 3 (upregulated during relapse). For each patient, we converted the FPKM of each gene included in each cluster into base-2 logarithms. Using these values, we calculated the relative expression change value, which was determined by subtracting the mean log value across samples from the original log value. On the basis of the relative expression change value, we calculated the mean and standard deviation (SD) for each sample.

To analyze the gene expression profile of each cluster, we used the Seurat package (v5.0.2) (22). A multimodal dataset of PBMCs, which was derived from single-cell RNA sequencing (scRNA-seq) combined with other modalities, such as CITE-seq, was loaded, followed by dimensionality reduction on the entire dataset (23). PCA generated a low-dimensional representation of the data, and uniform manifold approximation and projection (UMAP) was used for visualization. Gene module scores were calculated using the AddModuleScore function to evaluate the relative expression of specific gene sets in each cell, adjusted by subtracting the aggregated expression of a control feature set. In addition, weighted nearest neighbor (WNN) analysis was applied to integrate information from multiple modalities and capture cell-type specific patterns more accurately. The integrated WNN-based representation was visualized using UMAP, and the result was referred to as ‘W UMAP’ to denote the use of the WNN algorithm in the dimensionality reduction process.

To separately isolate monocytes, T cells and B cells, frozen PBMCs were thawed and stained with antibodies against surface antigens, after the treatment with FcR blocking reagent (Miltenyi Biotec, Tokyo, Japan). The antibodies used were as follows: CD4-Brilliant Violet 605 (BioLegend), CD20-PE/Cyanine7 (BioLegend), and CD3-Alexa Fluor 488 (BD Biosciences), CD14-Pacific Blue (BD Biosciences), CD19-APC-Cy7 (BD Biosciences). The samples were analyzed on a FACSAria Fusion (BD Biosciences), and monocytes, T cells and B cells were defined as CD3-CD19-CD20-CD14+ cells, CD3+CD4+ cells and CD19+CD20+ cells, respectively.

For the preparation of qPCR analysis, samples were dissolved in ISOGEN II (NIPPON GENE, Tokyo, Japan) kit. The total RNA was reverse transcribed with SuperScript IV VILO Master Mix (Thermo Fisher Scientific, Waltham, MA, USA) according to the manufacturer’s instructions. The qPCR analysis was performed by mixing the TaqMan Fast Advanced Master Mix (Thermo Fisher Scientific, Waltham, MA, USA), TaqMan Gene Expression Assay (Thermo Fisher Scientific, Waltham, MA, USA) and cDNA of each sample. The following TaqMan Gene Expression Assays were used: CX3CR1 (Assay ID: Hs01922583_s1), CCR2 (Assay ID: Hs00704702_s1), CCR5 (Assay ID: Hs99999149_s1), HBA2, HBA1 (Assay ID: Hs00361191_g1), CD180 (Assay ID: Hs00194403_m1), NFKBIA (Assay ID: Hs00355671_g1), DUSP1 (Assay ID: Hs00610256_g1), THBS1 (Assay ID: Hs00962908_m1), TNFAIP3 (Assay ID: Hs00234713_m1), NFIL3 (Assay ID: Hs00705412_s1), CXCL8 (Assay ID: Hs00174103_m1), EGR1 (Assay ID: Hs00152928_m1), SERPINB2 (Assay ID: Hs01010736_m1), and GAPDH (Assay ID: Hs02786624). The results were analyzed with QuantStudio 7 Flex Real-Time PCR System (Applied Biosystems, Waltham, MA, USA). To confirm the significant differences among each sample, Tukey’s multiple comparisons test was applied.

To determine the proportions of immune cell subsets of bulk RNA-seq data, we employed the CIBERSORTx deconvolution algorithm (https://cibersortx.stanford.edu/) (24). We configured the submitted job type to “Impute Cell Fractions”. The LM22 Signature Matrix file was used as a reference. The relative proportions of immune cells in each dataset were calculated using 1000 permutations, applying B mode batch correction.

To determine whether clinical relapse would affect the gene signature of PBMCs in NMOSD, we compared transcriptome data at several timepoints in a patient in our cohort who underwent relapse during steroid tapering.

On the third timepoint of four samples collected, the patient showed relapse and received immunoadsorption therapy. When PCA was performed, the samples obtained at the time of remission and after relapse (in the latter case, both immediately after relapse and after acute-phase treatment) were separated by PC1, indicating 54.1% variance ( Figure 1A ).

Samples obtained during the remission and relapse phases were further analyzed by K-means clustering. Each of the four clusters showed distinct changes in gene signatures before and after relapse. In general, genes comprising clusters 1 and 3 showed enhanced expression at relapse, while those comprising clusters 2 and 4 showed decreased expression ( Figure 1B ).

To further elucidate the types of immune cells in which the aforementioned gene signature changes took place, each cluster delineated using the K-means method was applied on single-cell resource of PBMCs ( Figure 1C ). While clusters 1 and 3 did not show specific distributions, both clusters 2 and 4 were predominantly observed in monocyte clusters ( Figure 1C ). In DEG analysis using samples in the relapse and remission phases, immune-related genes such as CX3CR1 and CCR2 were upregulated in relapse phases ( Figure 1D ). To validate these results of DEG analysis, qPCR was performed for five upregulated and downregulated genes in the same samples as those used for RNA seq, which were separately aliquoted from the same PBMCs. As a result, we were able to confirm that gene expression patterns were consistent with those observed with DEG analysis ( Figure 2 ). It is also important that cluster 2, which was downregulated at relapse, was characterized by both interleukin 10 (IL-10) signaling and glucocorticoid receptor pathways ( Figure 3 ). Moreover, cluster 3, which was characterized by enhanced gene expression at relapse, showed an association with the interferon pathway ( Figure 3 ). Similar results were observed for cluster 4, which was also downregulated at relapse (data not shown).

To further comprehend the sample variance depending on the status of steroid tapering, all the samples collected from 9 individual patients at different steroid dosage were analyzed for PCA analysis. The PC1 and PC2 variances were 20.4% and 14.1%, respectively, and the samples tended to be grouped according to patient identity rather than steroid dosage ( Figure 4A ).

To determine how gene signatures were related to steroid dosage, transcriptome data of samples collected from each patient were separately analyzed by the K-means method. Although the cohort as a whole did not exhibit a consistent pattern, three clustering types could be defined on the basis of the gene clustering patterns in monocytes. Specifically, in Patients 1, 2, and 5 (Group 1), the gene sets reduced by steroid tapering were highly expressed in monocytes ( Figure 4B ). By contrast, in Patients 7 and 9 (Group 2), those upregulated by steroid tapering were also highly expressed in monocytes ( Figure 4B ). Other clusters were randomly distributed across non-specific cell populations ( Supplementary Figure S1 ). Of note, genes with upregulated expression in Group 2 were associated with inflammatory pathways, especially interferon signaling, and those with downregulated expression were enriched in anti-inflammatory and glucocorticoid receptor pathways, as represented by IL-10 signaling ( Supplementary Figures S2 , S3 ). These varied observations among patients were not unexpected given their differing immunological statuses. The important thing is that patients who showed similar downregulation of genes associated with anti-inflammatory and glucocorticoid receptor pathways might have similar transcriptomic changes in the relapse phase, and represent the preconditioned status of NMOSD that is at the edge of reactivation while steroid is being tapered.

To validate whether the genes in the clusters obtained by the K-means method were actually expressed predominantly in a specific cell subtype, we have separately isolated CD14+ monocytes, CD4+ CD3+ T cells and CD20+ CD19+ B cells by flow cytometry and performed qPCR analysis to compare the gene expression of the interest. Although due to the limited amounts of residual samples which were aliquoted from the same PBMCs as used for RNA-seq analysis we were not able to perform the validations for the whole results, we were able to confirm that genes included in the cluster 2 of Patient 5 which were assigned to monocytes were significantly highly expressed in monocytes by qPCR analysis ( Figure 5 ).

Moreover, by applying deconvolution techniques of CIBERSORT we confirmed specific cell subtypes which were identified by gene set scoring in single-cell RNA-seq dataset. The deconvolution results of cluster 4 of patient 9, identified by K-means method, were shown to be as follows: the percentages of monocyte, macrophage, and dendritic cell populations were more predominantly expanded in the sample such as Patient9_5T, where steroid dosage was most reduced ( Supplementary Figure S4 ).

In addition, as control cohort, we utilized bulk RNA-seq data of PBMCs derived from three patients positive for anti-acetylcholine receptor antibody. K-means clustering and gene set scoring in single-cell RNA-seq dataset were performed following the same processes as those for NMOSD samples ( Supplementary Figure S5 ). In control samples, the genes in each cluster were predominantly distributed either to B cells or to NK cells, which pattern was remarkably distinct from those of NMOSD samples, where monocyte distribution has been their characteristic results.

Since a relapsed patient exhibited distinct gene clusters showing either upregulated or downregulated gene expression ( Figure 1B ), we sought to clarify the chronological landscape of each cluster in patients undergoing steroid tapering. Certain trend was observed with cluster 2, which showed downregulation at relapse, in patient 1, 2, 5 and 6 ( Supplementary Figure S6 ). By contrast, in Patients 5 and 6, cluster 3 exhibited enhanced expression at relapse, and this expression was positively correlated with steroid tapering ( Supplementary Figure S7 ). These observations suggest that gene signatures observed during flares of disease activity might also be present in the PBMCs of patients undergoing steroid tapering.

As gene expression was not significantly consistent between patients undergoing steroid tapering, as mentioned above, we used WGCNA to identify unbiased gene modules related to steroid dosing. WGCNA revealed 23 modules among the nine patients other than the one who experienced relapse (Patient 10) ( Figure 6A ).

Only the royalblue module was significantly correlated with steroid dosage, with a correlation coefficient of 0.4 ( Figure 6B ). The correlation coefficient value was negative, indicating that the expression levels of genes in the royalblue module increase when the steroid dosage decreases.

After identifying this steroid-correlated module, we examined how its associated genes changed in each patient as the steroid dosage was reduced ( Figure 7A ). Importantly, genes in the royalblue module showed remarkable enhancement at relapse in Patient 10. Given that this module was identified in the group of patients excluding Patient 10, it is interesting that the gene signature associated with steroid tapering was apparent in this case of relapse.

As for the other nine patients undergoing steroid tapering, Patients 5, 6, 7, and 9 showed enhanced expression of royalblue module genes. Royalblue module genes were expressed predominantly in the monocyte subset, although this preferential expression was rather subtle ( Figure 7B ).

To confirm specific cell subtypes which were identified by gene set scoring in single-cell RNA-seq dataset, we have validated the results by applying deconvolution techniques of CIBERSORT. The results of the deconvolution for WGCNA royalblue module were assigned to monocyte, macrophage, and dendritic cell populations. For example, in patient 5, the percentage of monocytes, macrophages and dendritic cells was shown to be more abundant in Patient5_4T sample, where steroid dosage was most reduced, compared to samples collected at other timepoints ( Supplementary Figure S8 ).

Importantly, pathway analysis showed that royalblue module genes were associated with the interferon pathway ( Figure 7C ). In the royalblue module, we identified genes important for the inflammatory response, including interferon-related genes and those involved in T cell and B cell differentiation ( Table 3 ).

All of these observations suggest the presence of a steroid-related gene module in PBMCs, and indicate that the innate pathway, as represented by the interferon signature, might represent fluctuating NMOSD activity during steroid tapering.