Identification of genes and pathways in the synovia of women with osteoarthritis by bioinformatics analysis.

Osteoarthritis (OA) has a high prevalence in female patients and sex may be a key factor affecting the progression of OA. The aim of the present study was to identify genetic signatures in the synovial membranes of female patients with OA and to elucidate the potential associated molecular mechanisms. The gene expression profiles of the GSE55457 and GSE55584 datasets were obtained from the Gene Expression Omnibus database. Data of two synovial membranes from normal female individuals (GSM1337306 and GSM1337310) and two synovial membranes from female patients affected by OA (GSM1337327 and GSM1337330) were obtained from the dataset GSE55457, and those of three synovial membranes from female patients affected by OA (GSM1339628, GSM1339629 and GSM1339632) were obtained from the dataset GSE55584. Differentially expressed genes (DEGs) were identified by using Morpheus software. Protein‑protein interaction (PPI) networks of the DEGs were constructed by using Cytoscape software. Subsequently, Gene Ontology (GO) function and Kyoto Encyclopaedia of Genes and Genomes (KEGG) pathway enrichment analyses of the top module of the PPI network were performed by using ClueGo. A total of 377 DEGs were identified in the synovial membranes of OA patients compared with those of normal individuals, including 164 upregulated and 213 downregulated genes. The top 10 hub genes were ubiquitin (UB)C, ribosomal protein (RP) L23A, mammalian target of rapamycin, heat shock protein 90 α family class A member 1, RPS28, RPL37A, RPS24, RPS4X, RPS18 and UBB. The results of the GO analysis indicated that the DEGs included in the top module of the PPI were mainly enriched in the terms 'nuclear‑transcribed mRNA catabolic process', 'nonsense mediated decay', and 'cytoplasmic translation and ribosomal small subunit biogenesis'. KEGG pathway analysis indicated that the DEGs included in the top one module were mainly enriched in the 'ribosome' pathway. The present study provides a systematic, molecular‑level understanding of the degeneration of the synovial membrane in the progression of OA in female patients. The hub genes and molecules associated with the synovial membrane may be used as biomarkers and therapeutic targets for the treatment of OA in female patients with OA.

Introduction

Osteoarthritis (OA) is a degenerative joint disease characterized by the erosion of articular cartilage and subchondral bone and inflammation of the synovial membrane (1,2). The etiology of OA is affected by factors including aging, injury, obesity and heredity. OA currently affects ~25% of the population aged >18 years and is expected to be the greatest cause of disability among patients aged >40 years (3,4).

Recent studies have increasingly recognized the importance of the synovial membrane in the progression of OA (5). The synovial membrane is significantly altered prior to the occurrence of visible cartilage degeneration (6). Prieto-Potin et al reported that the histological pattern of the synovial membrane in OA patients was characterized by hyperplasia of synovial lining and stromal vascularization (7). Inflammatory cells (including macrophages, as well as B- and T-lymphocytes) and cytokines [including interleukin (IL)-6, IL-1β and tumor necrosis factor-α] have been detected in the synovial membranes of OA patients (8–10). It has been demonstrated that the synovium is associated with pain and poor functioning of the knee (11). However, the detailed pathogenic mechanisms in the synovia in the progression of OA have remained poorly understood and there is currently no available intervention strategy to slow down disease progression in patients with OA.

Microarray technology has markedly expanded the capacity of investigators to examine the pathogenic processes of various diseases (12,13). This technology is an important tool for current functional genomic research. Numerous studies have used microarray technology to identify certain OA-specific proteins, including activating transcription factor 4, G protein-coupled receptor 18, S100 calcium-binding protein A9/A8, prostaglandin D2 synthase and stromal cell-derived factor 1 (14,15). However, the sex of a patient may have a key role in the development of OA: Epidemiological studies have revealed a higher incidence of OA in females than in males (16,17). Furthermore, the use of certain medicines, including non-steroidal anti-inflammatory drugs or methotrexate, may also affect gene expression in synovial membranes affected by OA. In the present study, to exclude the effects of sex and medication use on the synovial membrane, genetic expression data of samples from female patients with OA who did not take any medicines and from healthy females were extracted. A bioinformatics analysis was performed to explore key genes whose functions are important in the synovial membranes of female patients with OA. Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed, with the aim of gaining a better understanding of the molecular mechanisms of the synovia in female patients with OA.

Materials and methods

Microarray data

The Gene Expression Ominibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) stores original submitter-supplied records as well as curated DataSets. Two gene expression profiles (GSE55457 and GSE55584) were obtained from the GEO database. To eliminate the potential effect of sex and medication, only data of synovial membranes from healthy females and female patients with OA who did not take any medicines were used. Data of two normal synovial membranes (GSM1337306 and GSM1337310) and two synovial membranes affected by OA (GSM1337327 and GSM1337330) were obtained from the GSE55457 dataset, while the data of three synovial membranes affected by OA (GSM1339628, GSM1339629 and GSM1339632) were obtained from the GSE55584 dataset.

Identification of differentially expressed genes (DEGs)

Analysis was performed by using the matrix visualization and analysis platform Morpheus (https://software.broadinstitute.org/morpheus/). The expression of mRNAs with a signal-to-noise ratio of >2 or <-2 were defined as DEGs (18).

Protein-protein interaction (PPI) network analysis

To understand the functional interactions between these DEGs, a PPI network was constructed by using the web-based tool STRING (http://www.string-db.org). Subsequently, the PPI network was visualized by using Cytoscape software (http://www.cytoscape.org/). The top ten nodes ranked by the degree of interactions in the PPI network were selected. Subsequently, the plug-in Molecular Complex Detection (MCODE) was used to screen modules of the PPI network in Cytoscape. The selection criteria were set as follows: MCODE scores ≥3 and number of nodes ≥3.

GO and pathway enrichment analysis

In order to analyse the functions and pathways that the DEGs in the top module from the PPI network are involved in, GO and KEGG pathway enrichment analyses were performed using ClueGo (19). A cut-off of 0.4 was set for kappa score and terms including at least 3 genes were retrieved.

Results

Identification of DEGs

A total of 5 OA samples and 2 normal synovial samples were analyzed from the GSE55457 and GSE55584 datasets. The gene expression profiles were analyzed with Morpheus software. This analysis identified 164 up- and 213 downregulated DEGs in the OA samples compared with the normal samples. An expression heat map of the top 30 up- and downregulated DEGs is presented in Fig. 1.

Figure 1.

Heat map of the top 60 differentially expressed genes in patients with OA compared with those in a normal cohort (30 upregulated and 30 downregulated genes). Red, up-regulation; blue, down-regulation. Max, maximum; min, minimum; OA, osteoarthritis.

PPI network construction and module analysis

Based on the information in the STRING and Cytoscape databases, the top 10 hub nodes with the highest degree of interaction were screened. These hub genes were ubiquitin (UB) C, ribosomal protein (RP) L23A, mammalian target of rapamycin (mTOR), heat shock protein 90 α family class A member 1 (HSP90AA1), RPS28, RPL37A, RPS24, RPS4X, RPS18 and UBB (Table I). Among these genes, UBC had the highest degree of interaction with 97 nodes. Six modules from the PPI network satisfied the criteria of MCODE scores ≥3 and number of nodes ≥3 (Table II). The results of the GO analysis indicated that in the category ‘biological process’ (BP), the DEGs included in the top module (Fig. 2) were significantly enriched in the GO terms ‘nuclear-transcribed mRNA catabolic process’, ‘nonsense mediated decay’ (GO:000184), ‘cytoplasmic translation’ (GO:0002181) and ‘ribosomal small subunit biogenesis’ (GO:0042274) (Fig. 3). In the GO category ‘molecular function’ (MF), the DEGs included in the top module were significantly enriched in ‘ribosomal (r)RNA binding’ (GO:0019843) (Table III). The KEGG pathway enrichment analysis of the DEGs included in the top module revealed that these genes were mainly associated with the ribosome pathway (KEGG:03010) (Table III).

Table I.

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

Gene	Degree of interaction
UBC	97
RPL23A	20
MTOR	20
HSP90AA1	20
RPS28	17
RPL37A	16
RPS24	15
RPS4X	15
RPS18	14
UBB	13

All of the top ten genes were downregulated in the protein-protein interaction network. UBC, ubiquitin C; RPL23A, ribosomal protein L23a; MTOR, mechanistic target of rapamycin; HSP90AA1, heat shock protein 90 α family class A member 1; RPS28, ribosomal protein S28; RPL37A, ribosomal protein 37a; RPS24, ribosomal protein S24; RPS4X, ribosomal protein S4 X-Linked; RPS18, ribosomal protein 18; UBB, ubiquitin B.

Table II.

Six modules from the protein-protein interaction network.

Cluster no.	Score	Nodes (n)	Edges (n)	Node IDs
1	16.5	17	132	RPL19, RPS17, RPL23A, EEF1A1, RPS24, UBC, RPLP1, RPL37A, RPL22, RPLP2, RPL14, RPS4X, RPS28, RPS18, EEF2, PABPC1, UBB
2	4	10	18	ATG14, CSNK1A1, CHEK1, BECN1, MTOR, CSNK1E, CFLAR, GAPDH, YWHAB, UVRAG
3	3	3	3	PCBP1, SNRPD1, HNRNPUL1
4	3	3	3	BUB3, CLASP2, KIF18A
5	3	3	3	SEC31B, COPB1, ARF4
6	3	3	3	FPR1, PTGER3, GRM3

UBC, ubiquitin C; UBB, ubiquitin B; RPLP1, ribosomal protein lateral stalk subunit P1; RPLP2, ribosomal protein lateral stalk subunit P2; PABPC1, polyadenylate-binding protein; RPS28, ribosomal protein 28; RPS4X, ribosomal protein S4 X-Linked; RPS24, ribosomal protein S24; RPS17, ribosomal protein S17; RPS18, ribosomal protein S18; RPL23A, ribosomal protein L23a; RPL14, ribosomal protein L14; RPL19, ribosomal protein L19; RPL37A, ribosomal protein L37a; RPL22, ribosomal protein L22; EEF1A1, elongation factor 1-α 1; EEF2, eukaryotic translation elongation factor 2; ATG14, beclin 1-associated autophagy-related key regulator; CSNK1A1, casein kinase 1 α 1; CHEK1, checkpoint kinase 1; BECN1, beclin 1; MTOR, mechanistic target of rapamycin; CSNK1E, casein kinase 1 epsilon; CFLAR, CASP8 and FADD like apoptosis regulator; YWHAB, tyrosine 3-monooxygenase/tyrosine 5-monooxygenase activation protein β; UVRAG, UV radiation resistance associated gene; PCBP1, poly(rC) binding protein 1; SNRPD1, small nuclear ribonucleoprotein D1 polypeptide; HNRNPUL1, heterogeneous nuclear ribonucleoprotein U like 1; BUB3, budding uninhibited by benzimidazoles 3; CLASP2, cytoplasmic linker associated protein 2; KIF18A, kinesin family member 18 A; COPB1, coatomer protein complex subunit β 1; ARF4, ADP ribosylation factor 4; FPR1, formyl peptide receptor 1; PTGER3, prostaglandin E receptor 3; GRM3, glutamate metabotropic receptor 3.

Figure 2.

Top module from the protein-protein interaction network. UBC, ubiquitin C; UBB, ubiquitin B; RPLP1, ribosomal protein lateral stalk subunit P1; RPLP2, ribosomal protein lateral stalk subunit P2; PABPC1, polyadenylate-binding protein; RPS28, ribosomal protein 28; RPS4X, ribosomal protein S4, X-Linked; RPS24, ribosomal protein S24; RPS17, ribosomal protein S17; RPS18, ribosomal protein S18; RPL23A, ribosomal protein L23a; RPL14, ribosomal protein L14; RPL19, ribosomal protein L19; RPL37A, ribosomal protein L37a; RPL22, ribosomal protein L22; EEF1A1, elongation factor 1-α 1; EEF2, eukaryotic translation elongation factor 2.

Figure 3.

Functional distribution of Gene Ontology terms in the category ‘biological processes’ for the top module from the protein-protein interaction network. **P<0.01. rRNA, ribosomal RNA; ER, endoplasmic reticulum; SRP, signal recognition particle; number, associated genes with the biological function.

Table III.

Analysis of GO terms in the category MF and KEGG pathways in the top module in the protein-protein interaction network.

Term	GO ID	Function	Associated genes
MF	GO:0019843	rRNA binding	EEF2, RPL19, RPL23A, RPS18, RPS4X
KEGG pathway	GO:0003010	Ribosome	RPL14, RPL19, RPL22, RPL23A, RPL37A, RPLP1, RPLP2, RPS17, RPS18, RPS24, RPS28, RPS4X

MF, molecular function; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; rRNA, ribosomal RNA; RPLP1, ribosomal protein lateral stalk subunit P1; RPLP2, ribosomal protein lateral stalk subunit P2; RPS28, ribosomal protein 28; RPS4X, ribosomal protein S4 X-Linked; RPS24, ribosomal protein S24; RPS17, ribosomal protein S17; RPS18, ribosomal protein S18; RPL23A, ribosomal protein L23a; RPL37A, ribosomal protein L37a; RPL14, ribosomal protein L14; RPL19, ribosomal protein L19; RPL22, ribosomal protein L22; EEF2, eukaryotic translation elongation factor 2.

Discussion

In the present study, the DEGs and associated PPIs in the synovial membranes of female patients with OA compared with those in healthy females were investigated. These results may contribute to a better understanding of the etiology of the synovial membrane in OA in female patients. In addition, the identified hub genes may act as potential candidate biomarkers for the diagnosis of OA. As presented in Table I, the top 10 hub genes in the PPI network included UBC, RPL23A, mTOR, HSP90AA1, RPS28, RPL37A, RPS24, RPS4X, RPS18 and UBB. All hub genes were downregulated in the synovia of female patients with OA.

UBC and UBB are genes encoding ubiquitin in the mammalian genome. UB is a small protein that is involved in the maintenance of chromatin structure and the regulation of gene expression. Based on previous studies reporting that UB is essential for the growth of cancer cells (20,21), the present study hypothesized that the upregulation of UBB and UBC may be a therapeutic strategy for maintaining the normal function of cells in the synovium.

mTOR is a major repressor of autophagy and it was previously demonstrated that mTOR is overexpressed in cartilage affected by OA (22). By inhibiting mTOR, authophagy is activated, which reduces cartilage damage and synovial inflammation in an experimental model of OA (22,23). The present bioinformatics analysis indicated that mTOR was downregulated in the synovium of OA patients, suggesting that autophagy was activated in the early stage of OA.

RPs are the major constituents of the ribosome complex and their functions are essential for cell growth, proliferation and homeostasis. RPL23A has been detected in the cytoplasm of synovial cells from patients with rheumatoid arthritis (RA) or OA as well as in healthy individuals, and that in RA patients, RPL23A was attached to T cells (24). Downregulation of RPL23A may also cause abnormal functioning of the synovium in OA patients. Tang and Wade (25) reported that, in zebra finches, the total number of cells expressing RPL37 was lower in females than in males, and declined with increasing age in females. These results may explain why OA often occurs in older females. RPS24 insufficiency inhibits the proliferation of cancer cells (26) and causes distinct cell cycle defects in diamond-blackfan anemia (27). The biological functions of RPS28, RPS18 and RPS4X in OA remain elusive at present. Anthony and Liebman (28) reported that alterations in RPS-18 may affect the translational accuracy in Saccharomyces cerevisiae. Furthermore, studies suggested that RPS28, RPS18 and RPS4X are associated with Diamond-Blackfan anemia (29–31). Based on these previous studies, it may be speculated that decreased levels of RPS28, RPS18 and RPS4X in synovial membranes may affect the functions of key proteins in rRNA processing, which may result in the degeneration of the synovium.

The HSP90AA1 gene encodes the Hsp90α protein. At present, the relevance of its expression to the progression of OA remains largely elusive. In a previous study, the upregulated expression of HSP90AA1 was demonstrated to be significantly associated with tumors and was critical for the stability of proteins that are vital for tumor progression (32). This result indicated that the level of HSP90AA1 downregulation may have an effect on cell viability in the synovium.

The ClueGo analysis provided the GO terms in the categories BP and MF, as well as KEGG pathway terms associated with the DEGs included in the top module. With regard to BP terms, the DEGs included in the top module of the PPI network were mainly enriched in ‘nuclear-transcribed mRNA catabolic process’, ‘nonsense mediated decay’, ‘cytoplasmic translation’ and ‘ribosomal small subunit biogenesis’, which may otherwise lead to the synthesis of dysfunctional proteins (33). The other two enriched BP terms, ‘cytoplasmic translation’ and ‘ribosomal small subunit biogenesis’, have important roles in protein formation. In the GO category MF, the DEGs from the top module of the PPI network were enriched in the term ‘rRNA binding’. KEGG pathway enrichment analysis indicated a significant role of ‘ribosome’ pathways in the degeneration of the synovial membrane.

As the association of most of the genes that were identified to be associated with OA in the present study were not previously reported, a further study will be performed to verify the gene expression levels of the screened genes in synovial samples from female patients with OA and healthy patients. Furthermore, cells extracted from those samples will be cultured to identify the molecular mechanisms of OA associated with these genes. In addition, the progression of synovial membrane degradation will be assessed in gene knockout rat models to further verify the functions of these genes.

The limitation of the present study is that only two normal synovial membrane samples and five OA synovial membrane samples were included in the analysis, and inclusion of other samples may change the present results. It is therefore necessary to collect more synovial membrane samples from female patients with OA to detect the expression levels of significantly DEGs.

In conclusion, the present study describes the molecular mechanism that may be involved in the degeneration of synovial membranes in female patients with OA. These mechanisms include ‘nuclear-transcribed mRNA catabolic process’, ‘cytoplasmic translation’ and ‘ribosomal small subunit biogenesis’. In addition, the present study reported the major hub genes that may be involved in the molecular mechanisms associated with the synovium in female patients with OA. Sex disparities should be taken into consideration when further molecular biological experiments on OA are performed.