Identification of Candidate Genes Related to Synovial Macrophages in Rheumatoid Arthritis by Bioinformatics Analysis

Abstract

Objective

Rheumatoid arthritis (RA) is one of the most prevalent inflammatory arthritis worldwide. However, the genes and pathways associated with macrophages from synovial fluids in RA patients still remain unclear. This study aims to screen and verify differentially expressed genes (DEGs) related to identifying candidate genes related to synovial macrophages in rheumatoid arthritis by bioinformatics analysis.

Methods

We searched the Gene Expression Omnibus (GEO) database, and GSE97779 and GSE10500 with synovial macrophages expression profiling from multiple RA microarray dataset were selected to conduct a systematic analysis. GSE97779 included nine macrophage samples from synovial fluids of RA patients and five macrophage samples from primary human blood of HC. GSE10500 included five macrophage samples from synovial fluids of RA patients and three macrophage samples from primary human blood of HC. Functional annotation of DEGs was performed, including Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Protein–protein interaction (PPI) network of DEGs was established using the STRING database. CytoHubba was used to identify hub genes. MCODE was used to determine gene clusters in the interactive network.

Results

There were 2638 DEGs (1425 upregulated genes and 1213 downregulated ones) and 889 DEGs (438 upregulated genes and 451 downregulated ones) selected from GSE97779 and GSE10500, respectively. Venn diagrams showed that 173 genes were upregulated and 106 downregulated in both two datasets. The top 10 hub genes, including FN1, VEGFA, HGF, SERPINA1, MMP9, PPBP, CD44, FPR2, IGF1, and ITGAM, were identified using the PPI network.

Conclusion

This study provides new insights for the potential biomarkers and the relevant molecular mechanisms in RA patients. Our findings suggest that the 10 candidate genes might be used in diagnosis, prognosis, and therapy of RA in the future. However, further studies are required to confirm the expression of these genes in synovial macrophages in RA and control specimen.

Keywords:

• rheumatoid arthritis
• synovitis macrophages
• key genes
• pathways
• bioinformatics analysis

Introduction

Rheumatoid arthritis (RA) is a chronic and aggressive autoimmune disorder mainly characterized by destructive synovitis, cartilage destruction, and systemic organ involvement, with a prevalence of about 0.5–1% in the global population.^1,2 Chronic pain, tenderness, swelling, deformation, and stiffness of the joints are the clinical manifestations of RA, which ultimately lead to a high incidence of disability and a heavy socio-economic burden.^3,4

In RA patients, multiple chronic inflammatory cells infiltrate the synovium, usually bringing about pain and reducing joint functions. As we all know, fibroblast-like synoviocytes, macrophage-like synoviocytes, and macrophages are the main types of cells in the synovium, and synovitis can stimulate macrophages in synovial tissues to produce angiogenic factors, resulting in abnormal vascular formation in the synovium.⁵ As we all know, macrophages are central to the pathophysiology of RA, and involved in the initiation and maintenance of inflammation, the adhesion and migration of leukocytes, the degradation of matrix and angiogenesis.^6,7 Macrophages express adhesion molecules, chemokine receptors (CR), and other surface antigens, and also secrete chemokines, cytokines, growth factors, proteases, and other mediators, and these inflammatory mediators are important pathological mechanisms of inflammation and angiogenesis in arthritis.^6,7 For example, in the early stage of RA, macrophages in the synovium mainly secrete chemokines CXCL4 and CXCL7 to recruit neutrophils and monocytes, and chemoattractants (CCL2 and IL-8) and alarmins S100A8/9 as well as tissue remodeling enzymes (MMP-3 and MMP-12) were produced throughout the active disease;⁸ in the active stage of RA, macrophage population is characterized by M1-like inflammatory cells, which produce significant increase levels of proinflammatory cytokines TNF-α, IL-6, and IL-1β, while reduced expression of the M2 macrophage associated marker CD209.⁹ In view of the fact that macrophages and their expression and secretion products play a key role in the pathogenesis of RA, they have become an ideal target for the treatment of RA.

A large number of studies have shown that genetic changes and environmental influences might promote the development of RA.¹⁰ However, it is still unclear how the expression of gene and protein in the macrophage-like synoviocytes is involved in the development of RA.¹¹ Therefore, it is of great significance to reveal the underlying regulatory mechanisms of macrophage infiltration and activation, which might provide novel insights for more effective therapeutic and preventive strategies of RA.

Currently, the identification of differentially expressed genes (DEGs) between abnormal and normal samples, as a high-throughput technique, can be used to determine various disease-related target genes, potential molecular mechanisms, diagnosis, and prognosis.¹² Although there have been studies using Gene Expression Omnibus (GEO) to explore the mRNA expression profiles of pancreatic islets in RA, there is still a lack of research exploring pathogenic genes related to RA macrophages based on bioinformatics analysis.^13,14 We used DAVID database to perform Gene Ontology (GO) functional annotation analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DEGs, then constructed protein–protein interaction (PPI) network, and finally used cytoHubba and Cytotype MCODE to identify significant modules and key genes.

Therefore, our study aimed to identify hub genes related to synovial macrophages to further explore the pathogenesis of RA by Gene Ontology (GO) functional annotation analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis, the construction of protein–protein interaction (PPI) network, and the analysis of cytoHubba and Cytotype MCODE, thereby providing key synovial macrophages-related biomarkers of RA, and new insights for the pathogenesis of RA, which might improve the treatment strategies, and diagnostic approaches in the future.

Materials and Methods

Microarray Expression Data

We downloaded the microarray expression data of GSE97779 and GSE10500 from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). GSE97779 included nine macrophage samples from synovial fluids of RA patients and five macrophage samples from primary human blood of HC. This GSE97779_series_matrix was established using the platform GPL570 (HG-U133_Plus_2) Affymetrix Human Genome U133 Plus 2.0 Array. GSE10500 included five macrophage samples from synovial fluids of RA patients and three macrophage samples from primary human blood of HC. This GSE10500_series_matrix was established using platform GPL8300 (HG-U95Av2) Affymetrix Human Genome U95 Version 2 Array.

Identification of DEGs

We used LIMMA package (http://www.bioconductor.org/packages/release/bioc/html/limma.html) of R language (https://www.r-project.org/; version 3.4.0) to screen out DEGs related to synovial macrophages between RA samples and control samples. Genes with more than 1 probe set were averaged, and genes with probe sets lacking gene symbols were eliminated. Adjusted P-value <0.05 and Log [Fold Change] (Log [FC])>1, as the screening criteria, were considered statistically significant. Venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/) was drawn to show intersections of the two data sets.

KEGG Pathway Enrichment Analysis and GO Functional Enrichment Analysis

We used DAVID (Database for Annotation, Visualization and Integrated Discovery, https://david.ncifcrf.gov/) to perform KEGG pathway enrichment analysis and GO functional enrichment analysis to explore the underlying functions of the DEGs and the main pathways involved. We analyzed the DEGs in the GSE97779 and GSE10500 datasets by DAVID online program and used the “ggplot” package of R language to illustrate the results of GO and KEGG analysis. P value <0.05 and counts ≥2 were considered statistically significant.

Construction of the PPI Network and Identification of Hub Genes

We used the STRING website (https://string-db.org) to construct the PPI network. The cutoff value was set to a confidence score ≥0.7. We analyzed the interaction relationship of DEGs using Cytoscape software (version 3.4.0; http://www.Cytoscape.org). The plug-in Molecular Complex Detection (MCODE) was used to obtain significant gene clusters and the scores of the clusters (selection criteria: degree cutoff=2; max depth=100; k-core=2; node score cutoff=0.2). The nodes were the genes with high degree which was classified by using CytoHubba. We used the algorithm Degree to select the top 10 hub genes that might have important roles in the pathogenesis of RA. P < 0.05 was considered to be statistically significant.

Results

Identification of DEGs

There were 2638 DEGs (1425 upregulated genes and 1213 downregulated genes) and 889 DEGs (438 upregulated genes and 451 downregulated genes) were selected from GSE97779 and GSE10500, respectively (Figure 1A and B). Venn diagrams showed that 173 upregulated genes and 106 down-regulated genes existed in the both data sets. (Figure 1C and D).

Figure 1 Differentially expressed genes (DEGs) expression Volcano plots and Venn diagram between rheumatoid patients and healthy controls. (A and B) Volcano plots of GSE97779 and GSE10500 Data points in red represent upregulated and in green represent downregulated genes. (C and D) Venn diagram show the Co-expression of the upregulated and downregulated genes of the GSE97779 and GSE10500 databases.

GO and KEGG Enrichment Analyses of the 173 Upregulated Genes and 106 Downregulated Genes

The results of GO functional analysis showed that DEGs were mainly enriched in inflammatory response, immune response, cell adhesion, positive regulation of gene expression, and extracellular matrix disassembly of biological process (BP), plasma membrane, extracellular exosome, extracellular space, extracellular region, and integral component of plasma membrane of CC (cellular component), and protein binding, protein homodimerization activity, identical protein binding, GTPase activator activity, and catalytic activity of MF (Molecular function) (Figure 2).

Figure 2 Gene Ontology (GO) analysis of DEGs in the biological process (BP), molecular function (MF) and cellular component (CC). Select important GO entries with P < 0.05.

The results of KEGG pathway analysis suggested that DEGs were mainly enriched in PI3K-AKT signaling pathway, focal adhesion, lysosomes, transcriptional misregulation in cancer, hematopoietic cell lineage, acute myeloid leukemia, apoptosis, bladder cancer, prostate cancer, and HIF-1 signaling pathway (Figure 3A). The important signaling pathways—PI3K/AKT pathway is presented in Figure 3B.

Figure 3 (A) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of DEGs. Select significant enrichment pathways with P < 0.05. (B) The details of the important signaling pathways—PI3K/AKT pathway was presented.

Hub Genes and PPI Network Analysis

We used the STRING website to establish a PPI network, and Cytoscape to visualize the results. Figure 4A shows that the PPI network had a total of 168 nodes and 391 edges. Figure 4B and C present the top two modules of the PPI network. Figure 5 illustrates the interactions with each other among the top 10 hub genes. The top 10 hub genes included FN1 (Fibronectin 1, degree: 24), VEGFA (Vascular Endothelial Growth Factor A, degree: 24), HGF (Hepatocyte Growth Factor, degree: 17), MMP9 (Matrix Metallopeptidase 9, degree: 17), PPBP (Pro-Platelet Basic Protein, degree: 16), CD44 (CD44 Molecule, Indian Blood Group, degree: 16), ITGAM (Integrin Subunit Alpha M, degree: 16), IGF1 (Insulin Like Growth Factor 1, degree: 15), SERPINA1 (Serpin Family A Member 1, degree: 15) and FPR2 (Formyl Peptide Receptor 2, degree: 15). Detailed information about the genes is listed in Table 1.

Table 1 Function and Connectivity of the 10 Hub Genes

Gene Symbol	Full Name	Function	Degree
FN1	Fibronectin 1	Fibronectins bind cell surfaces and various compounds including collagen, fibrin, heparin, DNA, and actin. Gene Ontology (GO) annotations related to this gene include heparin binding and protease binding.	24
VEGFA	Vascular Endothelial Growth Factor A	Growth factor active in angiogenesis, vasculogenesis and endothelial cell growth. Gene Ontology (GO) annotations related to this gene include protein homodimerization activity and protein heterodimerization activity.	24
HGF	Hepatocyte Growth Factor	Potent mitogen for mature parenchymal hepatocyte cells, seems to be a hepatotrophic factor, and acts as a growth factor for a broad spectrum of tissues and cell types. Gene Ontology (GO) annotations related to this gene include identical protein binding and serine-type endopeptidase activity.	17
MMP9	Matrix Metallopeptidase 9	May play an essential role in local proteolysis of the extracellular matrix and in leukocyte migration. Could play a role in bone osteoclastic resorption. Gene Ontology (GO) annotations related to this gene include identical protein binding and metalloendo peptidase activity.	17
PPBP	Pro-Platelet Basic Protein	LA-PF4 stimulates DNA synthesis, mitosis, glycolysis, intracellular cAMP accumulation, prostaglandin E2 secretion, and synthesis of hyaluronic acid and sulfated glycosaminoglycan. Gene Ontology (GO) annotations related to this gene include growth factor activity and glucose transmembrane transporter activity.	16
CD44	CD44 Molecule (Indian Blood Group)	Cell-surface receptor that plays a role in cell-cell interactions, cell adhesion and migration, helping them to sense and respond to changes in the tissue microenvironment. Gene Ontology (GO) annotations related to this gene include transmembrane signaling receptor activity and cytokine receptor activity.	16
ITGAM	Integrin Subunit Alpha M	Integrin ITGAM/ITGB2 is implicated in various adhesive interactions of monocytes, macrophages and granulocytes as well as in mediating the uptake of complement-coated particles and pathogens. Gene Ontology (GO) annotations related to this gene include protein heterodimerization activity.	16
IGF1	Insulin Like Growth Factor 1	The insulin-like growth factors, isolated from plasma, are structurally and functionally related to insulin but have a much higher growth-promoting activity. Gene Ontology (GO) annotations related to this gene include growth factor activity and integrin binding.	15
SERPINA1	Serpin Family A Member 1	Diseases associated with SERPINA1 include Alpha-1-Antitrypsin Deficiency and Hemorrhagic Disease Due To Alpha-1-Antitrypsin Pittsburgh Mutation. Among its related pathways are Lung fibrosis and Response to elevated platelet cytosolic Ca2+. Gene Ontology (GO) annotations related to this gene include identical protein binding and protease binding.	15
FPR2	Formyl Peptide Receptor 2	Low affinity receptor for N-formyl-methionyl peptides, which are powerful neutrophil chemotactic factors. Gene Ontology (GO) annotations related to this gene include G protein-coupled receptor activity and N-formyl peptide receptor activity.	15

Figure 4 PPI network of DEGs and four cluster modules extracted by MCODE. (A) The interaction network between proteins coded by DEGs was comprised of 188 nodes and 391 edges. Red nodes represent upregulated genes, blue nodes represent downregulated. Two significant modules identified from the PPI network using the MCODE with a score>5.0. Cluster 1 (B) had the highest cluster score (score = 8.125, 17 nodes and 65 edges), followed by cluster 2 (C) (score = 5.778, 19 nodes and 52 edges). The nodes represented as ellipse (green) and edges as lines (gray).

Figure 5 Top ten hub genes with higher degree of connectivity from PPI network.

Discussion

In terms of pathogenesis, rheumatoid arthritis is related to chronic inflammation, and its severity is correlated with the levels of the pro-inflammatory cytokines derived from the macrophages in the synovium.⁹ Increasing evidence shows that the infiltration of macrophages in the synovium is related to the degree of joint erosion, and macrophages can mediate arthritis together with other cells.^15,16

The present study might provide a new insight into the phenotype of synovial macrophages and molecular mechanisms involved in RA. There were 279 DEGs in these two datasets, of which 173 were upregulated and 106 were downregulated. In order to further understand the interactions among these DEGs, we performed GO functional enrichment analysis and KEGG pathway enrichment analysis. GO analysis showed that DEGs were mostly enriched in inflammatory response, immune response, cell adhesion, positive regulation of gene expression, and extracellular matrix disassembly of biological process (BP). KEGG analysis suggested that DEGs participated in pathways related to immune regulation, such as PI3K-Akt pathway, focal adhesion, and lysosomes. These results were in accordance with previous findings in other research that immune abnormalities and inflammatory responses could contribute to autoimmune diseases like RA. One of the significantly enriched pathways, PI3K-Akt signaling pathway, was proved to be strongly related to the breakdown of synovium macrophages in RA.

In the pathogenesis of RA, the activation of PI3K can induce the formation of a cancer-like environment, which promotes the survival of multiple cells in the immune system, including macrophages/monocytes, FLS, and dendritic cells in the joint space.¹⁷ The Akt pathway regulates cell survival, and the activation of Akt pathway in the macrophages of RA patients can inhibit apoptosis via upregulating MCL1 (myeloid cell leukemia sequence 1).¹⁸ In addition, PI3K/Akt pathway can regulate the proliferation, survival, and migration of macrophages, and participate in the response to various inflammatory and metabolic signals in macrophages.^19,20 Therefore, activating the PI3K-Akt signaling pathway can further promote cell proliferation, inhibit cell apoptosis, and elevate the expression levels of cytokines, leading to synovial hyperplasia in RA. Nonetheless, in-vitro and in-vivo research proved that different isoforms of PI3K/Akt had different functions in macrophage activation.²¹ In addition, it was reported that inhibiting PI3K-Akt pathway not only played a role in ameliorating the injuries to the joints, but also led to apoptosis and death of synovial macrophages in RA patients, suggesting that PI3K/Akt pathway served as a potential treatment target for RA,²² and the PI3K inhibitor ZSTK474 could inhibit articular inflammation in vivo.²³

Furthermore, in our study, we performed bioinformatics analysis to identify the top 10 genes highly associated with RA, including FN1, VEGFA, HGF, MMP9, PPBP, CD44, ITGAM, IGF1, SERPINA1 and FPR2, and the top 10 hub genes that might play important roles in the diagnosis, treatment, and prognosis of RA in the future.

MMP9 can degrade type IV collagen in the extracellular matrix and basement membrane. MMP9 can stimulate fibroblast-mediated inflammation, pannus formation, and cartilage degradation. In addition, it can aggravate the erosion of RA synovium and participate in the destruction of joints in RA.²⁴ In MMP-9 knockout mice, arthritis induced by antibody was ameliorated, implying that MMP-9 could promote arthritis in vivo.²⁵

CD44 can regulate cAMP-modulated signaling pathway, and under the stimulation of pro-inflammatory cytokines (TNF-α and IL-1), the level of CD44 in serum, plasma, and synovial fluid could be elevated.²⁶ Previous studies reported that CD44 was upregulated and accumulated on the surface of stimulated macrophages in RA.³ Hyaluronic acid (HA) that can bind to CD44 is used as a target for RA treatment.²⁷ In addition, researchers found that monoclonal antibodies (mAbs) against CD44 could significantly ameliorate collagen-induced arthritis (CIA) in vivo models.²⁸

HGF can affect PI3K-AKT pathway in focal adhesion and ECM receptor pathways and compared with normal controls or patients with osteoarthritis, the expression of HGF in serum, synovium, and articular adipose significantly increases in RA patients.²⁹ In addition, HGF levels in the serum can be used to predict the severity of joint damage in RA patients.³⁰ Studies have shown that HGF can cause angiogenesis in articulation and pannus formation in patients with RA. Inhibition of HGF/c-Met pathway can significantly inhibit inflammation and bone loss in articulation of RA mice.³¹

PPBP, a member of the CXC chemokine family, is a mitochondrial gene that plays a role in regulating the susceptibility to chronic inflammation-related diseases such as RA.³² Beta-thromboglobulin (beta-TG), an alias for PPBP, was found to be markedly upregulated in female RA patients, compared with the control group.³³

IGF-1 is an important mediator in bone development and remodeling, which can promote the synthesis process and inhibit degradation activities of extracellular matrix.³⁴ Some research suggested that IGF1 could be used as a marker for RA-SF macrophages.⁹ It has been reported that inflammation in hippocampus could interfere with the signal transduction of IGF1 receptors, resulting in neural complications of RA.³⁵ In clinical practice, the cachexia and critical condition of RA patients were related to the downregulation of IGF1,³⁶ while ameliorating inflammation by using anti-rheumatic therapies such as TNF inhibitors and corticosteroids rarely resulted in the recovery of IGF1 system.³⁷

VEGFA is a dominant endothelial growth factor, a core member of vascular formation, and can be stimulated by pro-inflammatory factors such as TNF-α, IL-1, and IL-6.³⁸ Studies have shown that compared with osteoarthritis or normal controls, the expression level of VEGFA increased significantly in synovial fluids in RA patients. In addition, VEGFA levels in the serum were correlated with the severity of RA, especially the number of swollen joints.³⁹ The VEGFA rs699947 C/A functional polymorphism may be related to the risk of elderly RA cases.⁴⁰

ALX (lipoxin A4 receptor), also known as FPR2, is a member of the Gi-protein coupled receptor (GiPCR) family and can be stimulated by Annexin A1-derived peptides, Annexin A1, lipoxin A4, or amyloid A, a reactant formed in the acute phase. Activated FPR can exert a powerful anti-inflammatory function in synovial tissues.⁴¹ FPRs have an important part in CIA and antigen-induced arthritis.⁴² An in-vivo experiment also indicated that the receptor had an inhibitory role in the FPR2-deficient mice with K/B × N serum transfer arthritis.⁴³

Although FN1, SERPINA1, and ITGAM have not been confirmed to be associated with RA, our results suggest that they might be potential biomarkers and are worthy of further investigation.

This research had some limitations. The sample size of this study was relatively small. The effects of comorbidity, gender, and medication (non-steroidal anti-inflammatory agents or methotrexate) were not eliminated, which might influence gene expression levels in the synovial macrophages in RA. In addition, macrophages extracted from RA samples should be used for in-vitro experiments to explore the molecular mechanisms related to gene expression. Preclinical in-vivo experiments should be conducted to explore the functions of the identified genes in the development of the disease. Future validation using polymerase chain reaction assays and Western blot and synovial macrophages in RA is required.

Conclusion

In summary, we identified 10 candidate genes associated with RA, including MMP9, VEGFA, HGF, PPBP, CD44, IGF1, FPR2, SERPINA1, FN1, and ITGAM, and explored the candidate genes and the relevant molecular mechanisms in RA patients. SERPINA1, FN1, and ITGAM, except for MMP9, VEGFA, HGF, PPBP, CD44, IGF1, and FPR2, were rarely mentioned in previous studies on RA. Our findings suggest that the 10 candidate genes might be used in diagnosis, prognosis, and therapy of RA in the future. Further studies are required to confirm the expression of these genes in synovial macrophages in RA and control specimen.

Highlights

● Integrated bioinformatics analysis of two profile datasets based on macrophage samples from RA patients and healthy controls was performed, and 173 upregulated and 106 downregulated common DEGs were identified.

● GO analysis showed that common DEGs were enriched in inflammatory response and immune response, and KEGG pathway illustrated that DEGs were mainly enriched in PI3K-Akt signaling pathway.

● FN1, VEGFA, HGF, SERPINA1, MMP9, PPBP, CD44, FPR2, IGF1, and ITGAM might be used in diagnosis, prognosis, and therapy of RA in the future.

Abbreviations

RA, rheumatoid arthritis; DEGs, differentially expressed genes; GEO, Gene Expression Omnibus; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; PPI, protein–protein interaction; HA, hyaluronic acid; mAbs, monoclonal antibodies; CIA, collagen-induced arthritis; GiPCR, Gi-protein coupled receptor.

Data Sharing Statement

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Ethics Approval and Consent to Participate

Not applicable because GEO belongs to public databases, the patients involved in the database have obtained ethical approval, users can download relevant data for free for research and publish relevant articles, and our study is based on open-source data, and Guangzhou University of Chinese Medicine and the First Affiliated Hospital of Guangzhou University of Chinese Medicine specifically do not require research using publicly available data to be submitted for review to their ethics committee, so there are no ethical issues and other conflicts of interest.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis, and interpretation, or in all these areas; took part in drafting, revising, or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Acknowledgments

Thanks to the support from the funding, and GEO database for providing their platforms and contributors for uploading their meaningful datasets.

Disclosure

The authors report no conflicts of interest in this work.

Additional information

Funding

This work was supported by grants from the National Natural Science Foundation of China (No. 81573850) and Young scientific research talents training project of The First Affiliated Hospital of Guangzhou University of Chinese Medicine (2017QN07). The Guangzhou Science Technology and Innovation Commission Technology Research Projects (No. 201904010336).