Bioinformatics Analysis and Identification of Genes and Molecular Pathways Involved in Synovial Inflammation in Rheumatoid Arthritis.

BACKGROUND Rheumatoid arthritis (RA) has a high prevalence in the elderly population. The genes and pathways in the inflamed synovium in patients with RA are poorly understood. This study aimed to identify differentially expressed genes (DEGs) linked to the progression of synovial inflammation in RA using bioinformatics analysis. MATERIAL AND METHODS Gene expression profiles of datasets GSE55235 and GSE55457 were acquired from the Gene Expression Omnibus (GEO) database. DEGs were identified using Morpheus software, and co-expressed DEGs were identified with Venn diagrams. Protein-protein interaction (PPI) networks were assembled with Cytoscape software and separated into subnetworks using the Molecular Complex Detection (MCODE) algorithm. The functions of the top module were assessed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID). The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed. RESULTS DEGs that were upregulated were significantly enhanced in protein binding, the cell cytosol, organization of the extracellular matrix (ECM), regulation of RNA transcription, and cell adhesion. DEGs that were downregulated were associated with control of the immune response, B-cell and T-cell receptor signaling pathway regulation. KEGG pathway analysis showed that upregulated DEGs enhanced pathways associated with the cell adherens junction, osteoclast differentiation, and hereditary cardiomyopathies. Downregulated DEGs were enriched in primary immunodeficiency, cell adhesion molecules (CAMs), cytokine-cytokine receptor interaction, and hematopoietic cell lineages. CONCLUSIONS The findings from this bioinformatics network analysis study identified molecular mechanisms and the key hub genes that may contribute to synovial inflammation in patients with RA.

Background

Rheumatoid arthritis (RA) is an autoimmune disease that includes multiple disease subsets defined by the several inflammatory cascades that result in a final shared pathway where persistent synovial inflammation leads to damage to articular cartilage and underlying bone [1-4]. The etiology of RA continues to be investigated [5], but genetic factors are found in 50% of patients who develop RA [6], and more than 30 genetic regions have been linked to RA [7,8]. RA currently affects 0.1–0.5% of adults and is approximately three times more common in women and increases with age, with the highest incidence in women who are more than 65 years [9,10].

Studies have shown the role of synovial inflammation in the pathogenesis of RA [11,12]. The two main groups of synovial cells include fibroblast-like and macrophage-like synoviocytes. Inflammatory cytokines that are upregulated in rheumatoid synovium lead to the development of macrophage-like synoviocytes. Fibroblast-like synoviocytes may contribute to local damage to the joint in patients with RA. Pap et al. reported the findings from studies in mouse models of RA that showed that co-implantation of fibroblast-like synoviocytes and cartilage was associated with invasion of fibroblasts into cartilage, resulting in destruction of the joint [13]. The inflamed synovium in RA is associated with pain and reduced joint function [14]. However, the mechanisms of gene and protein expression in the synovium associated with the pathogenesis of RA remain poorly understood but are necessary to improve the prevention, diagnosis, and treatment of RA.

Microarray techniques have increased the capability to explore the pathogenic processes of several diseases and represents an important technology for functional genomic studies [15,16]. Microarrays have been used detect certain RA-specific proteins, including anacardic acid, histone acetyltransferase, nuclear factor-kappa B, and prostaglandin D2 synthase [17,18]. The data show that proteins produced by differentially expressed genes (DEGs) found within various cell structures and with a range of molecular functions (MFs) were linked to altered biological processes (BPs) in the progression of RA.

Therefore, this study aimed to identify DEGs linked to the progression of synovial inflammation in RA using bioinformatics analysis using Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis.

Material and Methods

Microarray data

The Gene Expression Ominibus (GEO) database () was used, which stored original submitted archives in addition to curated datasets. The gene expression profiles, GSE55235 and GSE55457, were acquired from the GEO database. Data from ten samples of synovium from normal individuals (GSM1332201, GSM1332202, GSM1332203, GSM1332204, GSM1332205, GSM1332206, GSM1332207, GSM1332208, GSM1332209, and GSM1332210), and synovium from ten patients with RA (GSM1332221, GSM1332222, GSM1332223, GSM1332224, GSM1332225, GSM1332226, GSM1332227, GSM1332228, GSM1332229, and GSM1332230) were obtained from the GSE55235 dataset. The data from synovium from ten patients with RA (GSM1337304, GSM1337305, GSM1337306, GSM1337307, GSM1337308, GSM1337309, GSM1337310, GSM1337311, GSM1337312, and GSM1337313) and 13 samples of synovium from patients with RA (GSM1337314, GSM1337315, GSM1337316, GSM1337317, GSM1337318, GSM1337319, GSM1337320, GSM1337321, GSM1337322, GSM1337323, GSM1337324, GSM1337325, and GSM1337326) were obtained from the GSE55457 dataset.

Identifying differentially expressed genes (DEGs)

The initial microarray data was assessed from the web-based tool, Morpheus () using heatmaps to view changes in gene expression. The data were then separated into normal and RA synovium groups for analysis. To determine the DEGs in RA synovium GEO2R () and GEO2R () were used. A DEG was defined by P<0.05 and a fold-change of log2 ≥1.

Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis

Annotation of cellular components (CCs), biological processes (BPs) and molecular functions (MFs) of DEGs were determined using Gene Ontology (GO) enrichment analysis. The route of gene cluster and related functions was determined using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (). This study used the Database for Annotation, Visualization and Integrated Discovery (DAVID) version 6.8 () to perform the KEGG pathway and GO enrichment analysis. A P-value <0.05 was considered to be statistically significant.

Protein-protein interaction (PPI) network analysis

To determine functional interactions between DEGs, protein-protein interaction (PPI) networks were created using the online Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) tool (). PPI networks were created with Cytoscape software (). The biological network, or collection of nodes that represent biological processes, had edges or vertices that represented PPIs. The ten top nodes were identified according to the level of interaction in the PPI network. Then, Molecular Complex Detection (MCODE) was performed to monitor PPI network modules with Cytoscape. The conditions for selection were: node score cutoff=0.2, degree cutoff=2, maximum depth=100, and k-core=2. The functional enhancement analysis of DEGs in the top module were performed using DAVID.

Results

Identification of differentially expressed genes (DEGs)

In total, 23 rheumatoid arthritis (RA) samples of synovium and 20 normal samples of synovium were analyzed. Gene expression profiles were assessed independently using Morpheus software. DEG expression heatmaps of the 30 most upregulated and down-regulated genes of the two expression groups are shown in Figure 1. There were 43 upregulated genes and 210 downregulated genes that were co-expressed in the two gene expression groups (Figure 2).

Figure 1

Heatmaps of the differentially expressed genes (DEGs) including the 30 most upregulated and 30 most down-regulated genes of the GSE55235 and GSE55457 expression groups. (A) Heatmap of the 60 most differentially expressed genes of GSE55235, including 30 upregulated (red) and 30 downregulated (blue) gene hubs. (B) Heatmap of the top 60 differentially expressed genes of GSE55457, including 30 upregulated (red) and 30 downregulated (blue) gene hubs.

Figure 2

Co-expression of the upregulated and downregulated genes of the GSE55235 and GSE55457 expression groups.

Gene Ontology (GO) enrichment

Gene Ontology (GO) showed that upregulated DEGs were significantly enhanced for cell adhesion, protein binding, cell cytosol, the organization of the extracellular matrix (ECM), and positive transcription regulation of RNA polymerase II promoter. Downregulated DEGs were enriched for the immune response, external plasma membrane, the T-cell receptor signaling pathway, regulation of immune response, and the B-cell receptor signaling pathway (Table 1).

Table 1

Gene Ontology (GO) analysis of upregulated and downregulated differentially expressed genes (DEGs) in the biological processes (BPs).

A. Upregulated
Term	Function	Count	P-value
GO: 0005515	Protein binding	427	1.56×10−10
GO: 0005829	Cytosol	186	1.57×10−9
GO: 0030198	Extracellular matrix organization	27	5.96×10−8
GO: 0007155	Cell adhesion	44	1.36×10−7
GO: 0045944	Positive regulation of transcription from RNA polymerase II promoter	72	4.78×10−7
B. Downregulated
Term	Function	Count	P-value
GO: 0006955	Immune response	36	3.68×10−23
GO: 0009897	External side of plasma membrane	21	8.26×10−15
GO: 0050853	B cell receptor signaling pathway	13	9.66×10−14
GO: 0050776	Regulation of immune response	15	1.76×10−9
GO: 0050852	T cell receptor signaling pathway	14	1.78×10−9

GO – Gene Ontology.

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis

The upregulated and down-regulated DEGs from the KEGG pathway analysis are shown in Table 2. The five most upregulated DEGs from the KEGG pathway analysis were those for the adherens junction, dilated cardiomyopathy, arrhythmogenic right ventricular cardiomyopathy (ARVC), osteoclast differentiation, and hypertrophic cardiomyopathy (HCM). The five most downregulated DEGs from the KEGG pathway analysis were for primary immunodeficiency, cell adhesion molecules (CAMs), the hematopoietic cell lineage, and cytokine-cytokine receptor interactions (Table 2).

Table 2

The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of upregulated and downregulated differentially expressed genes (DEGs), including the top five terms selected according to the P-value when more than five enriched terms were identified in each category.

A. Upregulated
Pathway ID	Name	Count	P-value	Genes
hsa04520	Adherens junction	15	7.34×10−6	EGFR, PTPRF, TGFBR2, LEF1, SMAD3, SNAI1, TCF7L2, IQGAP1, CTNNA2, IGF1R, CSNK2A1, SORBS1, AFDN, WASL, INSR
hsa05414	Dilated cardiomyopathy	15	5.39×10−5	ADCY1, TNF, ADCY2, ITGB4, CACNB2, CACNB3, ITGB3, TPM2, TPM1, TGFB2, DES, ITGA7, SGCD, PRKACA, SGCA
hsa04380	Osteoclast differentiation	19	7.00×10−5	IL1R1, TNF, FOSL2, SOCS3, CSF1, TGFBR2, NOX1, NFKB2, ITGB3, JUNB, TGFB2, TNFRSF11B, LILRA2, JUN, MAPK14, MAPK8, IKBKB, MAP2K7, AKT2
hsa05412	Arrhythmogenic right ventricular cardiomyopathy (ARVC)	12	3.89×10−4	DES, ITGA7, ITGB4, SGCD, CACNB2, LEF1, CACNB3, CDH2, ITGB3, TCF7L2, SGCA, CTNNA2
hsa05410	Hypertrophic cardiomyopathy (HCM)	13	3.98×10−4	IL6, TNF, ITGB4, CACNB2, CACNB3, ITGB3, TPM2, TPM1, TGFB2, DES, ITGA7, SGCD, SGCA
B. Downregulated
Pathway ID	Name	Count	P-value	Genes
hsa05340	Primary immunodeficiency	9	5.44×10−9	PTPRC, CD19, CD3D, CD8A, LCK, IL2RG, CD79A, IL7R, BLNK
hsa04060	Cytokine-cytokine receptor interaction	17	3.20×10−8	IL21R, CXCL9, TNFRSF17, IL15, CXCL11, CCL5, IL7R, CCL18, CXCL10, TNFSF10, CCR5, CXCL13, CCR2, CXCR6, CSF2RB, IL2RG, CD27
hsa04514	Cell adhesion molecules (CAMs)	10	5.74×10−5	PTPRC, SDC1, CD8A, CD2, HLA-DPA1, ITGA4, HLA-DMB, CD6, HLA-DOB, PDCD1LG2
hsa04062	Chemokine signaling pathway	11	8.83×10−5	ITK, CCR5, CXCL13, CCR2, CXCR6, CXCL9, CCL5, CXCL11, STAT1, CCL18, CXCL10
hsa04640	Hematopoietic cell lineage	8	9.01×10−5	CD38, CD19, CD3D, CD8A, MS4A1, CD2, ITGA4, IL7R

KEGG – Kyoto Encyclopedia of Genes and Genomes.

Protein-protein interaction (PPI) network construction and module analysis

Using data from the Cytoscape and STRING databases, the ten hub nodes with the greatest degree of network connection were determined. The top ten hub genes identified were PTPRC, CD19, LCK, CD2, CCL5, CCR5, CXCL10, STAT1, GZMB, and CXCL9 (Table 3). The gene with the greatest number of node networks (n=33) was the protein tyrosine phosphatase receptor type C (PTPRC) gene. Six PPI network modules fulfilled the degree cutoff (2), the maximum depth (100), the k-core (2), and the node score cutoff (0.2) (Table 4).

Table 3

Top ten genes with the highest degree of interaction in the protein-protein interaction (PPI) network.

Gene ID	Degree
PTPRC	33
CD19	29
LCK	27
CD2	27
CCL5	24
CCR5	23
CXCL10	22
STAT1	22
GZMB	19
CXCL9	18

All top ten genes were downregulated in the protein-protein interaction (PPI) network. PTPRC – protein tyrosine phosphatase, receptor type, C; LCK – lymphocyte-specific protein tyrosine kinase; CCL5 – chemokine (C-C motif) ligand 5; CCR5 – chemokine (C-C motif) receptor 5; CXCL10 – chemokine (C-X-C motif) ligand 10; STAT1 – signal transducer and activator of transcription 1; GZMB – granzyme B; CXCL9 – chemokine (C-X-C motif) ligand 9.

Table 4

Six modules from the protein-protein interaction (PPI) network satisfied the criteria of MCODE scores ≥4 and number of nodes >4.

Cluster	Score	Nodes	Edges	Node IDs
1	12	12	66	RAD51AP1, HMMR, RAD51, MAD2L1, DLGAP5, CDK1, BUB1, RRM2, CDC20, TYMS, KIAA0101, KIF11
2	6.889	19	62	SDC1, SLAMF1, CD79A, TLR7, TRAT1, CD19, CCR5, GZMK, CD3D, ITK, CD247, IL7R, GZMB, PTPRC, CD48, NKG7, ITGA4, CD2, GZMH
3	6.364	12	35	CXCL9, LCK, CD27, CCR2, CXCL13, CCL5, CXCR6, IL15, CXCL11, PNOC, CD38, CXCL10
4	4	5	8	TNFSF10, DDX60, RTP4, STAT1, GBP1
5	4	4	6	IGJ, TNFRSF17, POU2AF1, MZB1
6	3	3	3	ACVR1B, TGFBR2, SMAD3

Score = density × no. of nodes. RAD51AP1 – RAD51 associated protein 1; HMMR – hyaluronan-mediated motility receptor; MAD2L1 – mitotic arrest deficient-like 1; DLGAP5 – discs, large (drosophila) homolog-associated protein 5; CDK1 – cyclin-dependent kinase 1; BUB1 – budding uninhibited by benzimidazoles 1; RRM2 – ribonucleotide reductase M2; CDC20 – cell division cycle 20; TYMS – thymidylate synthetase; KIF11 – kinesin family member 11; SDC1 – syndecan 1; SLAMF1 – signaling lymphocytic activation molecule family member 1; TLR7 – toll-like receptor 7; TRAT1 – cell receptor associated transmembrane adaptor 1; GZMK – granzyme K; ITK – IL2-inducible T-cell kinase; PTPRC – protein tyrosine phosphatase, receptor type, C; LCK – lymphocyte-specific protein tyrosine kinase; CCL5 – chemokine (C-C motif) ligand 5; CCR5 – chemokine (C-C motif) receptor 5; CXCL10 – chemokine (C-X-C motif) ligand 10; NKG7 – natural killer cell group 7 sequence; ITGA4 – integrin, alpha 4; STAT1 – signal transducer and activator of transcription 1; GZMB – granzyme B; CXCL9 – chemokine (C-X-C motif) ligand 9; PNOC – prepronociceptin; TNFSF10 – tumor necrosis factor (ligand) superfamily, member 10; DDX60 – DEAD (Asp-Glu-Ala-Asp) box polypeptide 60; RTP4 – receptor (chemosensory) transporter protein 4; GBP1 – guanylate binding protein 1; IGJ – immunoglobulin J polypeptide; TNFRSF17 – tumor necrosis factor receptor superfamily, member 17; POU2AF1 – POU class 2 associating factor 1; MZB1 – marginal zone B and B1 cell-specific protein; ACVR1B – activin A receptor, type 1B.

GO and KEGG pathway gene enhancement were incorporated into the principal module (Figure 3), which showed that these genes were mainly linked to the T-cell receptor signaling pathway, protein binding, the immune response, the inflammatory response, immune cell receptor activity, transmembrane signaling receptor activity, cytokine-cytokine receptor interactions, primary immunodeficiency, and the human T-lymphotropic virus-1 (HTLV-I) infection (Table 5).

Figure 3

The top module from the protein-protein interaction (PPI) network.

Table 5

Functional and pathway enrichment analysis of the genes in the module, with the top three terms selected according to the P-value when more than three enriched terms were identified in each category.

A. Biological processes (BPs)
Term	Name	Count	P-value	Genes
GO: 0006955	Immune response	23	2.76×10−16	LY75, CD8A, GZMA, CXCL9, IL15, HLA-DMB, CXCL11, CCL5, IL7R, PDCD1LG2, CXCL10, NCR3, IGF1R, TNFSF10, FCGR2B, CCR5, CCR2, CTSC, HLA-DPA1, IL2RG, SEMA4D, HLA-DOB, CD27
GO: 0050852	T-cell receptor signaling pathway	11	5.23×10−9	PSMB10, ITK, PTPRC, CD3D, GATA3, PLCG2, CD247, LCK, HLA-DPA1, TRAT1, PSMB9
GO: 0006954	Inflammatory response	14	8.36×10−8	LY75, CXCL9, IL15, CXCL11, CCL5, TLR7, NCR3, CXCL10, SDC1, CCR5, CCR2, CXCR6, CD27, BLNK
B. Molecular functions (MFs)
Term	Name	Count	P-value	Genes
GO: 0005515	Protein binding	66	2.96×10−6	CD8A, LDLR, IL15, CXCL11, CXCL10, CD48, ACVR1B, SH2D1A, GATA3, DDX60, MS4A1, CSF2RB, LOXL2, CDK1, POU2AF1, MYO6, GZMA, LEF1, GZMB, RAD51, PRKCB, SDC1, MAD2L1, CCR5, JUN, RRM2, LCK, CTSC, SEMA4D, GAP43, GBP1, PSMB10, WNT5A, CD247, KIAA0101, CXCL9, MAP4K1, BTN2A2, CD72, CCL5, IL7R, HMMR, IGF1R, BUB1, CD2, IL2RG, CD27, BLNK, ITK, PTPRC, DLGAP5, TGFBR2, CDC20, FZD2, ACACB, ITGA4, STAT1, BIRC3, SLAMF1, PSMB9, TNFSF10, CD19, FCGR2B, PLCG2, IRF4, NKG7
GO: 0004888	Transmembrane signaling receptor activity	9	1.22×10−5	LY75, FCGR2B, CD3D, CD247, SEMA4D, CD72, TLR7, SLAMF1, CD27
GO: 0004872	Receptor activity	9	1.35×10−5	CD48, LY75, LDLR, CD2, CSF2RB, SEMA4D, TLR7, SLAMF1, NCR3
C. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways
Term	Name	Count	P-value	Genes
hsa04060	Cytokine-cytokine receptor interaction	16	7.26×10−9	TGFBR2, IL21R, CXCL9, IL15, CXCL11, CCL5, IL7R, CXCL10, ACVR1B, TNFSF10, CCR5, CCR2, CXCR6, CSF2RB, IL2RG, CD27
hsa05340	Primary immunodeficiency	8	2.64×10−8	PTPRC, CD19, CD3D, CD8A, LCK, IL2RG, IL7R, BLNK
hsa05166	HTLV-I infection	13	4.86×10−6	WNT5A, MAD2L1, CD3D, JUN, LCK, TGFBR2, IL2RG, CDC20, HLA-DPA1, IL15, FZD2, HLA-DMB, HLA-DOB

KEGG – Kyoto Encyclopedia of Genes and Genome.

Discussion

Currently, many factors are believed to contribute to the pathogenesis of rheumatoid arthritis (RA), including genetics, age, gender, and autoimmunity [19]. The management of RA is mainly directed to the treatment of symptoms associated with joint disease. The aims of this study were to identify differentially expressed genes (DEGs) and protein-protein interaction (PPI) networks linked to the progression of synovial inflammation in RA using bioinformatics analysis. The data from this study may provide greater insight into the etiology of RA that involves the synovium. Also, the hub genes defined in this study may lead to the development of potential diagnostic, prognostic, or therapeutic biomarkers for RA. Table 3 summarizes the top ten PPI network hub genes that were identified in the synovium of patients with RA, which included PTPRC, CD19, LCK, CD2, CCL5, CCR5, CXCL10, STAT1, GZMB, and CXCL9.

The PTPRC gene encodes protein tyrosine phosphatase receptor type C and is a gene that acts as a T-cell coactivation regulator that binds to dipeptidyl peptidase-4 (DPP4), also known as CD26. Protein tyrosine-protein phosphatase is needed for activation of T-cells through the T-cell antigen receptor. The early PTPase domain possesses enzymatic activity and T-cell activation results in the recruitment and dephosphorylation of Src kinase-associated phosphoprotein 1 (SKAP1) and the non-receptor tyrosine-protein kinase, FYN [20,21]. Previous studies have reported that PTPRC is involved in the progression of glioma and ovarian cancer cells [20,21]. It may be hypothesized that the upregulation of the PTPRC gene may preserve the normal function of the synovium.

The CD19 gene encodes for CD19, a molecule that is involved in binding with the antigen receptor of B-lymphocytes and is necessary for the development of a cellular immune response. CD19 binds to CD2 expressed by T-cells and interacts with lymphocyte function-associated antigen-3 (LFA-3) and CD48/BCM1 to facilitate T-cell adhesion to different cell types. CD19 is also involved in the activation and signaling of B-cells. In studies of fibroblast-like synoviocytes (FLS) from patients with RA grown in vitro, inhibition of CD2 and CD19 has been shown to reduce the recruitment of immune cells [22,23]. This findings from the present study showed both CD19 and CD2 were in the top ten PPI network hub genes identified in the synovium of patients with RA.

The CCL5 gene encodes for the CCL5 chemoattractant for memory T-helper cells, blood monocytes, and eosinophils, which causes histamine release from basophils, activates eosinophils, and several chemokine receptors, including CCR1, CCR3, CCR4, and CCR5 [24]. Downregulation of CCL5 can induce reduced inflammation of the synovium in patients with RA [24]. Agere et al. [25] reported that CCL5 induced degradation of collagen by activation of matrix metalloproteinase-1 (MMP-1) and MMP-13 in human synovial fibroblasts from patients with RA. Currently, the biological roles of CCL5 and CCR5 in RA remain unknown. Rump et al. [26] reported that CCL5 and CCR5 enhanced homing and differentiation of circulating osteoclast progenitor cells in rheumatoid synovium. Other studies have shown that CCL5 is secreted by circulating T-cells and contributes to the chemotactic activity of T-cells in RA [27-29]. Given these findings, it is possible that reduced levels of CCL5 and CCR5 in the synovium could alter inflammation in the synovium and the progression of RA.

The LCK gene encodes for a non-receptor tyrosine-protein kinase that has a vital role in T-cell maturation and selection in the thymus and in T-cell antigen receptor (TCR)-linked signal transduction pathways. Kumar et al. [30] showed that LCK is composed of three domains, SH3 linked to the SH2 domain at the amino terminal, and a kinase domain at the carboxy terminal. Physiologically, LCK assists with the development, differentiation, and function of T-cells controls neurite outgrowth and preserves long-term neuronal synaptic plasticity. Modification in the activity and expression of LCK is linked to several diseases, including asthma, cancer, diabetes, psoriasis, atherosclerosis, inflammatory bowel disease, and RA [30]. Farag et al. [31] showed that potent and selective LCK inhibitors could be used in inflammatory diseases, including RA. These findings suggest that the LCK downregulation might reduce inflammation in rheumatoid synovium.

The CXCL10 and CXCL9 genes encode for chemotactic agents for T-lymphocytes and monocytes and can bind to CXCR3, which may facilitate the proliferation, angiogenic activity, and survival of human mesangial cells within a heterotrimeric G-protein signaling pathway. Gao et al. showed that inhibition of CXCR3 (a CXCL10 receptor) significantly reduced chronic pulmonary inflammation by reducing inflammatory cell recruitment [32]. Studies have also reported that CXCL10 is necessary for cancer cell growth [33,34]. The findings of the present study have shown that CXCL10 and CXCL9 gene upregulation could be a therapeutic target for the maintenance of the normal functioning of cells in the synovium.

The STAT1 gene encodes for a transcription activator and signal transducer that facilitates cellular responses to cytokines (including KITLG/SCF), interferons (IFNs) and growth factors. Following the attachment of type I IFNs to cell surface receptors, signaling via protein kinases results in the stimulation of JAK kinases (JAK1 and TYK2) and to tyrosine phosphorylation of STAT 1 and 2. Phosphorylated STAT dimerizes and binds to ISGF3G/IRF-9 to create a complex called ISGF3 transcription factor that is transported to the cell nucleus.

Expression of the GZMB gene is required for target cell lysis in cell-mediated immune responses and it is associated with an activation cascade of caspases (aspartate-specific cysteine proteases) required for cell apoptosis. Bao et al. suggested that GZMB gene silencing can inhibit MAPK signaling pathway through regulation of the expression of inflammatory factors, apoptosis factors (caspase and Bcl-2), in addition to factors associated with angiogenesis (bFGF and VEGF), which can reduce hyperplasia of synovial tissue and injury to articular cartilage tissue in RA [35]. Therefore, the GZMB gene could be an effective therapeutic target for the management of RA.

In the present study, Gene Ontology (GO) terms of the biological processes (BPs) and molecular functions (MFs), in addition to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway terms that were linked to the DEGs identified in the top module were determined. Regarding the BPs, the DEGs in the top module of the PPI networks were mainly enhanced for T-cell receptor signaling pathways, the immune response, and the inflammatory response, which could lead to the development of autoantibodies. In the MF group, DEGs were enhanced for protein binding, receptor activity, and transmembrane signaling receptor activity. The KEGG pathway enrichment analysis identified cytokine-cytokine receptor interaction, primary immunodeficiency, and the human T-lymphotropic virus-1 (HTLV-I) infection, which suggested an important role for the breakdown in immune pathways in the rheumatoid synovium.

Because many of the genes identified in this study had not been previously associated with RA, further studies are needed to validate the expression and levels expression of these genes in rheumatoid synovium and healthy control synovial tissue samples. Also, cells isolated from synovial tissue samples should be cultured in vitro to determine the molecular mechanisms associated with the expression of these genes. Preclinical animal models using gene knockout may also identify the functions of the genes identified and assess their role in the progression of damage to the synovium.

As with most bioinformatic network analysis studies of human disease, this study had several limitations. The study sample size was relatively small and the sampling method did not eliminate the effect of gender, co-morbidity, and the use of certain medicines that could alter gene expression in rheumatoid synovium, including methotrexate or non-steroidal anti-inflammatory drugs. The genes investigated and their pathways were not confirmed through in vitro studies or other functional studies, but this would be an area for further investigation. Also, the findings from this bioinformatics network analysis would benefit from future validation by Western blot and polymerase chain reaction (PCR) assays using samples of rheumatoid synovial tissue.

Conclusions

The findings from this bioinformatics network analysis study identified molecular mechanisms and the key hub genes that may contribute to synovial inflammation in patients with rheumatoid arthritis (RA), which require further investigation for their potential role as future diagnostic, prognostic, or therapeutic biomarkers.