Bioinformatics Analysis of Some Novel Potential Biomarkers Associated With Osteoarthritis and Nervous System

Zhongkui Guo Nankai University School of Medicine Yong Qin Second A liated Hospital of Harbin Medical University Yang Chen Chinese PLA General Hospital Yi Li Chinese PLA General Hospital Ya Gu Chinese PLA General Hospital Ming Chen Chinese PLA General Hospital Pengbin Yin Chinese PLA General Hospital Yanpeng Zhao Chinese PLA General Hospital Licheng Zhang Chinese PLA General Hospital Peifu Tang (  pftang301@163.com ) General Hospital of Chinese PLA https://orcid.org/0000-0003-4279-1704


Introduction
Osteoarthritis (OA), mostly affects the knee joints, is an age-related chronic in ammatory and degenerative changes that usually accompanied by unbearable pain and dysfunction of joints, which brings great puzzles for individuals and the society [1]. Data shows that about 240 million people globally suffered this burdensome syndrome, which has been the 4th leading cause of years disability worldwide in 2020 [2]. Meanwhile, osteoarthritis can affect knees, hands, feet, hips, and spine, but knee osteoarthritis occupied approximately 85% part of all the osteoarthritis in the world [3]. Osteoarthritis is a complex and comprehensive disorder characterized by the degeneration of cartilage chondrocytes, subchondral bone remodeling, and synovium in ammation, and the loss of cartilage chondrocytes has been proved to be the primary pathological change [4].
Besides, a number of studies have illuminated that the nervous system, includes the peripheral nervous system and the central nervous system, also play vital parts in the ethiopathogenesis of osteoarthritis [5,6]. And as we all know, pain, the major reason that drives patients to seek medical advice, is the dominant symptom of osteoarthritis and is an important manifestation of neuromodulation of osteoarthritis as well. Nervous system could not only innervate the nociceptive pain but also innervate synovium, joint capsule and subchondral bone [7,8], and with the progression of osteoarthritis, the higher degree of damaged joints, the more nervous system were engaged in. Therefore, there must be some factors that has the capacity of affecting the nervous system as well as osteoarthritis.
Nowadays, bioinformatics technology has been a novel method to extract genetic information from many diseases.
In this work, we analysis three gene expression pro les, GSE114007, GSE51588, and GSE55457, that downloaded from the Gene Expression Omnibus database (GEO), by bioinformatics technology, and nally, three genes, HES1, JUN, and IER2, that may help to diagnosis and cure osteoarthritis in the future were identi ed, which were preliminary con rmed by using qRT-PCR.

Gene Expression Database
Gene expression pro les GSE114007, GSE51588, and GSE55457 were downloaded from the Gene Expression Omnibus (GEO, www.ncbi.nlm.nih.gov/geo/). Dataset GSE114007 includes 20 osteoarthritis and 18 normal mRNAs expression data of cartilage tissues. Dataset GSE51588 includes 20 osteoarthritis and 5 normal mRNA expression data of cartilage subchondral bone. Dataset GSE55457 includes 26 osteoarthritis and 20 normal mRNA expression data of synovial membrane.

Identi cation of Differentially Expressed Genes (DEGs)
The Limma package was used in R software to identify the DEGs. The corresponding P value of the gene symbols after t test was used, and genes with adjusted P < 0.05 and |logFC|>1.5 were de ned as DEGs. Then the pheatmap package was used in R software to draw a heatmap and a volcano plot according to the DEGs.

Protein-protein Interaction (PPI)
Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org) was used to illustrate the potential interactive relationships between DEGs. And then the Cytoscape software was used to draw the PPI network.

Key Genes Identi cation
To identify the key genes, rstly, biological process (BP) part of the GO analysis results was used to screen the neuro-related DEGs of the osteoarthritis. Then, the neuro-related DEGs list, together with the DEGs lists of GSE51588 and GSE55457 were submitted to Venn Diagrams (bioinformatics.psb.ugent.be/webtools/Venn/). Finally, genes in the intersection were regarded as our key genes.

Quantitative Real-time PCR (qRT-PCR)
Human chondrocyte cell line C28/I2 was chosen to con rm the results of bioinformatics analysis. Total RNA from normal group (normal C28/I2 human chondrocyte cells) and osteoarthritis group (IL-1beta induced C28/I2 human chondrocyte cell) was isolated using RNA Extraction Kit (Solarbio, Beijing, China). Then, a Reverse Transcriptase Kit (Solarbio, Beijing, China) was used to generate cDNA. qRT-PCR was conducted on the Bio-Rad Real-

Study Design
The study design was schematically depicted in Fig. 1. Dataset GSE114007 was used for DEGs screening (Step 1).
Then, the neuro-related genes of dataset GSE114007 together with the DEGs of dataset GSE51588 and GSE55457 were submitted to the Venn diagram to obtain the intersection to con rm the key genes (Step 5). Finally, qRT-PCR was used to further identify the accuracy of the key genes (Step 6).

Identi cation the DEGs in GSE114007
A total of 878 DEGs were identi ed with dataset GSE114007 (P < 0.05 and |logFC|>1.5), consisting of 495 upregulated genes and 383 down-regulated genes between the osteoarthritis and normal cartilage tissues. The top 10 up-regulated genes and the top 10 down-regulated genes were listed in Table 1 and Table 2 respectively. Heatmap of the DEGs were obtained in Fig. 2A, and Volcano plot of the DEGs were obtained in Fig. 2B.

GO function and KEGG pathway Enrichment Analysis of the DEGs in GSE114007
Gene ontology (GO) enrichment analysis, includes biological process (BP), cellular component (CC) and molecular function (MF), was used to explore the potential biological function of the 878 DEGs in dataset GSE114007. The top 5 functions for BP were extracellular matrix organization, collagen catabolic process, angiogenisis, collagen bril organization, and cell adhesion (Fig. 3A), the top 5 functions for CC were proteinaceous extracellular matrix, extracellular region, extracellular matrix, extracellular space, and plasma membrane (Fig. 3B), and the top 5 functions for MF were extracellular matrix structural constituent, transcription factor activity, RNA polymerase II core promoter proximal region sequence-speci c binding, transcriptional activator activity, RNA polymerase II core promoter proximal region sequence-speci c binding, scavenger receptor activity, and metalloendopeptidase activity ( Fig. 3C). KEGG pathway analysis revealed the mainly related pathway of the 878 DEGs in dataset GSE114007, which mainly were related with focal adhesion, PI3K-Akt signaling pathway, ECM-receptor interaction, staphylococcus aureus infection, protein digestion and absorption, and osteoclast differentiation (Fig. 3D).

PPI Network Construction of the DEGs in GSE114007 and Hub Gene Identi cation
In order to better understanding of the 878 DEGs in dataset GSE114007, STRING was used to analysis the potential interactive relationship between DEGs. Then, Cytoscape software was used to draw the PPI network (Fig. 4A).

Neuro-related Genes Identi cation
Previous studies have demonstrated that OA was correlated with the nervous system [2,4], therefore, in order to obtain the neuro-related genes that have the capacity of regulating osteoarthritis, the biological process (BP) functions of GO enrichment analysis results was used for further identi cation. The result showed that the neuro- . A total of 76 neuro-related differential expression genes (DEGs) were obtained as neuro-related genes with p value < 0.05 (Table 3).

Key Genes Identi cation
In order to detect the genes that associated with nervous system and osteoarthritis, dataset GSE99662 and GSE51588 were used for further veri cation. Dataset GSE51588 contains the differentially expressed genes of subchondral bones in osteoarthritis (obtains 1010 DEGs) and dataset GSE55457 contains the differentially expressed genes of synovial tissues in osteoarthritis (obtains 704 DEGs). To obtain the most relevant neuro-related genes, the 76 neuro-related DEGs together with the DEGs of dataset GSE51588 and GSE55457 were submitted to Venn diagram to get the intersection, and the genes in the intersection were regarded as the key genes (Fig. 5).
According to the Venn diagram, 3 key genes (HES1, JUN, and IER2) were identi ed at last, which may be associated with the nervous system as well as the pathophysiological process of OA.

qRT-PCR for Key Genes Validation
In order to validate the different expression levels of HES1, JUN, and IER2 in OA and normal chondrocytes, we chose the il-1 beta induced human chondrocytes cell line C28/I2 as OA group and the non-treated human chondrocytes cell line C28/I2 as normal group. By using qRT-PCR, we found that the expression of HES1, JUN, and IER2 genes were signi cantly changed in OA chondrocytes compared to normal condrocytes (Fig. 6), which were in full accord with our detection.

Discussion
Osteoarthritis is a leading cause of disability in the elderly with heavy burden on individuals, families and society. Today, for most non-terminally symptomatic patients, turning to medication for pain relief, but not a cure, is the main method for osteoarthritis treatment, and for terminal joint osteoarthritis patients, the only effective cure is the total joint replacement [2,4]. Therefore, due to the unclear mechanism induced treatment resistent and high burden of osteoarthritis, illuminate the relationship between the nervous system and osteoarthritis may provide us a new method to solve the puzzle that osteoarthritis brings us.
Evidence from animal to clinical models have demonstrated that osteoarthritis-like joint damage is not only limited to the destruction of cartilage chondrocytes, synovial tissues, subchondral bone, but also associated with the disorder of the peripheral and central nervous system [9,10].
It's well known that joint pain evoking is mainly correlated with the sensitization of peripheral nociceptive neurons and hyperexcitability of nociceptive neurons in central nervous system, which is closely associated with the synovial in ammation. Lei et al. [11] suggested that SDIM1, whose overexpression exhibited a role of improving survival of neuro-progenitor cells after injury, showed an signi cantly increased expression in synovial tissues of osteoarthritis patients with high pain. Bullock et al. [12] revealed CGRP release might play an important role in the peripheral sensitization during joint degeneration in osteoarthritis. Brutus et al. [13] summarized four known functional genetic variants (SCN9A, COMT, TRPV1 and P2X7) in pain candidate genes that in uence pain experiences in patients with osteoarthritis.
In this paper, bioinformatics technology was used to identify the DEGs of osteoarthritis, and at last, HES1, JUN, and IER2 were identi ed as the key genes that may take vital roles in the regulation of nervous system as well as osteoarthritis.
HES1, also known as Hes family bHLH transcription factor 1, is a downstream effector of Notch signaling. Activated Notch signaling can induce the activation of HES1 and hence inducing the expression of MMP13, which has the potential to promote the degradation of chondrogenic ECM, and resulted in the degeneration of cartilage [14]. Ni et al. [15] revealed that OSM is up-regulated in the synovial tissue of knee osteoarthritis and OSM-treated MC3T3-E1 cells showed a down-regulated HES1 in a time-dependent manner. Sugita et al. [16] reported that Hes1 induced Adamts5 and Mmp13, which are catabolic enzymes that break down cartilage matrix. Additionally, Matsuzaki, et al. [17] identi ed that Hes1 expression in mature neurons in the adult mouse brain is required for normal behaviors. Harris,et al. [18] demonstrated that HES1 has the potential of promoting the quiescence and proliferation of adult neural stem cells.
JUN, one of the transcription factor for activator protein 1, can mediate catabolic transcription, cell apoptosis and cell death. Chen et al. [19] reported that JUN was down-regulated in osteoarthritis knee cartilage compared to normal knee cartilage. Cai, et al. [20] demonstrated that osteoarthritis synovial tissues exhibit a decreased expression of JUN. While, studies also showed that blocking the JUN could prevent the degradation and apoptosis of chondrocytes [21], and the inhibition of JUN transcriptional activity protects against osteoarthritis cartilage destruction [22]. Moreover, studies also showed that JUN could participate in the regulation of antioxidant responses in neurons [23] and may have the potential of regulating the peripheral myelin development [24] and neuronal polarization during brain development [25].
IER2, short for immediate early response 2, is a protein function as a potential transcriptional factor or transcriptional coactivator which seems to play a pivotal role in regulating tumor cells. Xu et al. [26] reported that IER2 promotes the migration and invasion of hepatocellular carcinoma cell adhesion and motility. Xu et al. [27] demonstrated that IER2 promotes the migration and invasion of hepatocellular carcinoma cells via regulating the activity of Rho GTPases. Moreover, Moriya et al. [28] reported that ier2 mRNA was distributed in the telencephalon, midbrain and the hypothalamus. Besides, previous study [20] also reported that IER2 was down-regulated in osteoarthritis knee cartilage compared to normal knee cartilage. While, there is no research about IER2 on the speci c mechanism of osteoarthritis and nervous system by now.
In this study, we demonstrated that HES1, JUN, and IER2 may take vital parts in the regulation of nervous system and osteoarthritis, and the result was validated by using qRT-PCR nally. Meanwhile, we also found that there are still somewhat controversial between the previous studies. Considering of this, there must be still have some unknown mechanism between these three genes, osteoarthritis as well as nervous system. So, it's necessary to undertake more research to con rm our speculation.

Conclusion
Bioinformatics analysis has been a novel tool to nd new targets for diagnosis and cure of disease. In this paper, we screened out three genes that may offer new targets to understand osteoarthritis. While, at the same time, more studies are required to further con rm our results. Schematic illustration of study design. A total of 878 DEGs were identi ed with dataset GSE114007 (Step 1), and then Gene ontology (GO) enrichment analysis was performed (Step 2). All the 878 DEGs were tested according to their potential as neuro-related genes with the biological process (BP) result of GO analysis, and 76 genes were screened out (Step 3). GSE51588 (includes 1010 DEGs) and GSE55457 (includes 704 DEGs) were chosen for further veri cation (Step 4). The Venn diagram was performed and 3 genes were identi ed as our key genes (Step 5). Finally, qRT-PCR was used to identify the accuracy of the key genes (Step 6).  Venn diagram of DEGs from three datasets for key genes identi cation. The genes in the intersection were regarded as our key genes, and at last, three genes were identi ed as our key genes.