Identification of genes and pathways associated with subchondral bone in osteoarthritis via bioinformatic analysis.

Osteoarthritis (OA) is a high prevalent musculoskeletal problem, which can cause severe pain, constitute a huge social and economic burden, and seriously damage the quality of life. This study was intended to identify genetic characteristics of subchondral bone in patients with OA and to elucidate the potential molecular mechanisms involved. Data of gene expression profiles (GSE51588), which contained 40 OA samples and 10 normal samples, was obtained from the Gene Expression Omnibus (GEO). The raw data were integrated to obtain differentially expressed genes (DEGs) and were further analyzed with bioinformatic analysis. The protein-protein interaction (PPI) networks were built and analyzed via Search Tool for the Retrieval of Interacting Genes (STRING). The significant modules and hub genes were identified via Cytoscape. Moreover, Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) enrichment analysis were performed. Totally 235 DEGs were differentially expressed in the subchondral bone from OA patients compared with those of normal individuals, of which 78 were upregulated and 157 were downregulated. Eight hub genes were identified, including DEFA4, ARG1, LTF, RETN, PGLYRP1, OLFM4, ORM1, and BPI. The enrichment analyses of the DEGs and significant modules indicated that DEGs were mainly involved in inflammatory response, extracellular space, RAGE receptor binding, and amoebiasis pathway. The present study provides a novel and in-depth understanding of pathogenesis of the OA subchondral bone at molecular level. DEFA4, ARG1, LTF, RETN, PGLYRP1, OLFM4, ORM1, and BPI may be the new candidate targets for diagnosis and therapies on patients with OA in the future.

Introduction

Osteoarthritis (OA) is a common disorder form of arthritis, which can cause severe pain and stiffness of joints, constitute a serious social and economic burden, and seriously damage the quality of life among older individuals.[ The etiology is still enigmatic and it is generally acknowledged that the occurrence of the disease is associated with genetic predisposition and environmental factors. Several key cellular and molecular mechanisms of pathological progress of OA have been clarified. In addition to cartilage degeneration, which is the main pathological changes in OA, it also involves the entire articular structure, including inflammation, chondrocyte apoptosis, pathological transformation of synovitis and meniscus, structural remodeling of subchondral bone.[

Stable joint structure of cartilage and bone ensures normal joint weight-bearing. The disruption of this physiological relationship accelerates the pathological progress of OA.[ Subchondral bone comprises subchondral sponge and subchondral plate. Accumulating evidence demonstrate that it makes considerable contributions to the occurrence and development of OA.[ The significant differences in the morphology of the subchondral bone plate were found between the OA and the normal and osteoporosis through electron microscopic analysis, which strongly suggested the abnormal cell activity.[ Kraus et al have proved that the texture of subchondral bone can be applied as a biomarker for predicting the progress of knee OA.[ An animal study indicated that subchondral bone properties do influence load-induced OA progression.[ In addition, osteocyte culture studies have demonstrated that osteocytes may affect cartilage metabolism by releasing factors.[ These studies offer a new perspective for the molecular mechanism of subchondral structural degeneration and articular cartilage changes in the progress of OA. Therefore, to identify the potential mechanisms of OA involved in subchondral bone and develop effective biomarkers is the key to the prevention and treatment of OA.

Over the past few decades, whole-genome microarray and bioinformatic analysis have become a common technology for screening genetic alterations and studying the behavior of the differentially expressed genes (DEGs) simultaneously. In this work, the DEGs expressed in subchondral bone between OA and the normal were obtained from a microarray dataset through the Gene Expression Omnibus (GEO). To explore the potential mechanisms associated with subchondral bone in the development of OA at a molecular level, the protein–protein interaction (PPI) networks analysis, Gene Ontology (GO) and Kyoto Encyclopaedia of Genes and Genomes (KEGG) enrichment analysis among DEGs were carried out in turn. Our study explores the mechanism of OA associated with subchondral bone and the resulted candidate genes may be used as the targets to prevent, delay even block the progress of OA beta-defensins (HBDs).

Materials and methods

Data source and data preprocessing

The GSE51588 gene expression profile dataset, which was sequenced on the GPL13497 Agilent-026652 Whole Human Genome Microarray 4x44K v2 (Probe Name version), was obtained from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). According to the annotation information in the platform, the probes are transformed into corresponding gene symbols. GSE51588 included 10 normal-subchondral bone samples from knee tibial plateau (5 lateral vs 5 medial) and 40 OA-subchondral bone samples from knee tibial plateau (20 lateral vs 20 medial).[ The diagnosis of knee OA was based on the criteria for knee OA of the American College of Rheumatology.[ The Characteristic of the samples was shown in Table 1. The normalization of datasets was conducted by the preprocessCore package in R.

Table 1

Characteristics of the samples.

DEGs screening and hierarchical clustering analysis

After the data normalization, the DEGs between OA samples and normal samples were screened using the Multiple Linear Regression limma in R. A probe set having no corresponding gene symbol or a gene having a plurality of probe sets is deleted or averaged, respectively. Volcano plot was used to display both total DEGs and top 30 upregulated and downregulated DEGs. That was generated with ggplot2 package of R language. DEGs were identified using classical t test and statistically significant DEGs were defined with the criteria of |Log2(FC)| > 2 and adj. P < .05. A bidirectional hierarchical clustering heatmap of top 30 upregulated and downregulated DEGs were constructed using pheatmap package of R language after extracting the expression values from the gene expression profile.

KEGG and GO enrichment analyses of DEGs

The Database for Annotation, Visualization, and Integrated Discovery (DAVID; http://david.ncifcrf.gov) (version 6.8), an original web-accessible program that integrates and analyzes biological data, allows researchers to unravel the biological significance behind genes.[ KEGG is a bioinformatic resource for revealing gene advanced functions and connecting genomic information with functional information.[ GO is a tool for annotating and analyzing the biological function of large list of genes.[ The DAVID carried out the GO and KEGG enrichment analysis of the DEGs. The bubble plots were performed to demonstrate the results of enrichment analysis via ggplot2 package of R languages.

PPI network and module analyses

Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org) (version 11.0), an web-based tool that analyzes the functional interactions among proteins, was used to build the PPI network.[ Cytoscape (version 3.6.1) is a software providing an open platform for the visualization of biological processes and molecular interaction networks and integration of these networks. Molecular Complex Detection (MCODE) (version 1.5.1), a plugin of Cytoscpae, is applied to discover densely connected regions and identify the significant modules through clustering a given network. “MCODE scores ≥ 5,” “k-score = 2,” “Max depth = 100,” “node score cut-off = 0.2,” and “degree cut-off = 2” were set as the criteria for selection of significant modules. Finally, GO and KEGG enrichment analyses of significant modules were conducted via DAVID.

Hub genes selection and analyses

To obtain balance between the core genes and avoiding missing the key gene, we extracted the hub genes using cytoHubba. Through the cytoHubba plugin, 12 topological analysis methods were obtained. The top 20 hub-forming genes (10%) were identified based on MCC, MNC, and DMNC, respectively. Finally, we found the common genes via Venn diagrams and the overlapping genes were regarded as the hub genes.

Results

Normalization of datasets

The data of gene expression was downloaded from GEO database (GEO accession no. GSE511588). Normalization of the data of gene expression was performed using the preprocessCore package in R and presented in a boxplot (Fig. 1). The black lines in Figure 1A and B are basically at the same level, indicating a high consistency.

Figure 1

Normalization of samples. (A) Before normalization of total DEGs and (B) following normalization of total DEGs. DEGs = differentially expressed genes, OA = osteoarthritis.

Identification of DEGs and hierarchical clustering analysis

According to the established criteria, a total of 235 DEGs were obtained, consisting of 157 downregulated genes and 78 upregulated genes between OA subchondral tissues and non-OA subchondral tissues. An expression volcano plots and heatmap of the top 30 upregulated and downregulated DEGs were indicated in Figures 2 and 3, which demonstrated that the expression of identified DEGs could correctly distinguish the two kinds of samples.

Figure 2

Volcano plot of DEGs and top 30 upregulated and downregulated DEGs. Grey represents not-significant change in expression. Red represents up-regulation. Blue represents down-regulation. DEGs = differentially expressed genes.

Figure 3

Heatmap of the top 30 differentially expressed genes in patients with OA compared with those in a normal cohort. Red represents up-regulation. Blue represents down-regulation. Max = maximum, min = minimum, OA = osteoarthritis.

Functional enrichment analyses were conducted via DAVID to evaluate the biological classification of DEGs, as illustrated in Figure 4. GO analysis demonstrated that in the biological process (BP) ontology, the DEGs were significantly involved in the regulation of inflammatory immune response items, including defense response to bacterium (15 genes), immune response (21 genes), defense response to fungus (7 genes), inflammatory response (19 genes), and response to lipopolysaccharide (12 genes). In the cellular component (CC) ontology, the DEGs were significantly involved in the extracellular space (63 genes), extracellular region (67 genes), and extracellular exosome (61 genes), which suggested that extracellular components were closely related. In the molecular function (MF) ontology, the DEGs were significantly involved in the binding-related items, including RAGE receptor binding (6 genes), heparin binding (13 genes), oxygen binding (6 genes), and heme binding (9 genes). Besides, the other enriched category was associated with oxygen transporter activity (5 genes). Furthermore, KEGG analysis indicated that the DEGs primarily involved in inflammatory immune response, including amoebiasis (10 genes), protein digestion and absorption (8 genes), nitrogen metabolism (3 genes), cytokine–cytokine receptor interaction (8 genes), and Asthma (3 genes).

Figure 4

Bubble plot of GO and KEGG pathway analysis of DEGs in subchondral bone associated with osteoarthritis. (A) Top 5 significantly enriched biological processes in DEGs. (B) Top 5 significantly enriched cell component in DEGs. (C) Top 5 significantly enriched KEGG pathway in DEGs. (D) Top 5 significantly enriched molecular function in DEGs. RAGE = receptor of advanced glycation endproducts.

The information in the STRING was imported into the Cytoscape and the PPI network of DEGs was set up. Five significant modules were selected according to MCODE scores ≥ 5 (Table 2) and the top module was obtained as shown in Figure 5. Finally, the functional enrichment analyses of the significant modules were performed via DAVID (Table 3).

Table 2

Five significant modules from the protein–protein interaction network.

Figure 5

The most significant module was obtained from the PPI network with 17 nodes and 136 edges. Low value of combined score to circle with small size or bright colors and line with small sizes.

Table 3

GO and KEGG pathway enrichment analysis of DEGs in significant modules.

Hub gene selection and analysis

After identification of the top 20 genes based on MCC, MNC, and DMNC, respectively, the overlapping 8 genes were selected as the hub genes (Fig. 6). The basic characteristics of these hub genes are shown in Table 4.

Figure 6

Venn diagram of the hub genes. Eight hub genes were selected based on the overlap of top 20 genes identified by MCC, MNC, and DMNC.

Table 4

Functional roles of 8 hub genes.

Discussion

OA has long been considered as a degenerative disease with high prevalence and financial burden.[ Moreover, current treatments for OA are usually palliative, mainly pain management with anti-inflammatory drugs and analgesics. Regulators have not yet approved any pharmacological therapies with convincing disease-improving effects.[ Therefore, it is quite necessary to develop new, appropriate and effective therapies that can relieve, suspend, or even reverse the occurrence and progression of OA.

The microarray technology allows us to identify genetic changes in the development of OA and is considered as an effective way to deepen the understanding of the pathological mechanisms and detect new markers for early diagnosis and precise treatment in OA. To date, all of the bioinformatic analyses in OA have been conducted on human synovium, meniscus and articular cartilage; nevertheless, there is no study on human subchondral bone tissue.

In this work, 235 DEGs were identified in subchondral bones from OA patients compared with those of normal subjects, of which 157 were downregulated and 78 were upregulated. Then, functional enrichment analyses were conducted and the interactions between DEGs were discussed. GO analysis showed that DEGs were primarily involved in defense response to bacterium, immune response, inflammatory response, response to lipopolysaccharide, extracellular space, extracellular region, RAGE receptor binding and heparin binding.

This is consistent with our knowledge that the regulation of inflammatory immune response is the major mechanism for the pathogenesis of OA.[ The involvement of the immune system is one of the key factors in the development and progression of the disease.[ According to the previous studies, lipopolysaccharide (LPS) can be transferred from the intestinal tract into blood circulation, resulting in low-grade inflammation. OA is a low-grade inflammatory disease, and the increase of LPS related to obesity and metabolic syndrome could cause OA by priming the proinflammatory innate immune response.[ Recent work has indicated that the spread of substances through the synovial fluid or extracellular space may contribute to the cell-to-cell communication among cartilage, bone, and the synovial membrane. Some factors secreted by synovial fibroblasts could accelerate the progression of joint disease.[ It is worth noting that the pathophysiological processes in OA are mainly mediated by exosomes from synovial fluid or extracellular space.[ Besides, related heparin binding has already been reported to contribute a lot to the catabolic–anabolic imbalance seen in OA cartilage.[ VEGF production and inflammatory responses could be induced via RAGE receptor binding in human synoviocytes.[ In a word, all these studies confirm our findings. Previous studies indicated that these pathways may be involved in immune mediated mechanisms activating inflammatory response in OA, which reflecting the pathogenesis of multifaceted origin.[

After identification of the top 20 genes based on MCC, MNC, and DMNC, respectively, we selected 8 DEGs as hub genes based on the intersection of three groups of genes. Among these hub genes, DEFA4, a member of the family of antimicrobial and cytotoxic peptides, is considered to mediate host defense and are abundant in neutrophils granules.[ However, studies have demonstrated that beta-defensins (HBDs) are multifunctional antimicrobial peptides that are capable of linking inflammation and host defense mechanisms with tissue-remodeling processes in articular cartilage, which suggested that HBDs are additional factors in the pathogenesis of OA.[ ARG1, is a critical regulator of innate and adaptive immune responses through arginine metabolism. Burleigh et al indicated that ARG1 may associated with the mechanosensing mechanisms of OA.[ Absence of CCL2 strongly suppressed ARG1 and the CCL2/CCR2 axis make great contributions to the progression of pain in murine OA.[ LTF, is an important component of the non-specific immune system, including anti-inflammatory activity and host defense against a broad range of microbial infections. It has anabolic, differentiating and anti-apoptotic effects on osteoblasts and can also inhibit the formation of osteoclasts, possibly playing a role in the regulation of bone growth.[ The significant decrease in LTF was observed in the murine mucopolysaccharidosis I model with early pathogenic events in extracellular matrix. Alteration of extracellular matrix is the basis of progressive bone and joint disease symptoms caused by biomechanical failure of chondro-osseous tissue.[ RETN, has been recognized as a significant link between OA, inflammation and obesity.[ A study enrolling 1002 participants with early symptomatic OA demonstrated that a positive correlation between plasma RETN and knee radiographic OA events.[ In De Boer et al, there was no correlation between cartilage damage and serum RETN level, but there was still a weak positive correlation between RETN and local synovial tissue inflammation,[ suggesting that RETN has a potential pro-inflammatory effect in the pathogenesis of OA. Besides, RETN can be secreted by the OA cartilage and is abundant in OA joints, which is positively correlated with the degree of articular cartilage damage and the severity of symptoms and imaging of OA. Other studies have showed that RETN in synovial fluid from patients with OA are positively associated with catabolic and inflammatory factors, including CTX-II, MMP-1/-3, and IL-6.[ In addition, a study indicated that recombinant RETN is capable of causing decrease in proteoglycan and promoting the synthesis of PGE2 in mouse cartilage, and similarly inhibiting production of proteoglycan in human.[ Recent studies indicated that RETN strongly promotes large amounts of sulfated glycosaminoglycan (sGAG) consumption in meniscal explants, which demonstrated that the RETN with the ability of catabolism and pro-inflammation in OA.[

The interaction among OA and hub genes PGLYRP1, OLFM4, ORM1, and BPI have not been reported yet. PGLYRP1 is a pattern receptor that has bactericidal activity against Gram-positive bacteria by binding to murein peptidoglycans (PGN) and interfering with the biosynthesis of peptidoglycan. It also has bacteriostatic activity against Gram-negative bacteria and plays an important role in innate immunity.[ OLFM4 encodes a member of the olfactomedin family. The encoded protein is an extracellular matrix glycoprotein that promotes cell adhesion and is an anti-apoptotic factor that accelerates tumor growth.[ ORM1 encodes a key acute phase plasma protein, Alpha-1-acid glycoprotein 1, which acts to modulate the activity of the immune system in the acute phase responses.[ BPI encodes a LPS binding protein, which has antimicrobial activity against gram-negative organisms. A study has shown that BPI in synovial fluid may play a role in the pathogenesis of arthropathies.[

This article also has some limitations. It is more frequent that the medial tibial plateau is OA than lateral side because of the effects of biological load axes in human, especially in varus situation. Due to the lack of detailed raw data, we have not distinguished the changes in biomechanics and molecular level caused by this factor, which may cause corresponding errors in the analysis results. Looking forward to more research in the future to make up for the shortcomings of this article. In summary, the present study is intended to screen out DEGs that might be associated with the degeneration of subchondral bone in OA patients and attempts to explore the possible molecular mechanism of the occurrence and progression in OA. Totally, there were 235 DEGs and 8 hub genes identified via bioinformatic analysis, which could provide a series of promising diagnostic biomarkers and therapeutic targets for OA. However, further studies are necessary to clarify the biological function of these genes in the pathogenesis of OA.

Author contributions

Conceptualization: Zhanyu Yang, Jiangdong Ni.

Data curation: Zhanyu Yang.

Investigation: Letian Kuang.

Methodology: Yongquan Gao.

Validation: Shibin Tao.

Writing – original draft: Shibin Tao.

Writing – review & editing: Zhanyu Yang.