Comprehensive Analysis of Sphingolipid Metabolism-Related Genes in Osteoarthritic Diagnosis and Synovial Immune Dysregulation

Introduction

Osteoarthritis (OA) is a chronic degenerative disease that affects joint function. Its global prevalence rate is estimated to be 6–32%, with the elderly population being most affected [1]. Various risk factors contribute to the development of OA, including aging, excessive physical labor, dyslipidemia, and metabolic syndrome [2,3]. OA is characterized by progressive cartilage degeneration, subchondral sclerosis, synovial inflammation, meniscus degeneration, and osteophyte formation [4–6]. Clinical manifestations of OA include knee dysfunction, stiffness, and pain, which impose physical, psychological, and socio-economic burdens, particularly on the elderly [7]. Previous studies have highlighted the significant role of the synovium in OA development, as joint degeneration occurs following apparent pathological alterations in synovial inflammatory responses [8,9]. Osteoarthritic synovial fibroblasts help support immune cell activation [10] and secrete joint-damaged cytokines such as IL-6, TNF, and MMP3 [11]. The inflammatory response aroused by synovitis triggers local infiltration of macrophages, T cells, and B cells [12], which in turn lead to articular cartilage destruction and accelerated joint degeneration through inflammatory cascades [13,14]. Reducing synovial inflammation can help alleviate joint pain, or even reverse articular structural remodeling [15]. There is an urgent need to explore the underlying cellular mechanism of synovitis and methods of protecting the OA synovium from damaged cytokines in people with OA [16].

There are currently several therapies for OA, including limitation of activity, alleviation of pain, and mitigation of cartilage injury in the early stage, for delaying OA progression [17]. As joint degeneration progresses, many patients eventually receive knee replacement [18]. Thus, adequate diagnosis of OA becomes very important. However, the confirmation of OA depends on imaging data and patient symptoms [19,20], and timely diagnosis of OA is difficult due to lack of specific clinical manifestations. The effective recognition of OA is key to preventing the loss of articular cartilage and OA progression in joints [17, 21], and numerous studies have unveiled early diagnostic biomarkers for OA [22,23], but the clinical utilization of these biomarkers remains a challenge, underscoring the need for further exploration of potential candidate markers for clinical application.

Lipid metabolism disorder has been demonstrated to be important in OA pathogenesis, with significant differences in lipid metabolites observed between the synovial fluid of people with OA versus those without OA [24,25]. Specifically, sphingolipids have been found to be upregulated in OA and play a protective role in inhibiting immune infiltration and attenuating synovitis [26,27]. Therefore, sphingolipid metabolism pathways (SMP) may be involved in regulating synovitis in OA. However, the exact impact of SMP gene expression on synovitis has not been reported, and it remains unclear whether SMP genes can be applied to diagnosis of OA. In this study, we identified 4 hub genes associated with SMP in OA synovium. These hub genes have the potential to serve as diagnostic markers for OA and play a crucial role in immune cell infiltration and inflammatory factors. A deeper understanding of the molecular function of SMP genes in osteoarthritic synovium holds promise for exploring novel diagnostic biomarkers for OA and guiding future research endeavors.

Material and Methods

PREPROCESSING OF RAW DATA SOURCE:

The workflow chart of this study is shown in Figure 1. Gene sets expression (GSE) profiles containing the gene expression information of healthy (CT) and OA synovium (GSE29746, GSE12021, GSE55235, GSE55457) were downloaded from the GEO online database with search term “osteoarthritis AND (synovitis OR synovium)”. Then, the batch effects were removed with the “sva” package [28], and we obtained an integrated dataset including 40 OA and 40 healthy samples. Gene annotation file profiles, which used for converting gene IDs within included transcriptome data into gene symbols, were downloaded from the GPL96 (GSE55235, GSE55457, and GSE12021) and GPL4133 (GSE29746) platforms. The genes of SMP were acquired from Reactome (https://reactome.org/). The demographic information of these patients is shown in Table 1.

IDENTIFICATION OF SMP-RELATED DIFFERENTIAL EXPRESSION GENES (DEGS) IN OA SYNOVIUM:

Using the “limma” package in R software, the differential expression analysis was used to explore the DEGs between OA and healthy synovium in the integrated dataset. Using Venn analysis, the SMP-related DEGs were identified by overlapping the genes of SMP and DEGs. The “heatmap” and “ggplot2” packages in R software were used to draw the heatmap and volcano plots, and the “ggpubr” package was used to create boxplots.

IDENTIFICATION OF SMP-RELATED HUB GENES IN OA SYNOVIUM:

Using the “randomForest” package, we conducted random forest classifiers consisting of 1000 decision trees and then validated the 10-fold cross. The SMP-related genes were analyzed with the random forest function in the “randomForest” package and the “forestplot” package was used to generate a forest plot. The significance of each SMP-related gene was evaluated by importance using random forest analysis as previously described [29], and the top 10 were reserved. SVM recursive feature elimination (SVM-RFE) was employed to discard the relatively low predictive feature in each iteration [30]. Hence, SMP genes were ranked from most to least important, and the top 10 were acquired for further analysis. Then, the top 10 SMP-related genes that were acquired separately by random forest and SVM-RFE analysis were applied to Venn analysis, and the SMP-related hub genes were further identified in the peripheral blood samples of OA (GSE63359).

CONSTRUCTION OF OA PREDICTION MODEL:

Multi-factor logistic regression analysis was used to evaluate the diagnostic model. We calculated the area under the receiver operator characteristic (ROC) curve to predict whether the identified hub genes have diagnostic value in OA. According to the results of multi-factor analysis and clinical information, a predictive nomogram was constructed using the “regplot” and “rms” packages in R software [31]. Calibration plots and decision curve analysis (DCA) were used to assess the applicability of the constructed nomogram.

CORRELATION ANALYSIS BETWEEN THE EXPRESSION OF SMP-RELATED HUB GENES AND THE IMMUNE MICROENVIRONMENT:

Using the “relative” and “absolute” methods available on the CIBERSORT (https://cibersortx.stanford.edu/), the local abundance of immune cells in OA synovium was analyzed by the CIBERSORT module. Gene expression deconvolution was performed on RNAseq TPM level data of clinical samples. Based on CIBERSORT and MCP-counter analysis, we derived the relationship between SMP-related hub genes and the abundance of immune-related cells.

CONSENSUS CLUSTERING AND WGCNA ANALYSIS:

“ConsensusClusterPlus” package was utilized to conduct the consensus cluster analysis. The median absolute deviation was used to measure the top 5000 most variable genes applied for samples clustered. The cluster of OA samples and identification of the optimal number of clusters was conducted through unsupervised consensus clustering. Among each 1000 resampling iterations in consensus clustering, we applied ward linkage as an agglomeration method and Euclidean distance as a distance measure to perform hierarchical clustering. The gene sets with similar mRNA transcription models across OA synovium samples were weighted by WGCN analysis [32]. We first normalized the expression of OA synovium genes. Using a power adjacency function, we obtained connection strengths between any 2 genes by converting this matrix into an adjacency matrix. The scale-free topology (SFT) criterion was used to choose the parameter (soft threshold) for adjacency, where the scale-free topology is necessary to construct the network. Based on the SFT criterion recommendation, the optimal threshold parameter value was accepted due to its model-fit saturation of more than 0.85. Then, the optimal topology as soft-threshold power was identified by the “pickSoftThreshold” function. The module was determined using the “tree” method with “deepSplit” to identify modules with no less than 10 genes. Each modular feature was represented by an eigengene generated using the DynamicTreeCut algorithm in the context of WGCNA, as previously described [33]. Default values of WGCNA were consistently applied throughout the process unless stated otherwise. The co-expression modules that enriched for “Biological Process” in Gene Ontology (GO) enrichment analysis were described using Gorilla [34], which is based on pre-inputted gene sets into WGCNA.

GO/KEGG ENRICHMENT ANALYSIS:

The genes that were most obviously different in the WGCNA model were utilized to perform GO and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis based on the DAVID database (https://david.ncifcrf.gov/tools.jsp). All pathways with statistical significance (P<0.05) were retained.

STATISTICAL ANALYSIS:

All data analysis and processing were performed in R (version 4.2.1). SPSS software (version 23.0) was used for statistical analysis, which began by thoroughly evaluating normality and homoscedasticity to determine the appropriate test method, whether parametric or non-parametric. Normally distributed data underwent the two-tailed unpaired t test, whereas non-normally distributed data or cases not meeting t test assumptions were subjected to the Wilcoxon’s rank sum test between the CT and OA groups. All P values were 2-sided, and P<0.05 was considered statistically significant.

Results

IDENTIFICATION OF 6 SMP-RELATED HUB GENES IN OA AYNOVIUM:

After removing the batch effects, we performed the differential expression analysis in the integrated database (Figure 2A). We found that the SMP-related genes were expressed at higher levels in the OA group than in the CT group (Figure 2B, P<0.05). The 11 314 DEGs were detected by differential expression analysis between OA and the healthy group. Then, 48 SMP-related DEGs were obtained through overlapping the 11 314 DEGs and 93 SMP genes. A volcano and a box plot presented the level of 48 SMP-related DEGs (Figure 2C). The heatmap presents the expression heterogeneity of 48 SMP-related DEGs between healthy and OA synovium (Figure 2C). The 48 SMP-related DEGs were ranked by SVM-RFE and random forest analysis (Figure 2D). Then, 6 candidate hub SMP-related genes in OA synovium were obtained by intersecting the top 10 genes ranked by SVM-RFE and random forest analysis (Figure 3A).

THE 6 SMP-RELATED HUB GENES CAN SERVE AS DIAGNOSTIC MARKERS FOR OA:

To verify whether the 6 candidate hub genes can be used as clinical diagnostic markers for OA, we performed a binary univariate/multivariable logistic regression analysis for the 6 SMP-related hub genes, and statistically significant variables were included to develop a nomogram line chart for predicting the onset of OA. The nomogram shows that the possibility of OA was determined by the total scores corresponding to the sum of 6 included hub genes (Figure 3B, nomogram for 1 of the random samples). Then, the ROC curve was drawn to verify the reliability of the constructed diagnostic model. The model has a good diagnostic value for OA (Figure 3C, AUC=0.976). The random bootstrap procedure drawn with replacement from the OA samples (Figure 3C, n=100 bootstraps) was applied to assess the reliability of the constructed model. The calibration curves revealed that the model predictions of the nomogram plots were close to those of the ideal model (Figure 3C). The decision curve analysis revealed that the diagnostic model that included the 6 hub genes worked better for OA patients than the model that included only 1 hub gene (Figure 3D). Figure 4 shows the range of AUC, sensitivity, and specificity in the diagnostic model.

FOUR HUB GENES WERE IDENTIFIED IN PERIPHERAL BLOOD SAMPLES IN OA:

To determine whether the diagnostic marker we identified have clinical practicality, we assessed the expression levels of the 6 hub genes. We acquired a dataset comprising transcriptome data from 26 healthy individuals and 46 OA patients’ peripheral blood samples. Upon comparing the expression levels of 4 key genes in the peripheral blood of healthy subjects and those with OA, we observed significant disparities in the expression of these genes (Figure 5). Thus, the following 4 nominated hub genes were considered hub genes of SMP in OA synovium: beta-1,3-N-acetylgalactosaminyltransferase1 (B3GALNT1), sphingosine-1-phosphate phosphatase 1 (SGPP1), sphingomyelin synthase 1 (SGMS1), and sphingosine kinase 1 (SPHK1).

HUB GENES REGULATED THE IMMUNE MICROENVIRONMENT IN OA SYNOVIUM:

CIBERSORT analysis [35] was utilized to deconvolve the gene expression profiles of hub genes. Subsequently, adhering to the operational guidelines outlined on the CIBERSORT website (https://cibersortx.stanford.edu/), an mRNA-based CIBERSORT algorithm incorporating 100 permutations was created to investigate the relationship between the expression levels of the 4 SMP hub genes, immune cell infiltration, and inflammatory factors (Figure 6). Our findings revealed significant associations for B3GALNT1 with 9 out of 22 immune cell infiltrations, 5 out of 9 immune cell populations, and 13 out of 24 inflammatory cytokine levels. Additionally, we observed clear correlations between the expression level of SGMS1 and 7 out of 22 immune cell infiltrations, 6 out of 24 inflammatory cytokine levels, and a positive relationship with B lineage abundance. In the case of SPHK1, we found significant correlations with 5 out of 22 immune cell infiltrations, 9 out of 24 inflammatory cytokine levels, and the abundance of 2 out of 9 immune cell populations. Furthermore, the expression level of SGPP1 was associated with 4 out of 22 immune cell infiltrations, and 4 out of 24 inflammatory cytokine levels, and showed a negative correlation with endothelial cell abundance.

COMPARISON OF CLINIC PARAMETERS AND IMMUNE CHARACTERISTICS AMONG CLUSTERS:

The connection between identified hub genes and clinical factors was explored in R software with “ggpubr” and “ConsensusClusterPlus” packages [36]. We identified the 3 subtypes of OA patents based on the optimal K value (Figure 7A). We found the expression of 4 hub genes among 3 clusters was significantly heterogeneous, as were the clinical factors containing age and sex in OA samples (Figure 7B). There was no significant association between numbers and age, while there was a significant difference in the probability of different age periods among the 3 clusters (Figure 7C, 7D). The expressions of B3GALNT1, SGPP1, CSF3, and IL2 were upregulated, while SPHK1, PDGFA, and IL1A were downregulated in cluster 2 (Figure 8A, 8B). In cluster 2, there were fewer monocytes, but there were more T cells, cytotoxic lymphocytes, myeloid dendritic cells, and endothelial cells (Figure 8C).

CONSTRUCTION OF CO-EXPRESSION PATTERNS:

No samples were removed after performing quality checks among the 40 OA samples. By employing a power function (beta=10; Figure 9A), an adjacency matrix was created from a pairwise gene correlation matrix, initially computed from individual gene residual values. Subsequently, this matrix was converted into a topological overlap matrix, facilitating module assignments through the utilization of a dynamic tree cutting technique. Using the WGCNA algorithm, we identified all co-expression modules with statistical significance using optimal dynamic tree cut (Figure 9B) and hierarchical clustering (Figure 9C). The MEred modules and cluster 2 had the highest absolute correlation values, and the most relevant molecular changes of cluster 2 were detected in MEred modules. Thus, enrichment analysis was performed using genes from the MEred modules (Table 2). We observed a negative association between the phospholipase C-activating G protein-coupled receptor signaling pathway, signaling receptor activator activity, receptor ligand activity, and cluster 1. Conversely, these factors showed a positive association with cluster 2. No statistically significant difference was observed in cluster 3 (Figure 9D).

Discussion

OA is the most common degenerative osteoarthritic disease; it causes irreversible cartilage destruction and synovial inflammation [1]. In recent years, there have been increasing numbers of studies showing that synovitis plays a vital role in the pathogenic process of OA [8]. In this study, a total of 11 314 DEGs were identified between osteoarthritis (OA) and healthy synovium tissue. Among these, 19 genes were found to overlap (15 upregulated, 4 downregulated) and were associated with OA and SMP. Notably, B3GALNT1, SGPP1, SGMS1, and SPHK1 were identified as pivotal hub genes. Utilizing these hub genes, an OA risk prediction model was constructed, demonstrating their potential utility as diagnostic markers for OA. Additionally, the expression levels of these hub genes showed correlations with immune infiltration and inflammatory factors. Through unsupervised clustering and WGCNA analysis, heterogeneity in GO/KEGG enrichment was observed across different clusters.

Previous studies have shown that the triacylceramide 3-β-N-acetylgalactosamine transferase encoded by B3GALNT1 can encode glycosides and P antigens elicit immune responses [37]. Qin et al found that B3GALNT1 gene expression can regulate immune response [38]. In the present study, we found that B3GALNT1 expression shows a strong positive correlation with B lineage, memory B cells, and plasma cells, while it was negatively correlated with the abundance of naive B cells in OA. The result of B3GALNT1 expression may indicate the increased consumption of naive B cells and the activation of B cells in OA synovium. This observation suggests that the upregulation of B3GALNT1 in OA synovium can induce synovitis by enhancing the transformation of naive B cells into memory B cells and plasma cells. Furthermore, B3GALNT1 upregulates the level of IL15 and aggravates bone destruction [39], which is consistent with our results. Thus, downregulated B3GALNT1 expression may improve OA synovitis, but the specific mechanisms need further exploration. For SGPP1, previous studies reported that SGPP1 expression promotes the infiltration of immune cells [40] and the release of pro-inflammation factors [41]. Some pro-inflammation factors, such as IL-6, IL-1A, and IL-1B, can aggravate the development of local inflammation response of OA synovium [42]. In the present study, we found that the level of SGPP1 was negatively correlated with the level of IL1A but was positively correlated with TGFB1, TGFB2, and CSF3. Thus, SGPP1 expression may regulate synovitis development due to releasing different cytokines. Interestingly, SGPP1 increased the proliferation of rheumatoid arthritis (RA) synovial endothelial cells [43], but decreased the number of OA synovial endothelial cells. This inconsistency may indicate the different pathogeneses of OA and RA. The level of SGMS1 involves the release of inflammatory cytokines [44]. In the present study, SGMS1 expression was negatively correlated with IL6. In addition, it has been reported that IL6 accelerates cartilage matrix destruction [42]. Thus, the level of SGMS1 may be protective against OA injury. Previous studies found that SPHK1 regulates various immune cell functions and inflammatory mediator levels, such as IL6 and TGFB2 [45, 46]. In present study, SPHK1 level was positively correlated with IL-6 level and endothelial cell amount, but was negatively correlated with TGFB2. These results may indicate that SPHK1 can decrease the IL6 level and reduce the number of synovial endothelial cells. Meanwhile, consistent with a previous report, TGFB2 can inhibit the inflammatory response and alleviate OA progression [46]. Notable, the large confidence interval observed in the multifactorial analysis of B3GALNT1, SGMS1, and SGPL1 may be attributed to the relatively small sample size included in this study and confounding factors or unmeasured variables that were not accounted for in our analysis.

Early diagnosis plays a crucial role in timely intervention for osteoarthritis. Utilizing transcriptome data from various pathological tissues within joints, numerous studies have identified a range of biomarkers (hub genes) for early OA detection. For example, Xinyue et al [47] employed diverse bioinformatics screening methods to pinpoint key genes showing distinct expression in OA synovial tissue, suggesting their potential as diagnostic markers for OA. Similarly, Hannah et al [48] used single-cell analysis to highlight the significant role of the gene ZEB1, selected from single-cell transcriptome data of joint cartilage and meniscus, in the pathogenesis of OA. Previous research from our team has also indicated the involvement of genes within the arachidonic acid metabolism pathway in OA diagnosis and local inflammation within synovial tissues [23]. In this study, we found that the expression of SMP-related hub genes in synovium can serve as robust early markers for OA and are correlated with local inflammatory cell infiltration and levels of inflammatory factors, suggesting the tight correlation between the SMP-related hub genes and immune disorder in OA synovium. This provides a treatment direction of synovial immune disorders for OA in the future.

Four hub genes involved in sphingolipid metabolism showed heterogeneity among these 3 clusters. We observed that B3GALNT1 and SGPP1 were specifically upregulated in cluster 2, and both of them can be checked in OA patients’ peripheral blood. Thus, the B3GALNT1 and SGPP1 genes may serve as accessible potential biomarkers for OA diagnosis, especially in the elderly population. Utilizing peripheral blood testing for diagnosing OA presents the advantages of convenience and cost-effectiveness, highlighting the importance of gathering additional peripheral blood samples for validation in future studies. The relative score of T cells, endothelial cells, myeloid dendritic cells, and cytotoxic lymphocytes were upregulated, while monocytic lineage was downregulated in cluster 2. This indicates that elderly patients with OA may show a higher level of synovial vascularization and immune cell disorder [49]. We also found that CSF3 and IL2 were upregulated in cluster 2 and PDGFA was downregulated. Downregulation of CSF3 and IL2 significantly reduces the production of inflammatory markers [50,51]. Low levels of CSF3 and IL2 can relieve intractable chronic pain [52–54], and high levels of PDGFA can improve the development of OA synovitis [55]. These results may indicate that downregulation of CSF3 and IL2 and upregulation of PDGFA are potential therapeutic targets for cluster 2 (elderly OA patients). This study offers a direction for personalized treatment of osteoarthritis patients in the future, emphasizing the importance of addressing local immune dysregulation in elderly individuals with OA.

WGCNA analysis showed that MEred module eigengenes were negatively correlated with cluster 1 and cluster 2. MEred module eigengenes mainly involve the “phospholipase C-activating G protein-coupled receptor signaling pathway” in the biological process (BP). The opposite process may indicate higher activation of phospholipase C in cluster 2 compared with cluster 1. Interestingly, the activation of phospholipase C was proved to be harmful in OA, and inhibition of phospholipase C-mediated autophagy may benefit OA cartilage protection [56]. Therefore, targeting the activation of phospholipase C may be a sensitive intervention in cluster 2 patients. Meanwhile, the treatment strategies for inhibiting the activity of phospholipase C in cluster 1 patients may be ineffective. Therefore, the “phospholipase C-activating G protein-coupled receptor signaling pathway” should receive more attention in future research on treatment of elderly OA patients.

The study has limitations that need to be acknowledged. The influence of additional variables affecting the transcription level of inflammatory genes in OA synovium, such as corticosteroids, NSAIDs, and disease duration, remains unexplored. While the study proposed that SMP genes could potentially function as diagnostic markers for OA and interact with immune-mediated inflammation in the synovium, further validation through animal and clinical studies is crucial.

Conclusions

We conducted an investigation into the impact of sphingolipid metabolism in osteoarthritis (OA) synovium and observed significant upregulation of SMP genes. Among these, 4 hub genes (B3GALNT1, SGPP1, SGMS1, and SPHK1) identified from peripheral blood exhibit potential as novel diagnostic markers in OA patients. These hub genes showed significant correlations with inflammatory factors and immune-related cells. Additionally, our findings revealed differential expression of SMP hub genes among different clusters of OA patients, along with variations in the infiltration of immune-related cells and the level of inflammatory cytokines across the 3 clusters. This study suggests that SMP may play a role in regulating the immune microenvironment in OA. The participation of SMP-related hub genes in local synovial tissue immune dysregulation could be a pathogenic mechanism of OA, and targeting these hub genes may offer a potential solution to alleviate joint inflammation. However, further research is needed to fully elucidate the underlying mechanisms involved.

Data Sharing Statement

The datasets used and/or analyzed during the current study are available from the corresponding authors upon reasonable request.