Identification of Hub Genes and Biological Pathways in Inclusion Body Myositis Using Bioinformatics Analysis

Introduction

Inclusion body myositis (IBM) is an acquired inflammatory myopathy. Disease onset is gradual and usually begins after the age of 50 years.¹ IBM is characterized by asymmetrical muscular weakness and atrophy, mostly affecting quadriceps femoris and deep finger flexors.² IBM is a rare sporadic disease with an average prevalence of 24.8–45.6/1,000,000.^1,3 IBM has late-onset, asymmetric muscle weakness and progresses slowly, compared with other inflammatory myopathies, such as dermatomyositis and polymyositis. Besides inflammatory infiltrates, the muscle pathology of IBM is associated with rimmed vacuoles and protein aggregation. Unlike patients with other inflammatory myopathy, glucocorticoid, azathioprine, methotrexate, and other immunosuppressants or immunomodulators rarely benefit IBM patients.⁴ Moreover, there is no effective treatment for IBM patients due to the unclear etiology and pathogenesis. Previous studies have demonstrated that inflammatory and degenerative pathways played important roles in IBM.² Inflammatory infiltrate, abnormal protein deposition, and mitochondrial dysfunction occur in the muscle of IBM patients. Therefore, the precise molecular mechanisms and IBM pathways should be assessed to develop effective IBM treatment strategies.

Tremendous advances have been achieved in the genetic analysis of IBM in recent years due to the development of chip techniques and sequencing approaches. The international Myositis Genetics Consortium (MYOGEN) assessed 252 IBM patients and showed that HLA-DRB1*03:01(p = 5.77×10⁻³⁴), HLA-DRB1*01:01(p = 1.57×10⁻¹⁶), and HLA-DRB1*13:01(p = 3.8×10⁻⁸) are associated with a higher risk of IBM.⁵ Although HLA-DRB1*03:01 is associated with dermatomyositis, polymyositis and polymyositis, it is not the specific risk factor for IBM. Therefore, it is necessary to identify other indicators for the diagnosis and prognostic assessment of IBM patients.

Microarray technology has been widely used to study various diseases since it can simultaneously monitor several gene expression profiles. However, it is difficult to identify reliable key genes and pathways using microarray technology due to the false-positive rates of independent microarray studies. Herein, two microarray datasets of healthy and IBM samples, GSE3112 and GSE128470 were downloaded from the Gene Expression Omnibus (GEO) database. Bioinformatics methods were used to identify differentially expressed genes (DEGs). The DEGs were subsequently analyzed using PPI (protein-protein interaction), GO (Gene Ontology), and KEGG (Kyoto Encyclopaedia of Genes and Genomes). This research provides a preliminary understanding of the molecular pathogenesis of IBM, and several potential therapeutic targets for IBM treatment.

Materials and Methods

Microarray Data

The publicly available microarray data were obtained from GSE3112 and GSE128470 in GEO (ncbi.nlm.nih.gov/geo/). The expression profile GSE3112 had 23 IBM samples and 11 healthy samples, while the GSE128470 profile contained 12 healthy samples, 26 IBM samples, 12 DM (dermatomyositis) samples, six NM (necrotizing myopathy) samples, 14 NS (non-specific myositis) samples, and seven PM (polymyositis) samples. The sample information of IBMs and normal muscle tissues in the two datasets were showed in Table S1 (GSE3112) and Table S2 (GSE128470).

Identification of DEGs and Enrichment Analysis

GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) was used to identify DEGs between normal muscle tissues and IBM. The screening criteria for DEG determination were: | log2 fold change (log2FC) | ≥1.5 and P-value <0.5. The intersection of DEGs in the two datasets was conducted for further investigation. KEGG, GO, and GSEA (Gene set enrichment analysis) enrichment were performed using clusterProfiler R package to explore the functions of the DEGs.

Protein-Protein Interaction (PPI) Network and Subcluster Analysis

The protein-protein interactions of DEGs were predicted using the STRING online platform (http://string-db.org/). Molecular Complex Detection (MCODE) within Cytoscape was used to explore the gene clusters in the PPI network.⁶

Expression and ROC Analysis

The expression profiles of IBM and healthy muscle tissues were used to identify 10 DEGs with the highest scores. The Receiver operating characteristic (ROC) curves were generated using pROC R package. The ROC curves were used to evaluate the predictive power of the identified biomarkers.

Sample Acquisition and Immunohistochemistry

One normal muscle tissue and three IBM muscle tissues were obtained from the Second Hospital of Hebei Medical University, which was approved by Research Ethics Committee of the Second Hospital of Hebei Medical University. IBM diagnosis was conducted based on the diagnostic clinicopathological features.

The muscle samples were frozen, 5-μm-thick cryosections were cut under −20 °C. A citrate buffer (pH = 6.0) was used for antigen retrieval. Then rinsed twice using PBS for 5 minutes each time. The sections were then incubated with 10% goat serum at room temperature for 30 minutes to block non-specific binding sites. The goat serum was removed, then the samples were incubated with CD4 (Abcam ab183685), ITGAM (Abcam ab287572), and CCR5 (Abcam ab287959) at 4 °C overnight. The samples were then washed thrice (3 minutes each) using PBS. The samples were then incubated with a secondary antibody at 37 °C for 1 h. The samples were washed thrice with PBS, then each section was stained with DAB for visualization using the microscope.

Immune Infiltration Analysis

The deconvolution approach CIBERSORT has been widely used to study various diseases, including ulcerative colitis,⁷ lupus,⁸ and rheumatoid arthritis, due to its validity and reliability.⁹ The white blood cell gene bio-signature matrix in CIBERSORT comprises 547 genes, known as LM22, which is used to differentiate the 22 immune cell types consisting of regulatory T cells (Tregs), naive CD4 T cells, macrophages M2, activated memory CD4 T cells, plasma cells, follicular helper T cells, gamma delta T cells, CD8 T cells, macrophages M1, memory B cells, activated NK cells, macrophages M0, resting NK cells, monocytes, resting mast cells, resting memory CD4 T cells, activated mast cells, resting dendritic cells, neutrophils, naive B cells, activated dendritic cells, and eosinophils.¹⁰ Herein, CIBERSORT was used to determine the proportions of the 22 immune cells in the IBM and normal samples of GSE3112.

Gene Set Variation and Enrichment Analysis

GSVA (Gene set variation analysis) is a non-parametric and unsupervised tool used to estimate certain pathways or biological functions via the sample expression data.¹¹ Herein GSVA analysis was performed using “h.all.v7.2.symbols” gene sets (obtained from the MSigDB (Molecular Signatures Databases)) to explore the differences in biological processes between IBM and normal muscle tissues. Gene set enrichment analysis (GSEA) was used to validate the characteristics and the potential activated pathways of IBM in GSVA results.

Connectivity Map Analysis

Connectivity map (CMap)¹² was used to explore potential drugs for IBM. Enrichment score (>0.8) and P-value (< 0.01) were considered the screening standard.

Results

DEGs Identification

A total of 466 DEGs were obtained from GSE3112 (nine down-regulated and 457 up-regulated) (Figure 1A), and 310 DEGs were obtained from GSE128470 (23 down-regulated and 287 up-regulated) (Figure 1B). The DEGs could clearly distinguish between healthy and IBM samples, as shown in the heatmaps obtained through cluster analysis (Figure 1D and E). A Venn diagram identified 219 overlapping genes between GSE3112 and GSE128470, which were used for further analysis (Figure 1C).

Figure 1 DEGs identification (A and B) Volcanic maps showing DEGs from GSE3112 and GSE128470. Red, blue, and gray dots represent significantly up-regulated, downregulated, and non-significantly expressed DEGs, respectively (C) DEGs between GSE3112 and GSE128470. The two datasets had an overlap of 219 genes. (D and E) Heatmaps showing DEGs from GSE3112 and GSE128470.

PPI Network and Pathway Enrichment Analysis

The PPI analysis results conducted using STRING tools are shown (Figure 2A). The MCODE method identified four subclusters for the PPI networks (Figure 2B–E). GO analysis showed that the DEGs were mainly enriched in leukocyte migration, regulation of immune effector process, endocytic vesicle membrane, external side of the plasma membrane, chemokine receptor binding, and peptide antigen binding (Figure 2F). KEGG analysis showed that the DEGs were significantly enriched in phagosome, chemokine, Natural killer cell-mediated cytotoxicity, and Th1, Th2, and Th17 cell differentiation signaling pathways (Figure 2G).

Figure 2 PPI network and enrichment analyses (A) The PPI network of DEGs constructed using Cytoscape. (B–E) Four clusters identified via MCODE plugin. (F) GO enrichment analysis of DEGs. (G) KEGG enrichment analysis of DEGs.

Moreover, CD4 (degree=79), PTPRC (degree=70), ITGAM (degree=58), TYROBP (degree=55), CCR5 (degree=53), IRF8 (degree=51), STAT1 (degree=50), ITGB2 (degree=49), HLA-DRA (degree=47), and CXCL10 (degree=47) were screened as hub genes of IBM.

Differential Gene Expression and ROC Analyses

In this study, the GSE3112 and GSE128470 gene data sets were used to detect the 10 hub genes The hub genes were significantly up-regulated in IBM (P<0.05) (Figure 3A and B). The ROC curves demonstrated that these hub genes are potential biomarkers for IBM diagnosis (Figure 3C and D). Further analysis suggested that the hub genes were positively correlated with each other (Figure 3E and F).

Figure 3 Expression and ROC analysis (A) Expression of 10 hub genes between IBM and normal tissues in GSE128470 and (B) GSE3112. (C) ROC analysis in GSE128470 and (D) GSE3112. (E and F) The hub genes positively correlated with each other in GSE128470 and GSE3112. Red and blue represent high low correlations.

The expression level of the hub genes in IBM and other idiopathic inflammatory myopathies, including dermatomyositis (DM), necrotizing myopathy (NM), nonspecific myositis (NS), and polymyositis (PM) were then determined based on the GSE128470 cohort. The expression of several hub genes (CCR5, HLA-DRA, and PTPRC) could distinguish IBM and DM, NM, NS. Moreover, the expression of the hub genes was not significantly different between IBM and PM (Figure 4).

Figure 4 Expression level of the 10 hub genes in IBM and other idiopathic inflammatory myopathies ns: no significance; *P < 0.05; **P <0.01; ***P < 0.001.

Immunohistochemistry

Previous studies have reported that CD4 and CCR5 play an important role in IBM.¹³ ITGAM is associated with the inflammatory response.¹⁴ Immunohistochemistry showed that CD4, CCR5, and ITGAM were highly expressed in IBM than in normal muscle tissues (Figure 5). It was in keeping with the results of the biomarker analysis.

Figure 5 Immunohistochemistry Immunohistology revealed the expression of CD4, CCR5, and ITGAM in IBM and normal muscle tissues.

Immune Infiltration Analysis

The CIBERSORT algorithm results are shown in Figure 6A. Immune infiltration analysis results of IBM patients and healthy control subjects are shown in Figure 6B. The expression levels of naive B cells, plasma cells, CD8+ T cells, follicular helper T cells, and M2 macrophages were higher in the IBM group than in the healthy control group. Meanwhile, the degree of infiltration of memory B cells, CD4+ naive T cells, and resting natural killer (NK) cells was lower in the IBM group than in the healthy control group. The relationship between the hub genes and the level of the immune cell types was then assessed. STAT1, ITGAM, and CCR5 were positively correlated with M2 macrophages; CXCL10 and CD4 were negatively correlated with the level of resting natural killer (NK) cells (Figure 6C). These results suggested that the hub genes can regulate the immune infiltration of IBM and thus could be potential targets for IBM treatment.

Figure 6 Immune landscape of IBM (A) The relative proportion of infiltrated immune cells between IBM patients and healthy control subjects. (B) Boxplot showing the level of the 22 infiltrating immune cells in IBM and healthy control subjects. (C) Heatmap showing the relationship between the 10 hub genes and the 22 infiltrating immune cells. ns: no significance; *P < 0.05; **P <0.01; ***P < 0.001.

Gene Set Variation and Enrichment Analysis (GSVA and GSEA)

GSVA is used to analyze differences in pathway activities per sample. Herein, GSVA was used to calculate the enrichment scores of Hallmark gene sets in GSE3112 and GSE128470. The enrichment scores of 50 Hallmark gene sets between IBM and normal samples are shown in Figure 7A and B. Multiple pathways and biological processes, including apoptosis, inflammatory response, interferon-alpha response, IL2-STAT5 signaling, and P53 pathways, were activated in IBM.

Figure 7 GSVA and GSEA analysis (A and B) Heatmap showing the GSVA scores of potential pathways and biological processes in GSE128470 and GSE3112. (C-F) GSEA showing the different pathways and biological processes enriched in IBM.

GSEA analysis was used to explore the potential pathways and biological processes involved in IBM. The GSEA results of GSE3112 showed that genes in the IBM samples were enriched in allograft rejection, apical junction, apoptosis, complement, epithelial-mesenchymal transition, IL2-STAT5 signaling, inflammatory response, interferon-gamma response, KRAS signaling up, and P53 pathway (Figure 7C and D). The GSEA results of GSE128470 showed that genes in the IBM samples were enriched in allograft rejection, inflammatory response, KRAS signaling up, TNFα signaling via NFKβ, complement, epithelial-mesenchymal transition, P53 pathway, interferon-gamma response, apical junction, and IL2-STAT5 signaling (Figure 7E and F). These results indicate that genes in the IBM samples were all enriched in allograft rejection, apical junction, complement, epithelial-mesenchymal transition, IL2-STAT5 signaling, inflammatory response, interferon-gamma response, KRAS signaling up, and P53 pathway. Moreover, the relationship between the hub genes and potential pathways was also explored. The hub genes were positively correlated with these pathways and biological processes (Figure 8A and B).

Figure 8 Relationship between the hub genes and potential pathways, biological processes of IBM (A and B) The heatmap showing the 10 hub genes positively correlated with the potential pathways and biological processes in GSE3112 and GSE128470. Red represents a positive correlation.

The DEGs were analyzed using CMap to identify the small molecule drugs, including nadolol, nifuroxazide, 3-acetamidocoumarin, prestwick-1103, and sulfaquinoxaline, with potential therapeutic effects in IBM (Table 1).

Table 1 The Connectivity Map (CMap) Analysis Results

Discussion

Inclusion body myositis (IBM) has unique characteristics among the idiopathic inflammatory myopathies. However, there are no effective treatment interventions for IBM. It is widely acknowledged that IBM has a complex pathogenesis and has not been fully clarified. It is unclear whether immune-mediated and degenerative processes are causally linked in IBM. Previously studies have reported that IBM pathogenesis involves clonally-restricted and antigen-driven infiltrating CD8-positive T cells. IBM pathogenesis is also associated with abnormal protein deposition and mitochondrial dysfunction. Herein, bioinformatics and data mining analyses were performed to identify the DEGs from muscle biopsy samples from IBM patients and healthy groups. Enrichment analysis was used to screen biochemical processes and pathways. GSVA and GSEA analysis were conducted to confirm the potential pathways of IBM from enrichment analysis. The relative levels of immune cells in IBM and normal muscle groups were also evaluated. In summary, the molecular and immunological characteristics of IBM were investigated, providing possible molecules and signaling pathways for the understanding and treatment of IBM.

Ten hub genes (CD4, PTPRC, ITGAM, TYROBP, CCR5, IRF8, STAT1, ITGB2, HLA-DRA, and CXCL10), closely related to IBM were identified based on PPI network. CD4 encodes a glycoprotein expressed on the surface of T-lymphocyte populations of different functions.¹⁵ CD4 expression is correlated with the specificity of T cells to class II or I major histocompatibility complex (MHC) molecules on target cells. Therefore, MHC molecules can modulate immune cells by targeting CD4. Previous research showed that Cytotoxic CD8+T cells, some macrophages, and CD4+T cells could invade the muscle fibers of IBM. PTPRC (CD45) encodes the protein tyrosine phosphatase (PTP) family, regulating T cell and B cell antigen receptor signaling. Moreover, human PTPRC is associated with immunodeficiency disease and viral susceptibility.¹⁶ ITGAM encodes integrin αM chains, which play a key role in various inflammatory responses (in monocytes and neutrophils). Some studies have also demonstrated that ITGAM can act as a complement 3 receptor, involved in inflammatory response.¹⁴ TYROBP is a transmembrane receptor widely found in macrophages/monocytes, natural killer cells (NK), and neutrophils. Macrophages have played a key role in IBM pathogenesis in recent years. Previous studies have shown that TYROBP is one of the key drivers of multiple inflammatory pathways.¹⁷ Meanwhile, IRF8 stabilizes the transcription mechanism of type I IFN promoters, thus enhancing IFN expression in dendritic cells and monocytes.¹⁸ These results suggest that immune response is involved in IBM pathogenesis.

CCR5 modulates T-cell migration, and it is associated with some autoimmune disorders. Studies have shown that CCR5 binds up-regulate pro-inflammatory chemokines in IBM,¹³ consistent with this study. HLA-DRA encodes the subunit HLA-DR of α MHC class II molecule. It is expressed in activated CD4+ T cells, B cells, macrophages, and dendritic cells. Previous studies have shown that MHC class II positivity has high specificity in inflammatory myopathy. Herein, HLA-DRA expression was significantly increased in IBM patients.¹⁹ Therefore, it is necessary to analyze the expression pattern and function of HLA-DRA. CXCL10 belongs to CXC chemokine family. CXCL10 induces chemotaxis, apoptosis, and cell growth by binding to CXCR3 receptors. CXCL10 is widely expressed on the surface of IBM myofibers.²⁰ The inflammatory cells invading and surrounding non-necrotic muscle fibers are strongly associated with CCXL10 in the cytotoxic immune response of IBM. Herein, several hub genes (CCR5, HLA-DRA, and PTPRC) could distinguish IBM from other idiopathic inflammatory myopathy, including DM, NM, and NS. In contrast, all the 10 hub genes could not distinguish IBM and PM. However, these results should be verified using large cohorts containing enough samples.

GO and KEGG analyses showed that the DEGs were mainly enriched in Th1 and Th2 cell differentiation, chemokine activity, chemokine signaling pathway, MHC protein complex, antigen processing and presentation, and cellular response to interferon-gamma. These results indicate that immune infiltration and immune response in the local microenvironment play a key role in IBM. Th1-cells can activate the cytotoxicity of macrophages and CD8 T cells, thus inducing cellular immunity. Th2 cells induce the differentiation of B cells into plasma cells or memory B cells, thus producing antibodies. These findings suggest that Th1-mediated immune response is crucial in idiopathic inflammatory myopathy. Moreover, most T cells in muscle tissues are chemokine receptors associated with Th1 lineage. Studies have reported that IFN-γ, a cellular cytokine strongly associated with Th1, is up-regulated in IBM, leading to significant upregulation of autophagy activity in skeletal muscle and amyloid deposition in muscle fibers.²⁰ We found similar results in our research.

Immune cell infiltrations also play a crucial part in IBM pathogenesis. Greenberg reported a group of KLRG1⁺ T cells that could aggressively intrude IBM myofiber. Therefore, targeting this cell type is a potential treatment for IBM patients.²¹ Herein, CD8 T cells, Tregs, and macrophages were crucial in IBM. T cells and macrophages are the primary infiltrated immune cells in IBM, similar to this study. Tregs can inhibit effector response mediated by CD8+ T cells, thus protecting muscle cells.²² Herein, the overall degree of Tregs was significantly down-regulated in the IBM group than in the healthy controls. Moreover, the results showed that the hub genes may be involved in the regulation of multiple immune cells in IBM. These findings suggest that targeting hub genes and immune cells might be an effective therapeutic strategy for IBM patients.

The results of GSEA were similar to the DEG functional analysis results. Besides, GSEA revealed various enrichment functions, including IL2-STAT5 signaling, KRAS signaling, and P53 pathway. These signaling pathways are typical mutated pathways in tumorigenesis. IBM has been considered a separate inflammatory myopathy, independent of the tumor, in the past few decades. Although no study has assessed the relationship between IBM and malignancies, some studies have shown the relationship between dermatomyositis and malignant tumor. A retrospective study exhibited a significantly increased incidence of cancer in patients with dermatomyositis.²³ In some reports, patients have been diagnosed with uterine cancer, breast cancer, papillary thyroid carcinoma,²⁴ and renal cell carcinoma²⁵ within five years before or after IBM diagnosis. For instance, Michael et al assessed a 64-year-old man with IBM and concurrent transitional cell carcinoma of the bladder who recovered strength and self-care ability.²⁶ However, some rare sporadic IBM may be due to a paraneoplastic presentation. Herein, interferon-gamma response was also identified. Previous studies have reported that the induction of IFN-gamma significantly influences the proliferative response of HIV-infected peripheral blood mononuclear cells in vitro.²⁷ Moreover, the imbalance between the pro-inflammatory and anti-inflammatory responses mediated by polymorphisms in the IFN-gamma may influence the outcome of chronic HCV infection.²⁸ Some researchers have found that IBM may be related to some virus infections, such as hepatitis C virus and HIV virus. A retrospective study investigated the prevalence of HCV infection in IBM patients, and found a markedly higher infection rate of HCV in IBM than that in the age-matched PM.²⁹ It was reported that HIV virus might be a deciding factor for the progress of IBM because over 10% of infiltrating CD8+ T cells displayed specificity for human leukocyte antigen‐A* 0201‐HIV‐gag. However, these views are majorly theoretical.³⁰ Further basic and clinical researches are needed to determine the exact relationship between IBM and these disorders.

There are some limitations in this study. Our study was mainly based on bioinformatics analyses. The in vivo and in vitro researches should be validated to elaborate the molecular and immune-related mechanisms of IBM.

Conclusions

In this study, potential pathways and hub genes (CD4, PTPRC, ITGAM, TYROBP, CCR5, IRF8, STAT1, ITGB2, HLA-DRA, and CXCL10) correlated with IBM were identified through bioinformatics analysis. Further analyses demonstrated that these hub genes can regulate immune infiltrations and multiple pathways in IBM. Finally, CMap database was used to predict five candidate drugs for IBM patients. Therefore, this research can provide new insights into mechanism and potential therapeutic targets for the treatment of IBM.