LncRNA and mRNA expression associated with myasthenia gravis in patients with thymoma.

OBJECTIVE: Pathological alterations of the thymus are observed in the majority of patients with myasthenia gravis (MG). To explore the potential mechanisms of these alterations, we performed a transcriptome analysis and measured co-expression of aberrant long non-coding RNAs (lncRNAs) and messenger RNAs (mRNAs).
METHODS: RNA was extracted from eight patients with thymoma, five of whom had MG. Transcriptome profiles were acquired through mRNA and lncRNA microarray analysis. Quantitative reverse transcription polymerase chain reaction was used to verify the results of the microarray analysis. LncRNAs co-expressed with mRNA were analyzed with Pearson's coefficient. Next, cis-regulated and trans-regulated target genes were predicted. The functions of aberrant lncRNAs were explored on the basis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of target mRNAs.
RESULTS: The comparative microarray analysis identified 4360 lncRNAs and 2545 mRNAs with significant differential expression. The most significant GO enrichment terms were phosphoric ester hydrolase activity, phosphatase activity, and hydrolase activity, which were assigned as molecular functions. Regulation of endosome size was the most significant GO enrichment term assigned as a biological process, and Golgi apparatus was the most significant GO enrichment term assigned as cellular component. The reliability prediction terms of KEGG included calcium signaling pathway, glycosphingolipid biosynthesis, and caffeine metabolism.
CONCLUSION: MG-positive thymoma is associated with overactive biological processes and molecular functions, especially dephosphorylation and hydrolysis, which may affect thymocyte survival during selection in the thymus.

INTRODUCTION

Myasthenia gravis (MG) is an autoimmune disease caused by antibodies to acetylcholine receptor neuromuscular transmission. Pathological alterations of the thymus are observed in the majority of patients with MG, of which thymomas constitute ~10%–15%. Thymomas are common anterior mediastinal tumors derived from thymic epithelial cells, and it has been reported that 20% of thymomas are accompanied by MG. , However, the pathogenesis of MG with thymoma is not clear. Increasing evidence suggests that thymoma tissue can generate and export large numbers of CD4+/CD8+ T cells, which potentially stimulate autoantibody production and are related to MG progression. , However, the pathogenesis is still not clear.

Long non‐coding RNAs (lncRNAs) are a type of RNA molecule that consists of more than 200 nucleotides, which contribute to various biological cascades. , LncRNAs are involved in certain autoimmune diseases, including MG. , , However, few studies have compared lncRNA expression in patients with thymoma with and without MG. To explore the potential mechanism of thymoma with MG, this study investigated transcriptome profiles in patients with thymoma and MG [MG(+)] compared with patients with thymoma without MG [MG(−)], and performed target prediction and enrichment Gene Ontology (GO) analyses by calculating the co‐expression relationship of aberrant lncRNAs and messenger RNAs (mRNAs).

METHODS

Ethics

The Institutional Review Board at BeijingTongren Hospital approved this study, and all experiments were performed in accordance with relevant guidelines and regulations.

Clinicopathological patient features

All patients were pathologically diagnosed with thymoma. Patients were grouped into different subtypes according to their pathological features (Table 1).

TABLE 1

Baseline characteristics

Sample	Age (years)	Sex	MG	WHO thymoma subtype	Past disease
1	49	Female	IIA	AB	None
2	37	Male	IIB	B2/B3	None
3	38	Female	I	AB	None
4	69	Female	IIB	AB	HTN/DM
5	69	Male	I	B1	HTN
6	52	Female	None	B2	None
7	63	Female	None	AB	HTN
8	64	Male	None	B2	None

Abbreviations: DM, diabetes mellitus; HTN, hypertension; MG, myasthenia gravis; WHO, World Health Organization.

RNA extraction and microarray analysis

According to the manufacturer's protocol, total RNA was extracted from thymoma tissue using TRIzol reagent (Invitrogen) and purified using the miRNeasy Mini Kit (QIAGEN). Approximately 200 ng of total RNA from each thymoma sample was used for lncRNA/mRNA microarray analysis. LncRNA and mRNA expression was analyzed by lncRNA + mRNA Human Gene Expression Microarray V4.0 (CapitalBio Technology), which contained 40 916 lncRNAs and 34 235 mRNAs collected from almost all authoritative databases, including Agilent_ncRNA, ncRNA Expression Database, H‐invDB, GencodeV13, NONCODE v3.0, RefSeq, LNCipedia, UCR and UCSC_lncRNAs Transcripts, and lncRNAdb. Microarray images were analyzed and produced raw data using Agilent Feature Extraction (Agilent Technologies). Further analysis of raw data was performed using Agilent GeneSpring (Agilent Technologies). LncRNA/mRNA microarray experiments were performed by CapitalBio Technology Corporation, Beijing, People's Republic of China.

LncRNA co‐expression analysis and gene function annotation

Hierarchical clustering was used to show the distinguishable RNA expression patterns among samples, especially between different groups. Significantly aberrant lncRNAs and mRNAs were identified by Volcano plot filtering. GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were performed to annotate the functions of differentially expressed mRNAs. With Pearson's coefficient, the co‐expression relationship between lncRNAs and mRNAs was shown by a coding‐non‐coding gene co‐expression network. Furthermore, cis‐regulated and trans‐regulated target genes were predicted based on the co‐expression analysis. Next, the functions of aberrant lncRNAs were explored on the basis of GO and KEGG pathway analyses of target mRNAs. These analyses were performed by CapitalBio Technology Corporation.

Statistical analysis

Statistical analysis was performed using SPSS version 20.0 (IBM). LncRNA expression in verified samples is presented as mean ± standard deviation. Differences in lncRNA expression between groups were compared using the paired‐samples t‐test. Correlations between verified lncRNAs and clinicopathological characteristics were determined using Spearman's correlation analysis. All p‐values were two‐sided, and a p‐value of <0.05 was considered statistically significant.

RESULTS

Clinical information of samples for microarray analysis

Eight thymoma samples were selected for microarray analysis. The baseline characteristics of these samples are detailed in Table 1.

Microarray analysis overview

The lncRNA and mRNA expression profiles were presented in cluster and principal component analyses (Figure 1). These data validated that the samples from different groups were representative and appropriate for the experiment. Furthermore, images from the lncRNA analysis showed better concentricity, whereas images from the mRNA analysis showed a degree of discreteness. Collectively, these findings indicated that the results of the experiment were reliable.

FIGURE 1

Plots of microarray analysis data. Hierarchical cluster analysis plots of the expression of (a) long non‐coding RNAs (lncRNAs) and (b) messenger RNAs (mRNAs) in thymomas of patients with and without myasthenia gravis. Red denotes high relative expression, and green denotes low relative expression. Each RNA is represented by a single row of colored boxes, and each sample is represented by a single column. The three‐dimensional principal component analysis plots for normalized data of (c) lncRNAs and (d) mRNAs, respectively, show good similarity among groups

The comparative microarray analysis identified 4360 lncRNAs and 2545 mRNAs with significant differential expression (fold change ≥2.0, p < 0.05). Compared with MG (−) specimens, 3699 lncRNAs and 1331 mRNAs were upregulated, whereas 661 lncRNAs and 1214 mRNAs were downregulated in MG(+) patients, respectively. The results are presented in cluster and scatter plots (Figure 1).

Microarray analysis validation

To validate the disruption of lncRNA and mRNA expression in MG(+) patients, we selected the top four upregulated and downregulated lncRNAs and mRNAs, as well as four terms with the closest relationship to phosphorylation/dephosphorylation, for quantitative reverse transcription polymerase chain reaction. The expression of four lncRNAs (XLOC_006297, ENSG000000218510.3, XLOC_007052, and XLOC_002588) and four mRNAs (PRRSS56, PPP3R1, DUSP26, and PLCG1) was consistent with the results of the microarray analysis (Figure 2).

FIGURE 2

Comparison of long non‐coding RNA (lncRNA) and messenger RNA (mRNA) expression in patients with thymoma with myasthenia gravis [MG(+)] and patients with thymoma without MG [MG(−)] by real‐time polymerase chain reaction (PCR). The selected lncRNAs and mRNAs were validated by real‐time PCR of RNAs extracted from five MG(+) thymomas and 3 MG(−) thymomas. Each sample was analyzed in triplicate. Column heights represent mean fold changes in expression in the MG(−) group. The real‐time PCR results are consistent with microarray data

GO and pathway analyses of mRNAs

The GO analysis results showed that the main terms were biological processes in which dysregulated mRNAs were involved. Some terms were in cellular components, and the fewest were in molecular functions (Figure 3). Among the involved biological processes (Figure 4(a)), lipid phosphorylation, 4‐hydroxyproline metabolism, regulation of cell projection organization, cellular component biogenesis, and cellular component biogenesis were the five most significant processes associated with dysregulated mRNAs. Among the cellular components (Figure 4(b)), the troponin complex, striated muscle thin filaments, and myofilaments were the three most significant cellular components associated with dysregulated mRNAs. Among molecular functions (Figure 4(c)), RNA polymerase II‐activating transcription factor binding, ferric iron binding 3‐galactosyl‐N‐acetylglucosaminide, and 4‐alpha‐L‐fucosyltransferase activity were the three most significant molecular functions enriched by dysregulated mRNAs.

FIGURE 3

Degree of mRNA enrichment by function. The blue, green, and red bars represent gene ontology terms of biological processes, cellular components, and molecular functions, respectively

FIGURE 4

Top 30 Gene Ontology terms for differential messenger RNA (mRNA) genes between myasthenia gravis (MG)‐positive thymomas and MG‐negative thymomas. (a) Statistics of biological process enrichment; (b) statistics of cellular component enrichment; (c) statistics of molecular function enrichment

In the KEGG pathway analysis, among the top 30 significantly enriched KEGG pathways, arginine biosynthesis, glycosphingolipid biosynthesis, and sphingolipid metabolism were the most significant pathways involved in the pathogenesis of thymoma triggering MG and were associated with dysregulated mRNAs. The pathways that enriched the most dysregulated mRNAs were the calcium signaling pathway and purine metabolism (Figure 5).

FIGURE 5

The top 30 significantly enriched Kyoto Encyclopedia of Genes And Genomes pathways of messenger mRNAs (mRNAs) according to p‐values and gene number. Arginine biosynthesis, glycosphingolipid biosynthesis, and sphingolipid metabolism were the most significant pathways associated with dysregulated mRNAs. Pathways that enriched the most dysregulated mRNAs were the calcium signaling pathway and purine metabolism

LncRNA function annotation

To further understand the function of dysregulated lncRNAs between MG(+) thymomas and MG(−) thymomas, a co‐expression analysis was used to predict target mRNAs. GO and KEGG pathway analyses were used to annotate lncRNA function. According to p‐values and enrichment, there were a total of 30 GO enrichment terms (Figure 6) and 30 KEGG pathway items (Figure 7).

FIGURE 6

Top 30 Gene Ontology (GO) terms for the difference in messenger RNA (mRNA) genes between myasthenia gravis (MG)‐positive thymomas and MG‐negative thymomas. The most significant GO enrichment terms were phosphoric ester hydrolase activity, phosphatase activity, and hydrolase activity. (a) Statistics of biological process enrichment; (b) statistics of cellular component enrichment; (c) statistics of molecular function enrichment

FIGURE 7

The top 30 significantly enriched Kyoto Encyclopedia Of Genes And Genomes pathways of long non‐coding RNAs (lncRNAs) according to p‐values and gene number. Calcium signaling was the most significant pathway associated with dysregulated lncRNA/messenger RNAs

The most significant GO enrichment terms were phosphoric ester hydrolase activity (GO:0042578), phosphatase activity (GO:0016791), and hydrolase activity (GO:0016788), which were assigned as molecular functions (Figure 6(a)). All of these terms were associated with catalytic activity (GO: 0003824). Regulation of endosome size (GO: 0051036) and Golgi apparatus (GO:0005794) were the most significant GO enrichment terms assigned as biological processes (Figure 6(b)) and cellular components (Figure 6(c)), respectively. The reliability prediction terms of KEGG pathways (Figure 7) included the calcium signaling pathway (hsa04020), glycosphingolipid biosynthesis (hsa00601), and caffeine metabolism (hsa00232). These results from the GO and KEGG pathway analyses showed the likelihood that lncRNAs regulate the balance between phosphorylation and dephosphorylation in MG(+) thymomas.

DISCUSSION

Thymoma is a common type of anterior mediastinal tumor. It is sometimes detected in patients with paraneoplastic syndromes, such as MG. Furthermore, most patients with MG have an abnormal thymus, with ~10%–15% having a thymoma. Over the years, it has been suggested that thymoma and MG are closely linked. Despite decades of gene studies having been conducted to elucidate the pathogenesis of MG triggered by thymoma, the results are inconsistent. This is largely because of differences in research samples. The associations between MG and thymus were mostly explored, by research on blood samples in the early stage, such as peripheral blood mononuclear cells (PBMCs). Camilla et al. first showed that the thymus exports excess numbers of CD4+ and CD8+ T cells to peripheral blood in patients with MG and thymoma. This research supported the hypothesis that thymoma tissue itself can generate and export mature, long‐lived T cells and that these T cells reflect the thymic pathology and are likely to be related to the associated autoimmune diseases.

Using microarray assays of PBMCs from patients with MG with or without thymoma, Luo et al. found that aberrant lncRNAs regulate the expression of inflammatory cytokines. Upregulated mRNAs, such as interleukin (IL)‐8, CXCL1, CXCL3, and CXCL5, were observed among the oebiotech_11658, oebiotech_12721, oebiotech_21831, oebiotech_11933, and oebiotech_22652 “cis” target genes. Zhang et al. showed that lncRNA ENSG00000267280.1 and lncRNA ENSG00000235138.1 are associated with regulation of the inflammatory response and IL‐1 production. The most elevated lncRNA (XLOC 003810) in MG(+) thymoma patients was proven to modulate thymic T helper 17 cell/regulatory T cell balance.

Because thymectomy is a major theme in the study of MG, there are an increasing number of studies on thymomas and thymuses in patients with MG. Comparing MG+ thymomas with MG+ thymuses, several oncogenes, including PLK5, HMGA2, and REG4, are notably upregulated. Chemokines are rare aberrant genes in patients with MG(+) thymoma. Scarpino et al. found that CCL17, CCL22, and CCL19 are expressed significantly less in MG(+) thymomas compared with MG(−) thymomas. Using immunohistochemistry techniques, Zhang et al. showed CXCL13 overexpression only in generalized MG patients with type AB thymoma compared with MG(−) thymomas. In our previous research, we found that inflammatory chemokines, such as PNISR, CCL25, NBPF14, PIK3IP1, and RTCA, were significantly upregulated in MG(+) thymomas compared with MG(−) thymomas.

Based on the results of existing research, the dysfunction caused by thymoma was one of the pathogeny of paraneoplastic MG for thymoma patients. The results of blood sampling mostly show abnormal chemokines and cytokines, which may represent potential biomarkers for the diagnosis of MG. Research comparing pathological thymoma with the normal thymus mostly indicates the cause of thymoma instead of the mechanism of paraneoplastic MG. In our opinion, genome sequencing of MG(+) thymomas in comparison with MG(−) thymomas may help explore the direct relationship between thymoma and MG.

A large number of significantly anomalous lncRNAs and mRNAs were identified in the present study. Interestingly, differences were observed compared with other studies, which illustrates the complexity of MG triggered by thymoma. To verify the reliability of the microarray analysis, we chose several terms for PCR verification, including mRNAs and lncRNAs with the most obvious upregulation and downregulation in the microarray analysis (PRRSS56, PPP3R1, XLOC_006297, and ENSG000000218510.3), as well as those with the closest relationship to phosphorylation/dephosphorylation in the GO analysis (DUSP26, PLCG1, and XLOC_007052,XLOC_002588). The trends captured were broadly portable, but the regulation range was somewhat reduced.

The balance between phosphorylation and dephosphorylation, which is exerted by protein kinase and phosphatase, plays an important role in signal transduction pathways associated with protein–protein interactions. The Golgi apparatus, which is the most important structure for protein transport and modification, was the most significant GO enrichment cellular component term. As previously mentioned, long‐lived CD4+ and CD8+ T cells exported from thymomas potentially stimulate autoantibody production and subsequent MG. Coincidentally, phosphatases and dephosphorylation are involved in the regulation of T cell signaling and activation. Zheng et al. proved that protein phosphatase 2A plays an essential role in promoting thymocyte survival during selection. In the present study, the genes involved in dephosphorylation, including PLCG, LKB1, and DUSP26, among others, were regarded as key components of positive thymocyte selection and are closely related with subsequent autoimmune disease. , , ,

In conclusion, the present findings confirm that lncRNA and mRNA expression in MG(+) thymomas differ significantly compared with MG(−) thymomas. The microarray analysis showed that MG(+) thymomas are associated with certain overactive biological processes and molecular functions, especially dephosphorylation and hydrolysis. These abnormalities, which eventually cause MG, may affect thymocyte survival during selection through paracrine function. Although this study only produced some primary results, it does outline a potential mechanism for thymoma‐associated MG.