Introduction

Osteoarthritis (OA) is a multifaceted degenerative disease that involves both inflammatory and metabolic factors. It is one of the common joint diseases that mostly affects the knees, hands, feet, hips and spine.1 OA is more prevalent in women, the elderly population and patients with obesity. According to the recent statistical data reported by the World Health Organization (WHO), the prevalence of OA has increased from 12.3% to 21.6%.2,3 The quality of life and motor functions of patients with end-stage OA are substantially deteriorated, which is a leading cause of disability and an important health concern worldwide.4 Moreover, early diagnosis of OA remains a major challenge for clinicians.5 Therefore, early and personalized detection, diagnosis, intervention and treatment are particularly crucial.

With the development of modern medicine, predictive, preventive and personalized medicine (PPPM) has exhibited a thriving trend for future medical strategies. PPPM aims to continuously integrate the development trends of contemporary biomedical technology to provide novel insights into the early detection, diagnosis and treatment of diseases.6,7 Unlike previous medical models, PPPM is committed to taking action before or in the early stages of a disease to suppress its occurrence and development.8 At present, PPPM is used for the monitoring, prevention and early diagnosis of various diseases, especially some degenerative and chronic diseases.7

Although the concepts of primary, secondary and tertiary prevention apply to the management of OA, the application of the PPPM framework is significantly limited in OA compared with cancer, diabetes, cardiovascular and cerebrovascular diseases.9 Studies have continuously attempted to identify diagnostic biomarkers for OA, and some significant findings have been reported.10,11 However, there may not exist a perfect single marker.6 Therefore, within the framework of PPPM, multi-omics and comprehensive diagnostic strategies are key to personalized and accurate treatment of OA, which can improve the prognosis of OA from the following aspects:7

-Distinguishing the inflammatory status of patients with OA;

-Providing personalized therapeutic strategies for immune inhibition in patients with OA.

Programmed cell death (PCD) is the process by which cells stop growing and dividing and enter a phase that culminates in their controlled death without the leakage of their contents into the surrounding environment. PCD maintains the balance between pro-inflammatory and anti-inflammatory pathways during inflammatory responses.12 The pathogenesis and progression of OA have been associated with several PCD mechanisms, including apoptosis, autophagy, necroptosis, ferroptosis, cuproptosis, NETosis and pyroptosis.13 Some PCD mechanisms (ferroptosis, pyroptosis, cuproptosis and NETosis) play a crucial role in promoting OA, whereas others (autophagy and apoptosis) inhibit the development of OA.11,14–16 Overall, the roles of PCD mechanisms in OA are diverse. Early promotion of anti-inflammatory PCD and inhibition of pro-inflammatory PCD can improve the prognosis of OA.

Disulfidptosis is a new type of PCD discovered by Liu et al, and its role in OA remains elusive.17 At the intracellular level, the primary mechanisms of disulfidptosis include the abnormal formation of disulphide bonds in actin cytoskeletal proteins and the ingestion of homocysteine.17 The actin cytoskeleton controls the phenotypic stability of chondrocytes to a certain extent, and the cartilage phenotype maintains a fine-tuned balance between anabolism and catabolism. The collapse of the actin network and the loss of cartilage phenotypic stability contribute to the destruction of the cartilage matrix of chondrocytes, resulting in irreversible chondrocyte damage and eventually inducing OA.18,19 The ingestion of abundant cystine may promote the accumulation of glutathione to inhibit inflammation in OA.20,21 Therefore, disulfidptosis is a novel research direction, and elucidating its mechanisms may help to understand the pathogenesis of OA and its therapeutic targets.

Disulfidptosis plays a crucial role in various diseases.17,22 However, a comprehensive investigation of disulfidptosis is lacking. With the rapid development of sequencing technology and intelligent medicine, bioinformatics and machine learning have facilitated the identification of biomarkers for early and accurate diagnosis of OA, thus enabling early prevention, accurate diagnosis and personalized treatment of OA within the framework of PPPM.6,23,24 In this study, training and validation datasets related to OA were downloaded from the Gene Expression Omnibus (GEO) database, and disulfidptosis regulators were obtained from a study by Liu et al.17 Significant disulfidptosis regulators were screened using various machine learning algorithms and via differential expression analysis. Subsequently, we examined the role of these regulators in distinguishing patients with early-stage OA from those with end-stage OA. Based on the regulators, different disulfidptosis-related subtypes of OA were generated via unsupervised clustering, and a diagnostic model was constructed for accurate and personalized diagnosis and treatment of patients with OA within the framework of PPPM. Additionally, the representative characteristics of different subtypes were examined, including immune infiltration and the expression of immune checkpoints. The expression of the significant regulators was verified through single-cell RNA-seq profile and clinical synovial samples. Figure 1 illustrates the general function of disulfidptosis in OA and the step-by-step protocol of this study.

Figure 1 General status of disulfidptosis in OA and step-by-step protocol of this study.

Materials and Methods Data Inclusion and Acquisition

Several datasets containing the data of patients with OA were extracted from GEO (https://www.ncbi.nlm.nih.gov/geo/), a well-known large-scale multi-omic public database.25 The inclusion criteria for datasets were as follows: (1) selection of Homo sapiens as the organism, (2) use of synovial joints as tissue samples, and (3) inclusion of corresponding normal (NM) samples with a comprehensive introduction of each sample. Detailed information regarding the GEO datasets included in this study was provided in Table S1. The following disulfidptosis regulators were obtained from a study by Liu et al: GYS1, LRPPRC, NCKAP1, NDUFA11, NDUFS1, NUBPL1, OXSM, RPN, SLC3A2 and SLC7A11.17 A co-expression network of these 10 regulators was constructed using the GeneMANIA website (http://genemania.org/).26

Integration and Construction of the Training Set

Due to the three data sets (GSE55235, GSE55457, and GSE55584) were contributed by the same authors (Woetzel et.al), we selected them to construct our training set. First, a package in R software, “inSilicoMerging”, was utilized to merge the previous three data sets.27 Subsequently, the “limma” (version 3.42.2) package in R software was used to remove the batch effect existing in the above three data sets.28 As a result, our training set was obtained with 20 NM and 26 OA samples.

Differential Expression Analysis in the Training Set

To investigate the differential expression of genes between NM and OA at the mRNA level, the “limma” package in R was used to identify differentially expressed genes (DEGs) in the integrated training set.28 A p-value of <0.05 was defined as the threshold for screening DEGs. The results of differential expression analysis were visualized on a volcano plot constructed using the “ggplot2” package in R.

Identification of Significant Disulfidptosis Regulators in OA

The DEGs and the 10 disulfidptosis regulators were intersected to ensure the primary regulators in OA. Several machine learning (ML) methods were used to assess the diagnostic efficiency of the regulators. Random forest (RF) is used to evaluate the significance of each feature in classification problems.29,30 In this study, features with importance values of >1 were considered significant. Least absolute shrinkage and selection operator-Cox (Lasso-Cox) regression is an algorithm used for evaluating the influence of binary outcomes via integration of some variables using the R package “glmnet”.31 Significant regulators were identified by intersecting the results of RF and Lasso-Cox regression.

Identification of Different Disulfidptosis-Related Subtypes of OA

All the 26 OA samples in the training set were divided into different disulfidptosis-related subtypes with the corresponding expression value of the above significant regulators in each sample using the ‘ConsensusClusterPlus’ package in R.32 Based on the most appropriate K-value, two subtypes were identified.

The Uniform Manifold Approximation and Projection (UMAP) algorithm was used to examine the heterogeneity between the two subtypes.33 The disulfidptosis score was calculated based on the results of the UMAP algorithm to examine the biological characteristics of the two subtypes: Disulfidptosis Score = Σi (UMAP1i + UMAP2i).

Functional Enrichment Analysis

For gene ontology (GO) enrichment analysis, we utilized the ClueGO plug-in in the Cytoscape software with the threshold (threshold, p < 0.05).34 The Molecular Complex Detection (MCODE) plug-in in Cytoscape was used to extract significant models in the gene-pathway interaction network of the results of ClueGO. The R software package “clusterProfiler” was used to perform enrichment analyses under the background of KEGG and Reactome gene sets.35

In addition, gene set enrichment analysis (GSEA) was performed to comprehensively examine the biological characteristics of the two disulfidptosis-related subtypes of OA with the standard Hallmark gene set, respectively.36

Construction of a Diagnostic Model Based on the Disulfidptosis-Related Subtypes of OA

Differential expression analysis was performed between the two disulfidptosis-related subtypes to obtain DEGs.10 To select more prominent genes to construct a diagnostic model, the criteria for identifying DEGs included the absolute value of log2 fold change (|log2 FC|) >1 and p-value <0.05.

Weighted gene co-expression network analysis (WGCNA) was performed to screen for significant gene clusters associated with OA.37 The ‘gplots’ package in R was used for hierarchical cluster analysis, and abnormal samples or values were deleted. The gene correlation between samples was calculated using the WGCNA algorithm. Subsequently, the best soft-thresholding power was determined, and a standard non-proportional network was established. The WGCNA model or network was related to the characteristics of external samples. Various modules named by different colours were constructed through hierarchical clustering of genes using a dynamic tree-cutting strategy. The DEGs identified between the two disulfidptosis-related subtypes, the DEGs identified between NM and OA samples and the genes in significant clusters of WGCNA were intersected to identify candidate genes for subsequent analysis.

The Search Tool for the Retrieval of Interacting Genes (STRING, https://cn.string-db.org/) was used to construct a protein-protein interaction (PPI) network of the candidate genes, and the Cytoscape software was used to visualize the network.38,39 The cyttoHubba plug-in in Cytoscape was used to identify hub candidate genes.40 Subsequently, three machine learning methods, namely, SVM-RFE, RF and Lasso-Cox regression, were used to construct a diagnostic model. Finally, multivariate logistic regression analysis was performed to establish the model: Diagnostic Score = Σi (Coefficienti * Expression level of featurei).

Immune Infiltration Analysis

Single sample gene set enrichment analysis (ssGSEA) is a novel type of gene set variation analysis which is a method of unsupervised clustering based on a specific gene set to evaluate the score of each sample. In our study, ssGSEA was used to calculate the infiltration score of 28 different scores in OA. We collected genes from previous studies related to 28 different types of immune cells and predicted the infiltration abundance of immune cells based on these genes (Table S2).41,42 Moreover, based on our expression matrix, we used the “IOBR” package in R software to select the ESTIMATE method to calculate the ESTIMATE Score of each sample in our training set.43,44

Single-Cell Analysis

The Study: Cross-tissue, single-cell stromal atlas identifies shared pathological fibroblast phenotypes in four chronic inflammatory diseases in the public website: Single Cell Portal (https://singlecell.broadinstitute.org/single_cell) was included in this study to examine the expression of significant disulfidptosis regulators in OA. Samples from NM, OA, and rheumatoid arthritis (RA) in the above-mentioned datasets were selected for analysis.

Acquisition of Clinical Synovial Tissue Samples from Patients with OA and Healthy Individuals

To validate the expression of the significant disulfidptosis regulators between NM and OA samples, we collected clinical synovial tissue samples from the Department of Orthopedics in the Second Affiliated Hospital of Nanchang University. For the OA group (n = 6), synovial tissue samples were obtained through removal of the synovial membrane during joint replacement surgery for OA. For the NM group (n = 6), synovial tissue samples were collected from patients with a traumatic meniscus tear undergoing synovial membrane cleaning to expose the surgical field of vision during meniscus repair surgery. The protocol for collecting human tissue samples was approved by the Ethics Committee of the Second Affiliated Hospital of Nanchang University (Jiangxi, China). Written informed consent was obtained from all participants before enrolment.

qRT-PCR

The total RNA of serum was extracted from fresh synovial tissues utilizing Trizol reagent (Ambion, Singapore), following the manufacturer s instructions. The quantity and purity of the total RNA were determined using the Nanodrop®ND1000 (TIANJIN). The total RNA was reverse transcribed into complementary DNA (cDNA) employing the Script cDNA synthesis kit (TianGEN). According to the manufacturer’s protocol. Quantitative PCR reaction was performed on an EDC-810 Real-Time PCR Detection Instrument (Eastwin life, China), using SYBR Premix ExTaq kit (Takara, China), All genes were analyzed Relatively and calibrated to the expression of control groups. The gene expression levels were normalized to GAPDH and the 2−ΔΔCt method was used to calculate the definite RNA expression values. Primer pairs are shown in Table S3.

Statistical Analysis

The R software (version 4.2.2) and its related software packages were used to process and analyse data. To estimate differences between two groups, an unpaired t-test was used, and a one-way ANOVA was used to estimate differences among more than two groups. All correlation coefficients were calculated via Spearman correlation analysis. Results with p-value <0.05 were considered significant. Nomograms were visualized using the “RMS” package in R, and receiver operating characteristic (ROC) curves were visualized using the Sangerbox website (http://vip.sangerbox.com/home.html).

Results Exploration of the General Status of Disulfidptosis Regulators in OA

Firstly, we exhibited the remarkable landscape of disulfidptosis in OA. Based on the expression profiles of nine disulfidptosis regulators obtained from a study by Liu et al in the integrated dataset, Spearman correlation analysis indicated a remarkably positive correlation among the regulators17 (Figure 2A). A co-expression network of the nine regulators was constructed using the public website GeneMANIA, with physical interactions of 55.32%, pathway of 30.23%, co-expression of 12.63% and co-localization of 1.82%. Enrichment analysis revealed that nine regulators were mainly associated with the modification, transportation and synthesis of amino acids, including “amino acid transmembrane transporter activity”, “carboxylic acid transmembrane transport” and “peptidyl-asparagine modification” (Figure 2B).

Figure 2 Identification of significant disulfidptosis regulators in OA. (A) Heat map demonstrating the correlation among 9 disulfidptosis regulators in OA samples. (B) Co-expression network of the 9 disulfidptosis regulators and their potential biological functions. (C) Volcano map demonstrating DEGs with p-value of <0.05 between NM and OA samples. (D) The intersection between the DEGs and the 9 disulfidptosis regulators. (E, F) Results of Lasso-Cox regression analysis. (G, H) Results of RF (screening threshold: importance > 1). (I) Intersection between the results of Lasso-Cox regression and RF. (J-l) ROC curves of (J) NCKAP1, (K) OXSM and (l) SLC3A2. (M) Nomogram demonstrating the efficacy of the three regulators in distinguishing early- and end-stage patients with OA. (N) Calibration curve of the nomogram (*, p < 0.05; **, p < 0.01; ***, p < 0.001).

To identify prominent disulfidptosis regulators in OA, differential expression analysis was performed between the NM and OA groups. Genes with p-value of <0.05 were considered as DEGs (Figure 2C). Intersection of these DEGs with the nine disulfidptosis regulators revealed three prominent regulators, namely, NCKAP1, OXSM and SLC3A2 (Figure 2D).

Assessment of the Diagnostic Efficiency of the Three Regulators in OA via ML

The diagnostic efficiency of the three disulfidptosis regulators was assessed using ML methods. Lasso-Cox regression analysis, with 10-fold cross-validation, was performed to verify the diagnostic efficiency of the three regulators (Figures 2E-F). Furthermore, RF was used to assess the importance of the three regulators in distinguishing the NM and OA samples. With the importance of >1 as the criteria, all three regulators were selected (Figures 2G-H). After intersecting the results of the two ML methods, the three regulators (NCKAP1, OXSM and SLC3A2) were identified as significant disulfidptosis regulators (Figure 2I).

ROC curves were plotted to demonstrate the diagnostic efficiency of the three regulators. The AUC values of NCKAP1, OXSM and SLC3A2 were 0.58, 0.67 and 0.64, respectively, which indicated their strong diagnostic efficiency (Figures 2J-l). In addition, based on the corresponding expression of the three regulators in the each sample of GSE32317, a nomogram was constructed to assess the ability of the regulators to identify patients with early- and end-stage OA (Figure 2M). As shown in Figure 2N, the calibration curve validated the accuracy of the above nomogram.

Identification of Two Disulfidptosis-Related Subtypes of OA

With the background of the expression value of the three regulators (NCKAP1, OXSM and SLC3A2) in OA samples, all the 26 OA samples in the integrated training dataset were divided into two disulfidptosis-related subtypes (C1 and C2) (Figures 3A-C). As exhibited in Figure 3D, the UMAP plot demonstrated differences between the two subtypes of OA. In addition, the three disulfidptosis regulators were significantly differentially expressed between the two subtypes, which indicated the heterogeneity between the subtypes (Figures 3E-H). The UMAP algorithm was used to calculate the disulfidptosis score for objective recognition of the two subtypes. The score of the C1 subtype was higher than that of the subtype (p < 0.05, Figure 3I).

Figure 3 Identification of two disulfidptosis-related subtypes of OA and establishment of a disulfidptosis scoring model for OA. (A) A cumulative distribution curve was plotted to determine the most appropriate K value. (B) The area under the cumulative distribution curve in each K value. (C) Heat map of two disulfidptosis-related subtypes of OA. (D) UMAP plot of the two subtypes. (E) Heat map demonstrating the normalised expression level in the two subtypes. (F) Differences in the disulfidptosis score between the two subtypes. (G–I) Differential expression of (G) NCKAP1, (H) OXSM and (I) SLC3A2 between the two subtypes (*, p < 0.05; ****, p < 0.0001).

Biological Characteristics of the Two Disulfidptosis-Related Subtypes

The above-mentioned results suggest a significant difference between the two disulfidptosis-related subtypes of OA. DEGs were identified by differential expression analysis between the two subtypes (C2 versus C1). A total of 10,172 DEGs were identified based on the screening criteria specified in the above (Figure 4A). All DEGs were down-regulated in the C2 subtype. Similarly, a total of 182 DEGs were identified between NM and OA groups, including 48 up-regulated and 134 down-regulated genes (Figure 4B). Moreover, to select some candidate features with potential diagnostic efficiency, we intersected the DEGs between two disulfidptosis-related subtypes and the DEGs between NM and OA group. As a result, 121 candidate genes with potential diagnostic efficiency were selected for constructing a diagnostic model (Figure 4C). To limit the number of genes, the WGCNA algorithm was used to screen several gene modules, which were significantly correlated with OA. All modules (p < 0.05) negatively or positively correlated with OA were included, in which 782 genes were positively correlated and 692 genes were negatively correlated with OA (Figures 4D-F). Subsequently, the 121 candidate genes were intersected with the genes in the significant modules of WGCNA to obtain a total of 62 genes for subsequent establishment of the diagnostic model (Figure 4G).

Figure 4 Identification of significant diagnostic genes in the two OA subtypes. (A, B) Volcano map of DEGs between (A) the two disulfidptosis-related subtypes of OA and (B) between the NM and OA groups. (C) Venn diagram demonstrating the intersection between the two DEG sets. (D, E) Selection of the soft-thresholding power for WGCNA. (F) Division of modules to identify genes correlated with OA. (G) Venn diagram demonstrating the intersection among the overlapped DEGs, genes in the negatively correlated modules and genes in the positively correlated modules.

To elucidate the biological characteristics of the two subtypes, the 62 candidate genes were subjected to GO and KEGG enrichment analyses. GO enrichment analysis revealed that the genes were enriched in pathways related to immunity, metabolism, cell proliferation and differentiation, especially “chemokine production” (45%) (Figure S1a). The functional clusters identified using the MCODE plug-in were mainly related to immune responses (Figure S1b-m). The results of KEGG and Reactome enrichment analyses demonstrated the crucial role of the immune system in the two disulfidptosis-related subtypes (Figure S2a-b). In addition, the results of GSEA emphasised the role of the candidate genes in transcription and translation (Figure S2c-h).

Construction of a Diagnostic Model Based on the Two Disulfidptosis-Related Clusters

Given the crucial role of disulfidptosis regulators in OA and the lack of effective diagnostic biomarkers, we constructed a novel diagnostic model for OA based on the candidate genes. The 62 genes were imported into the STRING website to establish a PPI network (Figure 5A). The cyttoHubba plug-in was used to obtain the top 15 genes with the highest interaction score in each of the seven algorithms (MCC, MNC, Degree, EPC, Closeness, Radiality and Stress). After intersecting theses 15 genes in each algorithm, a total of 12 genes identified in all seven algorithms were selected for subsequent analysis (Figure 5B).

Figure 5 Construction and verification of the disulfidptosis-related diagnostic model of OA. (A) PPI network of significant genes in the two OA subtypes. (B) UpSet plot demonstrating the overlapped top 15 genes identified using seven algorithms of cyttoHubba. (C, D) Results of SVM-RFE revealed 9 significant genes. (E, F) Results of RF revealed 7 significant genes. (G, H) Results of Lasso-Cox regression analysis revealed 5 significant genes. (I) Venn diagram demonstrating the intersection of the results of the abovementioned three machine learning methods. (J) Calibration curve of logistic analysis to ensure the accuracy of the diagnostic model. (K, M, O) Differences in the diagnostic score between the NM and OA groups in the (K) training set, (M) GSE12021 dataset and (O) GSE1919 dataset. (L, N, P) ROC curve of the diagnostic score in the (L) training set, (N) GSE12021 dataset and (P) GSE1919 dataset (**, p < 0.01; ****, p < 0.0001).

Furthermore, three ML methods were used to select significant genes. SVM-RFE, a classical ML method, was used to screen for genes with more accurate diagnostic efficiency. A total of 9 genes were selected from the above-mentioned 12 genes (Figures 5C-D). Subsequently, RF was used to evaluate the importance of the 12 genes in recognizing OA and NM samples. As shown in Figures 5E-F, the lollipop chart demonstrates the diagnostic efficiency of the 12 genes in detail. Of the 12 genes, 7 genes with importance values of >1 were considered significant. Furthermore, five genes with more accurate diagnostic efficiency were identified via Lasso-Cox regression analysis. A crucial role may be played by these genes in the diagnosis of OA. (Figures 5G-H). Finally, the results of the three ML methods were intersected, and two genes (GADD45B and VEGFA) were selected for constructing the diagnostic model (Figure 5I).

Multivariate logistic regression analysis was performed to construct the diagnostic model for OA, and the formula to calculate the diagnostic score was as follows: (−0.00276107233300387 × GADD45B) + (−0.00190432657530599 × VEGFA). The calibration curve revealed no difference between the actual and predicted diagnoses (p = 0.338, Figure 5J). As a result, our diagnostic model revealed a statistically significant difference between the NM and OA samples in the training set (p < 0.0001, Figure 5K). The AUC of the model was 0.98 in the training set, which demonstrated the remarkable diagnostic efficiency of the model (Figure 5L). Two external validation sets (GSE12021 and GSE1919) were used to verify the reliability of the model. The diagnostic scores were significantly different between the NM and OA groups in the GSE12021 and GSE1919 datasets (Figures 5M and O). The AUC value was 0.93 and 1.00 in the GSE12021 and GSE1919 datasets, respectively, which indicated the predictive accuracy of the model (Figures 5N and P).

Investigation of the Immune Landscape of the Two Disulfidptosis-Related Subtypes

Both disulfidptosis-related subtypes of OA may be affected by the immune microenvironment. Two algorithms, namely, ESTIMATE and ssGSEA, were used to comprehensively explore the immune microenvironment of OA. The ESTIMATE score was lower in the C1 subtype than in the C2 subtype, indicating that immune responses and immune cell infiltration in the synovial tissue were different between these two disulfidptosis-related subtypes (Figure 6A). Subsequently, ssGSEA was used to visualize the superficial landscape of immune cell infiltration in all samples of the training set, which demonstrated substantial differences in the immune microenvironment between NM samples and the two OA clusters (Figure 6B). And between NM and OA, abundant immune cells exhibited significantly different infiltration status (Figure 6C). Most immune cells had significantly different infiltration scores between C1 and C2, which was consistent with the results of the ESTIMATE algorithm (Figure 6D). Using Spearman correlation analysis, the correlation coefficients of the three significant disulfidptosis regulators were calculated. The correlation between the immune microenvironment and the expression of the three regulators was stronger in C2 subtype than in the C1 subtype, which supported the results of the ESTIMATE algorithm (Figures 6E-J).

Figure 6 Investigation of the immune microenvironment and expression of ICKs in the two disulfidptosis-related subtypes of OA. (A) Differences in the ESTIMATE score between the two subtypes. (B) Stacked column chart demonstrating the infiltration of 28 types of immune cells in each sample of the training set. (C, D) Differences in the infiltration scores of each immune cell between (c) the NM and OA groups and (D) between the two subtypes of OA. (E–G) Lollipop chart demonstrating the correlation between the infiltration scores and the expression of (E) NCKAP1, (F) OXSM and (G) SLC3A2 in the C1 subtype and their corresponding p-values. (H–J) Lollipop chart demonstrating the correlation between the infiltration scores and the expression of (H) NCKAP1, (I) OXSM and (J) SLC3A2 in the C2 subtype and their corresponding p-values. (k–x) Differences in the expression of (K) CD2, (L) CD47, (M) CD96, (N) CD200, (O) CD226, (P) CTLA4, (Q) HHLA2, (R) KIR3DL1, (S) KLRD1, (T) LAG3, (U) PDCD1, (V) PDCD1LG2, (W) SIGIRR and (X) SIGLEC15 between the two subtypes of OA. (Y) Potential relationship among disulfidptosis, immune microenvironment and ICKs in OA. (-, p ≥ 0.05; *, p < 0.05; **, p < 0.01; ***, p < 0.001; ****, p < 0.0001).

The above-mentioned results suggest that C1 and C2 subtypes have different immune microenvironments. In light of the important role immune checkpoints (ICKs) play in regulating immune response, we examined their role in OA. The unpaired t-test was used to compare the expression of 14 classical ICKs between C1 and C2 subtypes. All ICKs except KIR3DL1 had significantly lower expression in the C2 subtype than in the C1 subtype (Figures 6K-X). As shown in Figure 6Y, the Sankey plot demonstrates the potential relationship between the two disulfidptosis-related subtypes and the immune microenvironment.

Overall, these two disulfidptosis-related subtypes of OA with different disulfidptosis score exhibited various inflammatory response and immune microenvironment. Specifically, OA patients with lower disulfidptosis score may experience higher infiltration of immune cells and exhibit more severe inflammatory response. Moreover, the expression of ICKs these OA patients with high disulfidptosis score probably can receive the immunotherapy for the pain and synovial inflammation, which may be regarded as a novel and personalized therapy for OA.

Assessment of the Expression of the Three Disulfidptosis Regulators at the Single-Cell Level and Validation of the Expression Trend in the Clinical Synovial Tissue

To examine the expression patterns of the three disulfidptosis regulators in OA, we selected a dataset including synovial fibroblasts in NM, OA and RA samples. The general landscape of individual synovial fibroblasts in different samples is demonstrated in Figure 7A. The expression of the three disulfidptosis regulators (NCKAP1, OXSM and SLC3A2) in the dataset was consistent with that in the training set, including the up-regulation of OXSM and down-regulation of NCKAP1 and SLC3A2 in OA. In addition, the expression of the three regulators was different between patients with OA and RA, indicating that the regulators are potential biomarkers that can help distinguish patients with OA from those with RA (Figures 7B-D). All three regulators exhibited similar expression patterns in the GSE12021 dataset (Figures 7E-G). Moreover, to further verify the expression trend of the three regulators between NM and OA, we collected clinical synovial tissue from healthy individuals and patients with OA. Then, qRT-PCR was performed to evaluate the relative concentration of mRNA. As a result, except for OXSM, other two regulators (NCKAP1 and SLC3A2) exhibited a significant down-regulation in OA, which was consistent with our bioinformatic analysis (Figures 7H-J). The contrary result of OXSM may result from the sample size and sample heterogeneity.

Figure 7 Expression of disulfidptosis regulators analysed via single-cell RNA sequencing. (A) General distribution of each cell in the NM, OA and RA groups. (B–D) Expression of (B) NCKAP1, (C) OXSM and (D) SLC3A2 as analysed via single-cell RNA-seq. (E–G) Expression of (E) NCKAP1, (F) OXSM and (G) SLC3A2 at the transcriptomic level in the NM, OA and RA groups. (H-J) Expression of (H) NCKAP1, (I) OXSM and (J) SLC3A2 at the mRNA level in the NM and OA patients in our in-house cohort via RT-qPCR. (-, p ≥ 0.05; *, p < 0.05; **, p < 0.01).

Discussion

In this study, the role of disulfidptosis regulators in the diagnosis, prevention and treatment of OA was comprehensively examined. Three significant disulfidptosis regulators, namely, NCKAP1, OXSM and SLC3A2were identified through differential expression analysis and machine learning. According to the transcriptomic profiles of these three regulators, clustering analysis was performed to investigate the biological characteristics and immune microenvironment of two disulfidptosis-related subtypes of OA and evualate the expression of ICKs in the subtypes. Finally, a disulfidptosis-mediated diagnostic model for OA was constructed, which demonstrated remarkable predictive efficiency.

Current Knowledge of the Relationship Between Disulfidptosis and OA

OA is a type of arthritis that causes joint pain, especially in the knee and hip joints. The pathology of OA is associated with cartilage ageing or overuse leading to the loss of chondrocyte homeostasis and death, resulting in friction between bones.45,46 The pathogenesis of OA is related to complex factors, including genetic susceptibility, ageing, inflammation, obesity and uneven joint arrangement.47 The relationship between OA and inflammation has become a major focus of research in recent years, and several studies have implicated inflammation in OA pathogenesis. In addition, anti-inflammatory strategies have become the primary means of treating OA.48,49 Disulfidptosis is a novel type of PCD. Disulfide stress occurs in glucose-deficient cells with high expression of SLC7A11, inducing abnormal disulfide cross-linking between actin cytoskeletal proteins and cytoskeletal contraction, eventually leading to actin network collapse and PCD.17 Patients with OA have lower levels of cytoskeletal proteins in their cartilage compared to healthy individuals.,50 indicating that the deterioration of OA is related to the destruction of actin. However, during disulfidptosis, the ingestion of abundant cystine may induce an increase in the concentration of intracellular glutathione and NADPH, which can lead to a decrease in the levels of reactive oxygen species (ROS).17 Consequently, the decreased levels of ROS play a crucial role in anti-inflammatory pathways, which can decrease the severity of OA and inhibits its progression.51,52 Maintaining pro- and anti-inflammatory responses in the internal environment of patients is particularly crucial for promoting early relief and inhibiting the progression of OA.53,54 Therefore, elucidating the relationship between OA and disulfidptosis is important for implementing PPPM in the diagnosis and treatment of patients with OA.

Potential Effects of the Findings of This Study on the Current PPPM Framework

It is promising and necessary to discuss the diagnosis and personalized treatment of OA based on the framework of PPPM. OA is a degenerative and multi-stage disease that should be effectively prevented and promptly treated using strategies based on PPPM.

First Stage of PPPM

The first stage of PPPM is vital for predicting patients with a high risk of OA and effectively alleviating the risk factors of OA, including age (>60 years), sex (female), obesity (BMI > 25), unbalanced and unhealthy diet and joint-related factors (joint injury and abnormal load).2,55 In addition, regular and healthy exercise and early prevention of pro-inflammatory factors are essential approaches.

Second Stage of PPPM

According to bone and joint experts worldwide, in order to improve the prognosis of patients with OA, early and timely diagnosis is crucial.56 Higher immune infiltration and stronger inflammatory responses have been significantly associated with the poor prognosis of patients with OA.55 Although predictive models have been widely used, their application in early diagnosis and determining immune/inflammatory cell infiltration in OA remains limited.10,57 Therefore, effective predictive strategies for distinguishing various inflammatory stages and for early diagnosis of OA should be developed. In this study, two disulfidptosis-related subtypes of OA were identified via consensus clustering. The disulfidptosis score was higher in the C1 subtype than in the C2 subtype, indicating that patients with the C1 subtype had more severe disulfidptosis with higher expression of the three regulators. Furthermore, there was a significant difference for the infiltration of the immune cells in the synovial tissue between these two disulfidptosis-related subtypes. Specifically, the infiltration score of Type 17 T helper cells were higher in C2, which emphasized that the inflammation response were more severe in the synovial tissue of OA patients with lower disulfidptosis score since the Type 17 T helper cells are regarded as the immune cell cluster playing a crucial role in inducing the inflammation response via releasing several cytokines like IL17, IL23, and so on. While other anti-inflammatory immune cells, including activated CD8 T cell, effector memory CD4 T cell, memory B cell, and monocyte, exhibited a higher infiltration in C1. The findings of this study suggest that disulfidptosis, which is considered as an anti-inflammatory type of PCD, is more inclined toward suppressing the inflammatory response of OA. On the other hand, the nomogram of GSE32317 and the disulfidptosis-mediated diagnostic model also emphasized the potential characteristics of disulfidptosis in early diagnosis of OA. Overall, within the framework of PPPM, we comprehensively investigate the role of disulfidptosis in OA, providing different management models for different subtypes, which may facilitate precise treatment of OA.

Third Stage of PPPM

After experiencing the pre-clinical stages, patients with end-stage OA or OA complicated with significant symptoms exhibit further deterioration of their physical condition and exercise ability. Given the anti-inflammatory role of disulfidptosis in OA and inflammatory reactions throughout OA, potential interventions based on the disulfidptosis regulators identified in this study may be beneficial for treatment.58

Comprehensive Investigation of the Disulfidptosis Regulators in OA SLC3A2

Member 2 of solute carrier series 3 (SLC3A2), also known as 4F2HC, encodes a multifunctional type II membrane glycoprotein involved in cell adhesion, fusion, transformation and amino acid transportation.59 The SLC3A2 protein, a chaperone protein of SLC7A11, forms a heterodimer with it, participating in the formation of the cysteine/glutamate transport system (Xc-system) and mediating cysteine uptake.60 Low expression of SLC3A2 may contribute to the development of inflammatory diseases.17 Relevant studies have demonstrated that SLC3A2 can inhibit cartilage degeneration in OA and its expression in the damaged area of the cartilage is down-regulated in OA. In addition, knockdown of SLC3A2 can promote OA-related cartilage degeneration.61 These findings are consistent with those of the present study, suggesting that the expression of SLC3A2 is low in patients with OA. The role of SLC3A2 in OA is the first indication that the pathogenesis of OA is heavily dependent on disulfidptosis.

NCKAP1

NCK-related protein 1 (NCKAP1), also known as Nap1, was first found to be significantly inhibited in the brain tissue of patients with sporadic Alzheimer’s disease (AD).62 The protein encoded by NCKAP1 and the WAVE, Sra, Abi and HSPC300 (also known as Brick1) proteins participate in the formation of the WAVE regulatory complex (WRC).63 NCKAP1 is the fourth on the list of genes required for disulfidptosis and can promote actin polymerisation and lamellipodia formation. Knock out of NCKAP1 inhibits disulfidptosis.17 In recent years, NCKAP1 has been widely used for tumor inhibition and as a novel prognostic biomarker for cancer. It can promote or inhibit tumor development in some types of cancers.64,65 As far as we know, this study is the first to demonstrate a relationship between NCKAP1 and OA.

OXSM

Mitochondrial 3-oxyl-ACP synthase (OXSM) is an enzyme required for the elongation of fatty acid chains in mitochondria, and its synthetic product is lipoic acid (LA),66 which is an important synthase involved in the mitochondrial fatty acid synthesis pathway II. Studies on the regulation of OXSM expression are limited. OXSM is a target of some genes and is involved in the regulation of ovarian cancer and glioma cells.66,67 The relationship between OXSM and OA warrants further investigation.

Feedback from This Study

Under the background of PPPM, integrating multi-stage and multi-omic analyses is key to diagnosing and treating OA. In this study, ML, bioinformatic tools and public transcriptomic data were predominantly used, whereas classical indices or biomarkers related to imaging omics, metabolomics and pathology were not included. Although two validation datasets were used, the sample size was relatively limited. Moreover, owing to ethical limitations, we could not obtain completely healthy human synovial tissue samples. More comprehensive in vitro and in vivo experiments should be conducted to validate the findings of this study. Disulfidptosis is a novel concept, and relevant genes associated with it should be further identified and investigated to expand the gene set.

Conclusion

Via integrating three independent datasets and utilizing a series of bioinformatics analysis and machine learning, this study demonstrated the anti-inflammatory role of disulfidptosis in OA and finally identified three crucial regulators of disulfidptosis in OA, including NCKAP1, OXSM, and SLC3A2. Moreover, the disulfidptosis score plays a crucial role in predicting the inflammatory responses and immune infiltration of OA, using which may predict the disease severity of OA and provide the personalized therapeutic strategies for OA patients. In addition, the disulfidptosis-mediated diagnostic model exhibited potential predictive efficiency. Furthermore, the transcriptomic level and single-cell level were merged to emphasize the crucial role of disulfidptosis, which was also validated in our in-house cohort. Altogether, disulfidptosis may be a new potential target for the diagnosis, stratification, targeted prevention and personalized treatment of OA in the future, which will promote the development of the PPPM framework for OA.

Data Sharing Statement

The original data used in this project can be downloaded in the public database GEO (https://www.ncbi.nlm.nih.gov/geo/). The data analyzed during the human experiments are available from the corresponding author on reasonable request.

Ethics Approval and Informed Consent

The protocol for collecting human tissue samples was approved by Ethics Committee of the the Second Affiliated Hospital of Nanchang University (No. Review [2022] No. (031) and 23 March 2022), and all procedures were performed in accordance with the 1964 Helsinki Declaration and its later amendments or comparable ethical standards. Written informed consent was provided by all participants before their inclusion in the study.

Acknowledgments

We sincerely acknowledge the contributions from the public platform Sangerbox (http://vip.sangerbox.com/home.html) accessed. And we thank Bullet Edits Limited for the linguistic editing and proofreading of the manuscript.

Author Contributions

All authors made a significant contribution to the work reported, whether that is in the conception, study design, execution, acquisition of data, analysis and interpretation, or in all these areas; took part in drafting, revising or critically reviewing the article; gave final approval of the version to be published; have agreed on the journal to which the article has been submitted; and agree to be accountable for all aspects of the work.

Funding

This research was funded by [Science and Technology project of Jiangxi Provincial Administration of Traditional Chinese Medicine, grant number 2023Z030], [the National College Students’ Innovative Entrepreneurial Training Plan Program, grant number 202210403001], and [the College Students’ Innovative Entrepreneurial Training Plan Program in Jiangxi Province, grant number S202210403002].

Disclosure

The authors have no relevant financial or non-financial interests to disclose.

References

1. Abramoff B, Caldera FE. Osteoarthritis: pathology, Diagnosis, and Treatment Options. Med Clin North Am. 2020;104(2):293–311. doi:10.1016/j.mcna.2019.10.007

2. Yao Q, Wu X, Tao C, et al. Osteoarthritis: pathogenic signaling pathways and therapeutic targets. Signal Transduct Target Ther. 2023;8:56.

3. Bernabei I, So A, Busso N, Nasi S. Cartilage calcification in osteoarthritis: mechanisms and clinical relevance. Nat Rev Rheumatol. 2023;19:10–27.

4. Hartnett DA, Milner JD, DeFroda SF. Osteoarthritis in the Upper Extremity. Am J Med. 2023;136(5):415–421. doi:10.1016/j.amjmed.2023.01.025

5. Mahmoudian A, Lohmander LS, Mobasheri A, Englund M, Luyten FP. Early-stage symptomatic osteoarthritis of the knee - time for action. Nat Rev Rheumatol. 2021;17:621–632.

6. Golubnitschaja O, Baban B, Boniolo G, et al. Medicine in the early twenty-first century: paradigm and anticipation - EPMA position paper 2016. EPMA J. 2016;7(1):23. doi:10.1186/s13167-016-0072-4

7. Golubnitschaja O, Liskova A, Koklesova L, et al. Caution, “normal” BMI: health risks associated with potentially masked individual underweight-EPMA Position Paper 2021. EPMA J. 2021;12:243–264. doi:10.1007/s13167-021-00251-4

8. Golubnitschaja O, Kinkorova J, Costigliola V. Predictive, Preventive and Personalised Medicine as the hardcore of ‘Horizon 2020’: EPMA position paper. EPMA J. 2014;5(1):6. doi:10.1186/1878-5085-5-6

9. Hewett TE, Myer GD, Ford KR, Paterno MV, Quatman CE. Mechanisms, prediction, and prevention of ACL injuries: cut risk with three sharpened and validated tools. J Orthop Res. 2016;34(11):1843–1855. doi:10.1002/jor.23414

10. Hao L, Shang X, Wu Y, Chen J, Chen S. Construction of a Diagnostic m(7)G Regulator-Mediated Scoring Model for Identifying the Characteristics and Immune Landscapes of Osteoarthritis. Biomolecules. 2023;14(1):13. doi:10.3390/biom14010013

11. Chang B, Hu Z, Chen L, Jin Z, Yang Y. Development and validation of cuproptosis-related genes in synovitis during osteoarthritis progress. Front Immunol. 2023;14:1090596. doi:10.3389/fimmu.2023.1090596

12. D’Arcy MS. Cell death: a review of the major forms of apoptosis, necrosis and autophagy. Cell Biol Int. 2019;43(6):582–592. doi:10.1002/cbin.11137

13. Liu S, Pan Y, Li T, et al. The Role of Regulated Programmed Cell Death in Osteoarthritis: from Pathogenesis to Therapy. Int J Mol Sci. 2023;25(1):24. doi:10.3390/ijms25010024

14. Medina CB, Mehrotra P, Arandjelovic S, et al. Metabolites released from apoptotic cells act as tissue messengers. Nature. 2020;580:130–135.

15. Xue JF, Shi ZM, Zou J, Li XL. Inhibition of PI3K/AKT/mTOR signaling pathway promotes autophagy of articular chondrocytes and attenuates inflammatory response in rats with osteoarthritis. Biomed Pharmacother. 2017;89:1252–1261. doi:10.1016/j.biopha.2017.01.130

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

17. Liu X, Nie L, Zhang Y, et al. Actin cytoskeleton vulnerability to disulfide stress mediates disulfidptosis. Nat Cell Biol. 2023;25:404–414.

18. Lauer JC, Selig M, Hart ML, Kurz B, Rolauffs B. Articular Chondrocyte Phenotype Regulation through the Cytoskeleton and the Signaling Processes That Originate from or Converge on the Cytoskeleton: towards a Novel Understanding of the Intersection between Actin Dynamics and Chondrogenic Function. Int J Mol Sci. 2021;23(1):22. doi:10.3390/ijms23010022

19. Leipzig ND, Eleswarapu SV, Athanasiou KA. The effects of TGF-β1 and IGF-I on the biomechanics and cytoskeleton of single chondrocytes. Osteoarthritis Cartilage. 2006;14(12):1227–1236. doi:10.1016/j.joca.2006.05.013

20. Sasaki E, Yamamoto H, Asari T, et al. Metabolomics with severity of radiographic knee osteoarthritis and early phase synovitis in middle-aged women from the Iwaki Health Promotion Project: a cross-sectional study. Arthritis Res Ther. 2022;24(1):145. doi:10.1186/s13075-022-02830-w

21. Liu W, Martinon-Torres M, Cai YJ, et al. The earliest unequivocally modern humans in southern China. Nature. 2015;526(7575):696–699. doi:10.1038/nature15696

22. Machesky LM. Deadly actin collapse by disulfidptosis. Nat Cell Biol. 2023;25(3):375–376. doi:10.1038/s41556-023-01100-4

23. Stranneheim H, Wedell A. Exome and genome sequencing: a revolution for the discovery and diagnosis of monogenic disorders. J Intern Med. 2016;279(1):3–15. doi:10.1111/joim.12399

24. Greener JG, Kandathil SM, Moffat L, Jones DT. A guide to machine learning for biologists. Nat Rev Mol Cell Biol. 2022;23:40–55.

25. Edgar R, Domrachev M, Lash AE. Gene Expression Omnibus: NCBI gene expression and hybridization array data repository. Nucleic Acids Res. 2002;30(1):207–210. doi:10.1093/nar/30.1.207

26. Franz M, Rodriguez H, Lopes C, et al. GeneMANIA update 2018. Nucleic Acids Res. 2018;46(W1):W60–W64. doi:10.1093/nar/gky311

27. Taminau J, Meganck S, Lazar C, et al. Unlocking the potential of publicly available microarray data using inSilicoDb and inSilicoMerging R/Bioconductor packages. BMC Bioinf. 2012;13(1):335. doi:10.1186/1471-2105-13-335

28. Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies