Osteoarthritis (OA) is the most common joint disease, mainly manifesting as pain, limited joint movement, and joint deformity. The risk factors of OA include trauma, aging, obesity, and heredity (Sandell, 2012; Chen et al., 2017a; Chang et al., 2021). During OA development, the entire joint are affected and undergo articular cartilage degeneration, osteophyte formation, subchondral sclerosis, synovitis, and meniscus degeneration, respectively, indicating the complicated and interactive OA pathogenic mechanisms (Chen et al., 2017a). The efficacies of the current treatments for OA in our clinics are limited. In recent years, exploration of disease-modifying osteoarthritis drugs (DMOADs) aiming at alleviating OA symptoms and/or prevent structural progression have drawn much attention. However, the DMOADs under research and development (R&D) and/or clinical trials mainly focus on one of the OA symptoms, such as cartilage degeneration, subchondral bone remodeling, local inflammation, or joint pain, and their potential downstream targets (Latourte et al., 2020). The possibility that newly explored-drug targets have heterogeneous expression profiles in different joint components raises uncertainty of the drug effectiveness. Besides, the R&D of joint component-specific drugs are also limited so far (Latourte et al., 2020).

Based on the rapid development of high-throughput genomics technologies, such as microarray and next-generation sequencing, bioinformatic analysis has been widely used to identify key genes, signaling pathways, and microRNAs (miRNAs) in various musculoskeletal disorders (Kang et al., 2021; Li et al., 2021; Umeno et al., 2021). So far, many studies have identified and explored potential therapeutic targets of OA based on bioinformatic screening. For example, upregulated arginase 2 (ARG2) in OA cartilage was screened out by microarray and further validated to facilitate cartilage destruction via upregulating matrix metalloproteinases (MMPs) (Choi et al., 2019). Activated osteochondral turnover, neurogenesis and inflammation in OA bone marrow lesions (BML) were also identified by microarray bioinformatically (Kuttapitiya et al., 2017). Besides, miRNA candidates that have potential as biomarkers and therapeutic targets in OA were identified and validated via comprehensively paired miRNA-messenger RNA (mRNA) analysis and functional enrichment analysis (Kung et al., 2018). In addition to identification of potential therapeutic targets, bioinformatics-based bulk sample analysis further helps classify potential OA subtypes for more precise diagnosis and personalized treatments. In recently years, several studies have stepped forward substantially in classifying potential OA subtypes based on bioinformatics (Soul et al., 2018; Yuan et al., 2020). Their surprising discoveries undoubtedly deepen our understandings on knee OA and will facilitate personalized treatments in the future. However, since previous bioinformatic studies mainly focus on one type of joint components as well, investigations on the overlapping DEGs (overlapping DEGs) in different joint components during OA development are still lacking. In 2016, about 5% overlapping DEGs were observed between the DEGs of the synovium and cartilage, while no further analysis was performed on these identified overlapping DEGs (Park and Ji, 2016). Recently, another study observed about 10% overlapping DEGs in OA cartilage and subchondral bone. They identified IL11 and CHADL as two potential therapeutic targets of OA by comparing their identified cartilage-subchondral bone overlapping DEGs with previously identified OA risk genes (Styrkarsdottir et al., 2018; Tuerlings et al., 2021).

Collectively, we believe it is meaningful to comprehensively analyze the key genes that are differentially expressed during OA development in the different joint components, including articular cartilage, subchondral bone, synovium, and meniscus, and their relevant pathways and miRNAs. Such approaches may provide clues to develop adequate treatments for OA by targeting at overlapping differentially expressed genes (DEGs) in different joint components and their relevant miRNAs and signaling pathways. The present study aims at identifying key genes, signaling pathways, and miRNAs in human knee OA by comparing the preexisting gene expression profiles derived from different joint components, including articular cartilage, synovium, subchondral bone, and meniscus. Specifically, the gene expression profiles (GSE) were obtained from the public available Gene Expression Omnibus database (GEO, http://www.ncbi.nlm.nih.gov/geo/). Gene Expression Omnibus 2R (GEO2R) was performed to identify the overlapping DEGs and followed by qRT-PCR validation. Furthermore, functional enrichment analysis, protein-protein interaction (PPI) analysis, and miRNA-mRNA interaction analysis were carried out to identify relevant signaling pathways and interactive miRNAs. This study may shed light on completer and undiscovered pathogenic mechanisms of knee OA development and pave the way toward the identification of new therapeutic targets for further R&D of effective therapies and clinical translation.

The present study aims to screen out key DEGs, their relevant signaling pathways, and interactive miRNAs in human knee OA based on bioinformatic analysis. The gene expression profiles derived from different joint tissues are covered and followed by identifying overlapping DEGs. To the best of our knowledge, this is the first time that overlapping DEGs in knee OA are identified from four different OA joint tissues, including articular cartilage, synovial membrane, subchondral bone, and meniscus. Gene expression profiles are all obtained from GEO database.

After identification of DEGs in different joint tissues, overlapping DEGs are then identified via Venn Diagram. As a result, 236 overlapping DEGs are identified, of which 160 are upregulated and 76 are downregulated. Those DEGs that are differentially expressed in at least three joint tissues are identified as the most overlapping DEGs. As a result, a total of 13 most overlapping DEGs are identified in our study, include AHR, HEY1, CXCL12, MMP9, COL8A1, GAP43, PTN, RUNX1, PIM1, OLFML2A, SLITRK6, RHBDL2, and ARL4C, while no overlapping genes are differentially expressed in all four components. Among above identified most overlapping DEGs, OLFML2A, SLITRK6, RHBDL2, and ARL4C are excluded from our constructed PPI network, suggesting their relatively limited values and will not be discussed in the following context.

MMP9 is an enzyme implicated in the cartilage destruction and has been identified as a hallmark of OA previously, the same as other MMPs like MMP1, MMP3 and MMP13. Previous studies have reported the upregulation of MMP9 in OA cartilage, synovial membrane, subchondral bone, and synovial fluid (Yang et al., 2010; Kim et al., 2014). COL8A1, which maintains cartilage stability through participating in collagen synthesis, is also reported to be highly expressed in OA cartilage (Fang et al., 2019). CXCL12, also known as SDF-1, is a chemokine that plays important role in angiogenesis, bone metabolism, cartilage homeostasis, and pro-inflammatory responses (De Klerck et al., 2005; Garcia-Cuesta et al., 2019). Upregulation of CXCL12 in OA cartilage, subchondral bone, synoviocytes, and synovial fluid has been reported by previous studies (Chen et al., 2017b; Bragg et al., 2019). Furthermore, inhibition of CXCL12/CXCR4 signaling was shown to prevent subchondral bone loss and significantly attenuate cartilage degradation (Dong et al., 2016; Chen et al., 2017b). GAP43 is a nervous tissue-specific cytoplasmic protein related to nerve regeneration. Previous study showed that the expression of GAP43 in pain-related sensory innervation of dorsal-root ganglia (DRG) was upregulated during OA progression (Kawarai et al., 2018). Thus, the upregulation of GAP43 in joint tissues may reflect joint sensory innervation, which is closed related to nociceptive sense and osteophyte formation (Wu et al., 2002; Orita et al., 2011). Pleiotrophin (PTN) is an 18-kDa heparin-binding neurite outgrowth-promoting growth factor. Previous studies have reported the upregulation of PTN in OA cartilage, synovial membrane, and synovial fluid. Notably, PTN is initially abundant in fetal or juvenile cartilage and then becomes absent in mature cartilage. During early stages of OA, PTN becomes re-expressed again (Pufe et al., 2003; Mentlein 2007). Besides, PTN is also proven to facilitate chondrocyte proliferation (Pufe et al., 2007). Interestingly, GAP43 and PTN are two nerve growth-related genes, suggesting the important role of nerve innervation and axon guidance during OA development in the entire joint. Their specific role in OA development needs to be further investigated. PIM1 is an enzyme that play important roles in cell cycle progression, apoptosis, and transcriptional activation (Bachmann and Moroy, 2005). However, there is no study reporting their relationships so far.

The rest of most overlapping DEGs, including HEY1, AHR, RUNX1, and HEY1, are all transcription factors. HEY1 is a basic helix-loop-helix protein (bHLH) transcription factor that belongs to HES/HEY family. Besides, HEY1 is also a direct target of canonical Notch signaling. Previous studies indicated that inhibition of Notch1 exacerbated experimental OA, while increased levels or activity of Notch2 contributed to the progression of OA (Hosaka et al., 2013; Lin et al., 2016). In addition, HEY1 is not only transcriptionally induced by Notch ligands, but also induced by BMP/TGF-β axis to exert as a transcription repressor. Its functions on repressing osteogenic differentiation, neuronal differentiation, and pro-inflammatory production have been reported before (Weber et al., 2014). Hes1, another transcription factor that belongs to HES/HEY family like HEY1, has been reported to be upregulated in OA cartilage and accelerate cartilage destruction via promoting MMP3, a disintegrin and metalloproteinase with thrombospondin motifs 5 (ADAMTS5), and interleukin-6 (IL-6) transcription (Sugita et al., 2015). Nonetheless, the specific role of HEY1 in OA progression remain to be investigated. AHR is a bHLH transcription factor that plays important role in development, immune system, and toxic response. Previous studies have demonstrated that AHR signaling activation significantly alleviated progression of rheumatoid arthritis (RA) through repressing C-reactive protein (CRP), NLR family pyrin domain containing 3 (NLRP3), tumor necrosis factor-alpha (TNF-α), and IL-6 expression (Jin et al., 2011; Huai et al., 2014; Liang et al., 2019; Piper et al., 2019), and enhancing nuclear factor erythroid 2-related factor 2 (NRF2) and IL-10 expression in B cells, macrophages, or hepatocytes (Tsuji et al., 2012; Piper et al., 2019; Rosser et al., 2020). However, several reports also indicated that AHR signaling activation exacerbated RA inflammation through activating cytokine-mediated inflammatory signaling in primary human fibroblast-like synoviocytes (Adachi et al., 2013; Lahoti et al., 2013). Ogando J et al. reported that the AHR signaling pathway was significantly more active in OA synovial tissues than in RA synovial tissues (Ogando et al., 2016). Furthermore, compared to resting chondrocytes, significantly upregulated AHR was also observed in hypertrophic chondrocytes (Cedervall et al., 2015). Nevertheless, the specific effects of AHR on OA are still in debate and remain to be investigated. RUNX1 was reported to be significantly downregulated in OA cartilage compared with normal control. Furthermore, its anabolic effect on chondrocytes contributes to the maintenance of cartilage homeostasis during OA development and that has also been recognized as a disease-modifying target of OA (Yano et al., 2013; Aini et al., 2016; Yano et al., 2019). In addition, the anti-angiogenic effect of RUNX1 through repressing vascular endothelial growth factor (VEGF) expression has also been reported (Ter Elst et al., 2011).

According to KEGG pathway analysis based on overlapping DEGs, we observe that overlapping DEGs are enriched in “PI3K-Akt signaling pathway,” “Wnt signaling pathway,” “Fluid shear stress,” “Nicotinate and nicotinamide metabolism,” “HIF-1 signaling pathway,” “MAPK signaling pathway,” “Cytokine-cytokine receptor interaction,” “PPAR signaling pathway,” “NOD-like receptor signaling pathway,” “Axon guidance,” and “TGF-β signaling pathway.” The above pathways are well consistent with existing research findings. For example, PI3K-Akt signaling and MAPK signaling are closely involved in inflammatory response to induce the production of catabolic markers such as MMPs, Adamts, IL-1β, and TNF-α in OA (Herrero-Beaumont et al., 2019; Chow and Chin, 2020). Nicotinate and nicotinamide metabolism play important roles in redox reaction. So far, several studies have reported the important role of nicotinate and nicotinamide metabolism in OA development (Yang et al., 2015; Junker et al., 2017). In 2015, Yang et al. (2015) demonstrated that nicotinamide phosphoribosyltransferase (NAMPT), a rate-limiting enzyme in the Nicotinamide adenine dinucleotide (NAD+) salvage pathway, acted as a crucial catabolic regulator of osteoarthritic cartilage destruction. In addition, aberrantly activated Wnt signaling or TGF-β signaling contributes to cartilage degradation, osteophyte formation, and formation of subchondral bone marrow osteoid islets (Zhen et al., 2013; van der Kraan 2017). Therapies targeting Wnt signaling have been undergoing clinical trials (ClinicalTrials.gov Identifier: NCT03928184). HIF-1α and HIF-2α were also reported to exert anabolic and catabolic effects on chondrocytes, respectively (Zhang et al., 2015). PPAR signaling, in particular PPARγ, which is significantly downregulated in OA cartilage, has been demonstrated to maintain articular cartilage homeostasis via regulating the mTOR pathway (Vasheghani et al., 2015; Zhu et al., 2019). NOD-like receptor signaling, which mediates innate immunity and participates in regulating inflammatory and apoptotic responses, mainly includes two subfamilies, NODs and NLRPs (Platnich and Muruve, 2019). Both NOD-dependent pathway and NLRP-dependent inflammasome pathway were validated to mediate OA development under external stimulus (Jin et al., 2011; Xu et al., 2015). Fluid shear stress (FSS) was known as one of the pathogenic mechanisms of OA and has also been used in the construction of an osteoarthritic cell model for a long time (Yang et al., 2020). In addition to above well-validated OA related-pathways, axon guidance pathway is lack of research up to now. Only several in vitro studies reported the effects of an axon guidance molecule, Semaphorin 3A (Sema3A), on osteoarthritic chondrocytes (Okubo et al., 2011; Sumi et al., 2018).

Through constructing the PPI network and miRNA-mRNA network in Cytoscape, a total of 25 hub genes, four significant modules, and 57 key miRNAs are identified. Among hub genes, MAPK14, IL2RB, and IL6R are all involved in cytokine-mediated inflammatory response during OA (Boileau et al., 2005, Liang et al., 2018, Shkhyan et al., 2018). A genome-wide association (GWAS) study has identified a single nucleotide polymorphism (SNP) of HLA-DQB1, rs7775228, associating with knee OA in Asian population (Valdes and Spector, 2011). Furthermore, a small molecule gp130 modulator (RCGD 423) was proven to improve chondrocyte proliferation and inhibit cartilage degradation via upregulating transcription factor MYC and suppressing IL-6-mediated inflammatory response (Shkhyan et al., 2018). For growth factor VEGFA, it is closely related to angiogenesis and inflammatory response (Gao et al., 2013). To validate above computational results, we performed qRT-PCR to validate the relative expression of identified key genes in OA cartilage clinically. Our results show that AHR and its downstream target CYP1A1, HEY1, MYC, and CXCL12 are indeed upregulated in degenerated cartilage compared to preserved cartilage (Figure 6).

According to module analysis, majority of hub genes are located in module one and module 2. In addition, MMP9, CXCL12, and HEY1 are not only the most overlapping DEGs, but also hub genes. Module analysis shows that Modules 1, 2, 3, and 5 are significant modules with the MCODE score ≥4. For Module 1, all of the genes within Module 1, except ANGPT2, are hub genes. MMP9 is located in Module 1. Enrichment analysis suggests that Module one plays important role in PI3K-Akt signaling pathway and cell adhesion molecules (CAMs) pathway. Regarding genes within Module 2, all of them are hub genes except ACP5 and FUT4. Genes in Module two play important roles in rheumatoid arthritis and cytokine-cytokine receptor interaction pathways. Furthermore, CXCL12 is located in Module 2. Module three is a group of genes associated with Wnt signaling pathway. No most overlapping DEGs or hub genes locate in Module 3. Module five includes COL8A1, COL13A1, and COL22A1, all of which are associated with collagen trimer pathway and have been demonstrated to be aberrantly expressed in degraded cartilage (Karlsson et al., 2010; Feng et al., 2019). A recent study reported that Lgr5+/Col22a1+ stem cells play important roles in differentiation toward articular chondrocytes and Col22a1-expressing cartilage superficial layer contributed to repair of cartilage defect (Feng et al., 2019).

Collectively, based on above identified overlapping DEGs, hub genes, pathways, and functional modules, the present study depicted the potential OA mechanisms covering the entire knee joint. We believe that our identified knee OA-related key genes, such as AHR, HEY1, MYC, GAP43, and PTN, and their relevant signaling pathways, including AHR signaling, Notch signaling and TGF-β signaling, C-MYC signaling, and axon guidance pathway may play important roles in knee OA development. Our predicted genes are worthy of being explored as novel targets of DMOADs in the future.

By comparing key miRNAs with OA-related miRNAs included in HMDD v3.2 database, 41 out of 57 miRNAs (about 71.9%) have been reported to be associated with OA, suggesting the considerable reliability of our miRNA-mRNA interactome predictions. For example, miR-29 family, including miR-29a, miR-29b, and miR-29c, was differentially expressed in OA cartilage and negatively regulated Smad, NF-κB, and canonical Wnt signaling (Le et al., 2016). In addition, upregulation of miR-101 significantly facilitated cartilage degradation and chondrocyte apoptosis (Dai et al., 2015; Lü et al., 2020). For miR-181a, Akihiro Nakamura et al. demonstrated that intra-articular injections of locked nucleic acid (LNA) miR-181a antisense oligonucleotides (ASO) significantly attenuated cartilage destruction in facet and knee joints in vivo (Nakamura et al., 2019). Another study reported that miR-9 promoted IL-6 expression and exacerbated cartilage degradation by targeting MCPIP1 expression (Makki et al., 2015). In 2019, a hydrogel-based drug delivery system equipped with locked nucleic acid (LNA) miR-221 inhibitor was constructed and enhanced cartilage regeneration significantly (Lolli et al., 2019). miR-199a* was also shown to inhibit IL-1β-induced Cyclooxygenase 2 (COX-2) expression as a direct regulator (Akhtar and Haqqi, 2012). Notably, a RNA sequencing-based miRNA-mRNA interactome study which confirmed a OA-specific miRNAs array showed considerable overlap with our identified key miRNAs as well, including miR-143, miR-155, let-7g, miR-7, miR-15, miR-101, miR-21, miR-19a, miR-16, miR-30a, miR-29, miR-1, miR-133a-3p, miR-20a, miR-320a, miR-31, miR-93, miR-221, and miR-335 (Coutinho de Almeida et al., 2019). Among our identified top 10 key miRNAs, in addition to those that have been reported before, miR-17, miR-21, miR-155, miR-185, and miR-1 are still lack of research and remain to be investigated in the future.

The present study also has several limitations. For example, although we originally intended to include as many data sets as possible, only four expression profiles met the criteria and were included into the data analysis. In addition, due to difficulties in obtaining appropriate specimens, only OA articular cartilage specimens were used to carry out clinical validation experiments. The relative expression of our identified knee OA-related key genes in other joint components need to be further validated.

In conclusion, the present study provides a comprehensive bioinformatics analysis of key genes, signaling pathways, and miRNAs in different joint tissues of knee OA patients. A total of 35 knee OA-related key genes and 57 key miRNAs were identified. Among them, key genes such as AHR, HEY1, MYC, GAP43, and PTN, and key miRNAs such as miR-17, miR-21, miR-155, miR-185, and miR-1 are lack of research so far. For key genes identified in the present study, their downstream mechanisms and specific effects on the different joint components also need to be explored, respectively. Through enrichment analysis, a number of OA-related pathways were identified, including PI3K-Akt signaling, Wnt signaling, fluid shear stress, nicotinate and nicotinamide metabolism, HIF-1 signaling, MAPK signaling, cytokine-cytokine receptor interaction, PPAR signaling, NOD-like receptor signaling, TGF-β signaling, and axon guidance pathways. Among them, many pathways have been well investigated and even under clinical trials, except for axon guidance pathway which is implicated in nerve innervation and axon guidance while still lack of research so far. Our study provides insight for the first time in identification of potential therapeutic targets of knee OA by comprehensively analyzing the overlapping genes differentially expressed in multiple joint components based on bioinformatics.