Bioinformatics analysis of differentially expressed genes in subchondral bone in early experimental osteoarthritis using microarray data

Abstract

Background

Osteoarthritis (OA) is the most common arthritic disease in humans, affecting the majority of individuals over 65 years of age. The aim of this study is to identify the gene expression profile specific to subchondral bone in OA by comparing the different expression profiles in experimental and sham-operation groups.

Methods

Gene expression profile GSE30322 was downloaded from the Gene Expression Omnibus (GEO) database. Differentially expressed genes (DEGs) were obtained by limma package. And Database for Annotation, Visualization and Integrated Discovery (DAVID) databases were further used to identify the potential gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. Furthermore, a protein–protein interaction (PPI) network was constructed and significant modules were extracted.

Results

Totally, 588 DEGs were identified including 199 upregulated DEGs and 389 downregulated DEGs screened in OA and sham-operation. GO showed that DEGs were significantly enhanced for ribosomal subunit export from nucleus and molting cycle. KEGG pathway analysis revealed that target genes were enriched in thiamine metabolism.

Conclusion

These key candidate DEGs that affect the progression of OA, and these genes might serve as potential therapeutic targets for OA.

Introduction

Osteoarthritis (OA) is a degenerative disease characterized by the gradual degeneration of articular cartilage, joint stiffness, and loss of function [1]. It was reported that over 27 million adults are affected by OA in the USA [2]. OA is a complex pathophysiological process involving inflammation, subchondral bone modification, and osteophyte formation. Subchondral bone alteration present to the cartilage degeneration and thus more studies should be focused on the subchondral bone alteration.

Subchondral bone consists tripartite: subchondral bone plate, trabecular bone, and bone marrow space [3]. It has been stated that most of the OA patients accompanied by the alterations of the subchondral bone [4]. Subchondral bone could transport nutrients or cytokines to the overlying cartilage. Meanwhile, subchondral bone cells contacted with chondrocyte and thus influence cartilage metabolism. A better understanding of the early molecular mechanism changes of subchondral bone in vivo may contribute to elucidating the pathogenesis of OA. Therefore, it is crucial to explore the differentially expressed genes (DEGs) in vivo and thus we could revealed new targets for OA [5].

Microarray technology has been used to obtain information on the genetic alteration that occurs during many diseases [6, 7]. Here, we downloaded the gene expression profile GSE30322 from the Gene Expression Omnibus database (GEO), including gene expression data for subchondral bone samples from five medial meniscectomy and medial collateral ligament transection group and five sham-operated group. Based upon this research, identifying DEGs and enriching their functions and signaling pathways may help reveal potential targets of early OA.

Materials and methods

Gene expression microarray data

The gene expression profile GSE30322 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE30322) was downloaded from the Gene Expression Omnibus (GEO, www.ncbi.nlm.nih.gov/geo/). GSE30322 was based on Agilent-014879 Whole Rat Genome Microarray 4x44K G4131F (Probe Name version) platform. GSE30322 dataset contained ten samples, including five bone 4 weeks post-surgery samples (E-group), and five sham-operated group (S-group) 4 weeks post-surgery samples.

DEGs in E-group and intact S-group samples

The raw data files were downloaded and then python scripts for matrix transformation were used. The analysis was carried out using Limma package from Bioconductor project. In this study, genes with P < .05 and [log fold change (FC)] > 2 were defined as DEGs. The DEGs data were then processed by R software (pheatmap package) to draw a heatmap and volcano plot.

GO and KEGG analysis of DEGs

Target genes list were submitted to the DAVID 6.8 (https://david.ncifcrf.gov/) to analyze candidate DEG functions and Kyoto Encyclopedia of Genes and Genomes (KEGG) of the overlapping genes. DEG functions, also named as Gene ontology (GO), mainly including biological process (BP), molecular function (MF), and cellular component (CC). P value less than 0.05 was considered as cut-off criterion [8,9,10].

Protein-protein interaction (PPI)

We used the online database STRING (Search Tool for the Retrieval of Interacting Genes, https://string-db.org/) to better illustrate the potential interactive relationships among the DEGs [11]. Then, the Cytoscape software was utilized for analyzing the interactions with a combined score > 0.4 (http://www.cytoscape.org/). Finally, the plug-in Molecular Complex Detection (MCODE) was used to filter the significant modules from the PPI network for the selection of hub genes (degree cut-off = 2, node score cut-off = 0.2, k-core = 2, and max. depth = 100) [12].

Results

Identification of DEGs

After analyzing, differentially expression gene profiles were obtained. Totally, 588 DEGs were identified including 199 upregulated DEGs and 389 downregulated DEGs screened in OA and sham-operation. Top 10 up-DEGs and down-DEGs were listed in Table 1 and Table 2, respectively. A box plot of the sample data is provided in Fig. 1. Volcano plot of the different genes can be obtained in Fig. 2. Moreover, we provided heatmap of the top 50 different genes between E-group and S-group (Fig. 3).

Table 1 The top 10 upregulated DEGs in early experimental osteoarthritis with P value < 0.05

Table 2 The top 10 downregulated DEGs in early experimental osteoarthritis with P value < 0.05

Fig. 1

Box plot for the sample data after normalization

Fig. 2

Volcano plot of the different genes in E-group and S-group

Fig. 3

Heat map of the top fifty different genes in E-group and S-group

GO term enrichment analysis of DEGs

Gene Ontology (GO) showed that up-DEGs were significantly enhanced for ribosomal subunit export from nucleus, ribosome localization, regulation of hemopoiesis, negative regulation of hemopoiesis, and rRNA-containing ribonucleoprotein complex export from nucleus. Downregulated DEGs were enriched for the molting cycle, hair cycle, molting cycle process, hair cycle process, and the skin epidermis development (Table 3 and Fig. 4).

Table 3 Gene ontology analysis of differentially expressed genes in E-group and S-group according to BP, CC and MF

Fig. 4

Gene ontology analysis of differentially expressed genes in E-group and S-group

KEGG pathway analysis of DEGs

The result of KEGG pathway analysis revealed that target genes were enriched in thiamine metabolism, regulation of lipolysis in adipocytes, central carbon metabolism in cancer, estrogen signaling pathway, collecting duct acid secretion, Rap1 signaling pathway, measles, sphingolipid metabolism, drug metabolism-other enzymes, and circadian rhythm. These key pathways were showed in Table 4 and Fig. 5.

Table 4 KEGG pathway enrichment analysis of differentially expressed genes in E-group and S-group

Fig. 5

The results of KEGG pathways enrichment analysis for DEGs based on clusterProfiler

Interaction network of DEGs and core genes in the interaction network

Using data from the Cytoscape and STRING databases, the 10 hub nodes with the greatest degree of network connection were determined. The top 10 hub genes identified were RASL2-9, PSMD6, CPOX, FTH1, PYGL, GNAI1, PTPN1, RIC8A, RAB7A, LOC680441, USP4, and HIST2H2BE (Fig. 6 and Table 3). We listed top three MCODE results in Fig. 7.

Fig. 6

Protein-protein interaction network of differentially expressed genes. Red nodes represent upregulated genes. Blue nodes represent downregulated genes. Green nodes represent the top 10 genes. Nodes > 10 was set as cut-off criteria

Fig. 7

The top 3 modules from the gene–gene interaction network. The squares represent the differentially expressed genes (DEGs) in modules, and the lines show the interaction between the DEGs

Discussion

Subchondral bone remodeling is regulated by the bone resorption and bone formation, which mainly regulated by osteoclast and osteoblast respectively. Despite advancements in the understanding of the mechanism of OA, an effective method for accelerating the process remains to be identified. In this analysis, we found that 588 DEGs between E-group and S-group.

Bone remodeling is regulated by the balanced processes of osteoclast-mediated bone resorption and osteoblast-mediated bone formation [13, 14]. Disequilibrium of this balance leads to dysregulated bone tissue remodeling and can result in excessive bone loss or extra bone formation and consequent skeletal disease [15]. We identified DEGs that play roles in osteoclast and osteoblast differentiation and function in the early stages of OA in this model. Previous studies identified that an altered phenotype of subchondral osteoblasts and osteoclasts contribute to OA progress. Lamas et al. [16] reported that COL10A1 downregulation seems to have a role in the establishment of a defective and/or unstable subchondral cartilage matrix in OA disease. Ren et al. [17] used gene expression profile GSE103416 to identify the different expression genes and potential pathways. Results show that Gna13/cGMP-PKG signaling pathway was identified as a potential research target for therapy and for further understanding the development of OA.

We further performed Kegg pathways to identify the potential pathways that involving in the progress of OA. We found that Rap1 signaling pathway was the most obvious different signaling pathways. For Rap1 signaling pathway, we found nine DEGs (GNAI3/RAPGEF5/CSF1/PIK3CD/MK1/LPAR2/P2RY1/GNAS/PRKD1/RASSF5). Another potential pathway was estrogen signaling pathway. Ren et al. [17] used gene expression profile GSE103416 and found that Gna13/cGMP-PKG signaling pathway was identified as a potential research target for therapy and for further understanding the development of OA. Feng et al. [18] revealed that PDGFRB, IFNG, EGR1, FASLG, and H3F3B may be the potential targets for OA diagnosis and treatment. Liang Y al [19]. found that estrogen deficiency is closely related to the development of menopausal arthritis including OA. Estrogen acts via ER and miR-140 to inhibit the catabolic activity of proteases within the chondrocyte extracellular matrix.

We also found that Rap1 signaling pathway participated into the pathological process of OA. Zhang et al. [20] performed a gene expression analyses of subchondral bone in early experimental osteoarthritis by microarray. Results found that Alp, Igf1, Tgf β1, Postn, Mmp3, Tnfsf11, Acp5, Bmp5, Aspn, and Ihh genes that involved in the pathological process of OA, and they also performed PCR to identify these DEGs in the OA and normal cartilage. Kovács et al. [21] also revealed that the Wnt and the OPG-RANKL-RANK signaling systems, as key mediators, interact in subchondral bone remodeling in OA development, which indicated that subchondral bone remodeling also affects OA. Zhou et al. [22] revealed that matrix metalloproteinase (MMP)1, MMP3, MMP13, and prostaglandin-endoperoxide synthase 2 (PTGS2) are associated with human developmental chondrogenesis.

Conclusion

In summary, totally, 588 DEGs were identified including 199 upregulated DEGs and 389 downregulated DEGs screened in OA and sham-operation. We identified 12 core genes, including RASL2-9, PSMD6, CPOX, FTH1, PYGL, GNAI1, PTPN1, RIC8A, RAB7A, LOC680441, USP4, and HIST2H2BE. Additional experimental studies will be needed to validate our findings.