Differential Transcription Profiling in Bone Marrow Mononuclear Cells between Myasthenia Gravis Patients with or without Thymoma

Purpose: Defective stem cells have been recognized as being associated with autoimmune diseases, such as systemic lupus erythematosus, rheumatoid arthritis, autoimmune cytopenias and Myasthenia Gravis (MG). However, the differential gene expression profile of Bone Marrow Mononuclear Cells (BMMCs) and the molecular mechanisms underlying MG pathogenesis have not been fully elucidated. Therefore, we investigated the abnormal expression and potential roles and mechanisms of mRNAs in BMMCs among patients with MG with or without thymoma. Methods: Transcription profiling of BMMCs in patients with MG without thymoma (M2) and patients with thymoma-associated MG (M1) was undertaken by using high-throughput RNA sequencing (RNA-Seq), and diseaserelated differentially expressed genes were validated by quantitative real-time polymerase chain reaction (qRT-PCR). Results: RNA-Seq demonstrated 60 significantly upregulated and 65 significantly downregulated genes in M2 compared with M1. Five disease-related differentially expressed genes were identified and validated by qRT-PCR analysis. Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses were performed to predict the functions of aberrantly expressed genes. Recombination activating 1 (RAG1), RAG2, BCL2-like 11, phosphatidylinositol 4,5-bisphosphate 3-kinase catalytic subunit alpha isoform and repressor element1-silencing transcription factor might play roles in MG pathogenesis involving the primary immunodeficiency signaling pathway, signaling pathways regulating pluripotency of stem cells and fork head box O signaling pathway. Conclusion: The aberrantly expressed genes of BMMCs in M1 or M2 patients demonstrate the underlying mechanisms governing the pathogenesis of MG.


INTRODUCTION
Myasthenia Gravis (MG) is a neuromuscular disorder characterized by muscle weakness and easy fatigability upon exertion. MG is mediated by autoantibodies against postsynaptic structural proteins of neuromuscular junctions: most antibodies target Acetylcholine Receptor (AChR), the less frequent Muscle-Specific Kinase (MuSK) or low-density lipoprotein receptor-related protein 4 (LPR4) [1]. Up to 60% to 80% of patients with MG have thymic hyperplasia and 21% have thymoma [2,3]. To the best of our knowledge, the thymus plays an initial role in triggering humoral immunity in MG because of either thymic hyperplasia or thymoma of lymphoproliferative origin [4][5][6][7]. However, the exact mechanisms leading to autoimmunity dysfunction and chronicity in patients with MG have not been elucidated. Despite improvements in treatment due to multiple therapies, including thymectomy, approximately 10% of patients with MG are considered treatmentrefractory [8], highlighting the need for further understanding of the pathophysiology and identification of organism biomarkers to characterize the pathogenesis of the disease and to monitor the therapeutic response.
Immune dysregulation may originate from defects in the proliferation and differentiation of bone marrow (BM)-derived Hematopoietic Stem Cells (HSCs) and/or hematopoietic support function of BM stromal cells, indicating that autoimmune diseases may primarily result from stem cell dysfunction [9][10][11]. Porta and Papadaki et al. compared the proliferation and differentiation of HSCs and hematopoietic support function of BM stromal cells between patients with Rheumatoid Arthritis (RA) and healthy controls; their results showed that the numbers of CD34+ BM cells, the proliferating capacity of BM myeloid precursors and the percentage of erythroid burst-forming units and colony-forming unit assays of granulocyte macrophage colonies significantly decreased, and the hematopoietic support function of BM stromal cells was abnormal in RA patients [12,13]. Furthermore, patients with autoimmune cytopenias secondary to systemic autoimmune disease exhibited increased apoptosis of BM progenitor cells, leading to low CD34+ cell numbers and abnormal hematopoietic supporting capacity of BM stromal cells due to aberrant, local or systemic inhibitory cytokine production or intricate interactions between hematopoietic and immune cells present within the BM microenvironment [14].
Nevertheless, the immune system emerging from the autologous HSCs of patients with MG is complex, and the pathogenesis of autoimmune dysfunction in the disease has not been fully elucidated. Over the past decade, Genome-Wide Association Studies (GWAS) and gene expression studies have provided novel insights into MG. Recently, gene expression profiles from specific cell subsets have been shown to be better determinants of immune status than whole-blood or Peripheral-Blood Mononuclear Cell (PBMC) profiles due to the diversity of leukocyte responses [15]. In this study, we investigated to explore the core transcriptional signature and crucial pathways of BM Mononuclear Cells (BMMCs) from patients with thymoma-associated MG (M1) and those with MG without thymoma (M2). A global landscape of dysregulated gene signatures was first obtained using high-throughput RNA sequencing (RNA-Seq), followed by pathway analysis and subsequent immune profiling. The findings obtained in this study may help to elucidate the heterogeneity of the disease and identify biomarkers and pathways driving disease pathogenesis.

Study patients
The cohort study was designed in two stages and included 36 patients from the outpatient department of the Department of Neurology of the First Affiliated Hospital of Guangxi Medical University. According to whether the patients had thymoma, all patients were divided into group M1 (n=18) and group M2 (n=18). For the RNA-Seq stage, we selected three patients from groups M1 and M2. For the quantitative real-time reverse transcriptionpolymerase chain reaction (qRT-PCR) step, all patients from the two groups were compared. MG was diagnosed based on the clinical evidence of muscle weakness and fatigue and at least one positive ancillary diagnostic test, including proof of pathological decrement in repetitive nerve stimulation or a positive response to a neostigmine test. Positive detection of autoantibodies against AChR in patients further supported the diagnosis. The Myasthenia Gravis Foundation of America (MGFA) [16] clinical classification was also employed to identify the MG subgroups. Regardless of the length of time over which the symptoms of MG appeared, all patients were coming to our clinic for the first time. All patients completed the collection of BM specimens prior to treatment with glucocorticoids, intravenous immunoglobulin, immunosuppression or thymectomy. The patients also had neither clinical symptoms of metabolic disease, infectious disease nor any autoimmune diseases, including thyroid disorders, Systemic Lupus Erythematosus (SLE), RA and Sjogren syndrome. Table 1 summarizes the demographic and baseline clinical data of the study subjects. All patients signed informed consent forms. The study was approved by the Medical Ethics Committees of the Institute of the First Affiliated Hospital of Guangxi Medical University.

RNA extraction and quality control
Total RNA was harvested for library construction and sequencing. In brief, total RNA was extracted from samples by using a TRIzol reagent kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer's protocol and subjected to an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and RNase free agarose gel electrophoresis for quality and purity tests, respectively. The RNA integrity number for each sample was higher than 7, indicating that the extracted RNAs possessed high quality and integrity. The ratios of Optical Density (OD) 260/280 in all samples were between 1.8 and 2.0, indicating high purity of extracted RNAs. The OD260/230 ratios in all samples ranged from 1.8 to 2.2, signifying an absence of organic extraction reagent contamination. The Q30 value of each sample was greater than 91%. Thus, the samples were determined to meet the quality and purity requirements for high-throughput sequencing.

mRNA-seq library construction and sequencing
After the total RNAs were extracted, ribosomal RNAs (rRNAs) were removed to retain mRNAs by using an oligo (dT) RiboZero Magnetic Gold Kit (Illumina, USA). The enriched mRNAs were fragmented and reverse-transcribed into cDNA. Next, the cDNA fragments were purified with a QiaQuick PCR extraction kit (Qiagen, Venlo, The Netherlands), end-repaired, subjected to poly (A) addition, and ligated to Illumina sequencing adapters. Next, uracil-N-glycosylase was utilized to digest second-strand cDNA. The digested products were size-selected using agarose gel electrophoresis, amplified by PCR, and sequenced using the Illumina HiSeq™ 4000 platform (Illumina, San Diego, CA, USA) by Gene Denovo Biotechnology Co. (Guangzhou, China).

Processing and mapping of mRNA-Seq reads
The Illumina HiSeq platform generated mRNA-Seq reads and removed the adaptors and low-quality reads (Q-value <20) bases at the 5′ and 3′ ends by FASTQ (version 0.18.0) [17]. The shortread alignment tool Bowtie2 (2.2.8) [18] was employed to map reads to the rRNA database. The rRNA mapped reads were subsequently removed. The remaining reads were further utilized in transcriptome assembly and analysis. The clean reads of each sample were then mapped to the reference genome by TopHat2 (version 2.1.1) [19].

Identification of differentially expressed mRNAs and enrichment analysis of Gene Ontology (GO) functions and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways
Scatter plots were generated, and hierarchical clustering was performed to assess and distinguish the variations in mRNA expression between the M1 and M2 groups. Aberrant mRNA expression was identified by volcano plot filtering and FC analysis. Differentially Expressed Gene (DEG) analysis was performed using the edgeR package. On the basis of this method, genes with |log2FC| ≥ 2 and a false discovery rate (FDR) <0.05 were classified as DEGs between the two groups. Next, differentially expressed coding RNAs were subjected to enrichment analysis of GO functions and KEGG pathways. An FDR <0.05 was considered to be the threshold for significant enrichment.  Table 2 shows the primers used for qRT-PCR analysis. The human glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene was employed as the reference gene for data normalization. The comparative threshold cycle (2-ΔΔCt) method of relative quantitation was employed to analyze the expression differences. Ct denotes the amplification cycles when the real-time fluorescence intensity of the reaction reached the set threshold, and when amplification was in logarithmic growth.

Statistical analysis
Data are presented as the mean ± Standard Deviation (SD) of the mean. The normal distribution of variables was investigated prior to statistical analysis using visual (histograms) and analytical methods (Kolmogorov-Smirnov/Shapiro-Wilk's test). Data were recorded and analyzed with Graph Pad Prism 8.0 (Graph Pad Software, La Jolla, CA, USA) and IBM SPSS Statistics version 24.0 (IBM Corporation, New York, USA). Differences in mRNA levels between groups were evaluated with one-way analysis of variance (ANOVA) and Pearson's chi-square test for categorical variables. P-values <0.05 were considered to be significant.     Figure 3).

GO function contains three terms: Cellular Component (CC), Molecular Function (MF) and Biological Process (BP). Figures 4 and 5 shows the GO analysis of upregulated and downregulated
DEGs, respectively. In the CC ontology, cell (GO: 0005623), cell part (GO: 0044464), organelle (GO: 0043226), and intracellular (GO: 0005662) were significantly enriched GO terms for upregulated ( Figure 4B) and downregulated genes ( Figure 5B). In the MF ontology, the most common GO terms for significantly upregulated ( Figure 4C) and significantly downregulated significant DEGs ( Figure 5C) were binding (GO: 0005488), ion binding (GO: 0043167) and metal ion binding (GO: 0046872). In the BP ontology, the upregulated genes ( Figure 4D KEGG pathway enrichment analysis indicated that the differentially expressed mRNAs were most significantly enriched in terms of signal transduction, immune system and cancer ( Figure 6).

Verification of DEGs
Five disease-related DEGs involving the primary immunodeficiency signaling pathway, signaling pathways regulating the pluripotency of stem cells and the FOXO signaling pathway were selected. These DEGs included four upregulated genes, namely, repressor element-1-silencing transcription factor (REST, also known neuronrestrictive silencer factor), BCL2-like 11 (BCL2L11, also known as BCL2 interacting mediator of cell death), Recombination Activating 1 (RAG1) and RAG2. One downregulated gene, phosphatidylinositol 4,5-bisphosphate 3-kinase catalytic subunit alpha isoform (PIK3CA), was selected for further validation by qRT-PCR analysis (Figure 7). Data regarding the expression of genes were consistent with the results obtained through RNA-Seq.

DISCUSSION
MG is a chronic autoimmune disorder with a high recurrence rate, and some cases can be refractory or life-threatening. Therefore, studying the genetic mechanism governing MG is important. Researchers have previously investigated the relationship between altered expression of mRNAs and the pathogenesis of MG.
To investigate a wide range of genes in BMMCs to dissect the physiological, metabolic and cellular processes of MG, we first analyzed the transcriptomic changes in BMMCs from thymomaassociated MG and MG without thymoma by conducting high-throughput sequencing and qRT-PCR analysis. RNA-Seq analysis showed that a significant number of mRNAs are altered in MG, indicating that they may play regulatory role in the pathophysiological processes of this disease. Precise annotation of mRNA functions remains complex. In this study, we interpreted the mRNA functions based on co-expression gene GO and KEGG pathway analyses by using pathway-related database [22] and GO databases. As shown in Table 3, RAG1 and RAG2 were associated with primary immunodeficiency and the PI3K-FOXO signaling pathway involving BPs, MFs and CCs. Recent studies [23] have highlighted the key role played by FOXO1, which has emerged as a critical target of the PI3K-FOXO pathway in the preservation of the HSC pool, regulation of progenitor commitment, and development and differentiation of B-cells by regulating its critical target genes, which include PIK3CA, early B-cell factor, RAG1, interleukin-7 receptor, activation-induced cytidine deaminase and L-selectin. RAG1 and RAG2 play important roles in the development and function of T and B cells and the immune system. In addition, the dysfunctions and hypomorphic mutations associated with these genes may result in malignancies of the hematopoietic system, primary immunodeficiency diseases and Severe Combined Immunodeficiency (SCID), Omenn Syndrome (OS) and autoimmunity [24][25][26][27]. Brauer et al. [28] investigated T-cell development and T-Cell Receptor (TCR) V(D)J recombination by using differentiation of human-induced pluripotent stem cells generated from patients (SCID versus OS) with different RAG1 mutations in vitro; these researches detected no TCR-α/δ excision circle (a marker of TCR -α/δ locus rearrangements) in SCID and OS-derived T-lineage cells, which was consistent with a pre-TCR block in T-cell development; this result elucidates important differences that explain the wide range of immunological Apoptosis plays an essential role in the development, function and homeostasis of the immune system. Experiments with transgenic and gene knockout models have shown that BCL2L11 is a critical regulator of the negative selection of B lymphocytes in BM HSCs and T lymphocytes in the thymus and periphery [29][30][31]. The phenotype of BCL2L11-deficient mice is an SLE-like autoimmune disease with an abnormal accumulation of self-reactive lymphocytes and HSCs and autoantibody production [32,33]. Furthermore, BCL2L11 knockout mice fail to delete autoreactive T cells in the thymus, which leads to enhanced self-reactive T effector cell function followed by exacerbation of kidney disease [31]. In our work, the BCL2L11 gene was significantly downregulated in M1 patients relative to M2 patients. In addition, the percentage of HSCs significantly increased in M1 patients compared with M2 patients. This result suggests that BCL2L11 may be an important regulator of homeostasis of the BM hematopoietic system, as described above. However, the correlation or coordinated regulation system of BCL2L11 with other apoptotic molecules has not been elucidated. Further studies are warranted to determine the coordinated aspects of the effect of BCL2L11 on the immune system.
REST is a zinc finger transcription factor that plays negatively regulated roles in Neuronal Stem and Progenitor Cell (NSPC) development and differentiation by recruiting transcriptional and epigenetic regulatory cofactors to repressor element-1 in the promoter regions of neuronal genes [34,35]. Recent studies [36] showed that in the T cell/neuron coculture model, activated CD4+ T cells reduced the activity of neurons and extracellular signalregulated kinase 1/2, causing REST mRNA upregulation and thereby mediating the downregulation of REST-mediated neuronal cell adhesion molecule L1, which exhibited an unexpected functions in the pathological crosstalk between the immune and nervous systems. Bergsland et al. [37] noted that blocking the function of REST led to loss of alterations in chromatin modifications and the suppression of nitric oxide-induced neuronal to glial switching. This result indicates that REST is an important factor in the NSPCspecific response to innate immunity. In our study, using RNA-Seq technology and qRT-PCR analysis, we determined that REST was upregulated in M2 patients, with the most enriched GO terms including 'negative regulation of nervous system development', 'regulation of mesenchymal stem cell differentiation', 'transcription factor binding' and 'muscle structure development' in the pathogenesis between M1 and M2 patients ( Figure 4). However, these findings suggest an unknown role of REST in the control of neuronal gene transcription in response to neuroimmunity. This unknown function may be investigated in future works involving gene transfection/knockdown, possibly helping to elucidate the pathogenesis of MG. Additional studies are required to confirm the biological function of these genes.
Our research features several limitations. First, the RNA-Seq data were limited to a small number of samples. The qRT-PCR results were acquired from a small sample consisting of 18 M1 patients and 18 M2 patients. Further studies employing a large sample size are necessary to evaluate and validate the results. Most of the conclusions of this study were also generated based on RNA-Seq data analysis. Furthermore, we annotated the functions of abnormally expressed mRNAs and were unable to validate the function of disease-related DEGs, which warrants further systematic research. In the future, we aim to further investigate and focus on the functional verification and molecular mechanisms of mRNAs in the pathological process of MG.

CONCLUSION
A total of 1876 DEGs were identified between M1 and M2 patients. RAG1, RAG2, BCL2L11, PIK3CA and REST were suggested to play important roles in the pathogenesis of MG. Nevertheless, the lack of biological experiments was a major limitation of the present study. Further research should be designed to support our results and investigate the molecular mechanisms underlying the pathological process of MG.