(1.1) Inclusion criteria: 1) aged 18–75 years old; 2) meeting the diagnostic criteria for NMOSD (Wingerchuk et al., 2015); 3) within 30 days of onset in the acute phase and not receiving relevant drug treatment; and 4) provision of signed informed consent.

Total RNA was extracted from whole blood stored in Paxgene blood RNA tubes by using a Paxgene Blood RNA Extraction Kit (PreAnalytiX, Qiagen). A Qubit fluorometer was used to quantify RNA at 260/280 nm, and an Agilent 2,100 Bioanalyzer was employed to determine the RNA integrity number (RIN) and 28S/18S values for quality control (RIN ≥ 6.0 and 28S/18S ≥ 0.7) and confirmed by electrophoresis on a denaturing agarose gel.

Analysis of m6A level was performed on the UPLC-QQQ-MS system, consisting of a Waters Acquity UPLC (Waters Corporation, Massachusetts, United States) and an AB SCIEX 5500 QQQ-MS instrument (AB Sciex LLC, Framingham, United States). As previously described (Su et al., 2014), we prepared a mixture containing all of the cofactors and enzymes needed for hydrolysis of the RNA into nucleosides. Then, 5 µl of the mixture was added to each sample of RNA, and the final volume was brought to 25 µl with RNase-free water. The samples were incubated at 37°C for 4 h, after which 75 µL of methanol was added to the system. Different dilutions were used for detection. UPLC separation was performed on an Acquity UPLC BEH amide column (1.7 µm, 2.1 mm × 100 mm, Waters Corporation, Massachusetts, United States) with a flow rate of 0.35 ml/min at 40°C. Formic acid in water (0.1%, v/v, solvent A, Aladdin, Shanghai) and methyl cyanides (v/v, solvent B, Sigma-Aldrich, Shanghai) were employed as the mobile phase. Mass spectrometry detection was performed in positive electrospray ionization mode, and multiple reaction monitoring (MRM) was used to monitor target nucleosides by using the following mass transitions: (precursor ions → product ions) of m⁶A (282.1 → 150.1), A (268 → 136), U (245 → 113), G (284 → 152), C (244 → 112), and I (269 →137). Quantification was performed using a standard curve originating from A and m⁶A. The calculated m⁶A/A ratio was used to indicate the m⁶A level. To achieve maximal detection sensitivity, we optimized the MRM parameters of all nucleosides.

MeRIP-seq was performed by Biotechnology Corporation (Shanghai, China). As previously described (Dominissini et al., 2012), GenSeq^® m6A MeRIP KIT (GenSeq, China) was performed, and the protocol involved the following four steps. Step 1: RNA fragmentation. The RNA was randomly divided into approximately 200-nt segments by using fragmentation buffer and stop buffer. The sizes and concentrations of the RNA fragments were detected by an Agilent Bioanalyzer and an Agilent RNA 6000 Pico kit (Bioptic Inc., Taiwan, China). Three micrograms of fragmented RNA was used as the input group and stored at −80°C. The remaining fragmented RNA was used for subsequent immunoprecipitation experiments. Step 2: Preparation of immunoprecipitated magnetic beads. The porcine gastric mucine (PGM) magnetic beads were gently blown with a pipette to promote full suspension, after which they were washed with 1 × IP buffer, and an m⁶A antibody was added. Step 3: Immunoprecipitation. The MeRIP reaction solution, which included 50 µl of fragmented RNA, 150 µl of nuclease-free water, and 50 µl of 5 × IP buffer, was prepared. Next, 250 μL of the reaction solution was added to the magnetic beads prepared in step 2. The magnetic beads were washed with 1 × IP buffer, LB buffer, and HS buffer separately and repeatedly. Step 4: RNA purification. The magnetic beads were resuspended in RLT buffer, washed with 75% ethanol, and then completely resuspended in 11.25 µl of nuclease-free water. The eluted RNA was transferred to a new centrifuge tube and then immediately used for subsequent experiments or stored at −80°C. Both the input samples without immunoprecipitation and the m⁶A IP samples were used for RNA-seq library generation. The library quality was evaluated with a Bioptic Qsep100 Analyser (Bioptic Inc., Taiwan, China). Library sequencing was performed on an Illumina NovaSeq 6,000 (LC-Bio Technology) instrument with 150 bp paired-end reads.

FastQC was applied to analyse the quality of the sequencing data by assessing the sequencing quality distribution, base content distribution, and proportion of repeated sequencing fragments. HISAT2 software was used to compare the filtered clean reads with the human reference genome (GRCh38/hg38) and thereby obtain unique mapped reads for further analysis. Peak calling was carried out with exomePeak software. After obtaining the peak, the gene structure locations of the peaks and the overall distribution characteristics were determined, and peak annotation analysis was performed with HOMER software (http://homer.ucsd.edu/homer/ngs/peakMotifs.html).GO analysis is used to assess the function of genes from various aspects and is divided into three main categories: biological process (BP), molecular function (MF) and cellular component (CC). The mRNAs of genes modified by m⁶A were evaluated based on GO annotations in the BP, MF, and CC categories in the database, and Fisher’s test was used to determine the significance level (p-value) of each category to screen for significant GO terms. KEGG analysis was based on the gene annotations, and selected genes with mRNA m⁶A modification were annotated according to their associated KEGG pathways. The significance level (p-value) was determined by Fisher’s test to screen significant pathway terms related to m⁶A gene enrichment.

MeRIP-qPCR was performed by ChoudSeq Biotech Inc. (Shanghai, China). Its analysis was used according to a previously reported method (Dominissini et al., 2012), as shown in the RNA MeRIP-seq library construction and sequencing section. m6A enrichment was determined by qPCR analysis. Reverse transcription was performed using a SuperScriptTM III Reverse Transcriptase kit (Invitrogen) for mRNA. For relative qRT-PCR, qPCR SYBR Green master mix (CloudSeq)was used to generate mRNA cDNA according to the manufacturer’s instructions. The reactions were performed on a QuantStudio five Real-Time PCR System (Thermo Fisher). The expression levels of mRNAs normalized to those of the endogenous control were calculated using the 2^−ΔΔCT method and are presented as the fold change relative to the control group. % (IP/Input) was used to calculate the differences.The sequences of the primers utilized in this study are listed in Additional file 1: Supplementary Table S1.

Data are presented as the mean ± SD. Significance differences between groups were determined by Student’s t-test using GraphPad Prism six sofeware. Significance was established at p < 0.05.

We enrolled three NMOSD patients and three HCs for MeRIP-seq, with mean ages of 47.33 ± 2.082 and 47.33 ± 3.055, respectively (p > 0.05). The female: male ratio was 3:0 in both groups (p > 0.05). The current disease duration at the time of recruitment was 5 ± 2 days (Table 1). We enrolled another three NMOSD patients and three HCs for MeRIP-qPCR validation, with mean ages of 59.00 ± 6.083 and 59.67 ± 6.807, respectively (p > 0.05). The female:male ratio was 3:0 in both groups (p > 0.05). The current disease duration at the time of recruitment was 5 ± 2 days (Table 1).

The levels of m⁶A differed significantly between the NMOSD patients and HCs (0.09531 ± 0.009259 (%), 0.1399 ± 0.02533 (%); p < 0.05) (Table 1; Figure 1A). Additionally, R ² was greater than 0.99, indicating that the linear relationship of the standard curve is good (Figure 1B).

Some joint and low-quality sequences were abtained in the original data, and adaptor and low-quality data were removed to abtain clean reads. In the MeRIP-seq library, 6,308,338,650, 7,779,927,600 and 6,841,024,050 total bases and 6,125,183,549, 7,628,452,162 and 6,635,546,539 valid total bases were abtained, and the effective bases accounted for 97.10, 98.10 and 97.00% of the totals, respectively. In the HC blood samples, 6,495,933,300, 6,730,312,650 and 6,736,672,200 total bases and 6,292,021,712, 6,453,732,389 and 6,604,071,380 valid total bases were obtained, with the effective bases accounting for 96.90, 95.90 and 98.00% of the totals, respectively. In the RNA-seq library, the effective bases accounted for more than 90% of the total bases in both the NMOSD and HC groups. The results are shown in Table 2.

We used HISAT2 to map reads to the GRCh38/hg38 genome with default parameters. Detailed statistical analyses were performed by comparing reads with reference sequences. In the MeRIP-seq library, the mapping ratios were 94.37, 87.19 and 91.55% in NMOSD patients, and 92.64, 93.87, 89.08% in HCs. In the RNA-seq library, the mapping ratios were 91.07, 89.26 and 86.16% in NMOSD patients and 89.31, 90.42, 88.04% in HCs. The ratios of unique mapped reads are shown in Table 3.

MeRIP-seq analysis of RNA derived from whole blood revealed 14,444 nonoverlapping m⁶A peaks within 7,034 coding transcripts (mRNAs) and 440 nonoverlapping m⁶A peaks within 330 lncRNAs in the NMOSD group; there were three biological replicates. In the HC group, there were 12,806 nonoverlapping m⁶A peaks within 6,480 coding transcripts (mRNAs) and 369 nonoverlapping m⁶A peaks within 265 lncRNAs, with three biological replicates. Of these, 5,568 coding transcripts (70.1%) within mRNAs (Figure 1C) and 199 long noncoding transcripts (50.3%) within lncRNAs (Figure 1D) overlapped between the NMOSD patients and HCs. Of these, there were 3,600 differential peaks (p < 0.05), including 3,371 differentially methylated peaks within mRNAs and 95 within lncRNAs (p < 0.05) (Figure 1I).

To analyse the distribution profiles of m⁶A peaks within mRNAs, the peaks were categorized into five transcript segments: the 5’UTR; start codon (400 nucleotides centred on the start codon); coding sequence (CDS); stop codon (400 nucleotides centred on the stop codon); and 3’UTR. In our study, m⁶A was most often detected in the CDS, with some sites being observed near the 3’UTR and stop codon (Figures 1E–G).

Motif analysis of peaks within mRNAs with the highest scores (p-value = 1e-141) obtained from three biological replicates revealed a consensus sequence (GGAC) as well as other motifs in the NMOSD and HC samples (Figure 1H), indicating the reliability of the data.

In total, we identified 3,371 differentially methylated m⁶A sites (DMMSs) within mRNAs, of which 68.5% (2,310/3,371) showed significantly increased methylation and 31.5% (1,061/3,371) significantly reduced methylation (NMOSD vs. HCs; Table 4). We also identified 95 DMMSs for lncRNAs, of which 68.5% (71/95) exhibited significant increased methylation levels and 31.5% (24/95) exhibited significantly decreased methylation levels (Table 4). Tables 5, 6 provide the top ten m⁶A sites displaying increased and reduced methylation levels with the highest fold change values (Tables 5, 6).

To elucidate the functions of m⁶A in NMOSD, we selected protein-coding genes containing DMMSs for GO and KEGG pathway analyses. In the BP category, genes with increased m⁶A site methylation were significantly (p < 0.05) enriched in RNA splicing, mRNA processing, and gene expression (Figure 2A); genes with decreased m⁶A methylation were highly enriched in the regulation of transcription, transcription, DNA templating, and viral processes (Figure 2B). In the MF category, poly(A) RNA binding and protein binding were enriched among mRNAs with increased m⁶A methylation (Figure 2A), and those with lower levels of methylation were enriched in protein binding and DNA binding (Figure 2B). With respect to the CC category, the nucleus and nucleolus showed enrichment for both increased (Figure 2A) and decreased m⁶A methylation (Figure 2B). We also used a GO bubble chart to illustrate the top twenty enriched GO terms, mainly involving transcription, DNA templating, regulation of transcription, and gene expression, among others (Figure 2C).

We performed KEGG pathway analysis of DMMSs and observed that genes with increased m⁶A methylation were significantly (p < 0.05) enriched in the mRNA surveillance pathway, neurotrophin signalling pathway, and Epstein-Barr virus infection, among others (Figure 3A). Genes with decreased m⁶A methylation were significantly (p < 0.05) enriched in the TNF signalling pathway, oestrogen signalling pathway, and focal adhesion (Figure 3B). We also used a KEGG bubble chart to show the top twenty enriched pathways, including Epstein-Barr virus infection and the mRNA surveillance pathway (Figure 3C). These results suggest that m⁶A may have many key roles in NMOSD.

We drew a four-quadrant map based on two sets of data: RNA-seq, all gene differential expression table (including non-significant results); and RNA modification, all differential peak table (including nonsignificant results). The top 10 genes (the largest absolute values of diff.log2. fc are shown in purple) were determined to be MT-RNR1, C60rf203, XK, CD22, SEMA3A, HECW2, TIGD5, VISIG4, HLA-DQB1 and NOC2L. (Figure 4A).

In order to verify the reliability of MeRIP-seq results, MeRIP-qPCR was used to verify the modification levels of modification sites on the target genes.We chose the top four hypomethylated genes and the top four hypermethylated genes as target genes (Figure 4A). Among the eight target genes, MeRIP-qPCR results of MT-RNR1, XK, CD22, SEMA3A, HECW2 and TIGD5 were consistent with the sequencing results, and there were significant differences (p < 0.05) (Figure 4B). MeRIP-qPCR results of VISIG4 and C60rf203 were opposite with the sequencing results, although there were significant differences (p < 0.05) (Figure 4B). Nevertheless, it suggested that our sequencing results were credible.