Comprehensive Analysis of N6-Methyladenosine (m6A) Methylation in Neuromyelitis Optica Spectrum Disorders

Background: N6-Methyladenosine (m6A) methylation is the most prevalent internal posttranscriptional modification on mammalian mRNA. But its role in neuromyelitis optica spectrum disorders (NMOSD) is not known. Aims: To explore the mechanism of m6A in NMOSD patients. Methods: This study assessed the m6A methylation levels in blood from two groups: NMOSD patients and healthy controls. Methylated RNA immunoprecipitation Sequencing (MeRIP-seq) and RNA-seq were performed to assess differences in m6A methylation between NMOSD patients and healthy controls. Ultra-high performance liquid chromatography coupled with triple quadruple mass spectrometry (UPLC-QQQ-MS) method was performed to check m6A level. Differential m6A methylation genes were validated by MeRIP-qPCR. Results: Compared with that in the control group, the total m6A level was decreased in the NMOSD group. Genes with upregulated methylation were primarily enriched in processes associated with RNA splicing, mRNA processing, and innate immune response, while genes with downregulated methylation were enriched in processes associated with the regulation of transcription, DNA-templating, and the positive regulation of I-kappa B kinase/NF-kappa B signalling. Conclusion: These findings demonstrate that differential m6A methylation may act on functional genes to regulate immune homeostasis in NMOSD.


INTRODUCTION
N 6 -Methyladenosine (m 6 A) is one of the most abundant internal modifications of eukaryotic messenger RNA (mRNA) and plays an important role in gene expression regulatory processes, including the maintenance of stability, splicing and translation (Chen et al., 2019a). Approximately 25% of mRNAs are estimated to contain at least one m 6 A nucleotide site, and most of them occur at the sixth position of adenosine. m 6 A was shown to be specifically catalysed at the consensus gene sequence DRACH (D A/G/U, R A/G, H A/C/U). In addition to m 6 A, other types of posttranscriptional modifications exist, including 5-methylcytosine (m 5 C) and pseudouridine (ψ). To detect regions containing m 6 A peaks, researchers have developed an immunoprecipitation-based method combined with high-throughput sequencing, termed methylated RNA immunoprecipitation sequencing(MeRIPseq) (Dominissini et al., 2012). To determine m 6 A level, researchers have developed the UPLC-QQQ-MS method, which has the advantages of short analysis time, good resolution and high sensitivity (Su et al., 2014). m 6 A modification is dynamic and reversible and regulated by m 6 A methyltransferases (writers), demethylases (erasers) and reading proteins (readers). Writers include Methyltransferase Like 3 (METTL3), methyltransferase-like 14 (METTL4) and Wilms tumour 1-associated protein (WTAP). METTL3 catalyses the production of m6A; METTL4 forms a complex with METTL3, participates in interactions with target mRNA and recruits RNA, and WTAP stabilizes the complex (Liu et al., 2014;Ping et al., 2014;Lin et al., 2016). Erasers, including fat and obesityassociated protein (FTO) and alk B homologue 5 (ALKBH5), catalyse the demethylation of bases that have been modified by m 6 A (Jia et al., 2011;Zheng et al., 2013). Readers, including YTH N6 methyladenosine RNA-binding protein 1-3 (YTHDF1-3) and insulin-like growth factor 2 mRNAbinding protein 1-3 (IGF2BP1-3), are also necessary for this process Wang et al., 2015;Hanniford et al., 2020). These proteins recognize sites of mRNA methylation and affect RNA metabolism directly or indirectly. Overall, m6A-modifying enzymes play important roles in tumorigenesis and in the development of the human central nervous system (Lin et al., 2017;Zhang et al., 2018;Heck et al., 2020;Xue et al., 2021). In recent years, studies have revealed that m 6 A participates in regulating the immune microenvironment (Chen et al., 2019b;Tang et al., 2020;Xu et al., 2020), but its involvement in NMOSD is unknown.
NMOSD is a kind of recurrent antibody-mediated neuroimmune disease, in which the optic nerve, spinal cord, and area postrema of the medulla are damaged, often causing severe symptoms (Wingerchuk et al., 2015). The specific diagnostic and pathogenic biomarker for NMOSD is aquaporin 4immunoglobulin IgG (AQP4-IgG), which binds to the foot processes of astrocytes, activating complement and thereby recruiting neutrophils and eosinophils; this leads to astrocyte death, followed by axonal degeneration and oligodendrocyte injury (Lennon et al., 2004). The prevalence and incidence of NMOSD are 0.52-4.4 and 0.05-0.40 per 100,000, respectively (Pandit et al., 2015). Approximately 60% of NMOSD patients relapse within 1 year and approximately 90% within 5 years. Furthermore, approximately 50% of NMOSD patients require a wheelchair within 5 years of their diagnosis due to limb paralysis, and 62% become blind, affecting their quality of life and imposing a heavy burden on patients, families and society (Pittock and Lucchinetti, 2016). Currently, methylprednisolone and plasma exchange used in the acute stage as well as rituximab, mycophenolate mofetil and azathioprine used in the remission stage are only partially effective and do not substantially reduce the high disability associated with the disease or its recurrence. Therefore, it is necessary to further understand the pathogenesis of NMOSDs and explore more effective drugs.
Regardless, the m 6 A status has not been assessed in NMOSD. To investigate differences in m 6 A modification patterns between healthy controls (HCs) and NMOSD patients, we performed UPLC-QQQ-MS to assess the m 6 A level and MeRIP-seq to identify the first known transcriptome-wide m 6 A sites in NMOSD patients. In general, the total m6A level was decreased in patients compared to HCs. We also detected 3,371 differentially methylated peaks within mRNAs and 95 within long noncoding RNAs (lncRNCs, p < 0.05). Intriguingly, mRNAs harbouring differentially methylated peaks were shown to be involved in many important biological pathways associated with immune homeostasis.

Blood Collection From NMOSD Patients and HCs
(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. (1.2) Exclusion criteria: 1) severe mental symptoms precluding cooperation with treatment; 2) no informed consent signed; 3) receiving immunosuppressive therapy within 30 days of new clinical symptoms; 4) other autoimmune diseases or tumours; 5) new clinical symptoms and not receiving relevant treatment but with clinical symptoms for more than 30 days before seeing a doctor. (1.3) On the second day after admission, 10 ml of blood was collected into a BD Paxgene blood RNA tube; the samples were stored at −4°C for 24 h and then at −80°C for preparation of RNA. (1.4) In total, six healthy individuals with no acute or chronic illness from the Health Care Centre of the First Rehabilitation Hospital of Shanghai were enrolled in the study as a control group.

RNA Preparation and Quality Control
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.

UPLC-QQQ-MS Method
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.

RNA MeRIP-Seq Library Construction and Sequencing
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 6 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 6 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 6 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 6 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 6 A gene enrichment.

Methylated RNA ImmunoprecipitationeRIP-qPCR
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 MeRIPseq 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 Analysis
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.

Demographic and Clinical Features of the Patients With NMOSD
We enrolled three NMOSD patients and three HCs for MeRIPseq, 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).

Levels of m 6 A in NMOSD Patients and HCs
The levels of m 6 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 2 was greater than 0.99, indicating that the linear relationship of the standard curve is good ( Figure 1B).

Sequencing Statistics and Quality Control
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 Table 3.

General Features of m 6 A Methylation in NMOSD Patients and HCs
MeRIP-seq analysis of RNA derived from whole blood revealed 14,444 nonoverlapping m 6 A peaks within 7,034 coding transcripts (mRNAs) and 440 nonoverlapping m 6 A peaks within 330 lncRNAs in the NMOSD group; there were three biological replicates. In the HC group, there were 12,806 nonoverlapping m 6 A peaks within 6,480 coding transcripts (mRNAs) and 369 nonoverlapping m 6 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 6 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 6 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.
Frontiers in Genetics | www.frontiersin.org November 2021 | Volume 12 | Article 735454 displaying increased and reduced methylation levels with the highest fold change values (Tables 5, 6).

Differentially Methylated RNAs Are Involved in Important Biological Pathways
To elucidate the functions of m 6 A in NMOSD, we selected protein-coding genes containing DMMSs for GO and KEGG pathway analyses. In the BP category, genes with increased m 6 A site methylation were significantly (p < 0.05) enriched in RNA splicing, mRNA processing, and gene expression ( Figure 2A); genes with decreased m 6 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 6 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 6 A methylation ( Figure 2B). We also used a GO bubble chart to illustrate the top twenty enriched GO terms, mainly involving transcription, DNA templating, The consensus motif for m 6 A modification is "GGAC".(I) Volcano map of differential m 6 A modification. Diff.log2 FC < 0 indicates hypomethylation, indicating that the site was demethylated under treatment conditions; Diff.log2 FC > 0 indicates hypermethylation, indicating hypermethylation of this site under treatment conditions. Diff.lg.fdr: log10 (FDR) of differential methylation analysis is the log10 conversion value of the FDR value obtained from the differential analysis. Red indicates hypermethylation, and blue indicates demethylation. Red and blue represent genes with a more than 2-fold difference in m 6 A modification.
We performed KEGG pathway analysis of DMMSs and observed that genes with increased m 6 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 6 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 6 A may have many key roles in NMOSD.

Association Analysis of Differential RNA Modification and Differential Gene Expression
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).

Validation of Differential RNA Modifications by MeRIP-qPCR
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. DISCUSSION m 6 A modification is a dynamic and reversible RNA modification occurring in eukaryotes that is carried out by "writers", "erasers", and "readers". In recent years, studies have found that m 6 A is involved in the regulation of the immune microenvironment (Chen et al., 2019b;Tang et al., 2020;Xu et al., 2020). NMOSD is a kind of autoimmune disease that is triggered by disruption of immune homeostasis. In this study, we assessed the m 6 A status in NMOSD patients and HCs, revealing differences between the groups that support the dynamic features of m 6 A modification. Indeed, we found that the m 6 A level was much lower in patients than in HCs, demonstrating that m 6 A modification is very important in the process of NMOSD. In addition, Zhang found that altering m 6 A levels might be a novel way for improving the efficacy of platinum in cancer treatment (Zhang et al., 2021a), and Zhao reported that m 6 A regulation could influence the progression of lung adenocarcinoma (Zhao and Xie, 2021). We hypothesized that increasing the level of m 6 A might relieve the symptoms of NMOSD, but further studies are needed to confirm this hypothesis.
In the present study, many important biological pathways were found to be related to differentially methylated mRNAs. Because of the construction of strand-specific libraries, we also predicted lncRNAs harbouring m 6 A peaks and identified DMMSs. While lncRNAs are known to play important roles in mediating transcriptional and posttranscriptional regulation (Cantile et al., 2021;Moison et al., 2021), their biological functions in regards to m 6 A modification remain unknown. In general, the m 6 A modification of lncRNAs is very common in cancers. For example, Wu found that m 6 A-induced lncRNA RP11 could cause the dissemination of colorectal cancer cells . In our study, we investigated the m 6 A modification of lncRNAs between NMOSD patients and HCs, which highlighted immune homeostasis. However, further analysis is needed to validate these results.
Evidence supports a strong relationship between immune homeostasis and m 6 A modification, with m 6 A modification playing important biological roles in NMOSD (Chen et al., 2020;Chong et al., 2021;Zhang et al., 2021b;Zhong et al., 2021). In our study, the GO and KEGG analyses of mRNAs harbouring DMMSs showed that those with increased methylation were mostly enriched in RNA splicing, mRNA processing and gene expression, which supports the importance of m 6 A modification in NMOSD. For example, the increased expression of phosphatase and tensin homolog deleted on chromosome Ten (PTEN) and an approximately twofold enhancement methylation were observed in NMOSD patients compared with the HCs, indicating PTEN serves as a key enzyme in NMOSD ( Figure 5A) (Han et al., 2019). Furthermore, PTEN is involved in inositol phosphate metabolism. Gene expression is important in NMOSD, and our results suggest a link between NMOSD and the regulation of gene expression, subsequently resulting the disruption of immune homeostasis. In contrast to m 6 A sites with increased methylation, those with decreased methylation were mostly enriched in regulation of transcription, DNA templating, and positive regulation of I-kappa B kinase/NF-kappa B signalling pathways, which include genes such as v-myc avian myelocytomatosis viral oncogene homologue (MYC) ( Figure 5B) and nucleotide-binding oligomerization domain containing 2 (NOD2) ( Figure 5C). "Readers", which include YTHDF1-3, recognize and bind to m 6 A modifications on mRNAs Li et al., 2017;Jiang et al., 2019), thereby playing an important role in regulating translation. According to recent studies, we hypothesize that the methylation of mRNAs is associated with their translation, as it may influence their expression, leading to subsequent changes in global translation.
In our association analysis of differential RNA modification and differential gene expression, the top 10 genes included SEMA3A and CD22 ( Figure 4A). Zhang also reported that SEMA3A may negatively regulate axonal regeneration in retinal ganglion cells , and CD22 is an important inhibitory molecule on the surface of B cells that negatively regulates B cell activation (Gross et al., 2009). All these findings are very beneficial for future studies.
Although few studies have explored the roles of m 6 A modification in NMOSD, there are still some limitations to our study. For example, we may need to expand the number of enrolled patients. In addition, data from the same patient should be obtained before and after treatment.

CONCLUSION
Our study demonstrates the differential m 6 A methylome in NMOSD patients compared to healthy controls, suggesting a strong association between m 6 A methylation and the regulation of immune homeostasis in NMOSD. The findings fundamentally contribute to future studies on immune homeostasis in NMOSD.

DATA AVAILABILITY STATEMENT
The raw data have been made publicly available under SRA accession number PRJNA737585.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by The Ethics Committee of Shanghai First