Screening of gene signatures for rheumatoid arthritis and osteoarthritis based on bioinformatics analysis.

The current study aimed to identify gene signatures during rheumatoid arthritis (RA) and osteoarthritis (OA), and used these to elucidate the underlying modular mechanisms. Using the Gene Expression Omnibus database, the present study obtained the GSE7669 mRNA expression microarray data from RA and OA synovial fibroblasts (n=6 each). The differentially expressed genes (DEGs) in RA synovial samples compared with OA samples were identified using the Linear Models for Microarray Analysis package. The Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses were performed using the Database for Annotation Visualization and Integrated Discovery. A protein‑protein interaction network was constructed and the modules were further analyzed using the Molecular Complex Detection plugin of Cytoscape. A total of 181 DEGs were identified by comparing RA and OA synovial samples (96 up‑ and 85 downregulated genes). The significant DEGs in module 1, including collagen, type I, α 1 (COL1A1), COL3A1, COL4A1 and COL11A1, were predominantly enriched in the extracellular matrix (ECM)‑receptor interaction and focal adhesion pathways. Additionally, significant DEGs in module 2, including radical S‑adenosyl methionine domain containing 2 (RSAD2), 2'‑5'‑oligoadenylate synthetase 2 (OAS2), myxovirus (influenza virus) resistance 1 (MX1) and ISG15 ubiquitin‑like modifier (ISG15), were predominantly associated with immune function pathways. In conclusion, the present study indicated that RSAD2, OAS2, MX1 and ISG15 may be notable gene signatures in RA development via regulation of the immune response. COL3A1, COL4A1, COL1A1 and COL11A1 may be important gene signatures in OA development via involvement in the pathways of ECM-receptor interactions and focal adhesions.

Introduction

Rheumatoid arthritis (RA) is a chronic inflammatory disease which can lead to progressive disability, early mortality and other systemic complications (1). RA is incurable, has high socioeconomic costs and severely reduces the quality of life of patients (2). Osteoarthritis (OA) is the most common joint disease; it causes pain, disability and the loss of joint function (3). OA is also a chronic and incurable illness that lacks effective treatment strategies (4). Thus, advances in the understanding of the molecular mechanisms underlying these two diseases may lead to the development of novel therapeutic strategies.

RA is characterized by autoimmune and synovial inflammation, the destruction of multiple joints and the formation of pannus (5). OA is characterized by an inflammatory response and progressive breakdown of the articular cartilage of the joint (6). Genetic factors have been implicated in the development of these diseases (7,8). The human leukocyte antigen molecules and their relevant immunological pathways have been associated with RA pathogenesis (9). Sun et al (10) demonstrated that paired immunoglobin-like type receptor α was associated with inflammatory cell infiltration and was elevated in the synovial tissue from mice with RA. In addition, a genome-wide association and functional study suggests that DOT1-like histone H3K79 methyltransferase is associated with cartilage thickness and hip OA (11). Valdes et al (12) also confirmed that genetic variation in the SMAD family member 3 gene may result in the progression of hip and knee OA. Together, these findings indicate the importance of genetic mechanisms in the pathogenesis of RA and OA. Despite previous progress, the gene signatures associated with the pathogenesis of RA and OA remain unknown, and reliable predictive biomarkers for prognosis and treatment are lacking.

Microarray analyses have been increasingly used to identify disease-associated genes and pathways for elucidation of the molecular mechanisms of RA and OA (13,14). In a previous study, the GSE7669 microarray data was used to analyze differentially expressed genes (DEGs) between RA and OA using a gene co-expression network (15), or to screen candidate genes associated with RA by investigating core and periphery interaction structures (16). By contrast, the current study used this microarray data and comprehensive bioinformatics methods to identify DEGs in synovial RA samples compared with OA samples. Additionally, the present study performed functional enrichment analysis for DEGs and functional module analysis of the protein-protein interaction (PPI) network. The current study aimed to identify important disease-associated genes and the molecular mechanisms involved in RA and OA.

Materials and methods

Affymetrix microarray data

The GSE7669 gene expression profile deposited by Pohlers et al (17) was downloaded from Gene Expression Omnibus (www.ncbi.nlm.nih.gov/geo/), which was based on the platform of Affymetrix Human Genome U95 Version 2 Array (Affymetrix, Inc., Santa Clara, CA, USA). This dataset included the gene expression profiles from the synovial fibroblasts of 6 patients with RA and 6 patients with OA.

Data preprocessing and DEG screening

All the raw expression data was preprocessed using the Affymetrix package (18) in R (cran.at.r-project.org) and Bioconductor (www.bioconductor.org), and the normalization was performed using the robust multiarray average algorithm (19). The gene expression matrix of samples was acquired.

DEGs in RA synovial samples compared with OA samples were identified using the Linear Models for Microarray Analysis (Limma; www.bioconductor.org/packages/release/bioc/html/limma.html) package (20) in R/Bioconductor. t-test in the Limma package was used to analyze the P-value of each gene symbol. Only DEGs with P<0.05 and |log2 fold change|>0.5 were considered to indicate a statistically significant difference.

Functional enrichment analysis of DEGs

Gene Ontology (GO; www.geneontology.org) (21) is widely used in biology for the collation of large-scale gene lists, including biological process (BP) ontology. The Kyoto Encyclopedia of Genes and Genomes (KEGG; www.genome.ad.jp/kegg) (22) is used for extracting the pathway information from molecular interaction networks. To understand the biological significance of DEGs, GO BP enrichment analysis and KEGG pathway analysis were performed using the Database for Annotation Visualization and Integrated Discovery (DAVID; david.abcc.ncifcrf.gov) online tool (23). The P<0.05 in the hyper-geometric test and gene count >2 were defined as the cut-off values.

PPI network construction

The Search Tool for the Retrieval of Interacting Genes (STRING) (24) is a database providing information on experimentally verified and predicted protein interactions by calculating their combined score. Based on the information of the STRING database, DEGs with the combined protein interaction score >0.4 were selected and used to construct a PPI network in the present study. Hub proteins (25) in the PPI network were identified based on connectivity degree analysis. The PPI network was visualized using Cytoscape software (www.cytoscape.org) (26).

Functional module analysis of PPI network

The functional modules of the PPI network were subsequently identified using the Molecular Complex Detection (MCODE) (27) plugin of Cytoscape. The parameters were set as follows: Degree cut-off, 2; node score cut-off, 0.2; K-core, 2; and max depth, 100. Functional enrichment analyses for DEGs in functional modules with higher degree and node scores were subsequently performed using the DAVID tool. P<0.05 and gene count >2 were considered to indicate a statistically significant difference.

Results

DEG screening

Comparing OA and RA samples, a total of 181 DEGs were obtained, in which 96 genes were upregulated in RA and downregulated in OA samples, and 85 genes were downregulated in RA and upregulated in OA samples.

Functional enrichment analyses

The GO BP terms and KEGG pathway analyses were conducted for functional annotation of the DEGs. The top 10 GO BP terms are presented in Table I. The results demonstrated that the upregulated genes were significantly enriched in functions associated with cell adhesion, biological adhesion and extracellular matrix (ECM) organization. Downregulated genes were predominantly associated with the regulation of cell death, and responses to hormone and endogenous stimuli.

Table I

The top 10 significantly enriched GO biological progress terms.

GO term	Description	No. enriched genes	P-value
Upregulated
0007155	Cell adhesion	24	1.35×10−11
0022610	Biological adhesion	24	1.39×10−11
0007166	Cell surface receptor linked signal transduction	19	2.74×10−2
0043062	Extracellular structure organization	17	1.99×10−15
0030198	Extracellular matrix organization	16	4.32×10−17
0001568	Blood vessel development	12	2.48×10−7
0001944	Vasculature development	12	3.16×10−7
0001501	Skeletal system development	12	3.32×10−6
0030199	Collagen fibril organization	10	6.38×10−14
0010033	Response to organic substance	10	3.01×10−2
Downregulated
0010941	Regulation of cell death	9	4.73×10−2
0009725	Response to hormone stimulus	7	9.50×10−3
0009719	Response to endogenous stimulus	7	1.49×10−2
0010942	Positive regulation of cell death	7	2.05×10−2
0001558	Regulation of cell growth	6	2.71×10−3
0040008	Regulation of growth	6	2.67×10−2
0009615	Response to virus	5	2.07×10−3
0045596	Negative regulation of cell differentiation	5	2.22×10−2
0032870	Cellular response to hormone stimulus	4	2.82×10−2
0006820	Anion transport	4	3.39×10−2

GO, gene ontology.

In addition, 3 and 2 KEGG pathways were significantly enriched by up- and downregulated DEGs, respectively (Table II). The enriched KEGG pathways for the upregulated genes were ECM-receptor interactions, focal adhesions and the TGF-β signaling pathway, whereas the pathways enriched by the downregulated genes included the TGF-β signaling pathway and vascular smooth muscle contraction (Table II).

Table II

The significantly enriched KEGG pathways.

KEGG pathway term	Description	No. enriched genes	P-value
Upregulated
hsa04512	Extracellular matrix-receptor interaction	11	3.77×10−10
hsa04510	Focal adhesion	13	1.60×10−8
hsa04350	Transforming growth factor-β signaling pathway	4	2.84×10−2
Downregulated
hsa04350	Transforming growth factor-β signaling pathway	4	1.28×10−2
hsa04270	Vascular smooth muscle contraction	4	2.51×10−2

KEGG, Kyoto Encyclopedia of Genes and Genomes.

PPI network analysis

Based on the information of the STRING database, a total of 343 protein interactions with combined scores >0.4 were included in the PPI network (Fig. 1). The top 20 hub proteins were identified according to connectivity degree (Table III), including collagen type I α 1 (COL1A1), COL1A2, COL3A1, COL4A2, integrin α 1, COL4A1, biglycan and COL11A1. Notably, in addition to interferon-induced protein with tetratricopeptide repeats 3, these hub nodes were upregulated and the majority of them are collagen proteins.

Figure 1

Protein-protein interaction network of DEGs. Red nodes represent upregulated genes in RA samples (downregulated in OA samples). Green nodes represent downregulated genes in RA samples. Color depth represents level of significance, with deeper colors indicating greater significance. DEG, differentially expressed gene; RA, rheumatoid arthritis; OA, osteoarthritis.

Table III

The top 20 differentially expressed genes with higher connectivity degree in the protein-protein interaction network.

Gene	Degree
Collagen, type I, α 1	28
Collagen, type I, α 2	27
Collagen, type III, α 1	27
Collagen, type IV, α 2	22
Integrin subunit α 1	17
Collagen, type IV, α 1	17
Biglycan	17
Collagen, type XI, α 1	16
Versican	15
Collagen, type IV, α 2	15
Secreted protein, acidic, cysteine-rich	15
Periostin, osteoblast specific factor	15
Elastin	13
Cadherin 11	12
Interferon induced protein
with tetratricopeptide repeats 1	12
Fibromodulin	12
Collagen, type XIV, α 1	12
Interferon induced protein
with tetratricopeptide repeats 3	11
Collagen, type XV, α 1	11
Collagen, type XVI, α 1	11

In order to improve analysis of the PPI network, 5 modules were detected using the MCODE plugin. Only 2 module scores were >10 and the scores of other modules were not >3.5. Module 1, with the highest node score (11.455), consisted of 12 nodes and 63 edges, in which upregulated genes, including COL3A1, COL4A1, COL1A1 and COL11A1 were included (Fig. 2A). Module 2 (node score=11) consisted of 11 nodes and 55 edges, in which the downregulated genes were included (Fig. 2B). The downregulated genes in module 2 included radical S-adenosyl methionine domain containing 2 (RSAD2), 2′-5′-oligoadenylate synthetase 2 (OAS2), myxovirus (influenza virus) resistance 1 (MX1) and ISG15 ubiquitin-like modifier (ISG15).

Figure 2

Protein-protein interaction network of DEGs in (A) Module 1 and (B) Module 2. Red nodes represent upregulated genes in RA samples (down-regulated genes in OA samples). Green nodes represent downregulated genes in RA samples. Color depth represents level of significance, with deeper colors indicating greater significance. DEG, differentially expressed gene; RA, rheumatoid arthritis; OA, osteoarthritis.

In addition, functional enrichment analyses for DEGs in the functional modules were performed. Module 1 was predominantly associated with the regulation of cell proliferation, response to wounding and wound healing (Table IV). Module 2 was predominantly associated with immune functions, including the response to viruses, defense response and immune response (Table IV). The significantly enriched pathways of module 1 were ECM-receptor interaction and focal adhesion, whereas, no pathways were significant for module 2 (Table IV).

Table IV

The significantly enriched GO biological progress terms for modules.

GO term	Description	No. enriched genes	P-value
Module 1
0042127	Regulation of cell proliferation	12	6.47×10−10
0009611	Response to wounding	10	4.98×10−9
0042060	Wound healing	7	1.89×10−7
0001666	Response to hypoxia	7	3.28×10−7
0070482	Response to oxygen levels	7	4.70×10−7
0009725	Response to hormone stimulus	9	4.75×10−7
0031099	Regeneration	6	8.12×10−7
0009719	Response to endogenous stimulus	9	1.15×10−6
0008284	Positive regulation of cell proliferation	8	1.20×10−6
0001817	Regulation of cytokine production	6	3.26×10−6
Module 2
0009615	Response to virus	5	5.91×10−8
0006952	Defense response	3	2.74×10−2
0006955	Immune response	3	3.40×10−2

GO, gene ontology.

Discussion

There is a lack of effective treatments for RA and OA, and exploring the gene signatures associated with the diseases may elucidate the molecular pathogenesis and provide opportunities for biomarker development. The results of the present study demonstrated that the significant DEGs of module 1, including COL3A1, COL4A1, COL1A1 and COL11A1, were predominantly enriched in the ECM-receptor interaction and focal adhesion pathways. Additionally, the DEGs of module 2, including RSAD2, OAS2, MX1 and ISG15, were predominantly associated with immune responses. All these DEGs and pathways may be important mechanisms in the development of RA and OA.

In a previous study, Olex et al (28) identified the ECM-receptor interaction pathway to be an important signaling and metabolic pathway during the progression of OA. Koelling et al (29) also demonstrated that several dysregulated genes in OA samples were associated with ECM-receptor interactions and focal adhesions. In the present study, DEGs, including COL3A1, COL4A1, COL1A1 and COL11A1, were predominantly enriched in the ECM-receptor interaction and focal adhesion pathways. It can, therefore, be speculated that the these pathways may contribute to OA progression. Additionally, COL3A1 is a gene important for cartilage function, and its expression was previously observed to be correlated with the radiographic severity of canine elbow OA (30). Cui et al (31) demonstrated that COL3A1 expression was enriched in the focal adhesion pathway, which may suggest a molecular mechanism of OA. Furthermore, COL4A1 was previously identified as an OA-associated gene involved in the development of the disease (32). Gene expression analysis also indicated that COL1A1 was dysregulated in TGF-β-stimulated OA samples (33). COL11A1 was demonstrated to be an OA susceptibility gene in human joint tissues and is important in the development of this degenerative musculoskeletal disease (34). Additionally, a previous investigation demonstrated that multiple collagen genes (COL1A1, COL2A1, COL3A1 and COL4A1) were associated with the progression of OA (28). Collagen derivatives are candidates for disease-modifying OA drugs and are marketed as having therapeutic effects on reducing the symptoms of OA (35). Therefore, it is speculated that these collagens may be important gene signatures of OA and have implications in the progression of this disease through effects on ECM-receptor interactions and focal adhesions.

Furthermore, the innate immune response is understood to cause inflammation and joint destruction in RA, and a genome-wide association study has demonstrated that immune regulatory factors underlie this disease (36). Additionally, the potential of Tank-binding kinase 1 as a therapeutic target in RA suggests that there is an association between the synoviocyte innate immune responses and RA development (37). McInnes and Schett (1) also investigated the importance of immune responses in the pathogenesis of RA. In the current study, GO functions associated with the immune response were enriched. Thus, the results of the present study support previous findings and suggest that immune responses may contribute to RA development.

Additionally, functional enrichment analysis performed in the present study demonstrated that DEGs in RA samples, including RSAD2, OAS2, MX1 and ISG15, were significantly enriched in the immune response pathway. RSAD2 is a type I interferon (IFN) response gene, and has been used in the clinic for the prediction of RA development (38,39). OAS2 is involved in the IFN β signaling pathway, and investigation of core and periphery interaction structures previously identified it as a candidate gene associated with RA (16). MX1 is also an IFN response gene and was previously demonstrated to be correlated with disease activity in fibroblast cells of RA synovial tissue (40). ISG15 sensitizes the IFN-activated JAK-STAT pathway, which is established to be important in RA development (41). Additionally, type I IFN has been demonstrated to enhance immune responses in vivo and acts as a signal linking innate and adaptive immunity (42). Smith et al (43) demonstrated that, as type I IFN signatures, RSAD2, OAS2 and MX1 may be important for predicting treatment response in RA. Therefore, the DEGs identified in the present study may be crucial gene signatures for elucidating the molecular pathogenesis of RA.

In conclusion, the results of the current study indicate that RSAD2, OAS2, MX1 and ISG15 may be RA gene signatures, and may be associated with RA development via effects on immune responses. COL3A1, COL4A1, COL1A1 and COL11A1 may be important gene signatures contributing to OA development via involvement in ECM-receptor interactions and focal adhesion. The present findings aid the clarification of the molecular mechanisms of RA and OA. However, the sample size used in the current study was small. Additional experiments, including reverse transcription-quantitative polymerase chain reaction and western blot analysis, were not performed to confirm the mRNA and protein expression levels of these gene signatures. Thus, further studies are required to investigate the potential clinical applications of these gene signatures.

Table V

Significantly enriched KEGG pathways for modules.

Term	Description	Count	P-value
Module 1
hsa04512	ECM-receptor interaction	8	2.61E-13
hsa04510	Focal adhesion	8	1.36E-10

No significantly enriched pathways were identified for Module 2. Description represents the name of KEGG pathways and count represents the number of enriched genes. KEGG, Kyoto Encyclopedia of Genes and Genomes.