Introduction

Osteoarthritis (OA) is a prevalent, progressive, and degenerative joint disease characterized by reduced chondrocytes and degradation of the extracellular matrix (ECM) that causes joint pain, stiffness and functional limitations, severely impacting patients’ quality of life by leading to impaired sleep, depression, and even disability.1,2 The incidence and prevalence of OA are increasing annually due to the aging population and rising obesity rates. It currently affects approximately 530 million people worldwide and is projected to pose a significant public health challenge in the years to come.3 The lack of early prediction and targeted treatment for OA are due to its complex and diverse pathological changes, coupled with an incomplete understanding of its pathogenesis. Thus, understanding the pathogenesis of OA is essential for enhancing clinical diagnosis, alleviating symptoms, and improving prognosis.

Bioinformatics analyses have emerged as efficient and promising tools for screening significantly aberrantly expressed genes and genetic pathways. With the rapid development of high-throughput sequencing technology, integrated analysis techniques based on biological networks have been developed to offer new insights into disease processes and to identify biomarkers.4,5 As research progresses, it has been discovered that numerous specific genes can participate in the development of OA.6,7 Recent studies have shown that PGAM5 is significantly associated with an increase in OA incidence compared to controls, and specific inhibition of PGAM5 could substantially attenuate OA symptoms.8 Meanwhile, researchers have identified promising OA biomarkers such as C-terminal telopeptide of collagen type II (CTX-II)9 and type I collagen cross-linked N-telopeptide (NTX),10 but their discrimination capacity remains insufficient for clinical applications. These findings suggest that certain functional genes play an important role in OA progression, and it is worth noting that weighted gene co-expression network analysis (WGCNA) is a widely used systems biology approach for elucidating the correlation between gene clusters (modules) with similar expression patterns and disease phenotypes in transcriptome sequencing data or microarray samples.11 WGCNA can be utilized to screen for disease molecular markers or potential therapeutic targets in various conditions such as atherosclerosis,12 sepsis13 and rheumatoid arthritis.14 However, there is a lack of research on studying pivotal genes for OA using WGCNA, indicating the need for further in-depth investigation into its biological mechanisms through bioinformatics methods.

In this study, we aimed to identify novel differentially expressed genes (DEGs) in OA. Two GEO datasets (GSE169077 and GSE114007) were analyzed to determine DEGs between OA and control specimens, based on which a weighted gene expression network to identify modules was constructed. A protein-protein interaction (PPI) network was also created to screen for hub genes associated with OA, which were validated using RT-qPCR on rat OA models. Collectively, results could contribute to identifying novel key genes and pathways involved in OA progression, as well as providing potential therapeutic targets.

Materials and Methods Data Collection and Preprocessing

We searched the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo) database using the keyword “Osteoarthritis (OA)”. Raw data from the GSE169077 and GSE114007 datasets were downloaded. The GSE169077 dataset comprised 5 normal and 6 OA pool samples, with each pool containing 5 individual samples. The GSE114007 dataset included 18 normal and 18 OA samples. Data processing was performed using the R software.

Differentially Expressed Genes (DEGs) Analysis

Differentially expressed genes (DEGs) was performed using limma software package of R software.15 Genes with adjusted false discovery rate (FDR) < 0.05 and absolute log2 (fold change) > 1 were considered significantly differentially expressed between the control and OA groups.

Gene Set Enrichment Analysis (GSEA)

GSEA was utilized to identify the associated signaling pathways in the progression of OA. Significant enrichment was determined based on a p-value < 0.05 and an FDR (q-value) < 0.05. The results were visualized using the normalized enrichment score (NES).

Weighted Gene Co-Expression Network Analysis (WGCNA) and Identification of Significant Modules

WGCNA was performed on the merged dataset to elucidate gene interactions and identify co-expression patterns of genes and modules.16 Briefly, the samples were clustered to detect any significant outliers. Then, a scale-free network was constructed by assessing the strength of connections between genes. The “PickSoftThreshold” function, with the networkType parameter set as “unsigned”, was employed to determine the value of β (a soft threshold power parameter). Co-expression modules were then formed based on genes with similar expression patterns. Module identification was achieved using a dynamic tree-cutting algorithm, with the minimum module size set to 30. Lastly, modules with high similarity scores were merged using a threshold value for each dataset, and module affiliation (MM) and gene significance (GS) were calculated for each module and utilized for hub gene selection. Potential pivotal genes associated with OA were assessed using criteria of MM > 0.8 and GS > 0.5. Hub genes were identified as those present in significant modules across both datasets.

PPI Network Construction and Functional Enrichment Analysis

To further analyze the biological functions of the genes, the search tool for retrieval of interacting genes/proteins (STRING) database was utilized to construct a PPI network, which was visualized using Cytoscape software. The hub genes were analyzed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis using the Metascape database. Enrichment analysis of GO terms and KEGG pathways was conducted using the “ClusterProfiler” package in R software to predict potential functions.

Animals and Treatments

Pathogen-free adult male Sprague–Dawley rats (8 weeks old, weighing 280–300 g) were obtained from the Animal Center of Chongqing Medical University, housed under standard conditions in a temperature-controlled room (22 ± 3 °C) with a 12-hour light/dark cycle, and were provided food and water ad libitum. All animal procedures to be employed in the project was approved by Institutional Animal Care and Use of Chongqing Medical University (IACUC-CQMU), Approval number: IACUC-COMU-2024-0473. A rat model of OA was induced by intra-articular injection of 4% papain (Biosharp, China) under 1% sodium pentobarbital anesthesia (5mL/kg). The surgical area was shaved, and the right knee was sterilized. Then, 0.2 mL of 4% papain solution was injected into the joint cavity, the control group was injected with an equal amount of saline. This procedure was repeated three times on days 1, 4, and 7, three days after the end of modeling, the rats were sacrificed, the synovium and cartilage of their right knee joints were removed for subsequent experiments. All experiments were conducted following the approved guidelines of the Animal Care and Use Committee of Chongqing Medical University, strictly adhering to all ethical principles.

Quantitative Real-Time PCR (RT-qPCR) Analysis

Total mRNA from tissue samples (Mixed proteins of synovial tissue and joint synovial fluid) (n = 6) was extracted using TRIzol (Thermo Fisher Research, USA) following the manufacturer’s standard protocol. After determining the concentration of the extracted RNA, cDNA synthesis was performed using 1μg total RNA with the HiScript® II Q Select RT SuperMix for qPCR kit (VAZYME, China) according to the manufacturer’s instructions. RT-qPCR was conducted using the CFX96 Real-Time PCR Detection System (Bio-Rad) with ChamQ SYBR qPCR Master Mix (VAZYME). The primer sequences are provided in Table S1. The expression of the target genes was normalized to the levels of GAPDH. The experiments were performed in triplicate.

Western Blotting

Total protein was extracted from rat tissues (Mixed proteins of synovial tissue and joint synovial fluid) using RIPA with phosphatase and protease inhibitors. The protein concentration was quantified using a BCA Protein Assay Kit (Beyotime, China). The extracted proteins were separated using SDS-polyacrylamide gels electrophoresed (SDS-PAGE) and transferred onto polyvinylidene difluoride (PVDF) membranes. After blocking with 5% nonfat dried milk diluted in tris-buffered saline (TBS), the membranes were incubated with the primary antibodies at 4 °C overnight.

The primary antibodies are as follows: anti-GADD45B (1:2000, Bioss), anti-CDKN1B (1:1000, Cell Signaling Technology), anti-HILPDA (1:1000, Santa Cruz), anti-CLDN5 (1:4000, Novus Biologicals), anti-GAPDH (1:1000, Cell Signaling Technology).

The membranes were washed and incubated with corresponding secondary antibodies. Reagents for Western blot detection by ECL substrate working solution (Abbkine, China) and images were captured using an imaging system (Bio-Rad). The density of bands was analyzed by ImageJ software.

Histological Analysis

Knee joints were harvested from rats and fixed with 4% paraformaldehyde. Briefly, all specimens were decalcified in a 10% EDTA solution at room temperature for 28 days. The knee joints were embedded in paraffin wax and serially sectioned at thicknesses of 4μm. The sections were stained with safranin O-fast green or hematoxylin and eosin (H&E) according to manufacturer’s instructions. A minimum of five random images from three samples were analyzed per group.

Statistical Analysis

Statistical analysis was conducted using GraphPad Prism version 9.0. Two-tailed Student’s t-tests were employed to analyze the differences between two groups. All data are expressed as means ± SEM, and P < 0.05 was considered statistically significant.

Results Identification of the DEGs in OA Patients

The workflow of the study design is shown in Figure 1. Two expression datasets containing both control and OA samples were utilized to explore the potential role of mRNAs in OA. The profiles of DEGs from GSE169077 and GSE114007 were analyzed separately between control and OA samples. In the GSE169077 dataset, a total of 675 DEGs were identified based on a q value < 0.05 and |log2FC| > 1. Among these, 462 genes were found to be upregulated, and 213 genes were downregulated (Figure 2A and B). Conversely, in the GSE114007 dataset, a total of 2736 DEGs were identified, with 1688 genes showing significant upregulation and 1048 genes displaying significant downregulation (Figure 2C and D).

Figure 1 The workflow of the study.

DEGs, differentially expressed genes; PPI, Protein–protein interaction; KEGG, Kyoto Encyclopedia of Genes and Genomes; GO, Gene Ontology.

Figure 2 Differentially expressed genes (DEGs) between control and OA samples in the GEO datasets. (A) Volcanic map for differential expression analysis of GSE169077. Blue represents down-regulated genes, red represents up-regulated genes and black represents undifferentiated genes. (B) Heat map for differential expression analysis of GSE169077. (C) Volcanic map for differential expression analysis of GSE114007. Blue represents down-regulated genes, red represents up-regulated genes and black represents undifferentiated genes. (D) Heat map for differential expression analysis of GSE114007.

Functional Enrichment Analysis of DEGs

GSEA was performed to identify distinct pathways associated with DEGs in the progression of OA. The results revealed that pathways related to endochondral ossification and collagen formation were significantly upregulated (Figure 3A). Conversely, various metabolic pathways, particularly those involved in glycogen metabolism, NADP metabolism, cysteine metabolism, and lactate metabolism, were significantly inhibited in the OA samples (Figure 3B). Additionally, high enrichment scores were observed for pathways related to the formation of skeletal muscle in OA.

Figure 3 GSEA enrichment analysis of the two datasets. (A) GSEA enrichment analysis of the upregulated genes.(B) GSEA enrichment analysis of the downregulated genes.The left panel, GSE169044; the right panel, gse114007.

Identification of Clinically Significant Modules

The WGCNA algorithm was utilized to construct the co-expression gene network, and a soft threshold (scale-free R2 = 0.9) was chosen for each dataset to construct a scale-free network, which closely resembles the real biological network state (Figure 4A and B). Then, gene modules were identified through average hierarchical clustering and dynamic tree cutting. The minimum module size was set to 30, with a height cut of 0.1 for both datasets to merge the modules. A total of 5 and 20 distinct modules, labeled with different colors, were identified using the GSE169077 and GSE114007 datasets, respectively (Figure 4C and D). Heatmaps illustrated the correlation between these modules in control and OA samples (Figure 4E and F). In the GSE169077 dataset, the highest correlation between MM and GS was observed in the blue module (r=−0.84, p=0.001). Similarly, in GSE114007, the magenta (r = 0.94, p = 2 e - 18), purple (r = 0.82, p = 2 e - 10), cyan (r = 0.86, p = 3 e - 12), pink (r = 0.84, p=6e-11), royalblue modules (r=−0.83, p=1e-10) are most significantly associated with OA. Therefore, modules exhibiting the highest association were selected as the clinically significant modules for subsequent analysis.

Figure 4 Construction of the co-expression network of the OA and control samples in the two datasets. (A and B) Soft-thresholding power analysis and mean connectivity was used to obtain the scale-free fit index of network topology in GSE169077 (A) and GSE114007 (B). (C and D) The hierarchical clustering of all genes in GSE169077 (C) and GSE114007 (D) dataset with corresponding color assignments. Each branch represents a gene, and each color represents a module in the constructed gene co-expression network by WGCNA. The color band shows the results obtained from the automatic single-block analysis. (E and F) The heatmap of the correlation between gene clustering module in GSE169077 (E) and GSE114007 (F) dataset. Each module contains the correlation coefficient and p value.

Identification of the Hub Genes

To further screen for the most promising candidate hub genes in OA, a total of 35 intersected genes from the gray module were identified based on previous analyses (Figure 5A) and further analyzed using GO and KEGG pathways. GO analysis revealed that these hub genes were primarily enriched in biological processes such as cellular response to hypoxia, regulation of cell proliferation and vascular permeability, and response to lipopolysaccharide. At the molecular level, these genes were found to be enriched for functions such as receptor binding and protein kinase A binding (Figure 5B). KEGG analysis indicated enrichment in pathways, including the PI3K-Akt signaling pathway, ECM-receptor interaction, HIF-1 signaling pathway, and FoxO signaling pathway (Figure 5C). Next, the STRING database was used to analyze the interaction relationships of intersected genes. Through PPI analysis, the top 10 hub genes, including GADD45B, CLDN5, CITED2, CDKN1B, ANGPTL4, ADM, SLC2A3, PFKFB3, MXI1, and HILPDA, were identified based on the highest degrees of connectivity using the cytoHubba plug-in in Cytoscape (Figure 5D). Box plot results demonstrated that these genes were all significantly downregulated in OA samples compared to controls (Figure 5E).

Figure 5 Screening of hub genes related with OA. (A) The Venn diagram demonstrates the intersection genes of module significantly associated with OA in the two datasets. (B) GO enrichment analysis of genes related to OA. The y-axis shows the number of genes enriched in this term. (C) KEGG enrichment analysis of genes related to OA. The abscissa is the rich factor, and the ordinate is the KEGG path name. The greater the rich factor, the more significant the enrichment level of DEGs in this pathway. The color of the circle represents the P value. The smaller of P value, the more reliable the enrichment significance of DEGs in this pathway. The size of the circle indicates the number of genes enriched in the pathway. (D) Establishment of PPI network on the STRING database. PPI, protein-protein interaction. (E) The expression level of the hub genes in the GSE169077 and GSE114007. HC, healthy control, two-tailed t-test. All data are presented as the means ± SEM. ns, not significant; *P < 0.05; **P < 0.01; ***P <0.001.

Validation of Hub Genes in OA Animal Models

To further validate the expression levels of these hub genes in the progression of OA, we successfully established OA models in rats (Figure 6A and B). The histopathological results demonstrated that the surface of cartilage in the NC group was smooth, and morphology of the chondrocytes was normal and arranged neatly and orderly. While rat in the OA group had defective cartilage surface and hypertrophied chondrocytes. Moreover, The OA group showed more pronounced synovial inflammation than that of the NC group with enlargement of the synovial lining cell layers and infiltration of inflammatory cells. As shown in Figure 6C, although these hub genes were expressed in the control group, the mRNA levels of GADD45B, CLDN5, HILPDA, CDKN1B, and ADM were consistently lower in the OA groups, which was largely consistent with the results of our previous screening. The receiver operating characteristic (ROC) curve, designed to evaluate diagnostic values, revealed that the area under the curve (AUC) values of GADD45B, CLDN5, HILPDA, and CDKN1B were >0.85, while the AUC value of ADM was 0.8333. (Figure 6D). Subsequently, we further validated the lower expression of GADD45B, CLDN5, HILPDA, and CDKN1B in OA rat model by Western blotting (Figure 6E).

Figure 6 Validation of the expression of hub genes in rat OA model. (A) Schematic strategy of the OA rat model. (B) Representative images of histopathological staining of cartilage and synovial sections from different groups, n = 3. Scale bar = 100 μm (HE staining), Scale bar = 50 µm (SO/FG staining).(C) RT-qPCR for the levels of the hub genes in rat OA model and control group, n = 6. NC, normal control. Two-tailed t-test. All data are presented as the means ± SEM. ns, not significant; *P < 0.05; **P < 0.01. (D) Determination of the diagnostic efficiency of pivotal genes by ROC curve analysis. AUC, area under the curve. (E) Western blotting and semi-quantitative analysis of protein, n = 3. Two-tailed t-test. All data are presented as the means ± SEM. *P < 0.05; **P < 0.01.

Discussion

With the global population gradually aging, OA, the most prevalent chronic joint disease, is being increasingly diagnosed in a significant number of elderly individuals worldwide, posing a potential threat to patient health in severe cases.17 Timely diagnosis of OA is essential for slowing disease progression, yet effective early diagnostic methods remain lacking. The search for specific biomarkers and key genes has become a focal point in OA research in recent years. The utilization of bioinformatics based on high-throughput sequencing data has played a pivotal role in various medical research fields, leading to improved disease diagnosis and personalized treatment.18 Identifying biomarkers associated with OA is instrumental in enhancing OA diagnosis and early detection. In this study, we utilized OA tissue expression data from the GEO database to screen for key genes influencing OA through differential expression analysis, WGCNA, and PPI analysis. Additionally, we employed a rat OA model to validate the expression levels of key genes in OA tissues.

In this study, DEGs were analyzed using data from two GEO databases. GSEA revealed significant upregulation of pathways related to endochondral ossification and collagen formation, alongside inhibition of various metabolic pathways in OA. Chondrocytes and osteoblasts, pivotal in skeletal development, have been implicated in the progression of conditions such as OA, bone dysplasia, and osteoporosis.19,20 Furthermore, aberrant metabolic activities may contribute significantly to cartilage degeneration and OA progression.21 Next, we constructed a gene co-expression network based on gene expression associated with clinical characteristics of OA, based on which 35 genes were identified. GO and KEGG enrichment analyses were used to explore biological functions and signaling pathways. GO analysis revealed enrichment primarily in biological functions related to cellular response to hypoxia. KEGG analysis showed enrichment in pathways, including the HIF-1 signaling pathway, ECM-receptor interaction, and FoxO signaling pathway. In OA, HIF-1α has been identified as an important regulator of chondrocyte differentiation and matrix synthesis, essential for maintaining normal chondrocyte function, as demonstrated in both in vitro and in vivo studies.22,23 Studies have also highlighted the pivotal role of HIF-1α degradation in cartilage degeneration in OA, with stabilization of HIF-1α shown to prevent OA progression.24 Moreover, OA progression is closely linked to ECM catabolism in cartilage.25,26 Increased ECM catabolism disrupts articular cartilage integrity and homeostasis.27,28 Additionally, the FoxO signaling pathway emerges as highly dysregulated in OA cartilage compared to normal tissues.29 Combined deletion of FoxO1, FoxO3, and FoxO4 in the myeloid lineage has been found to promote osteoclast formation and bone resorption, resulting in decreased trabecular and cortical bone mass,30 indicating that FoxOs play an important role in maintaining articular cartilage homeostasis.31 These findings corroborate the significant role of the OA-related genes identified by WGCNA in our study.

In our pursuit to further identify pivotal genes influencing the progression of OA, we conducted PPI analysis and ultimately obtained 10 hub genes: GADD45B, CLDN5, CITED2, CDKN1B, ANGPTL4, ADM, SLC2A3, PFKFB3, MXI1, and HILPDA. Among these, GADD45B, CLDN5, HILPDA, CDKN1B, and ADM exhibited consistently differential expression by RT-qPCR based on rat OA samples. In addition, GADD45B, CLDN5, HILPDA and CDKN1B had AUC values > 0.85 by ROC curve analysis, suggesting that they may have good diagnostic performance. Further, through Western blotting, we verified the low expression of GADD45B, CLDN5, HILPDA, CDKN1B in the OA compared to control group. GADD45B, a small molecular protein, regulates cell cycle progression and cellular stress response.32 Previous research has demonstrated its regulatory role in various inflammatory diseases.33 In our study, we observed a reduction in GADD45B expression and validated this through RT-qPCR in OA tissues. Given the decreased expression of GADD45B, its inhibitory effect on apoptosis may be attenuated. We hypothesize that GADD45B may be involved in the pathogenesis of OA by promoting the apoptosis of chondrocytes and accelerating the destruction of cartilage.34,35 Additionally, studies have reported upregulation of HILPDA expression during chondrogenic differentiation and association with promoting chondrogenic differentiation and regulating the expression of chondrogenic-related factors.36 Under physiological conditions, chondrocytes are in a relatively hypoxic environment, which is essential for maintaining chondrocyte homeostasis. However, more research focuses on the HIF family, and there are few studies on HILPDA in OA, which may be worthy of further study to help us understand the regulatory mechanism of hypoxia response in the occurrence and development of OA more comprehensively. Furthermore, increasing evidence suggests that CDKN1B is a key regulator of cell cycle progression and serves as a valuable marker in conditions such as osteoporosis,37 atherosclerosis,38 and Alzheimer’s disease.39 This further indicates an imbalance between cell proliferation and apoptosis of OA chondrocytes. An in-depth understanding of the mechanism may help to reveal the regulatory network of OA and provide new insights for the prediction and treatment of OA. The expression of ADM is regulated by TNF-α, can impact adipogenesis and osteogenesis, and has been shown to be significantly downregulated in OA. Moreover, there is evidence suggesting that ADM downregulation may be involved in the development of OA by facilitating synovial cell apoptosis and chondrocyte dedifferentiation in inflammatory arthritis by increasing oxidative stress products and pro-inflammatory cytokines.40,41 However, the role of CLDN5 in the occurrence and development of OA has not yet been elucidated. As a tight junction protein, we speculate that the decrease of CLDN5 may lead to the increase of endothelial permeability, which in turn affects the metabolism of chondrocytes. However, the specific regulatory mechanism still needs to be further investigation. These findings underscore the significant roles played by pivotal genes in the progression of OA. However, this requires further validation in an independent cohort, which is a limitation of the study.

In conclusion, we identified four key genes (GADD45B, CLDN5, HILPDA and CDKN1B) associated with OA through integrated bioinformatics analysis. Enrichment analysis revealed that these genes affected the occurrence and development of OA in multiple pathways. Taken together, these genes hold potential as biomarkers for OA diagnosis and treatment. Further investigations through in vitro and in vivo experiments are warranted to elucidate the precise mechanisms by which these key genes affect OA.

Ethics Statement

This study was reviewed and approved by Institutional Animal Care and Use of Chongqing Medical University (IACUC-CQMU), with the approval number: IACUC-CQMU-2024-0473, dated: 2024-06-07.

Acknowledgments

We acknowledge support from the Graduate Student Research and Innovation Projects in Chongqing (Grant No. CYS23376); Natural Science Foundation of Chongqing, China (Grant No. CSTB2022NSCQ-MSX0134); X9449 Xie Youhong Demonstration of Internet+ Medical and Nursing Base Model Innovation (Grant No. 41021300060299).

Disclosure

The author(s) report no conflicts of interest in this work.

References

1. Gu Y, Hu Y, Zhang H, et al. Single-cell RNA sequencing in osteoarthritis. Cell Prolif. 2023;56(12):e13517. doi:10.1111/cpr.13517

2. Sharma L, Solomon CG. Osteoarthritis of the Knee. N Engl J Med. 2021;384(1):51–59. doi:10.1056/NEJMcp1903768

3. Englund M. Osteoarthritis, part of life or a curable disease? A bird’s-eye view. J Intern Med. 2023;293(6):681–693. doi:10.1111/joim.13634

4. Cecchini MJ, Hosein K, Howlett CJ, et al. Comprehensive gene expression profiling identifies distinct and overlapping transcriptional profiles in non-specific interstitial pneumonia and idiopathic pulmonary fibrosis. Respir Res. 2018;19(1):153. doi:10.1186/s12931-018-0857-1

5. Samy BA, Raman K, Velayutham S, et al. Natural product extract fractions as potential arthritis treatments: a detailed analysis using in-silico, in-vivo, and in-vitro methods. Int Immunopharmacol. 2025;144:113595. doi:10.1016/j.intimp.2024.113595

6. Jin P, Liu H, Chen X, et al. From bench to bedside: the role of extracellular vesicles in cartilage injury treatment. Biomater Res. 2024;28:0110. doi:10.34133/bmr.0110

7. Su S, Tian R, Jiao Y, et al. Ubiquitination and deubiquitination: implications for the pathogenesis and treatment of osteoarthritis. J Orthopaedic Transl. 2024;49:156–166. doi:10.1016/j.jot.2024.09.011

8. Liu Y, Hao R, Lv J, et al. Targeted knockdown of PGAM5 in synovial macrophages efficiently alleviates osteoarthritis. Bone Res. 2024;12(1):15. doi:10.1038/s41413-024-00318-8

9. Ali N, Turkiewicz A, Hughes V, et al. Proteomics profiling of human synovial fluid suggests increased protein interplay in early-osteoarthritis (OA) that is lost in late-stage OA. mol Cell Proteomics. 2022;21(3):100200. doi:10.1016/j.mcpro.2022.100200

10. Mobasheri A, Thudium CS, Bay-Jensen A-C, et al. Biomarkers for osteoarthritis: current status and future prospects. Best Pract Res Clin Rheumatol. 2023;37(2):101852. doi:10.1016/j.berh.2023.101852

11. Nangraj AS, Selvaraj G, Kaliamurthi S, et al. Integrated PPI- and WGCNA-retrieval of hub gene signatures shared between barrett’s esophagus and esophageal adenocarcinoma. Front Pharmacol. 2020;11:881. doi:10.3389/fphar.2020.00881

12. Wen J, Ren T, Zheng J, et al. Identification and verification of pivotal genes promoting the progression of atherosclerosis based on WGCNA. Artif Cells Nanomed Biotechnol. 2023;51(1):276–285. doi:10.1080/21691401.2023.2203185

13. Gao XM, Zhou X-H, Jia M-W, et al. Identification of key genes in sepsis by WGCNA. Prev Med. 2023;172:107540. doi:10.1016/j.ypmed.2023.107540

14. Chen Y, Liao R, Yao Y, et al. Machine learning to identify immune-related biomarkers of rheumatoid arthritis based on WGCNA network. Clin Rheumatol. 2022;41(4):1057–1068. doi:10.1007/s10067-021-05960-9

15. Yu G, Wang L-G, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. 2012;16(5):284–287. doi:10.1089/omi.2011.0118

16. Ravasz E. Hierarchical organization of modularity in metabolic networks. Science. 2002;297(5586):1551–1555. doi:10.1126/science.1073374

17. Sampath SJP, Venkatesan V, Ghosh S, et al. Obesity, metabolic syndrome, and osteoarthritis-an updated review. Curr Obes Rep. 2023;12(3):308–331. doi:10.1007/s13679-023-00520-5

18. Akalin PK. Introduction to bioinformatics. mol Nutr Food Res. 2006;50(7):610–619. doi:10.1002/mnfr.200500273

19. Chen Y, Mehmood K, Chang Y-F, et al. The molecular mechanisms of glycosaminoglycan biosynthesis regulating chondrogenesis and endochondral ossification. Life Sci. 2023;335:122243. doi:10.1016/j.lfs.2023.122243

20. Anwar A, Sapra L, Gupta N, et al. Fine-tuning osteoclastogenesis: an insight into the cellular and molecular regulation of osteoclastogenesis. J Cell Physiol. 2023;238(7):1431–1464. doi:10.1002/jcp.31036

21. Zheng L, Zhang Z, Sheng P, et al. The role of metabolism in chondrocyte dysfunction and the progression of osteoarthritis. Ageing Res Rev. 2021;66:101249. doi:10.1016/j.arr.2020.101249

22. Hu S, Zhang C, Ni L, et al. Stabilization of HIF-1α alleviates osteoarthritis via enhancing mitophagy. Cell Death Dis. 2020;11(6):481. doi:10.1038/s41419-020-2680-0

23. Zhang C, Yang F, Cornelia R, et al. Hypoxia-inducible factor-1 is a positive regulator of Sox9 activity in femoral head osteonecrosis. Bone. 2011;48(3):507–513. doi:10.1016/j.bone.2010.10.006

24. Zhang H, Wang L, Cui J, et al. Maintaining hypoxia environment of subchondral bone alleviates osteoarthritis progression. Sci Adv. 2023;9(14):eabo7868. doi:10.1126/sciadv.abo7868

25. Lin AC, Seeto BL, Bartoszko JM, et al. Modulating hedgehog signaling can attenuate the severity of osteoarthritis. Nat Med. 2009;15(12):1421–1425. doi:10.1038/nm.2055

26. Miao Y, Chen Y, Xue F, et al. Contribution of ferroptosis and GPX4’s dual functions to osteoarthritis progression. EBioMedicine. 2022;76:103847. doi:10.1016/j.ebiom.2022.103847

27. Rahmati M, Nalesso G, Mobasheri A, et al. Aging and osteoarthritis: central role of the extracellular matrix. Ageing Res Rev. 2017;40:20–30. doi:10.1016/j.arr.2017.07.004

28. Hodgkinson T, Kelly DC, Curtin CM, et al. Mechanosignalling in cartilage: an emerging target for the treatment of osteoarthritis. Nat Rev Rheumatol. 2022;18(2):67–84. doi:10.1038/s41584-021-00724-w

29. Fisch KM, Gamini R, Alvarez-Garcia O, et al. Identification of transcription factors responsible for dysregulated networks in human osteoarthritis cartilage by global gene expression analysis. Osteoarthritis Cartilage. 2018;26(11):1531–1538. doi:10.1016/j.joca.2018.07.012

30. Bartell SM, Kim H-N, Ambrogini E, et al. FoxO proteins restrain osteoclastogenesis and bone resorption by attenuating H2O2 accumulation. Nat Commun. 2014;5:3773. doi:10.1038/ncomms4773

31. Almeida M, Porter RM. Sirtuins and FoxOs in osteoporosis and osteoarthritis. Bone. 2019;121:284–292. doi:10.1016/j.bone.2019.01.018

32. De Smaele E, Zazzeroni F, Papa S, et al. Induction of gadd45beta by NF-kappaB downregulates pro-apoptotic JNK signalling. Nature. 2001;414(6861):308–313. doi:10.1038/35104560

33. Shen XY, Shi S-H, Li H, et al. The role of Gadd45b in neurologic and neuropsychiatric disorders: an overview. Front Mol Neurosci. 2022;15:1021207. doi:10.3389/fnmol.2022.1021207

34. Zhang X, Bu Y, Zhu B, et al. Global transcriptome analysis to identify critical genes involved in the pathology of osteoarthritis. Bone Joint Res. 2018;7(4):298–307. doi:10.1302/2046-3758.74.Bjr-2017-0245.R1

35. Cao H, Fu Y, Zhang Z, et al. Unbiased transcriptome mapping and modeling identify candidate genes and compounds of osteoarthritis. Front Pharmacol. 2022;13:888533. doi:10.3389/fphar.2022.888533

36. Zhao Y, Wang X, Nie K. IRF1 promotes the chondrogenesis of human adipose-derived stem cells through regulating HILPDA. Tissue Cell. 2023;82:102046. doi:10.1016/j.tice.2023.102046

37. Xu Y, Sun L, Hu J, et al. Knockdown of hsa_circ_0001275 reverses dexamethasone-induced osteoblast growth inhibition via mediation of miR-377/CDKN1B axis. PLoS One. 2021;16(5):e0252126. doi:10.1371/journal.pone.0252126

38. Yamada S, Senokuchi T, Matsumura T, et al. Inhibition of local macrophage growth ameliorates focal inflammation and suppresses atherosclerosis. Arterioscler Thromb Vasc Biol. 2018;38(5):994–1006. doi:10.1161/atvbaha.117.310320

39. Ogawa O, Lee H-G, Zhu X, et al. Increased p27, an essential component of cell cycle control, in alzheimer’s disease. Aging Cell. 2003;2(2):105–110. doi:10.1046/j.1474-9728.2003.00042.x

40. Qi X, Yu F, Wen Y, et al. Integration of transcriptome-wide association study and messenger RNA expression profile to identify genes associated with osteoarthritis. Bone Joint Res. 2020;9(3):130–138. doi:10.1302/2046-3758.93.Bjr-2019-0137.R1

41. Li C, Wei P, Wang L, et al. Integrated analysis of transcriptome changes in osteoarthritis: gene expression, pathways and alternative splicing. Cartilage. 2023;14(2):235–246. doi:10.1177/19476035231154511