Introduction

Osteoarthritis (OA), regarded as a degenerative and chronic joint disease, typically presents with joint pain that is exacerbated by use and alleviated with rest,¹ has troubled more than 240 million people globally and has been a leading cause of disability in the elderly.^2,3 Specifically, the prevalence of hip OA has reached 0.85%, and that of knee OA has reached 3.8% worldwide.⁴ There are numerous risk factors linked to OA, such as age, obesity, sex, previous joint injury, ethnicity, and malalignment of joints.^2–4 The risk of OA typically tends to increase with age, and women are more susceptible to OA than men are.⁵ Hence, more diagnostic methods and novel treatments need to be invented, and more exploration of pathogenesis should be carried out.

Currently, many techniques are available for diagnosing OA. In terms of patient history, X-ray, ultrasound, magnetic resonance imaging (MRI), micro-CT, and arthroscopy have been widely used in the clinic.^6,7 However, numerous current diagnostic methods have limitations, triggering prolongation of the illness course. Early pathological changes cannot be detected by X-ray. In addition, the symptoms of patients may be inconsistent with the severity examined by X-ray, partly owing to the invisibility of endogenous pain sources and difficulty in assessing exogenous pain features. Ultrasound is unable to detect subchondral bone lesions, and outcomes may differ because of the proficiency of the operators. MRI has not been fully developed and is not completely harmless to the human body, increasing the possibility of causing nephrogenic fibrosis due to intravenous administration.^6–8 In terms of OA treatment, lifestyle changes, the use of anti-inflammatory agents and surgery⁹ have been used to cure OA. Nevertheless, none of the treatments can eliminate symptoms or delay progression. We still face great challenges in the diagnosis and treatment of OA.

To date, an increasing number of signal transcription pathways associated with OA have been discovered. Canonical Wnt signalling, which is transduced by 19 Wnt genes, plays a pivotal role in the homeostasis of joints. The interaction between Wnt and the β-catenin destruction complex is initiated by the binding of Wnt and Frizzled receptors. β-catenin then moves into the nucleus and associates with T-cell factor/lymphoid-enhancer factor family members, activating the expression of Wnt target genes whose overexpression leads to cartilage hypertrophy, while inactivation is responsible for cartilage degradation, indicating that Wnt signalling should be managed carefully to prevent OA.^10,11 The classical NF-κB pathway is also linked with the development of OA. Numerous immune mediators, such as IL-1β, TNF-α, and TLRs, activate the IκΒ-kinase (IKK) complex (IKKα, ΙΚΚβ, ΙΚΚγ), triggering the phosphorylation and subsequent degradation of IKB. Matrix-degrading enzymes such as matrix metalloproteinases (MMPs) and the ADAMTS family are ultimately expressed via the transcription of target genes stimulated by IKB, promoting osteoarthritic cartilage erosion.^11–13 In addition to key pathways, genes relevant to OA genesis have also been explored. Yao et al¹⁴ reported that chondrocyte ferroptosis promotes the progression of OA. Liu et al¹⁵ reported that the SLC3A2 gene inhibited ferroptosis and suppressed the deterioration of OA. The expression of SLC3A2 is downregulated in people with OA compared with healthy controls.¹⁵ Moreover, there is also some evidence that the expression of the BCL6 gene is positively correlated with the severity of OA.¹⁶ The transcription level of BCL6 is upregulated in OA, which indicates that the BCL6 gene may be a predictive biomarker for OA. Mesenchymal stem cells (MSCs) have become a promising therapeutic approach for osteoarthritis (OA) because of their capacity to differentiate into chondrocytes and promote cartilage regeneration, effectively alleviating the detrimental effects of a macrophage-induced inflammatory microenvironment.¹⁷

The diagnostic dilemma and many molecular research foundations have encouraged further studies on genes and biomarkers associated with OA. In this study, we performed synovial tissue sequencing on samples from 35 patients with OA and 32 healthy controls. Six key genes were identified to construct models for accurately predicting OA. Furthermore, we verified the expression of the KDELR3 gene in clinical samples and discovered markedly upregulated expression of the KDELR3 gene in patients with OA.

Methods

Patient Summary

We collected a total of 35 OA patients and 32 healthy controls from four GEO cohorts. The filter standard was the cohort containing both synovial tissue samples from OA patients and healthy controls, including GSE1919, GSE55235, GSE55457, and GSE82107, the basic information of these cohorts is displayed in Table 1. To validate the identified biomarkers, we also collected 8 real-world samples from our institute between October 1st, 2023, and January 1st, 2024, four from OA patients and four from healthy controls. In the OA group, samples were collected from patients undergoing total knee arthroplasty, whereas in the control group, samples were obtained from patients without a history of OA who had sustained knee injuries due to trauma. To ensure the integrity of the samples, they were promptly fixed in formalin within half an hour postsurgery. This study was conducted in compliance with the principles of the Declaration of Helsinki and was approved by the Ethics Committee of the Second Affiliated Hospital of Anhui Medical University [approval number: YX2023-137(F1)].

Table 1 Basic Information of the Enrolled Cohorts

Elimination of Batch Effects

Batch effects represent the nonbiological disparities between two or more datasets. To nullify the bias engendered by batch effects in this investigation and harmonize the transcription profiles across the three GEO cohorts, the ComBat algorithms from the “sva” package were deployed to eradicate the batch effects between these GEO-derived cohorts.

Calculation of the Scores of Signalling Pathways

Gene set enrichment analysis (GSEA) is a computational approach for ascertaining whether a designated gene set exhibits statistically significant deviations between two groups. We implemented GSEA to initially contrast the various activated signalling pathways between OA patients and healthy controls. The backdrop file of molecular signature gene sets was procured from MSigDB, C5: Biological Process, comprising a total of 7530 gene sets.^18,19

Diagnostic Model Derived from Integrative Approaches Utilizing Machine Learning Methodologies

To develop a consensus model with high accuracy and stable performance for distinguishing OA patients from normal patients, we integrated 10 machine learning algorithms and 101 algorithm combinations. The integrative algorithms included elastic network (Enet), Lasso, Ridge, Stepglm[both], Stepglm[backward], glmBoost, Latent Dirichlet Allocation (LDA), NaiveBayes, plsRglm, Random Forest (RF), Support Vector Machine (SVM). The signature generation procedure was as follows: (a) Identifying the most activated pathways in OA patients from all four GEO cohorts; (b) performing 101 algorithm combinations on the selected genes from the most activated pathways; (c) training all the models in GSE55235 and validation in the other three cohorts, which were not used for the filtration of pathways; and (d) calculating the AUC value across all the enrolled cohorts.

Infiltration of Immune Cells

Single-sample gene set enrichment analysis (ssGSEA) was used to quantify and assess the infiltration of immune cells via transcriptional data. The gene sets encompassing 28 immune cells were procured from the study conducted by Charoentong et al.²⁰

Immunohistochemical Staining

To validate the diagnostic value of KDELR3, we performed IHC staining of normal and OA tissues from clinical samples. FFPE samples were sectioned into 4.5 μm thick slices. These slices were subsequently hydrated. Antigen retrieval was then conducted with citrate buffer solution. The tissue sections were blocked for one hour and then incubated with a primary anti-KDELR3 antibody (1:1000, PA6579; Abmart, Shanghai, China) overnight at 4°C. The following day, secondary antibodies labelled with goat anti-rabbit IgG-HRP (1:1000, M21002; Abmart, Shanghai, China) were applied, and the samples were stained with DAB.

Statistics

All the statistical computations were executed via R (Version: 4.2.2). A Pearson correlation analysis was applied for continuous data. Pearson correlation analysis was used for continuous data. The chosen pathways were illustrated in a heatmap via the R package “pheatmap”. The kappa statistic is a statistical measure employed for juxtaposing predictive and actual subtypes. The Kruskal‒Wallis test was used for statistical comparisons among more than two groups, whereas the Wilcoxon test was used for pairwise comparisons.²¹ A two-tailed P value <0.05 was considered statistically significant.

Results

In the present study, we enrolled a total of 67 samples from four cohorts, 35 OA patients and 32 healthy controls, and we also collected clinical samples from four OA patients and four controls to validate the protein levels of the identified key genes. The process of the current research is displayed in Figure 1, which mainly includes cohort collection, pathway identification, model construction and real-world validation.

Figure 1 Flow chart summarizing the current research.

Preparation of the Enrolled Cohorts

To keep the findings more stable, we first removed the batch effects via the “sva” package. The first component variance was 49.8%, and the second component variance was 18.1% for the original difference in the gene expression matrix of the four GEO cohorts. After the bath effects were removed, the variance of the first component decreased to 9.8%, and that of the second component decreased to 8.3% (Figure 2A). These results indicate that the batch effects were removed. We subsequently used the “GSVA” package to calculate the pathway activation scores of 7530 gene sets among all the enrolled patients, and heatmaps were generated to display the distribution of these pathways in each cohort (Figure 2B).

Figure 2 Preprocessing of the input data for subsequent analysis. (A) The combat algorithm eliminates the batch effect. (B) Heatmaps showing the activation of signalling pathways in each cohort.

Screening for Activated Signalling Pathways in Osteoarthritis Patients

With the preset threshold value, p value less than 0.05 and log2fc value greater than 0, we screened the activated signalling pathways in OA patients and selected 897 pathways from GSE55235, 545 pathways from GSE55457, 478 pathways from GSE1919 and 754 pathways from GSE82107 (Figure 3A). In addition, with the help of UpSet plots and Venn diagrams, 16 pathways were identified for subsequent analysis, with a total of 478 nonduplicated genes. The 16 pathways refer to the cellular respiration chain, protein transport, lysosomal and endocytosis, bone remodelling and skeletal muscle regeneration (Figure 3B and Table S1).

Figure 3 Screening for activated signalling pathways in osteoarthritis patients. (A) Different signalling pathways between OA patients and healthy controls in the four GEO cohorts. (B) UpSet plot and Venn diagram delineating the 16 mutually activated signalling pathways across the four cohorts.

Machine Learning-Based Integration Develops an Osteoarthritis Diagnostic Model

Among the 478 selected genes, 298 were expressed in all four GEO cohorts; therefore, they were used as the input data. As mentioned in the methods, 101 types of algorithms were considered for the construction of the diagnostic model, and the results of 58 combined approaches were successfully output because the other approaches are not suitable for the current research, possibly because of the small amount of input data. The AUC values of each approach are listed in Figure 4A, with the two-step model of glmBoost plus RF showing the highest average AUC value of 0.95. This model was generated with six genes, including LY86, SORL1, KDELR3, CSK, PTGS1, and PTGS2. Furthermore, we double checked the predictive value of this model via kappa analysis and observed preferable consistency in GSE55235 (Kappa = 1.000, P < 0.001), GSE55457 (Kappa = 0.700, P < 0.001), GSE82107 (Kappa = 0.643, P = 0.004) and GSE1919 (Kappa = 1.000, P < 0.001) (Figure 4B).

Figure 4 Machine learning-based integration has led to the development of an osteoarthritis diagnostic model. (A) A total of 101 kinds of prediction models were constructed in GSE55235 and further validated across the other three validation cohorts via AUC values. (B) Evaluating the consistency of the prediction and actual subgroups via kappa statistics.

Selection of Significant Biomarkers for Real-World Validation

For the six selected genes in the diagnostic model, we compared their expression between 35 OA patients overall and 32 healthy controls. We observed higher levels of LY86, SORL1, KDELR3, CSK and PTGS1 in OA patients than in controls, whereas PTGS2 showed the opposite trend (all P < 0.05, Figure 5A and Table S2). We observed that KDELR3 was expressed at higher levels in OA patients than were the other genes; therefore, we validated its protein level via IHC staining of tissues collected from our own institute, four from OA patients and four from healthy controls. We observed that synovial tissues expressed very low levels of KDELR3, whereas OA patients presented high levels of KDELR3 in synovial tissues (Figure 5B).

Figure 5 Differential expression of the six genes selected for the osteoarthritis diagnostic model. (A) The expression of six genes among normal samples and OA patients from the four cohorts. (B) Immunohistochemical staining showing the protein level of KDELR3 in real-world samples from normal controls and OA patients.

Activated KEGG Signalling Pathways in OA Patients

To further evaluate the activated pathways, we also performed KEGG pathway analysis. When we compared the KEGG pathway activation between 35 OA patients and 32 healthy controls, we observed the most activated signalling pathways in OA, including the lysosome (NES = 1.936, P = 0.0038), riboflavin metabolism (NES = 1.900, P < 0.001), glycan biosynthesis (NES = 1.822, P = 0.010), and glycosaminoglycan degradation (NES = 1.694, P = 0.0040) pathways (Figure 6), which are tightly associated with cellular respiration and cell proliferation, which is consistent with the findings in Figure 3. In addition, we also observed altered pathways related to base excision repair (NES = 1.717, P = 0.018) and DNA replication (NES = 1.660, P = 0.047) (Figure 6), which indicated that genetic alterations in the process of OA might also stimulate the immune microenvironment.

Figure 6 The activation pathways were investigated with all the enrolled samples.

Status of Immune Cells in the Synovial Microenvironment of OA Patients

Immunocytes play pivotal roles in the onset and progression of OA.²² We compared infiltration of the 28 immunocytes and observed that infiltration of most types, including macrophage, natural killer T cell, regulatory T cell, type 1 T helper cell, activated B cell, activated CD8 T cell, effector memory CD4 T cell, gamma delta T cell, immature B-cell and immature dendritic cell, was increased in OA patients, whereas infiltration of type 17 T cells, type 2 T cells, and eosinophils was increased in normal controls (Figure 7A). These results indicate that OA is an immune-hot disease. We further directly compared the correlation between KDELR3 and immunocytes in OA tissues and observed that KDELR3 might be involved in the activation and infiltration of most types of immunocytes (Figure 7B), especially for effector memory CD4 T cells (R_pearson = 0.58, P < 0.001, Figure 7C) and natural killer T cells (R_pearson = 0.53, P < 0.001, Figure 7D).

Figure 7 Infiltration of immunocytes among OA patients and correlation with KERLR3. (A) Differentially infiltrated immunocytes among OA patients and normal controls. (B) Correlations between KDELR3 and all types of immunocytes. (C) KDELR3 was positively associated with the infiltration level of effector memory CD4 T cells. (D) KDELR3 was positively associated with the infiltration level of natural killer T cells. *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.

Discussion

OA is a degenerative joint disease characterized by the progressive deterioration of cartilage and chronic inflammation of the synovial membrane, serving as a principal source of incapacitation among the global geriatric population.²³ While the exact etiology contributing to osteoarthritis remains elusive, it is increasingly evident that inflammation plays an instrumental role in joint tissue degradation subsequent to injury.²⁴ The synovium is a specialized form of connective tissue that envelopes diarthrodial joints, enwraps tendons, and constitutes the linings of bursae and adipose pads. In synovial joints, this tissue separates the synovial cavity and its fluid from adjacent tissues and, via the synovial fluid, facilitates the nourishment of chondrocytes.²⁵ The histopathological architecture of the synovium in OA patients is characterized by hyperplasia of the synovial lining, fibrosis in the sublining, and vascularization of the stroma.²⁶ A considerable influx of leukocytes from the vascular compartment, orchestrated by cytokines and cell adhesion molecules, has been observed, with various studies pinpointing macrophages and T-cell lymphocytes as the predominant immune cells within the OA synovium.²⁷

While conventionally perceived as predominantly a malady of hyaline cartilage with concomitant bone implications attributed to overload or excessive utilization, the pathophysiological construct underlying the development of OA is now understood to be highly complex. An expanding corpus of evidence suggests that synovitis and the subsequent proinflammatory mediators are important in the pathogenesis of OA, exerting effects on articular cartilage.²⁸ Macrophage pyroptosis instigates the induction of IL-1β, IL-18, and HMGB1, potentially facilitating the translocation of HMGB1 from the fibroblast nucleus and its subsequent engagement with RAGE, thus activating the TGF-β1/SMAD3 signalling cascade and influencing synovial fibrosis.²⁹ Hypoxia in the subchondral bone constitutes a crucial element in the interaction between chondrocytes and osteoblasts, promoting the transition of chondrocytes toward an OA-like phenotype, likely through the activation of the Wnt/β-catenin signalling pathway in osteoblasts.³⁰ Synovitis manifests early in the disease trajectory of OA, a stage where radiographic evidence of cartilage damage remains minimal. Consequently, it is vital to identify biomarkers associated with altered synovial tissue.

In the present investigation, we endeavoured to establish an OA diagnostic model utilizing machine learning, which was founded on 35 OA patients and 32 healthy controls from four diverse clinical cohorts. The batch effects were initially eliminated to maintain a uniform level of the gene expression matrix. To construct a holistic model for diagnosing OA, we initially computed the enrichment scores of signalling pathways for all samples, followed by the identification of overlapping pathways impacting OA across all cohorts. Eventually, 16 signalling pathways, including the cellular respiration chain, protein transport, lysosomal and endocytosis, bone remodelling, and skeletal muscle regeneration, were selected. Upon eliminating redundant genes, a total of 298 genes were incorporated into the analysis. We considered 101 types of algorithms for the construction of the diagnostic model and ultimately successfully outputted the results of 58 combined approaches. We selected the model constructed via the glmBoost plus RF algorithm, which employs six genes and significantly differentiates OA patients from healthy controls; these genes include LY86, SORL1, KDELR3, CSK, PTGS1, and PTGS2.

LY86, also known as Myeloid Differentiation 1, is a secreted glycoprotein that forms a complex with radioprotective 105 that is expressed on DCs, B cells and macrophages and plays an important role as a regulator of the TLR4 pathway.³¹ The SORL1 gene, located on chromosome 11q23.2-q24.2, encodes a 250 kDa protein dubbed SorLA, a functional sorting receptor for the amyloid-β precursor protein. An increasing body of evidence underscores that variants in the SORL1 gene are risk factors linked with Alzheimer’s disease³² and act as biomarkers of colorectal cancer³³ and ankylosing spondylitis.³⁴ Upon binding to the heterotrimeric signalling G protein G(q/11), the KDEL (Lys-Asp-Glu-Leu) receptor family instigates transport via the Golgi complex, subsequently activating the proteins.³⁵ CSK is the gene encoding the C-Src tyrosine kinase, one of the key kinases controlling cartilage and bone tissue remodelling.³⁶ The development of inflammation during OA is linked with C-Src tyrosine kinase, resulting in elevated production of IL-1β and other inflammatory modulators that expedite the catabolism of cartilage and bone tissue.^37,38 Prostaglandin-endoperoxide synthase (PTGS), alternatively termed cyclooxygenase (COX), is the principal enzyme facilitating prostaglandin biosynthesis. Two isozymes of COX, encoded by different gene products, exist, a constitutive COX-1 and an inducible COX-2, which exhibit divergence in their regulatory expression and tissue distribution. COX-1 and COX-2 constitute vital targets through which nonsteroidal anti-inflammatory drugs exert their effects and play crucial roles in the clinical treatment of OA.³⁹ The knockout of COX-2 in osteocytes or the administration of a COX-2 inhibitor efficaciously salvages the structure of subchondral bone and mitigates cartilage degradation in a spontaneous OA mouse model.⁴⁰

Along with clinical samples from patients who underwent surgery at our institute, we validated the increased protein level of KDELR3 in OA patients compared with healthy controls. KDELR3 has exhibited significant aberrations in an array of diseases, including type 2 diabetes, uveal melanoma,⁴¹ and renal cell carcinoma.⁴² The endoplasmic reticulum (ER) lumen is crucial for lipid metabolism, calcium storage, and protein folding and trafficking. Luminal ER-resident proteins perform these functions and are retained in the ER under homeostasis by carboxy-terminal ER retention sequences (ERSs), which are thought to interact with KDEL receptors. Located in the cis-Golgi, KDEL receptors retrogradely transport ERS-containing proteins back to the ER via COPI vesicle-mediated transport.⁴³ KDELR2 and KDELR3, but not KDELR1, are upregulated during ER calcium depletion as part of the unfolded protein response (UPR). This adaptation aims to retain ERS-containing proteins and maintain ER homeostasis, indicating isoform-specific roles of KDELR3 in managing ER stress and protein retention.⁴⁴ The E96V mutation in the Kdelr3 gene influences beta cell function and contributes to type 2 diabetes (T2D) susceptibility, which disrupts insulin granule maturation and proinsulin levels in beta cells, linking KDELR3 to altered pancreatic beta cell function and increased ER stress in T2D.⁴⁵ In addition, KDELR3 plays a role in the development of extreme myopia (EM) through its involvement in scleral fibroblast function and extracellular matrix organization, and mutations in KDELR3 were found to be associated with elongated ocular axial length.⁴⁶

The findings from our study have potential clinical implications, particularly in the context of early diagnosis and targeted therapeutic strategies for osteoarthritis (OA). The machine learning-based diagnostic model, which leverages the expression profiles of six genes, represents a novel approach for distinguishing OA patients from healthy individuals. Given the early onset of synovitis in OA, the model’s ability to detect subtle changes in synovial tissue could allow for earlier intervention, potentially slowing disease progression before significant cartilage degradation occurs. Moreover, the high level of the identified biomarkers, such as KDELR3, was validated in OA patients via IHC staining of real-world samples. These findings could inform the development of personalized medicine approaches, ultimately improving patient outcomes.

However, the study has several limitations. First, the relatively small sample size of 35 OA patients and 32 healthy controls may limit the generalizability of our findings to broader populations. Future studies should validate the model in larger and more diverse cohorts to ensure robustness and reproducibility. Additionally, while we identified key pathways and biomarkers, functional validation through in vitro or in vivo studies is needed to confirm their roles in the pathogenesis of OA.

Conclusion

We constructed and validated a diagnostic model for OA patients via machine learning, and the six-gene diagnostic model can accurately identify patients with osteoarthritis. KDELR3 has emerged as a novel biomarker for OA, demonstrating promising prospects for clinical application.