This article is mentioned in:

Introduction

Osteoarthritis (OA) is the most common arthritic disease in humans, affecting the majority of individuals over 65 years of age (1). OA most affects the joint including hands, knees, spine and hips, and is a leading musculoskeletal disease with declining joint functions (1). It takes a worldwide toll in the light of decreased physical ability, increased morbidity, and causes a substantial economic burden as well as serious socioeconomic consequences (2). Treatment of OA includes pain alleviation, the functional capacities maintenance and quality-of-life improvement (3). However, there are no effective interventions to decelerate the progression of OA since the precise mechanisms that involved in the pathogenesis of OA remain largely unknown.

OA has been recently considered as a multifactorial whole-joint disease that can affect the whole joint, including structural defects, cellular changes, and dysfunction of all compartments of the joint, such as cartilage, bone and synovium (4). Evidence has demonstrated that degradation and loss of articular cartilage, hypertrophic changes in bone, subchondral bone remodeling, and inflammation of the synovium are the main characteristics of this disease (5). Recently, Wang et al had showed that proteins of the complement system were differentially expressed in osteoarthritic synovial fluids compared to those from healthy individuals (6). Vance et al observed that highly expressed genes in the synovial fluid from patients with early meniscal injury having no OA symptoms were associated with OA and inflammation, such as Leukocyte-Associated Immunoglobulin-Like Receptor 1 (LAIR1) and Chemokine (C-C Motif) Receptor 6 (CCR6) (7). Although progresses have been gained about the pathogenesis of OA, concrete genetic mechanisms of OA remain to be elucidated.

In the present study, we downloaded the microarray data of GSE32317 from a public database. With the same microarray data, Zhu et al demonstrated that Tachykinin, Precursor 1 (TAC1) and the G protein-coupled receptor pathway might play significant roles in the progression of the early and late stages of OA (8). In addition, the work of Ma et al showed that the genes related immune response, cartilage development, and genes involve in the Toll-like receptor (TLR) signaling pathway and Wnt signaling pathway might be the potential target genes for the OA treatment (9). In this study, we analyzed the differentially expressed genes (DEGs) in samples of synovial membrane from patients with early stage of knee OA and late stage of knee OA compared with healthy specimens. Comprehensive bioinformatics analysis was applied to analyze the significant functions and pathways that were enriched by the common DEGs identified both in early-stage knee OA samples and late-stage knee OA samples. Besides, a protein-protein interaction (PPI) network was constructed and significant modules were extracted. Moreover, transcription factors (TFs) that could regulate genes in the significant modules were identified. We sought to gain more insights into the molecular circuitry in OA and identify more potential target genes for OA treatment.

Materials and methods

Microarray data

The microarray data of GSE32317, which was deposited by Wang et al (6), was downloaded from the NCBI Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE32317). The platform information is GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array. This dataset consisted of 26 samples of synovial membrane from the suprapatellar pouch of 10 early-stage knee OA patients who were undergoing arthroscopic procedures for degenerative meniscal tears (with documented cartilage degeneration but no full-thickness cartilage loss, Kellgren Lawrence score ≤2), and of 9 patients with end-stage knee OA (diffuse full thickness cartilage erosion) who were undergoing total knee joint replacement, and 7 healthy individual samples which were run on the same platform and array (from GEO accession number GSE12021) as our OA samples. In this study, the samples from early-stage knee OA patients, end-stage knee OA patients and healthy individuals were named as OA_Early, OA_End and Healthy, respectively.

Data preprocessing and screening of DEGs

The preprocessed microarray data that was conducted by robust multi-array average (RMA) algorithm (10) were obtained. Subsequently, probe ID was converted into gene symbol using the R/Bioconductor platform annotation packages hgu133plus2.db (11), org.Hs.eg.db (12), annotate (13). When several probes were mapped to a same gene symbol, the average value of these probes was calculated as the expression value of this gene. Then, 20389 gene expression matrix were obtained.

DEGs in OA_Early vs. Healthy and OA_End vs. Healthy were analyzed. Unpaired t test in limma package (14) was used to calculate P-values to determine the significance of DEGs. False discovery rate (FDR) (15) was applied to carry out the correction of multiple testing using Benjamini and Hochberg (BH) method (16). In this study, |log2fold change (FC) |≥1 and FDR <0.05 were selected as the threshold for DEGs screening. Moreover, the common DEGs in OA_Early vs. Healthy and OA_End vs. Healthy with the same gene change were considered to be potentially associated with the mechanism of OA. In this study, the common DEGs with same gene change in these two groups were identified and were our focus. Venn diagram was constructed using Venny 2.0 (http://bioinfogp.cnb.csic.es/tools/venny/).

Function and pathway enrichment analysis

ToppGene Suite (http://toppgene.cchmc.org) which is free and open to all users is a one-stop portal for gene list enrichment analysis, and candidate gene prioritization (17). ToppFun application from Toppgene suit can provide gene enrichment analysis in many biological categories, such as Gene Ontology (GO) terms and biological pathways (18). Specially, in the pathway analysis, ToppFun can use a comprehensive collection of pathways from several major databases such as KEGG, Reactome, and BioCarta (18). In this study, we undertook gene enrichment analysis of common upregulated genes and common downregulated genes using the ToppGene Suite ToppFun software (17). GO term or pathway was identified as significant under a BH multiple correction at a FDR value cut-off of 0.05.

PPI network construction

In this study, the PPI network was retrieved from the Search Tool for the Retrieval of Interacting Genes (STRING) database (19), which is a database providing ease of access to known and predicted protein interactions. Proteins included in the PPI network were all common DEGs with the same gene change in OA_Early vs. Healthy and OA_End vs. Healthy. Interaction pairs with a confidence score >0.7 were input to construct the PPI network. Additionally, Cytoscape software (6) was used for network visualizations, where nodes indicated proteins and edges represented interactions between any two proteins. Moreover, the degree of each node in this PPI network was calculated and the nodes with higher degree were considered as hub proteins.

Besides, significant modules in the PPI network were analyzed using Clustering with Overlapping Neighborhood Expansion (ClusterONE) tool, which is designed to identify densely connected regions (20). P<0.0001 was set as the cut-off value. In addition, pathway enrichment analysis was performed on the significant modules selected.

Search for candidate TFs for genes in significant modules

iRegulon (21), available as a Cytoscape plugin (6), implements a genome-wide ranking-and-recovery approach to detect enriched TF motifs and their optimal sets of direct target genes. In this study, we executed iRegulon and screened TFs for genes in the significant modules identified. The parameters were set as: Minimum identity between orthologous genes, 0.05; maximum false discovery rate on motif similarity, 0.001. The default Normalized Enrichment Score (NES) cut-off in iRegulon is set at 3.0, corresponding to FDR between 3% and 9% (21). In the present study, we selected the TF-gene pairs with NES above 5 for the follow-up analysis.

Results

DEGs screening

Compared with the Healthy samples, a total of 1,207 and 1,575 DEGs were identified in OA_Early and OA_End samples, respectively. The Venn diagrams representing upregulated genes and downregulated genes identified in OA_Early and OA_End were shown in Fig. 1A. From the results, we found that total 740 genes were upregulated in both OA_Early and OA_End samples. While a total of 308 common downregulated genes were identified in OA_Early and OA_End samples. The heat map of these common upregulated and downregulated genes was represented in Fig. 1B, showing that these common DEGs could significantly distinguish knee OA samples from Healthy samples.

Figure 1. Differentially expressed genes (DEGs) screening. (A) Venn diagrams representing upregulated genes and downregulated genes identified in OA_Early and OA_End samples. (B) heat map of common DEGs. Blue indicates decreased gene expression and sandybrown represents increased gene expression.

GO and KEGG pathway enrichment analysis

The top 5 over-represented GO terms, including molecular function (MF), biological processes (BP), and cellular component (CC) categories of the common up and downregulated genes were summarized in Tables I and II. Results showed that the common upregulated genes were significantly associated with immune response, receptor binding, and cell activation. Besides, the common downregulated genes were significantly enriched in different GO terms, such as structural constituent of muscle, structural molecule activity, and neuropeptide binding. On the other hand, the results of pathway enrichment analysis were shown in Table III (only most significant pathways were shown). In special, the common upregulated genes were found to be enriched in distinct pathways, such as rheumatoid arthritis. The downregulated genes were enriched in only two pathways, namely, GPCRs, Class A Rhodopsin-like and Striated Muscle Contraction.

Table I. Top 5 over-represented GO terms mapped to each category of common upregulated genes.

Table I.

Top 5 over-represented GO terms mapped to each category of common upregulated genes.

Category	GO ID	Name	FDR	Count
MF	GO:0038024	cargo receptor activity	2.13E-02	11
	GO:0004715	non-membrane spanning protein tyrosine kinase activity	2.13E-02	9
	GO:0005102	receptor binding	4.90E-02	79
BP	GO:0006955	immune response	2.22E-19	131
	GO:0050776	regulation of immune response	3.49E-15	85
	GO:0002682	regulation of immune system process	1.22E-14	108
	GO:0002684	positive regulation of immune system process	1.22E-14	79
	GO:0001775	cell activation	1.78E-14	90
CC	GO:0005764	lysosome	4.58E-10	52
	GO:0000323	lytic vacuole	4.58E-10	52
	GO:0005773	vacuole	5.51E-09	53
	GO:0005765	lysosomal membrane	7.78E-06	28
	GO:0005774	vacuolar membrane	2.57E-05	29

[i] Only three over-represented GO terms of MF category were enriched. GO, gene ontology; MF, molecular function; BP, biological process; CC, cell component; FDR, false discovery rate.

Table II. Top 5 over-represented GO terms mapped to each category of common downregulated genes.

Table II.

Top 5 over-represented GO terms mapped to each category of common downregulated genes.

Category	GO ID	Name	FDR	Count
MF	GO:0008307	structural constituent of muscle	2.68E-03	6
	GO:0005198	structural molecule activity	1.78E-02	19
	GO:0042923	neuropeptide binding	1.89E-02	3
	GO:0003779	actin binding	3.97E-02	13
	GO:0051373	FATZ binding	3.97E-02	2
BP	GO:0003012	muscle system process	7.23E-05	19
	GO:0055002	striated muscle cell development	1.94E-03	11
	GO:0055001	muscle cell development	2.73E-03	11
	GO:0006936	muscle contraction	5.84E-03	14
	GO:0018149	peptide cross-linking	9.77E-03	5
CC	GO:0030016	myofibril	1.31E-06	15
	GO:0030017	sarcomere	1.31E-06	14
	GO:0043292	contractile fiber	1.31E-06	15
	GO:0044449	contractile fiber part	2.22E-06	14
	GO:0031674	I band	6.58E-05	10

[i] GO, gene ontology; MF, molecular function; BP, biological process; CC, cell component; FDR, false discovery rate.

Table III. The most significant pathways enriched by common upregulated genes and common downregulated genes.

Table III.

The most significant pathways enriched by common upregulated genes and common downregulated genes.

Gene change	Database	ID	Name	P-value	Count
Up	BioSystems: KEGG	172846	Staphylococcus aureus infection	4.79E-09	18
	BioSystems: KEGG	200309	Rheumatoid arthritis	1.68E-07	20
	BioSystems: KEGG	213780	Tuberculosis	1.21E-05	25
	BioSystems: KEGG	128760	Intestinal immune network for IgA production	1.35E-04	12
	BioSystems: KEGG	99052	Lysosome	2.47E-04	18
	BioSystems: KEGG	83051	Cytokine-cytokine receptor interaction	7.82E-04	27
	MSigDB C2: BioCarta	M917	Complement pathway	9.05E-04	7
	BioSystems: KEGG	144181	Leishmaniasis	1.04E-03	13
	BioSystems: REACTOME	366160	Adaptive immune system	1.04E-03	48
	MSigDB C2: BioCarta	M7146	Classical complement pathway	1.16E-03	6
	BioSystems: KEGG	842771	Inflammatory bowel disease (IBD)	1.47E-03	12
	BioSystems: WikiPathways	198823	Complement activation, classical pathway	3.62E-03	6
Down	BioSystems: WikiPathways	198886	GPCRs, Class A Rhodopsin-like	1.24E-02	12
	BioSystems: WikiPathways	198903	Striated muscle contraction	1.84E-02	5

PPI network analysis

The PPI based on the common DEGs was constructed as shown in Fig. 2, consisting of 378 nodes (proteins) and 895 interactions (edges). Nodes with higher degree were Angiotensinogen (Serpin Peptidase Inhibitor, Clade A, Member 8) (AGT) (degree=45), Chemokine (C-X-C Motif) Ligand 12 (CXCL12) (degree=33), Jun Proto-Oncogene (JUN) (degree=32), Phosphoinositide-3-Kinase, Regulatory Subunit 1 (Alpha) (PIK3R1) (degree=28), Lysophosphatidic Acid Receptor 5 (LPAR5) (degree=27), Neuromedin U Receptor 2 (NMUR2) (degree=25), Regulator Of G-Protein Signaling 18 (RGS18) (degree=23), RGS19 (degree=23), Vav 1 Guanine Nucleotide Exchange Factor (VAV1) (degree=23).

Figure 2. Protein-protein interaction (PPI) network. Circular nodes indicate the common upregulated genes. Square nodes indicate the common downregulated genes. Bigger node represents node with higher degree.

Extraction of significant modules from the PPI network and TFs identification

Three significant modules were obtained. Besides, the TFs that regulated genes in these 3 modules were identified and the regulation network was constructed correspondingly (Fig. 3). From the results, we found that a total of 7, 8 and 14 TFs could regulate genes in module 1, module 2, and module 3, respectively. Thereinto, TFs identified to regulate genes in module 1 were not DEGs (Fig. 3A). While 2 DEGs were found to be TFs that regulated genes in module 2, including Lysine (K)-Specific Demethylase 2B (KDM2B, a common upregulated gene) and SRY (Sex Determining Region Y)-Box 8 (SOX8, a common upregulated gene) (Fig. 3B). Two DEGs were found to be TFs that regulated genes in module 3, including KDM2B and V-Maf Avian Musculoaponeurotic Fibrosarcoma Oncogene Homolog F (MAFF, a common downregulated gene) (Fig. 3C).

Figure 3. Regulation network construction and the enriched pathways of genes in modules (A) 1, (B) 2 and (C) 3. Triangle represents transcription factors. Rounded nodes in the regulation network indicate overlapping upregulated DEGs and square nodes indicate the overlapping downregulated DEGs.

Discussion

Gene expression profiling analysis revealed genes with abnormal expression concerned with OA, enabling the identification of targets for therapeutic options of OA. In the current study, 1,207 and 1,575 DEGs were identified in knee OA_Early and knee OA_End samples compared with healthy samples, respectively. Total 740 genes were upregulated and 308 genes were downregulated both in OA_Early and OA_End samples. These common DEGs were enriched in different GO terms and pathways. The upregulated genes AGT and CXCL12 were identified to be hub proteins in the PPI network or in the selected module 1. Besides, the upregulated DEG KDM2B was identified to be a TF that could regulate genes in the significant module 2 and module 3.

The risk factors for OA include biological factors and mechanical injury factors, such as aging, obesity and joint trauma (22). Meniscal tear is one of the risk factors for OA and is strongly associated with the development and progression of OA (23). In this study, synovial membranes of early knee OA patients were from those patients with meniscal tear. While the OA_End samples were derived from end-stage knee OA patients (diffuse full thickness cartilage erosion) who were undergoing total knee joint replacement regardless of the specific cause. To eliminate the influence of risk factor, we focused on the common DEGs both identified in the OA_Early and OA_End samples. Moreover, we chose 3 DEGs (AGT, CXCL12, and KDM2B) with higher node degree or having the function of TF to be our further focus.

The protein encoded by AGT, pre-angiotensinogen or angiotensinogen precursor, is cleaved by the renin enzyme in answer to lowered blood pressure (24). Then, the resulting product, angiotensin I, is cleaved by angiotensin converting enzyme (ACE) to generate the angiotensin II, a physiologically active enzyme (24). Evidence had demonstrated that local renin-angiotensin system components expressed particularly in hypertrophic chondrocytes and did not in hyaline chondrocytes (25). Recently, Yamagishi et al revealed that activation of the renin-angiotensin system may introduce OA (26). Evidence had shown that an ACE inhibitor (Ramipril) could improve vascular function in patients with rheumatoid arthritis (27). However, few study directly showed the effect of ACE inhibitor on OA. In our study, we found that AGT was an upregulated DEG both in the early stage of knee OA and late stage of knee OA. Besides, AGT was a hub protein in the PPI network. In line with the previous studies, it is reasonable to conclude that that AGT may play a key role in the development and progression of knee OA, which needs further validations.

CXCL12 functions as the natural ligand for the G-protein coupled receptor, chemokine (C-X-C motif) receptor 4 (CXCR4) and plays a role in many diverse cellular functions, including inflammation response (28). The work of Zhu et al showed that the G protein-coupled receptor pathway might play significant roles in the progression of the early and late stages of OA (8). Recently, He et al found that CXCL12 levels in the plasma and synovial fluid might serve as effective biomarkers for the severity of OA (29). In the present study, we found that CXCL12 was an upregulated DEG both in the early stage of knee OA and late stage of knee OA compared with controls. Besides, CXCL12 was another hub protein in the PPI network. Collectively, we suggest that CXCL12 may be essential in the pathogenesis of knee OA and it may be used as a biomarker of the disease severity as well as a potential drug target (30).

Furthermore, we found that the DEG KDM2B was identified to be a TF that could regulate genes in the significant module 2 and module 3 in our study. KDM2B encodes a member of the F-box protein family which function in phosphorylation-dependent ubiquitination (31). Farcas et al had demonstrated that an direct link between recognition of CpG islands by KDM2B and targeting of the polycomb repressive system (32). DNA methylation occurs naturally and mostly at positions where cytosine is bonded to guanine to form a CpG dinucleotide (32). Moreover, the involvement of DNA methylation has been most studied in the context of OA in musculoskeletal diseases (33). Taken together, we suggested that KDM2B might play a critical role in the pathogenesis of knee OA and KDM2B may be used as a novel drug target. However, experimental verifications are needed to confirm this finding.

In conclusion, AGT, CXCL12, and KDM2B were identified to be key genes associated with the pathogenesis of knee OA. Though the results are promising, the sample size in our study is small and this study needs further experimental verifications. Thus, further investigations with larger samples in this direction may provide a more conclusive result and may emphasize the potential of these identified genes as targets for therapeutic strategy.

Acknowledgements

This study was supported by the Institutional Research Funds for The First Affiliated Hospital of Xi'an Jiaotong University (no. 2015YK27).

References