Blood Transcriptome Profiling in Myasthenia Gravis Patients to Assess Disease Activity: A Pilot RNA-seq Study.

Myasthenia gravis (MG) is an antibody-mediated autoimmune disease characterized by exertional weakness. There is no biomarker to reflect disease activity and guide treatment decision. Here, we reported a pilot blood transcriptome study using RNA sequencing (RNA-seq) that identified differences of 5 samples in active status and 5 in remission from 8 different patients and 2 patients provided samples for both active and remission phase. We found a total of 28 differentially expressed genes (DEGs) possibly related to disease activity (23 up-regulated and 5 down-regulated). The DEGs were enriched for the cell motion and cell migration processes in which included were ICAM1, CCL3, S100P and GAB2. The apoptosis and cell death pathway was also significantly enriched, which includes NFKBIA, ZC3H12A, TNFAIP3, and PPP1R15A. Our result suggests that transcript abundance profiles of the genes involved in cell trafficking and apoptosis may be a molecular signature of the disease activity in MG patients.

INTRODUCTION

Myasthenia gravis (MG) is an antibody-mediated autoimmune disease characterized by exertional skeletal muscle weakness, caused by auto-antibodies to components of muscle membrane at neuromuscular junction (NMJ) [1]. Acetylcholine receptor (AChR) antibodies are detected in about 85% of MG patients, and antibodies to muscle-specific tyrosine kinase (MuSK) are positive in approximately half of the generalized seronegative MG patients [2]. The list of pathogenic auto-antibodies in MG has been recently expanded with discovery of antibodies to low-density lipoprotein receptor-related protein 4 (LRP4) and agrin in some of the double-seronegative MG patients [3]. These highly specific auto-antibodies in MG have a crucial role for the diagnosis of disease, and MuSK antibody in particular seem to be useful in predicting adverse effects to acetylcholinesterase inhibitors and refractoriness to conventional immunosuppressive treatment [4]. Although the levels of these auto-antibodies tend to decrease with immunosuppressive treatment, they are highly variable between patients, and do not correlate with disease severity [5]. There is an unmet clinical need to develop biomarkers that can reflect disease activity and/or severity, and to guide treatment decisions.

Gene expression profiling on a genome-wide scale has been increasingly used to investigate pathogenesis, and also to develop potential biomarkers in various human diseases. Advances of high-throughput techniques such as DNA microarrays and recently RNA sequencing (RNA-seq) allow us to analyze gene expression in an unbiased and comprehensive way [6]. In MG, distinct gene expression signatures were found in the thymus of AChR antibody positive patients using microarrays, which include dysregulation of chemokines such as CCL21, CXCL13, CXCL10 and CXCR3 [789]. Interferon type I overexpression together with dysregulated expression of dsRNA signaling molecules was found in thymoma-associated MG patients [10]. Of note, previous gene expression studies mostly have used thymic tissues which are not readily accessible. Furthermore, it is not feasible to sample repeatedly in order to monitor disease activity and response to treatment. Since the status of human immune system may be best monitored by the changes of the composition and transcript abundance of circulating immune cells, blood transcriptome could be a good alternative target of investigation. Indeed, recent studies in autoimmune disease such as systemic lupus erythematosus, multiple sclerosis, and psoriasis have demonstrated perturbation of blood transcriptome, and identified molecular signatures to predict clinical relapse and response to treatment [111213]. Here, we report the results of blood whole transcriptome study using RNA-seq in 8 patients with MG who were in either remission or active disease states. We found 28 genes that are differentially expressed (DEGs) according to disease activity. Functional analysis suggested that genes expression profiles related to immune cell trafficking and apoptosis might be a molecular signature of the disease activity in MG.

MATERIALS AND METHODS

Patients

This study was approved by the Institutional Review Board of the Seoul Metropolitan Government Boramae Medical Center (IRB no. 16-2013-116). Written informed consent was obtained from the participants who were diagnosed as having MG and on the status of either active or remission. The active group consisted of 5 patients with de novo or refractory MG (moderate to severe symptoms despite long-term immunosuppressive treatment). Disease severity was graded according to the Myasthenia Gravis Foundation of America (MGFA) Clinical Classification [14]. Remission was defined by the MGFA post-intervention status, and included complete stable remission (CSR), pharmacologic remission (PR), and minimal manifestation (MM) (Table 1). Two patients provided samples at different time points, one during active disease status and the other in remission state. There were not statistically differences of mean age (p=0.69), disease duration (p=0.31), and AChR antibody titer (p=0.69) between 2 groups.

Table 1

Demographics of study population

*Patient 1 and 6 are same patients.

†Patient 2 and 7 are same patients.

Abbreviations: MG: myasthenia gravis; AChR Ab: acetylcholine receptor antibody; MGFA: Myasthenia Gravis Foundation of America; M: male; F: female; NL: normal; Pd: prednisolone; MM: minimal manifestation; PR: pharmacologic remission.

PBMC isolation and RNA purification

For isolation of peripheral blood mononuclear cells (PBMC), the Lymphoprep™ was used according to the manufacturer's protocol (Axis-shield, Oslo, Norway). Medium was placed in the tube, and then blood sample diluted with saline with 1:1 was added. After centrifugation for 20 minutes at 600×g , sedimented PBMCs were harvested. RNA purification was performed with the RNeasy Mini kit with the isolated PBMC sample (Qiagen, Seoul, Korea). The cell pellet was mixed with RLT buffer and 70% ethanol. The lysate was then loaded onto the RNeasy Mini spin column to facilitate the binding of RNA to the column and for the removal of contaminants. DNase was added to remove residual DNA efficiently.

RNA-Seq

The mRNA-Seq sample was obtained using Illumina TruSeq™ RNA Sample Preparation Kit (Illumina, Inc., San Diego, CA, USA). In brief, purifying the poly-A containing mRNA molecules with poly-T oligo-attached magnetic beads was the first step, followed by thermal mRNA fragmentation. The RNA fragments were then transcribed into first strand cDNA using reverse transcriptase and random primers. The cDNA was synthesized to second strand cDNA using DNA Polymerase I and RNase H. After the end repair process, single 'A' bases were added to the fragments and adapters were then ligated, preparing cDNA for hybridization onto a flow cell. Finally, the products were purified and enriched with PCR to create the cDNA library (Macrogen, Seoul, Korea).

Aligning RNA-Seq reads and abundance estimation

Fragmented cDNAs were aligned using TopHat v.2.0.11 [15] and subsequently aligned with sequences obtained from the human genome (UCSC version hg19) using the Bowtie 2.1.0 algorithm [16]. Abundance of aligned reads were estimated by Cufflinks v.2.1.1 [17], which accepted aligned reads and assembled the alignments into a simple and clear set of transcripts. Next, RNA-seq fragment counts were measured by the unit of fragments per kilobase of exon per million fragments mapped (FPKM) [18]. DESeq, another tool for DEG analysis, was used to compare the results with Cuffdiff analysis. Cuffdiff determines differential expression using t-test from FPKM values and is based on beta negative binomial model [19], while DESeq uses exact test based on negative binomial model [20]. We compared the results from Cuffdiff and DESeq analyses, and took the intersection of them for downstream pathway analysis.

Statistical analysis

For DEG analysis, the values of log2 (FPKM+1) were calculated, and these were normalized by quantile normalization. p-values were obtained by t-test between the active and remission groups, and fold changes were calculated with the mean log2 (FPKM+1) values, gene by gene. All data analysis of DEG was conducted using R 2.14.1 (http://www.r-proj ect.org). To segregate the samples according to the disease activity, a multi-dimensional scaling (MDS) analysis was done.

Pathway analysis using DAVID and IPA

For functional enrichment analysis using gene ontology (GO), the Database for Annotation, Visualization and Integrated Discovery (DAVID v.6.7) was used. The list of commonly detected genes both in Cuffdiff and DESeq analysis was uploaded via the web interface (http://david.abcc.ncifcrf.gov), and the background was designated as Homo sapiens [21]. Functional annotation clusters were selected using the cutoff threshold enrichment score 1.3 which is equivalent to a non-log scale score of 0.05. Gene interaction was visualized by Gene Multiple Association Network Integration Algorithm (GeneMANIA) (http://www.genemania.org/). Another functional analysis was performed with the Ingenuity Pathway Analysis (IPA) software (licensed use of Ingenuity Systems, http://www.qiagen.com/ingenuity).

RESULTS

Differentially expressed genes (DEG)

In total, there were 48,385 transcripts, and we excluded any transcripts with an FPKM value of '0,' leaving 10,640 transcripts to be analyzed. MDS analysis showed that our MG samples could be differentiated according to the disease activity along the first plot dimension (Fig. 1).

Fig. 1

Multidimensional Scaling plot of gene expression profiles. The first plot dimension roughly corresponds to the disease activity. Paired samples are linked by dotted arrows.

The level of expression of 98 genes was significantly different between the two groups (fold change≥2, p-value<0.05) with 63 genes up-regulated and 35 genes down-regulated in the remission group (Supplementary Table 1, Fig. 2). The number of DEGs derived with DESeq analysis was greater with 165 up-regulated and 127 down-regulated in the remission group (Supplementary Table 2). Twenty-eight genes were common in the results of both Cuffdiff and DEseq analyses (23 up-regulated and 5 down-regulated genes, Table 2, Fig. 3). When only the two paired samples were analyzed separately, we found 18 up-regulated genes (fold change≥2, p-value<0.05), with no down-regulated genes (Supplementary Table 3). All the up-regulated genes from the paired sample analysis except for TPGS1 were contained in the list of DEGs analyzed with the whole study samples.

Fig. 2

Differentially expressed genes by Cuffdiff. (A) Volcano plot illustrating the differential expression levels of genes of the active and remission group. Genes that are significantly up- and down-regulated in the active compared to remission group are shown in red and green dots, respectively. (B) Heat map of the hierarchical clustering based on 98 differentially expressed genes (fold change≥2, p-value<0.05).

Table 2

List of differentially expressed genes (from both Cuffdiff and DESeq analyses)

Fig. 3

Venn diagram of DEGs from DESeq and Cuffdiff. (A) Up-regulated genes in the remission group. (B) Down-regulated genes in the remission group.

DAVID analysis with up-regulated genes revealed two enriched GO functional annotation clusters (Table 3). The first cluster consisted of biological functions of cell motion and cell migration. The second cluster consisted of apoptosis and cell death. Number of down-regulated genes was not sufficient to analyze. Gene interaction of up-regulated genes was shown in Fig. 4.

Table 3

Functional annotation cluster analysis using DAVID

*Gene Ontology.

Fig. 4

Gene interactions of up-regulated genes from GeneMANIA. Blue line indicates co-localization of the genes and predicted functional relationships between genes are indicated in orange. Circles filled with black signify the commonly detected genes by Cuffdiff and DESeq. Circles filled with gray indicate their interaction and added by GeneMania.

According to the downstream effect analysis performed by using IPA, top 5 functions affected were hematologic system development and function, immune cell trafficking, cellular movement, inflammatory response and cell-to-cell signaling and interaction pathways. Most of the genes related to these categories were down-regulated in active group.

DISCUSSION

In the present study, whole transcriptome RNA sequencing was performed in 10 PBMC samples obtained from eight MG patients to investigate systemic changes of gene expression possibly relevant to disease activity. Because there is no general consensus regarding the best method for the differential expression analysis of RNA-seq data, two different analysis tools (Cuffdiff and DESeq) were used separately with the intersection of the results being taken for downstream pathway analysis. The results demonstrated that among 10,640 transcripts only 28 genes were differentially expressed with 23 being up-regulated and 5 down-regulated in remission compared to active status of the disease. Interestingly, major functional themes of the differentially expressed genes include differentiation, trafficking and apoptosis of immune cells, inflammatory response and cell-to-cell signaling.

S100B (S100 calcium-binding protein B), which was the most down-regulated gene in the remission group, is a ligand of the receptor for advanced glycation end products (RAGE). In the experimental autoimmune myasthenia gravis (EAMG) model, S100B levels were significantly higher [22]. Mean value of S100B of MG was higher than the normal controls, but it was not statistically significant [23]. CTTN (cortactin) was another down-regulated gene in the remission group, which is needed for the formation of the AChR cluster and antibody against cortactin was reported in MG patients [24]. ABL1 (Abelson murine leukemia viral oncogene homolog 1) is up-regulated in the remission group, which is critical mediator of postsynaptic assembly at the NMJ via providing a specific tyrosine kinase activity downstream of the MuSK receptor that is required for agrin-induced AChR clustering [25].

One of the significant functional annotation cluster in DAVID analysis was the apoptosis and programmed cell death. Immunosuppressive agents are the mainstay of treatment for MG and their mechanisms are related with these pathways. Corticosteroids induce T lymphocyte apoptosis, and other agents such as azathioprine also inhibit cell proliferation [26]. All of the remission group patients received prednisolone and three patients also received tacrolimus, another T cell apoptosis inducer [27]. CXCL13 and CCL21 were overexpressed in hyperplastic thymus [828] and CXCR5 was overexpressed on T cells of MG patients [29]. Most of the genes associated with chemotaxis, however, were down-regulated in the active group of the present study. This discrepancy needs to be confirmed in further studies, but might be explained by the followings. First, the previous studies had compared between MG patients and healthy controls. Second, difference of the target samples (thymus or neuromuscular junction in the previous studies) might account for the discrepancy. In a previous study of rheumatoid arthritis, chemokine receptor profile of PBMCs was increased with treatment, which the authors explained as systemic compensation for the changes in the inflamed tissue [30].

There are several limitations to be acknowledged in this pilot study. Sample size was too small to make a robust conclusion. Immunosuppressive treatment varied among patients, which might affect the results. Despite these limitations, to our knowledge, this is the first study investigating the systemic changes of blood whole transcriptome according to the disease activity in MG patients. Genes that are involved in the biological process of immune cell trafficking and apoptosis might be promising candidates for biomarkers, and warrant further studies.