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

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.

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 plugin 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 cutoff = 0.2, k-core = 2, and max. depth = 100) [12].

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).

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).

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 metabolismother enzymes, and circadian rhythm. These key pathways were showed in Table 4 and Fig. 5.  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.

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 osteoblastmediated 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]   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.