Epigenome‐wide DNA methylation analysis of myasthenia gravis

Myasthenia gravis (MG) is a common neuromuscular junction disorder and autoimmune disease mediated by several antibodies. Several studies have shown that genetic factors play an important role in MG pathogenesis. To gain insight into the epigenetic factors affecting MG, we report here genome‐scale DNA methylation profiles of MG. DNA was extracted from eight MG patients and four healthy controls for genome‐wide DNA methylation analysis using the Illumina HumanMethylation 850K BeadChip. Verification of pyrosequencing was conducted based on differential methylation positions. Subsequently, C2C12 and HT22 cell lines (derived from mouse) were treated with demethylation drugs. Transcribed mRNA of the screened differential genes was detected using quantitative real‐time PCR. The control and MG group were compared, and two key probe positions were selected. The corresponding genes were CAMK1D and CREB5 (P < 0.05). Similarly, the myasthenic crisis (MC) and non‐MC group were compared and four key probe positions were selected. The corresponding genes were SAV1, STK3, YAP1, and WWTR1 (P < 0.05). Subsequently, pyrosequencing was performed for verification, revealing that hypomethylation of CAMK1D was significantly different between the MG and control group (P < 0.001). Moreover, transcription of CREB5, PKD, YAP1, and STK3 genes in the C2C12 cells was downregulated (P < 0.05) after drug treatment, but only YAP1 mRNA was downregulated in HT22 cells (P < 0.05). This is the first study to investigate genome‐scale DNA methylation profiles of MG using 850 K BeadChip. The identified molecular markers of methylation may aid in the prevention, diagnosis, treatment, and prognosis of MG.

Myasthenia gravis (MG) is a common neuromuscular junction disorder and autoimmune disease mediated by several antibodies. Several studies have shown that genetic factors play an important role in MG pathogenesis. To gain insight into the epigenetic factors affecting MG, we report here genome-scale DNA methylation profiles of MG. DNA was extracted from eight MG patients and four healthy controls for genome-wide DNA methylation analysis using the Illumina HumanMethylation 850K BeadChip. Verification of pyrosequencing was conducted based on differential methylation positions. Subsequently, C2C12 and HT22 cell lines (derived from mouse) were treated with demethylation drugs. Transcribed mRNA of the screened differential genes was detected using quantitative real-time PCR. The control and MG group were compared, and two key probe positions were selected. The corresponding genes were CAMK1D and CREB5 (P < 0.05). Similarly, the myasthenic crisis (MC) and non-MC group were compared and four key probe positions were selected. The corresponding genes were SAV1, STK3, YAP1, and WWTR1 (P < 0.05). Subsequently, pyrosequencing was performed for verification, revealing that hypomethylation of CAMK1D was significantly different between the MG and control group (P < 0.001). Moreover, transcription of CREB5, PKD, YAP1, and STK3 genes in the C2C12 cells was downregulated (P < 0.05) after drug treatment, but only YAP1 mRNA was downregulated in HT22 cells (P < 0.05). This is the first study to investigate genome-scale DNA methylation profiles of MG using 850 K BeadChip. The identified molecular markers of methylation may aid in the prevention, diagnosis, treatment, and prognosis of MG.
Myasthenia gravis (MG), a common neuromuscular junction (NMJ) disorder, is an autoimmune disease mediated by several antibodies. Clinical features of MG mainly include fluctuating muscle weakness and fatigability after an activity. MG often affects the extraocular muscles, bulbar muscles, and the proximal limb muscles. Around 15-25% of MG patients will have myasthenic crisis (MC), which causes respiratory muscle weakness and respiratory failure, increasing the risk of death. According to reports, the annual incidence of MG in southern China is estimated to be 1.55-3.66 per 100 000, and the prevalence is estimated to be 2. .07 per 100 000 [1]. Current studies have shown that genetic susceptibility, autoimmune environment, and autoantibodies participate in the pathogenesis of NMJ disorder. Moreover, the above factors affect the production of endplate potential and muscle contraction causing MG [2]. Many studies have shown that autoantibodies against normal body antigens, including acetylcholine receptors (AChRs), musclespecific kinase, and lipoprotein-related protein 4, are related to the development of MG [3]. Genetic factors are also one of the important factors for MG development. Moreover, cytotoxic T lymphocyte correlated antigen-4 (CTLA4) gene, the human leukocyte antigen (HLA) gene, the cholinergic receptor nicotinic alpha 1 (CHRNA1) gene, and the protein tyrosine phosphatase nonreceptor type 22 (PTPN22) gene participate in MG development [4][5][6]. With the development of detection methods, the role of epigenetics in diseases has gradually been valued.
DNA methylation is one of the most common phenomena of epigenetic alterations. DNA methylation is a form of chemical modification of DNA, which alters RNA transcription without changing the DNA sequence. Abnormal DNA methylation is one of the most extensive epigenetic changes in cancer development and can alter gene expression [7,8]. Changes in DNA methylation are heritable and affect the active state of the X chromosome, and gene expression. Mamrut et al. compared the transcriptome and methylome of peripheral monocytes in female twins with MG. The results showed that alteration of gene expression or methylation level may jointly cause MG [9], suggesting that genetic susceptibility may contribute more to the pathogenesis of MG than previously thought.
In this research, the epigenetic DNA methylation pattern in MG patients and healthy controls, and the differentially methylated positions (DMPs) were identified. The new potential biological functions of genes are proposed through bioinformatics tools analysis.

Ethical considerations
This study was conducted strictly in accordance with the Helsinki Declaration and was approved by the permission Ethics Committee of the Second Affiliated Hospital of Wenzhou Medical University and Yuying Children's Hospital (LCKY2018-12). In addition, patients provided informed consent to participate in this study in writing.

Participants and design
The study participants received treatment at the Department of Neurology, the Second Affiliated Hospital of Wenzhou Medical University and Yuying Children's Hospital. The patients were diagnosed by two neurologists at least at the deputy director level in accordance with the 2015 China Myasthenia Gravis Diagnosis and Treatment Guidelines. The diagnostic criteria were as follows: (a) with muscle weakness of some specific striated muscles, with mild symptoms in the morning but severe in the afternoon; (b) positive for the neostigmine test; (c) positive for the repetitive nerve stimulation; and (d) with AChR antibodies in the peripheral blood of most patients with generalized MG. Basing on the typical clinical characteristics MG, has pharmacological and/or neurophysiological characteristics, and can be diagnosed in the clinic. Qualified units can detect anti-AChR antibodies in the blood of patients for further verification of the diagnosis. In summary, the diagnosis was done if the above first point and one of 2, 3, 4 points on selection criteria were met. Patients with other diseases, such as Lambert-Eaton myasthenic syndrome, progressive bulbar palsy, botulism, and metabolic myopathy, were excluded. The MC diagnostic criteria were as follows: sudden increase in muscle weakness, progressive weakness, or paralysis of the respiratory and swallowing muscles, all of which are life-threatening.

DNA methylation experiment
DNA was extracted from the peripheral blood of the patients using the DNeasy Blood and Tissue Kit (Qiagen, Beijing, China). The purity and concentration of DNA were measured using Nanodrop 2000 (Thermo, Beijing, China). A concentration of 500 ng DNA of each sample was converted by bisulfite using EZ DNA Methylation Kits (Zymo Research, Tustin, CA, USA). The converted products were loaded into 850K BeadChips in accordance with the manufacturer's guide and protocol (Illumina, Hayward, CA, USA). The chip detects more than 850 000 sites in each sample, comprehensively covering CpG islands, promoters, coding regions, etc. The whole processing steps related to chips, such as amplification, labeling, cleaning, hybridization, and scanning, were performed by Shanghai Pudong Decode Life Science Research Institute.

Data analysis
The data were analyzed using CHAMP [10] package in R. The β value was used to represent the proportion of each methylated CpG position. Firstly, probes with detection Pvalue < 0.01, probes with < 3 beads in at least 10% of samples, non-CpG probes, multihit probes, and probes located on chromosomes X and Y and (SNP-related probes) were filtered. The final probes for subsequent analysis were 820 000. Afterward, the β value matrix was normalized using BMIQ [11] for adjusting type I and type II probe bias. Next, singular value decomposition analysis (SVD) [12] was then performed to the batch effect caused by BeadChip Slide and Array. COMBAT [12] was applied to correct the batch effect. All CpG sites were annotated using EPICAN-NO.ILM10B4.HG19 [13]. Differential methylated CpG positions (probes) were identified using champ. DMP function initiates the limma package to calculate the P-value for differential methylation using a linear model. The DMP function and the adj.P values were computed using the Benjamini-Hochberg method [14]. CpGs having |Δβ| ≥ 0.20 and adjusted adj.P-value ≤ 0.05 were considered as DMPs. Differential methylated regions (DMRs) were calculated using PROBE LASSO [15]. All CpGs were for DMR calling, and minSigProbesLasso, lassoRadius, and boundary were the default values. Gene Ontology (GO) [16] and Kyoto Encyclopedia of Genes and Genomes (KEGG) [17] enrichment analysis was conducted using the R package to evaluate the biological processes, cellular components, molecular mechanism, and pathway of DMPs and DMRs. Finally, functional epigenetic modules (FEM) were applied to infer differential methylated gene modules in protein-to-protein interaction network (PPI) using the FEM [18] package. The copy number variation was also analyzed using the CONU-MEE [19] package.

Pyrosequencing verification
Pyrosequencing is gradually replacing bisulfite sequencing polymerase chain reaction and has become the gold standard for allele quantification in methylation research [20]. Therefore, pyrosequencing was preferred in this study to confirm the methylation level of the 850K array. Ten samples of MG patients were also included for pyrosequencing, besides those used for the 850K array.  Table S1. Gene expression was calculated as 2 ÀΔΔCt , with GAPDH mRNA as the internal control.

Western blot
Cells were first harvested from the culture medium. Subsequent protein extraction was performed on ice for 30 min using RIPA lysis buffer (P0013B; Beyotime, Beijing, China) and protease and phosphatase inhibitor cocktail (P1046; Beyotime). The protein sample was at 13 400 g for 30 min at 4°C. The proteins were concentrated using a BCA kit (MA0082, Meilunbio, Dalian, China), and western blotting was performed as previously described. Primary antibodies against DNMT1 were purchased from Santa Cruz Biotechnology (Santa Cruz, CA, USA). The β-actin was purchased from Affinity Biosciences (Affinity, Cincinnati, OH, USA). The goat anti-rabbit IgG and anti-mouse IgG secondary antibodies were purchased from Affinity Biosciences (Affinity).

Patient characteristics
Blood samples from eight MG patients and four healthy controls were collected for Infinium Methyla-tionEPIC BeadChip. The proportion of male-to-female patients was 37.5% and 62.5%, respectively. The average age of the study participants was 50 AE 12.9 years old. Three of the eight MG patients previously had MC, while three had MC combined with thymoma (Table 1). After analysis of the sequencing results, 10 MG patients were selected for pyrosequencing verification to confirm the universality of the selected gene methylation changes (Table 1).

Quality control of genome-wide methylation detection
After quality control of the samples and chips, the median, P25, and P75 quantile values were used to describe the data. The results showed that the values of the experimental samples were evenly distributed, indicating that the data obtained by chip detection were symmetrical (Fig. 1A). CpGs that interfered with subsequent analysis were filtered out using the SVD [21] for quality control on the obtained raw data. The interbatch differences using were corrected using the combat algorithm [12]. Sample cluster analysis mainly showed the degree of similarity among all samples. In the tree graph, the samples were mainly divided into two clusters, with the upper group being MG group. When the cutting height further decreased, the samples were further divided into MC group and the non-MC group, with the lower group being the control group (Fig. 1B). The principal component analysis (PCA) identified the clusters of samples based on methylation level, which showed the difference between the MG and the control group (Fig. 1C). Similarly, the sample correlation heat map demonstrated the consistency and difference in DNA methylation between those two groups (Fig. 1D). After the homogenization process above, the possible influence of the methylation level variation was eliminated.

Differential methylation position (DMP) analysis
The  Table 2). Table 2 showed that DMPs located in the body region accounted for a large proportion. Hence, the main focus was the analysis of loci located in the promoter region of functional genes. The comparison between MG and the control group showed that DMPs are widely distributed on chromosomes 1-22 (Fig. 2). Furthermore, the comparison between the MG and the control group showed that the top 30 GO terms of GO enrichment were enriched in the MG group, including homophilic cell adhesion via plasma membrane adhesion molecules, calcium ion binding, and neuron projection (Fig. 3A). Among the 30 GO terms, eight were related to synapses, and two were related to nerves. This study showed that differentially methylated sites are related to NMJ dysfunction, consistent with MG. The GO enrichment analysis results between MC and non-MC groups are shown in Fig. 3B. The top 30 GO enrichment items included multicellular organism development and homophilic cell adhesion via plasma membrane adhesion molecules. Of these, four are related to actin and six were related to cell junctions. Actin is one of the structural proteins in muscles. Actin plays an important role in muscle movement and is also consistent with symptoms in MC patients.
In addition, 30 enriched pathways were identified after the KEGG analysis of DMP among the MG group compared with the control group. Some of the enriched pathways included aldosterone synthesis and secretion, morphine addiction, and calcium signaling (Fig. 4A). Among the 30 KEGG-enriched pathways, multiple pathways related to calcium ion channels were found, providing a basis for the subsequent screening of functional genes of abnormal methylation sites in the MG group. The KEGG analysis of the MC versus non-MC group showed that hippo signaling pathway, morphine addiction, axon guidance, etc. were enriched in the top 30 enrichment pathways (Fig. 4B). Among the top 30 enriched pathways, the hippo signaling pathway showed a significant difference, and this pathway participates in cell death, differentiation, and inhibition of cell proliferation [22]. Clinically, patients with MG often have abnormal thymus, raising the need for subsequent research in this field. While comparing the MG and control group, 20 key DMPs related to MG were identified, and two key probe sites were selected based on the KEGG results. The corresponding genes were CAMK1D and CREB5 (P < 0.05), which were located in the key position of the calcium ion signal pathway. The MC and non-MC groups were also compared, and 20 important DMPs related to MG were identified, and four key probe sites were selected based on the KEGG results. The corresponding genes were SAV1, STK3, YAP1, and WWTR1 (P < 0.05), which were also located in the key position in the hippo signaling pathway.

Differential methylation region analysis
The MG and control group were compared, and a total of 409 DMRs were filtered in the whole MG patients, such as HLA-DQB1, and other genes of the The sample correlation heat map for two groups of differential methylation sites. The color ranges from red to blue, indicating that the correlation between the two samples is from low to high. The higher the sample correlation, the more similar the overall methylation pattern. Besides, sample G1-G8 represents MG patients in the non-MC group, W1-W3 represents MG patients in the MC group, and C1-C4 represents the healthy control group. HLA family were identified in the MG group. Among these 409 DMRs, 116 (28.4%) were located on chromosome 6. According to several studies, chromosome 6 contains the HLA gene, which is associated with autoimmune diseases. The GO enrichment analysis of DMR showed that the GO terms were significantly enriched in neurotransmitter transporter activity, neurotransmitter: sodium symporter and organic acid: sodium symporter activity, which was related to the pathogenesis of MG (Fig. 5A). In addition, we also compared the MC and the non-MC groups were compared, and 343 DMRs were identified. In addition, plasma lipoprotein particle assembly, neurotransmitter transporter activity, neurotransmitter: sodium symporter, and organic acid: sodium symporter activity also have significant enrichment in the GO enrichment analysis, consistent with the phenomena between the MG and control groups (Fig. 5B).

Functional epigenetic modules analysis
Epigenetics plays an important regulatory role in cell development, cell differentiation, complex disease formation, and tumorigenesis. FEM analysis is a function supervision algorithm that helps identify some genes or proteins with important functions. As shown in Fig. 6, the seed gene of the module was DTNA. In addition, SYNM, SNTG1, SNTG2, etc. were all positively correlated. The negatively significant genes included FMN1, MAGEE1, and ADNP2. DTNA, a seed gene, demonstrated its importance in the whole network. According to FEM analysis, the differential modular seed genes and related interacting genes were screened out in the MG and the control group, as well as the comparison between the MC group and non-MC group.

Pyrosequencing confirms differences in CAMK1D
The current research believes that based on the genome-wide level of DNA methylation studies, although the sensitivity and breadth have been improved compared with previous techniques, the results obtained still have false-positive results. During validation, the potential MG differential methylated genes were revalidated using the gold standard for methylation level detection-pyrosequencing technology. The MG and the control groups, and added 10 other samples of the MG group were analyzed based on the existing data. Subsequently, the CAMK1D (cg02323098) and CREB5 (cg16989544) were verified by pyrosequencing. The experimental results showed that the methylation level of CAMK1D (cg02323098) in the peripheral blood of MG patients decreased by 52% compared with the control group. The statistical results were significantly different (P < 0.001) (Fig. 7).  6. Function epigenetic module display diagram. Each node is a gene, and the color of the node is whether the calculated difference statistic is significant, blue represents positive significance, yellow represents negative significance, and gray represents no difference. The line between nodes is the degree of connectivity between genes, and the thickness of the line represents the level of correlation.

Detection of target gene expression in C2C12 cell line and HT22 cell line after demethylating drugs
To further verify the above results, procaine hydrochloride was used to verify demethylation in the C2C12 cell. The mRNA levels of the selected part of the target genes changed significantly, suggesting that the drug treatment was effective. The mRNA levels of CREB5 and PKD were downregulated after procaine hydrochloride [23] treatment compared with the control group (P < 0.05) (Fig. 8A). Similarly, compared with the control group, mRNA levels in YAP1 and STK3 were also decreased after procaine hydrochloride treatment (P < 0.05) (Fig. 8B). No significant difference was observed between CAMK1D and CALM. 5-Aza-2 0 -deoxycytidine [24] and procaine hydrochloride were used for demethylation analysis in HT22 cells. The expression of DNMT1 protein decreased after treatment with the demethylation drugs (Fig. 9A). The experiment showed the expression of DNMT1 was significantly inhibited by DAC (Fig. 9B), indicating that the demethylation effect could be achieved after treatment with demethylation drugs. After 3 days of procaine hydrochloride treatment, only the YAP1 mRNA transcription was lower after demethylation drugs treatment than in the control group (P < 0.05) (Fig. 10B). Conversely, no significant difference was observed in CAMK1D, CALM, CREB5, PKD, and STK3 (Fig. 10A,B).

Discussion
In this research, the genome-wide methylation of MG and healthy controls was analyzed using the Infinium MethylationEPIC BeadChip, and the potential methylation sites associated with MG were screened out. Compared with the control group, 1722 DMPs and 409 DMRs were screened out in the MG group. MG is a disorder in the NMJ, CAMK1D and CREB5 were identified, and both were key nodes in the calcium ion signaling pathway using functional annotation. The probes were further validated through pyrosequencing, which showed that the expression of the cg02323098 probe, named CAMK1D, was significantly downregulated compared with MG and control groups (P < 0.05). Related studies have found that CAMK1D is associated with lung adenocarcinoma cells, type 2 diabetes mellitus, and essential hypertension [25][26][27]. In addition, some studies have found that CAMK1D was a key regulator of tumor innate immune resistance  [28]. As an autoimmune disease mediated by autoantibodies, MG can be observed in patients with immunomodulatory defects [29]. However, the ability of CAMK1D to regulate this immune response remains to be validated.
The MC and the non-MC group were compared and screened out 2924 DMPs and 343 DMR. Similarly, GO enrichment and KEGG analysis were performed, and SAV1, STK3, YAP1, and WWTR1, which were located in the key position of the hippo signaling pathway, were screened out. The hippo signaling pathway is a conservative signal network that regulates tissue homeostasis and organ development. Studies have suggested that the hippo signaling pathway promotes cell death and differentiation, and inhibits cell proliferation [22]. Moreover, the hippo signaling pathway is related to a variety of human diseases, including cancer, autoimmune diseases, and developmental abnormalities [30]. Most AChR + MG patients have thymic changes, 10% have a thymoma, and 80% of the patients with early-onset MG have thymic follicular hyperplasia [31,32]. However, whether the changes in the thymus of MG patients are related to the hippo signaling pathway remains to be validated.
In the results, some DMPs/DMRs (such as the presence of low-level methylation in HLA-DQB1) have been previously studied and found to be related to the pathogenesis of MG. The past 30 years of research show that the HLA genome contains a large number of genes related to the human immune system. Hence, it is a strong risk factor related to the development of MG [33]. The specific HLA alleles are also thought to be related to the susceptibility of MG [34]. This study is more advantageous because Illumina 850K Beadchips, instead of the Illumina 450K, were used to cover doubled methylation sites. Therefore, more extensive and detailed information was obtained. In the current study, some analyses did not give satisfying results due to the limitation of our small sample size. For example, morphine addiction was observed in the two groups of analysis for the significantly different pathways after KEGG analysis. The above observation lacks a solid explanation for MG patients. Increasing the sample size could validate the results of this study.
Due to the limited number of samples, cell culture experiments were performed further to verify the above findings at the cellular level. The results showed no significant difference in the expression of CAMK1D and CALM in the C2C12 cell line and the HT22 cell line. However, the expression of the target gene was partially upregulated under treatment with demethylating drugs. Ca 2+ /calmodulin-dependent protein kinases (CaM-Ks) regulate various biological events mediated by intracellular Ca 2+ , including muscle contraction, neurotransmitter release, and gene expression [35]. Studies have established Xenopus neuromuscular cocultures and found that activation of CaMK is essential for neurotrophic factors to regulate NMJ development [36]. Other studies have found that CaMK mediates the inhibitory effect of BDNF on NMJ maturation, and Ca 2+ is released from the cell through the IP3 receptor or Ryanodine receptor, regulating the neurotrophic effect of NMJ maturation [37]. Another differential gene, CREB5, was significantly low in C2C12 cell line (P < 0.001). The level of PDK, an upstream gene was also significantly low (P < 0.05). Although no significant difference was observed in HT22 cells, its expression still decreased. As a transcription regulator, CREB5 is a transcription factor involved in cell survival, cell proliferation, cell adaptation, and differentiation. CREB plays an important role in the immune system by regulating the expression of several inflammatory mediators in white blood cells and macrophages [38][39][40]. Due to limited conditions, the pyrosequencing of YAP1 and STK3 genes could not be verified in this study. However, following the interest in these genes and their corresponding pathway, the mRNA levels of YAP1 and STK3 genes were analyzed at the cellular level. In the C2C12 cell line, the mRNA levels of YAP1 and STK3 genes were significantly lower compared with the control group (P < 0.05). Similarly, the YAP1 mRNA level in the HT22 cell line after 3 days of procaine treatment was lower compared with control group (P < 0.05). These results are contrary to the expectation of this study. Considering that hippo signaling pathway is a complex signaling pathway, and each key step is affected by multiple factors, the pathway could be regulated by inhibition and promotion factors. However, the nonspecific demethylation drugs used in this study resulted in extensive knockdown of hypomethylation levels, which could not exclude the influence of other genes on the expression of the target gene. However, bioinformatic analysis revealed changes in methylation levels of several key genes of hippo signaling pathway. The hippo signaling pathway is associated with the pathogenesis of several human diseases, including cancer, autoimmune diseases, and developmental abnormalities [30]. The hippo signaling pathway is still a hotspot research area for future studies.
Epigenetic marks are considered one of the important factors related to the development of MG. The screening of methylation molecular markers provides new research ideas and clues for the pathogenesis of MG, and the occurrence and development of the disease.

Conclusion
850K BeadChips sequencing was performed on the collected MG patients and healthy individuals. The differences in the aldosterone synthesis, secretion, and hippo signaling pathway between MG patients and healthy individuals were determined using GO enrichment and KEGG enrichment analysis. Pyrosequencing analysis confirmed the methylation changes of CAMK1D, which is worthy of further studies. This is the first study to investigate the genome-scale DNA methylation profiles of MG disease using the Illumina Human-Methylation 850K BeadChip. DNA methylation may be involved in the occurrence and development of MG. It is hoped that through future research, abnormal methylation positions could be used as markers for the clinical screening of MG, screening of related populations, and monitoring MG progression. Thus, methylation analysis has future clinical application potential.