SLC38A1 and STX11 are mitochondria-related biomarkers associated with immune infiltration in osteoarthritis

Background: Mitochondrial dynamics and mitophagy play crucial roles in osteoarthritis (OA); however, the specific contributions of mitochondrial dynamics-related genes (MD-RGs) and mitophagy-related genes (MP-RGs) remain unclear. This study aimed to elucidate the precise mechanisms linking these genes in the context of OA.

Methods: OA-related transcriptome datasets and single-cell RNA sequencing (scRNA-seq) dataset incorporating MD-RGs and MP-RGs were utilized in this study. Hub genes were identified through differential expression analysis, weighted gene co-expression network analysis (WGCNA), and machine learning. A nomogram was then constructed based on the hub genes. Enrichment and immune infiltration analyses were performed on the hub genes, and key cell types were identified based on hub gene expression. Finally, the expression of the hub genes was validated using reverse transcription-quantitative polymerase chain reaction (RT-qPCR).

Results: SLC38A1 and STX11 were identified as hub genes linked to mitochondrial dynamics and mitophagy in OA. These genes enabled the construction of a reliable nomogram for predicting OA risk. Enrichment analysis revealed that the top biological processes converged on the ECM–receptor interaction, underscoring its critical role in OA pathogenesis. Immune infiltration analysis uncovered significant disparities in 10 immune cell types, including activated CD4 T cells and central memory CD4 T cells, between OA patients and healthy controls. The levels of these immune cells were strongly correlated with the expression of SLC38A1 and STX11. Additionally, endothelial cells, monocytes, and T cells emerged as key cellular players in OA. RT-qPCR validation showed that SLC38A1 was significantly downregulated in OA samples, and STX11 exhibited a similar trend, suggesting their potential roles in OA progression.

Conclusion: This study identified SLC38A1 and STX11 as key genes linked to mitochondrial dynamics and mitophagy in OA. These findings provide a theoretical basis and valuable reference for the diagnosis and treatment of OA.

1 Introduction

Osteoarthritis (OA) is a chronic joint disease characterized by degenerative changes in the articular cartilage, synovium, and subchondral bone (Steinmetz et al., 2023). It is the most common joint disorder, with a prevalence that increases with age, affecting approximately 595 million people worldwide in 2020 (Steinmetz et al., 2023). OA has a series of complex pathological changes, including articular cartilage abrasion, synovial inflammation, subchondral bone remodeling, and osteophyte formation, which influence the development and progression of OA (Chu et al., 2020; Chilelli et al., 2024). OA is an irreversible degenerative disease, presenting with late-stage symptoms such as joint pain and deformity, which severely impact the daily activities of the patients (Duruöz et al., 2023). The causes of OA include age, gender, genetics, metabolism, and joint injury (Tibor and Ganz, 2022). Chondrocyte death, including autophagy, ferroptosis, apoptosis, and pyroptosis, contributes to the progression of OA (Guan et al., 2024; Zhu et al., 2023). Therefore, early diagnosis and timely intervention are extremely important as they can help delay disease progression and alleviate symptoms. In recent years, an increasing number of studies have confirmed that various dysregulated genes can serve as important diagnostic markers and therapeutic targets (Hadzic and Beier, 2023).

Mitochondria are essential organelles within cells that carry out and coordinate various metabolic processes and play a significant role in the development of OA (Blanco et al., 2004). Studies have shown that ERK1/2 is a key factor in promoting IL-1-induced mitochondrial fission and apoptosis in chondrocytes (Ansari et al., 2022). The mitochondrial network in normal chondrocytes remains intact, whereas several chondrocytes in OA cartilage exhibit excessive mitochondrial fragmentation (Blanco et al., 2004; Ansari et al., 2022). Mitochondrial dysfunction is also associated with OA (Blanco et al., 2004). Mitophagy is a process that selectively removes damaged or dysfunctional mitochondria through autophagy, thereby maintaining mitochondrial quality control and homeostasis. It has been reported that HIF-1α can alleviate cell apoptosis and senescence in chondrocytes under hypoxic conditions through mitophagy, thereby ameliorating cartilage degeneration in surgically induced OA mouse models (Hu et al., 2020). Additionally, curcumin exerts chondroprotective effects in osteoarthritis by promoting AMPK/PINK1/Parkin-mediated mitophagy (Jin et al., 2022). Therefore, genes related to mitochondrial dynamics and autophagy may serve as biomarkers for OA patients.

In summary, based on the single-cell and bulk transcriptome analyses, we have explored the potential of mitochondrial dynamics and autophagy as biomarkers for OA patients and their underlying molecular mechanisms, providing new insights for the early clinical diagnosis and treatment of OA patients.

2 Materials and methods

2.1 Data collection

OA-related transcriptome datasets (GSE57218 and GSE117999) and the single-cell RNA sequencing (scRNA-seq) dataset (GSE152805) (GPL20301) were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). The GSE57218 dataset (GPL6947), containing 7 control and 33 OA cartilage tissue samples, was considered the training set (Ramos et al., 2014). The GSE117999 dataset (GPL20844), consisting of 12 OA and 12 control cartilage tissue samples, was considered the validation set. Detailed information on the samples in the training and validation sets is presented in Supplementary Table S1. The GSE152805 dataset included three OA cartilage tissue samples for scRNA-seq. A total of 23 mitochondrial dynamics-related genes (MD-RGs) were obtained from the published literature (Zhang et al., 2024). In total, 29 mitophagy-related genes (MP-RGs) were acquired from the Reactome database (Yang et al., 2022).

2.2 Single-cell analysis

Seurat objects from the scRNA-seq data were created using the Seurat package (version 4.1.0) (Tan et al., 2023). In the GSE152805 dataset, cells with fewer than 200 detected genes and genes expressed in fewer than three cells were filtered out. Subsequently, cells meeting the following criteria were retained: (1) the number of detected features (nFeature_RNA) was greater than 1,000 and less than 5,000; (2) total RNA counts (nCount_RNA) were below 30,000; and (3) the proportion of mitochondrial gene expression was less than 5%. Following data normalization using the NormalizeData function in the Seurat package, the FindVariableFeatures method was employed to select high-variable genes. Subsequently, principal component analysis (PCA) was performed, and the p-value of each principal component (PC) was calculated using the JackStraw function in the Seurat package. The significance of PCs was assessed using the ScoreJackStraw function, and PCs with statistically significant differences (p < 0.05) were selected for subsequent analysis. The scree plot of PCs was drafted using the ElbowPlot function in the Seurat package. Subsequently, unsupervised cluster analysis of cells was performed using the FindNeighbors and FindClusters functions of the Seurat package to identify cell clusters with a resolution of 0.3. Furthermore, cell clusters were annotated using the SingleR package (version 1.831) (Zhang et al., 2022), with the HumanPrimaryCellAtlasData from the celldex package used as the reference gene set, and cell subpopulations were subsequently identified. Moreover, marker genes for the different cell subpopulations were identified using the FindAllMarkers function in the Seurat package, and the scRNA-seq of differentially expressed genes (scRNA-seq DEGs) were selected using the Wilcoxon test, with a threshold of |log₂fold change (FC)| > 0.25.

2.3 Differential expression analysis

The limma package (version 3.54.0) (Li et al., 2022) was used to identify bulk differentially expressed genes (bulk DEGs) between the OA and control groups in the GSE57218 dataset, with conditions of adj.p < 0.05 (FDR correction was used) and |log₂FC| > 0.5. A volcano map of bulk DEGs was generated using the ggplot2 package, and the top 10 up- and downregulated bulk DEGs (based on log₂FC) were displayed in a heatmap, which was created using the ComplexHeatmap package (version 2.15.1) (Gu, 2022).

2.4 Weighted gene co-expression network analysis

The MD-RGs and MP-RGs scores of samples in the GSE57218 dataset were calculated using the single-sample gene set enrichment analysis (ssGSEA) algorithm in the GSVA package (version 1.42.0) (Hänzelmann et al., 2013). Then, the difference in the MD-RGs and MP-RGs ssGSEA scores was compared between the OA and control groups using the Wilcoxon test (p < 0.05). To identify genes related to the ssGSEA scores of MD-RGs and MP-RGs, weighted gene co-expression network analysis (WGCNA) was performed using the WGCNA package (version 1.70–3) (Langfelder and Horvath, 2008). First, samples in the GSE57218 dataset were clustered to determine whether any outlier samples needed to be removed. Subsequently, a soft threshold (β) selection was performed to construct a co-expression network by selecting a scale-free R² value greater than 0.9 and a mean connectivity value close to 0, ensuring that the constructed network corresponded more closely to a scale-free topology. The dynamic tree cutting algorithm was employed to partition modules, with a minimum requirement of 300 genes per module. Finally, the ssGSEA scores were integrated with the WGCNA module eigengenes, and the Pearson correlation between the modules and the phenotypic traits (MD-RGs and MP-RGs) was calculated using the cor function (|cor| > 0.3; p < 0.05). Genes in key modules were identified as the key module genes.

2.5 Identification and analysis of candidate genes

Candidate genes were obtained by taking the intersection of scRNA-seq DEGs, bulk DEGs, and key module genes. The results were visualized using a Venn diagram, which was created using the VennDiagram package (version 1.7.1) (Chen and Boutros, 2011). Furthermore, biological functions and pathways involved with the candidate genes were explored by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. The above enrichment analyses were processed using the clusterProfiler package (version 4.2.2) (Yu et al., 2012) with p < 0.05.

2.6 Identification of hub genes

Based on candidate genes, least absolute shrinkage and selection operator (LASSO) (conducted a 10-fold cross-validation) and XGBoost (the importance of genes is quantified using the gain indicator) machine-learning algorithms were applied using the glmnet package (version 4.1–2) (Shi et al., 2023) and XGBoost package (version 1.6.2.1) (Kui et al., 2022) to identify signature genes. Then, the intersections of signature genes identified using two machine-learning algorithms were identified as the candidate hub genes. Furthermore, the expression of candidate hub genes was compared between the OA and control groups using the Wilcoxon test (p < 0.05). The genes that showed significant differences between the OA and control groups, with consistent expression trends in the GSE57218 and GSE117999 datasets, were identified as hub genes.

2.7 Construction of the nomogram

According to hub genes, a nomogram was established to predict the risk of OA using the rms package (version 6.5.0) (Pan et al., 2021). The receiver operating characteristic (ROC) curve was generated using the pROC package (version 1.18.0) (Robin et al., 2011) to explore the ability of the nomogram to distinguish between OA and control samples. Decision curve analysis (DCA) was performed to further evaluate the efficacy of the nomogram.

2.8 Functional analysis of hub genes

According to the expression of hub genes, OA samples in the GSE57218 dataset were divided into high and low expression groups. Gene set enrichment analysis (GSEA) was performed to find pathways using the clusterProfiler package. The reference gene set was c2.cp.kegg.v7.5.1.symbols.gmt, and the enriched threshold was | normalized enrichment score (NES)| > 1, NOM P < 0.05, and q < 0.25. GeneMANIA (https://genemania.org/) was used to identify genes that were functionally related to hub genes, and a gene–gene interaction (GGI) network was constructed. To explore tissue-specific expression of hub genes, BioGPS (http://biogps.org/) was utilized to predict the expression of hub genes in various tissues (species: human).

2.9 Immune infiltration analysis

To explore the correlation of hub genes and immune cells, the enrichment scores of 28 types of immune cells in the training set were calculated using the ssGSEA algorithm (Bindea et al., 2013). Subsequently, the Wilcoxon test was employed to compare the differences between immune cell enrichment scores between the OA and control groups (p < 0.05). After that, correlations between immune cells and between immune cells and hub genes were analyzed using the corrplot package (version 0.92) (Liu et al., 2022).

2.10 Regulatory network and drug prediction analyses

MicroRNA (miRNA)-targeting hub genes were predicted using the miRDB (https://mirdb.org/) and TargetScan databases (https://www.targetscan.org/vert_80/). The intersections of miRNAs derived from two databases were selected as key miRNAs. Then, long non-coding RNAs (lncRNAs) targeting key miRNAs were obtained from the starBase database (http://starbase.sysu.edu.cn/). To find potential drugs for the treatment of OA, Comparative Toxicogenomics Database (CTD, http://ctdbase.org) analysis was carried out. The lncRNA–miRNA–mRNA regulatory network and the drug–hub gene network were constructed using Cytoscape software (version 3.8.2) (Yue et al., 2021).

2.11 Cell communication and pseudo-time analysis

Cell communication of cell subpopulations was evaluated using the CellChat package (version 1.6.1) (Fang et al., 2022). After creating the CellChat object, importing objects to the ligand–-receptor database (CellChatDB.human), and performing preprocessing, the cell communication network was inferred. A bubble diagram was drafted to show the communication probability of ligand–receptor pair regulation from some cell groups to other cell groups. To identify key cells in the GSE152805 dataset, the expression of hub genes across cells was analyzed, and cells exhibiting high and widespread expression of the hub genes were defined as key cells. Cell trajectory differentiation of key cells was analyzed using the monocle package (version 1.0.0) (Zhai et al., 2020). In addition, the expression of hub genes in key cells differentiation was explored. Subsequently, enrichment of transcription factors (TFs) was evaluated using single-cell regulatory network inference and clustering (SCENIC, https://github.com/aertslab/SCENIC) based on single-cell data.

2.12 Reverse transcription-quantitative polymerase chain reaction

To investigate the expression levels of hub genes, an RT-qPCR experiment was performed. Initially, total RNA was extracted from the tissue samples using TRIzol reagent (Ambion, Austin, USA), followed by determination of RNA concentration. Subsequently, cDNA synthesis was carried out via reverse transcription using the SweScript First-Strand cDNA Synthesis Kit (ServiceBio, Wuhan, China). Finally, quantitative analysis was performed using the Universal Blue SYBR Green qPCR Master Mix (ServiceBio, Wuhan, China), and the gene’s relative expression level was calculated using the 2^–△△Ct method, with GAPDH used as an internal reference. Additionally, the primer sequences for the relevant genes are shown in Supplementary Table S2.

2.13 Statistical analysis

R software (version 4.2.2) was used to process and analyze data. Statistical significance between the two groups was assessed using the Wilcoxon rank-sum test. A p-value < 0.05 was considered statistically significant. The overall analysis process of this study is shown in Figure 1.

Figure 1

Figure 1. Analysis flowchart.

3 Results

3.1 Single-cell profiling revealed OA-related immune and endothelial subpopulations

After quality control of the scRNA-seq data from GSE152805, 10,032 cells and 19,050 genes remained following filtering (Supplementary Figure S1). Out of 2,000 highly variable genes, the p-value of the top 50 PCs was far less than 0.05. Then, the first 20 PCs were chosen for subsequent analysis through PCA (Figures 2A–C). As shown in Figure 2D, 11 cell clusters were obtained through uniform manifold approximation and projection (UMAP) cluster analysis. Subsequently, 6 cell subpopulations were annotated based on 11 cell clusters (chondrocytes, endothelial cells, macrophages, monocytes, T cells, and tissue stem cells) (Figure 2E). A total of 1,860 scRNA-seq DEGs were further selected among different cell subpopulations, and the expression of marker genes that relied on |log₂FC| of each cell subpopulation was determined (Figure 2F).

Figure 2

Figure 2. Single-cell sequencing analysis of osteoarthritis cartilage tissues. (A) Acquisition of highly variable genes. (B) Selection of meaningful PCs. (C) Principal component and standard deviation distribution plot. (D) UMAP plots of cartilage-associated cells colored by cluster. (E) Annotation of different cell clusters. (F) Expression levels of marker genes in various cell clusters.

3.2 Bulk DEGs associated with OA were screened out

Through differential expression analysis, 1,124 bulk DEGs (509 upregulated and 615 downregulated genes) were screened out between the OA and control groups (Figures 3A, B). In the OA group, the ssGSEA score of MD-RGs and MP-RGs was significantly lower than that in the control groups (Figures 3C, D).

Figure 3

Figure 3. Bulk RNA sequencing analysis of osteoarthritis and control cartilage tissues. (A) Volcano plot of differentially expressed genes in bulk RNA-seq. (B) Circle heatmap of differentially expressed genes. (C) Comparison of mitochondrial dynamics_single-sample gene set enrichment analysis (MD_ssGSEA) scores between the OA and control groups. (D) Comparison of mitophagy_ssGSEA (MP_ssGSEA) scores between the OA and control groups. *p < 0.05 and **p < 0.01.

3.3 SLC38A1 and STX11 emerge as mitochondrial hub genes via integrative analysis

Cluster analysis was performed on all samples in the GSE57218 dataset, and no outlier sample was found. Thus, all samples were used for subsequent analysis (Supplementary Figure S2). By setting R² = 0.9 as the threshold, a soft threshold of 4 was selected, and data with mean connectivity close to 0 were filtered out (Figure 4A). Consequently, 9 co-expression modules were obtained through similarity analysis and by setting a minimum of 300 genes for each gene module (Figure 4B). Notably, as shown in Figure 4C, MEblue and MEturquoise modules, respectively, exhibited a strong positive (cor = 0.78 and 0.76; p < 0.05) and negative correlation (cor = −0.65 and −0.69, p < 0.05) with MD-RGs and MP-RGs ssGSEA scores, respectively. By combining 2,866 genes in the MEturquoise module and 2,420 genes in the MEblue module, 5,286 key module genes were finally acquired.

Figure 4

Figure 4. WGCNA. (A) Selection of the optimal soft threshold (β) value. (B) Dendrogram of gene clusters. (C) Correlation between modules and MD-RGs and MP-RGs.

3.4 SLC38A1 and STX11 were selected as hub genes

Intersections of 1,860 scRNA-seq DEGs, 1,124 bulk DEGs, and 5,286 key module genes yielded 12 candidate genes (Figure 5A). Enrichment analysis of candidate genes showed that 236 GO terms and 7 KEGG pathways were enriched. Regulation of leukocyte-mediated cytotoxicity, solute: sodium symporter activity, and phagocytic vesicle membrane-related functions were GO terms significantly enriched by the candidate genes (Figure 5B). For KEGG pathways, viral myocarditis, type-I diabetes mellitus, SNARE interactions in vesicular transport, cell adhesion molecules, etc., were the main pathways (Figure 5C). In total, eight signature genes (SLC7A8, SOD3, STX11, TSHZ1, NPTX2, SLC38A1, EFHD1, and ARID5B) were selected by LASSO analysis when lambda was 0.0232298 (Figure 5D). Through XGBoost analysis, six signature genes, namely, SOD3, SLC7A8, NPTX2, SLC38A1, NPDC1, and STX11, were selected as gene features (Figure 5E). The intersection of eight signature genes in LASSO analysis and six signature genes in XGBoost analysis yielded five candidate hub genes, including SLC7A8, SOD3, and STX11 (Figure 5F). Expression analysis of candidate hub genes showed that SLC38A1 and STX11 were obviously less expressed in the OA group in the GSE57218 and GSE117999 datasets (Figure 5G). Therefore, SLC38A1 and STX11 were identified as hub genes for further analysis.

Figure 5

Figure 5. Screening of candidate genes with machine learning analysis. (A) Intersection of scRNA-seq differentially expressed genes (scRNA-seq DEGs), bulk DEGs, and key module genes. (B) GO enrichment of candidate genes. (C) KEGG pathway enrichment of candidate genes. (D) Lasso regression. (E) XGBoost analysis. (F) Intersection of candidate genes in Lasso and XGBoost analysis. (G) Violin plot of relative expression levels of candidate genes in GSE57218 and GSE117999 datasets. *p < 0.05, **p < 0.01, ***p < 0.001, and ****p < 0.0001.

3.5 Functional analysis of hub genes

Based on the hub genes, a nomogram was constructed to predict the risk of OA (Figure 6A). The area under the curve (AUC) value of the ROC curve was 0.944, which demonstrated that the nomogram had a good ability to diagnose OA and control samples (Figure 6B). Furthermore, DCA showed that the nomogram and hub genes had a higher net benefit than the extreme curve and single factor (Figure 6C). To find biological functions involved in hub genes, GSEA was performed. The results showed that the ECM–receptor interaction was consistently among the top five enriched pathways based on enrichment scores of the hub genes (Figures 6D, E). Through GeneMANIA, we identified 20 genes that are functionally related to hub genes. A GGI network was constructed, showing interactions such as GLUL–SLC38A1, GLS2–SLC38A1, and GLS–SLC38A1. This analysis revealed that hub genes were involved in multiple vital biological processes, including amino acid transport, carboxylic acid transport, and regulation of neurotransmitter levels (Figure 6F).

Figure 6

Figure 6. Construction of a nomogram and functional analysis of hub genes. (A) A nomogram is used to predict the risk of OA. (B) Area under the curve value of the ROC curve. (C) Model evaluation curves show that the model containing the two identified hub genes has a higher net benefit. (D) Gene set enrichment analysis of the SLC38A1 gene. (E) Gene set enrichment analysis of the STX11 gene. (F) GGI network using GeneMANIA.

3.6 Hub genes correlate with NK cell infiltration in OA cartilage

Among 28 immune cell types, 10 types of immune cells had significant differences between the OA and control groups, including activated CD4 T cells, activated dendritic cells, and central memory CD4 T cells (Figure 7A). Correlation analysis showed that activated B cells had the strongest negative correlation with natural killer (NK) T cells (cor = −0.43; p < 0.05), while regulatory T cells had the strongest positive correlation with macrophages (cor = 0.88; p < 0.05) (Figure 7B). Moreover, hub genes had a negative correlation with NK cells (cor < −0.3; p < 0.05), and SLC38A1 was positively correlated with plasmacytoid dendritic cells (cor >0.3; p < 0.05) (Figure 7C).

Figure 7

Figure 7. Analysis of immune cell components in single-cell sequencing. (A) Cell fraction difference of 28 immune cell types between the OA and control groups. (B) Correlation analysis of all immune cells. (C) Correlation between hub genes and all immune cells. *p < 0.05, **p < 0.01, and ***p < 0.001.

3.7 Comprehensive regulatory network and drug prediction of hub genes

To understand the regulatory relationship of hub genes, an lncRNA–miRNA–mRNA regulatory network was constructed with SLC38A1, 3 miRNAs, and 18 lncRNAs. In this network, all lncRNAs regulated SLC38A1 through hsa-miR-23a-3p, hsa-miR-23b-3p, and hsa-miR-23c (Figure 8A). Through computer simulation assumptions, SLC38A1 predicted a total of 14 drugs, and STX11 predicted a total of 18 drugs. Among these drugs, methyl methanesulfonate was predicted by SLC38A1 and STX11, quercetin was predicted by SLC38A1, and oxazolone was predicted by STX11 (Figure 8B). Furthermore, the top 10 tissues with the highest expression scores for the two hub genes commonly included the testis interstitium, ovary, dorsal root ganglion, ciliary ganglion, and fetal lung (Figure 8C).

Figure 8

Figure 8. Analysis of the lncRNA–miRNA–mRNA regulatory network and drugs prediction of two hub genes. (A) LncRNA–miRNA–mRNA regulatory network of SLC38A1. (B) Predicted drugs for SLC38A1 and STX11. (C) Top 10 tissues with the highest expression scores for the two hub genes.

3.8 Endothelial cells, monocytes, and T cells were identified as key cells

Through cell communication analysis, chondrocytes were shown to have a strong correlation with other cell subpopulations (Figure 9A). Following that, monocytes had the strongest correlation with monocytes, and the ligand–receptor pairs were IL-1B–IL-1R2 (Supplementary Figure S3). To select key cells, the expression of hub genes was analyzed, and the results showed that hub genes exhibited higher and widespread expression in endothelial cells, monocytes, and T cells. Thus, these three cell types were identified as key cells for subsequent analysis (Figure 9B). There were five stages of differentiation for endothelial cells, arranged from top to bottom (Supplementary Figure S4A). For monocytes and T cells, there were eight and three differentiation stages, respectively, arranged from left to right (Supplementary Figure S4B, C). In endothelial cells, the expression of SLC38A1 initially increased in early differentiation stages and then decreased before gradually increasing again. In contrast, STX11 showed a gradual decrease initially, followed by a slight increase and a subsequent rapid increase during differentiation. For monocytes, SLC38A1 expression increased gradually throughout differentiation, while STX11 initially increased slowly, followed by a plateau phase. In T cells, SLC38A1 expression initially increased and then declined rapidly, whereas STX11 showed a slow initial decrease, followed by a gradual increase and a subsequent decline during differentiation (Supplementary Figure S5). According to scRNA-seq data, 195 TFs were screened out, and the top 30 TFs included TBX21, SPIB, and MAF (Supplementary Figure S4D).

Figure 9

Figure 9. Communication between different cell types. (A) Interaction number and strength between different cell clusters. (B) Violin plots of the expression levels of two hub genes in various cell clusters.

3.9 Differential expression of SLC38A1 in the OA and control groups

RT-qPCR is a molecular biology technique that is widely used in clinical diagnostics and drug development to measure the expression levels of specific genes in samples (Falini and Dillon, 2024). In this study, RT-qPCR was utilized to assess the mRNA expression levels of hub genes to investigate their potential differential expression in OA samples. The results showed differential expression of SLC38A1 between the control and OA groups, with lower expression levels detected in the OA group (p < 0.01) (Figure 10A). The result consisted of the expression trend of SLC38A1 in the OA group in the GSE57218 and GSE117999 datasets. However, the expression level of STX11 in the OA group exhibited a decrease, but the change did not reach statistical significance (p > 0.05) (Figure 10B). This might be due to the small sample size used, and it will need to be confirmed with a larger sample size in future studies.

Figure 10

Figure 10. Detection of relative expression levels of SLC38A1 and STX11. (A) Relative expression level of SLC38A1 in OA and control tissues. (B) Relative expression level of STX11 in OA and control tissues. **p < 0.01; ns represents no statistical significance.

4 Discussion

In this study, we aimed to unravel the intricate mechanisms linking MD-RGs and MP-RGs in OA. The findings significantly contribute to the understanding of the pathogenesis of OA at the molecular level and offer new perspectives for its diagnosis and treatment.

The pathways enriched by candidate genes—such as phagocytic vesicle membrane, type-I diabetes mellitus, and cell adhesion molecules—have been previously reported in the literature to be associated with OA (Rios-Arce et al., 2022; Njoto et al., 2019; Cheng et al., 2020; Chen et al., 2020). According to research, TMEM230 promotes antigen processing, transport, and presentation by regulating the membrane-bound organelles. It also plays a crucial role in regulating the secretion of metalloproteinase and heparanase, which are required for tissue remodeling in OA and rheumatoid arthritis (RA) (Abeni et al., 2024). Type-I diabetes mellitus may suppress the development of post-traumatic osteoarthritis (Rios-Arce et al., 2022). Neural cell adhesion molecule (NCAM) could inhibit hypertrophic chondrocyte differentiation of mesenchymal stem cells by suppressing ERK signaling, thereby reducing chondrocyte hypertrophy in OA models. These findings suggest that the potential utility of NCAM as a novel therapeutic target for OA (Cheng et al., 2020).

We found that SLC38A1 was downregulated and that STX11 showed a downregulated expression trend in OA, which might serve as potential biomarkers associated with mitochondrial dynamics and mitophagy in OA through bioinformatics analysis. SLC38A1 (solute carrier family 38, member 1) is a sodium-coupled neutral amino acid transporter, which is primarily responsible for the transport of amino acids across the cell membrane. Although there are currently no reports on the functional mechanism of the SLC38A1 gene in osteoarthritis, many other members of the solute carrier family have been reported to play a role in OA. For example, capsiate inhibited the expression of HIF-1α by activating SLC2A1, thereby reducing the progression of ferroptosis-related osteoarthritis (Guan et al., 2023). Additionally, miR-19b-3p in exosomes from osteoarthritic fibroblast-like synoviocytes enhanced chondrocyte ferroptosis and damage in OA by sponging off SLC7A11 (Kong et al., 2023). STX11 (Syntaxin 11) is mainly expressed in immune cells, is involved in cytotoxic granule exocytosis, and is crucial for the cytotoxic functions of CD8⁺ T cells and NK cells (D'Orlando et al., 2013). In a bioinformatics study of rheumatoid arthritis and osteoarthritis, STX11 was found to be highly enriched in protein–protein interaction networks (Qiu et al., 2021). Additionally, STX11 may modulate OA-associated inflammation and OA progression by affecting the activity of macrophages, T cells, and NK cells (D'Orlando et al., 2013; Offenhäuser et al., 2011; Dabrazhynetskaya et al., 2012). Thus, STX11 may be vital in OA pathogenesis via immune cell function and cytokine regulation, further suggesting its potential relevance in OA and RA. Future research should use knockdown or overexpression techniques in cell models and gene agonists or inhibitors in animal models to further clarify the specific molecular mechanisms of SLC38A1 and STX11 in OA progression.

Single-cell analysis showed high SLC38A1 and STX11 expression in monocytes and immune cells, and the proportion of NK cells was found to be significantly higher in the OA group than in the control group, according to immune infiltration analysis. Monocytes, when exposed to various stimuli, express ligands for the immune receptor found on NK cells. These ligands activate NK cells and enhance their cytotoxicity (Getz et al., 2007). Additionally, both monocytes and NK cells secrete cytokines that can influence each other’s functions. Monocytes produce IL-15, which boosts NK cells’ cytotoxicity and IFN-γ production (Cooper et al., 2001). In turn, IFN-γ from NK cells activates monocytes and enhances their antigen presentation ability (Cooper et al., 2001). These immune cells’ interactions may influence OA progression. Moreover, the correlation analysis of SLC38A1 and STX11 with all immune cells showed a negative correlation with NK cells, further suggesting that SLC38A1 and STX11 may influence the development of osteoarthritis through the immune infiltration of NK cells.

The construction of the lncRNA–-miRNA–mRNA regulatory network for SLC38A1 provided insights into the post-transcriptional regulation of this hub gene. The prediction of drugs associated with SLC38A1 and STX11, such as methyl methanesulfonate, quercetin, and oxazolone, offers potential therapeutic targets for OA treatment. Quercetin alleviates the progression of osteoarthritis by regulating inflammatory cascades and chondrocyte apoptosis (Hu et al., 2019; Wang et al., 2023), and SLC38A1 may serve as an intermediate through which it exerts these effects. These predicted drugs could potentially modulate the functions of the hub genes and, thus, alleviate the pathological symptoms of OA.

However, there are limitations to this study, including limited in vitro and in vivo validation of the functions of SLC38A1 and STX11, the use of RT-qPCR to validate the expression of hub genes, and the need for more versatile experiments, such as gene knockout or overexpression in cell or animal models, to fully understand their role in OA. Second, the small sample size of RT-qPCR and scRNA-seq data, along with the lack of control samples in the scRNA-seq dataset, may affect the accuracy of the result. Only SLC38A1 showed significant differential expression in RT-qPCR. In addition, STX11 did not reach statistical significance, which might be due to the small sample size. However, the STX11 expression trend matched the bioinformatics results in the GSE57218 and GSE117999 datasets, and it was widely expressed in key cell subsets such as endothelial and mononuclear cells in single-cell data, leading us to designate it as a hub gene. Additionally, while the study primarily focused on the relationship between mitochondrial-related genes and OA, other contributing factors were not fully explored. Moreover, the lncRNA–miRNA–mRNA regulatory network and drug predictions were based on database predictions. The actual regulatory roles need further exploration of their molecular mechanisms through molecular experiments (such as siRNA and mimic) and cellular models (such as drug induction). Future studies should expand the scope to deepen our understanding of OA and develop more comprehensive diagnostic and therapeutic strategies.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Ethics statement

These studies involving humans were approved by the Ethical Committee of the Affiliated Hospital of Shandong University of Traditional Chinese Medicine. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

Author contributions

WL: data curation, formal analysis, writing – original draft. MY: data curation, formal analysis, writing – original draft. WY: visualization writing – review and editing. YC: conceptualization, supervision, writing – review and editing. WW: conceptualization, supervision, writing – review and editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work was supported by the Project of China Association for the Promotion of Traditional Chinese Medicine Research (2021-18) and the Shandong Provincial Archives Science and Technology Project (2024-42).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2025.1585775/full#supplementary-material

References

Abeni, E., Cocola, C., Croci, S., Martino, V., Piscitelli, E., Gualtierotti, R., et al. (2024). Single-cell transcriptomic analysis to identify endomembrane regulation of metalloproteins and motor proteins in autoimmunity. Adv. Protein Chem. Struct. Biol. 141, 299–329. doi:10.1016/bs.apcsb.2024.03.007

CrossRef Full Text | Google Scholar

Ansari, M. Y., Novak, K., and Haqqi, T. M. (2022). ERK1/2-mediated activation of DRP1 regulates mitochondrial dynamics and apoptosis in chondrocytes. Osteoarthr. Cartil./OARS, Osteoarthr. Res. Soc. 30 (2), 315–328. doi:10.1016/j.joca.2021.11.003

CrossRef Full Text | Google Scholar

Bindea, G., Mlecnik, B., Tosolini, M., Kirilovsky, A., Waldner, M., Obenauf, A. C., et al. (2013). Spatiotemporal dynamics of intratumoral immune cells reveal the immune landscape in human cancer. Immunity 39 (4), 782–795. doi:10.1016/j.immuni.2013.10.003

CrossRef Full Text | Google Scholar

Blanco, F. J., López-Armada, M. J., and Maneiro, E. (2004). Mitochondrial dysfunction in osteoarthritis. Mitochondrion 4 (5-6), 715–728. doi:10.1016/j.mito.2004.07.022

CrossRef Full Text | Google Scholar

Chen, H., and Boutros, P. C. (2011). VennDiagram: a package for the generation of highly-customizable venn and euler diagrams in R. BMC Bioinforma. 12, 35. doi:10.1186/1471-2105-12-35

CrossRef Full Text | Google Scholar

Chen, W. C., Lin, C. Y., Kuo, S. J., Liu, S. C., Lu, Y. C., Chen, Y. L., et al. (2020). Resistin enhances VCAM-1 expression and monocyte adhesion in human osteoarthritis synovial fibroblasts by inhibiting MiR-381 expression through the PKC, p38, and JNK signaling pathways. Cells 9 (6), 1369. doi:10.3390/cells9061369

CrossRef Full Text | Google Scholar

Cheng, B. F., Lian, J. J., Yang, H. J., Wang, L., Yu, H. H., Bi, J. J., et al. (2020). Neural cell adhesion molecule regulates chondrocyte hypertrophy in chondrogenic differentiation and experimental osteoarthritis. Stem Cells Transl. Med. 9 (2), 273–283. doi:10.1002/sctm.19-0190

CrossRef Full Text | Google Scholar

Chilelli, B., Matthews, G., Minas, T., and Gomoll, A. H. (2024). “Articular cartilage and subchondral bone pathology,” in Biologic knee reconstruction (Oxford, United Kingdom: CRC Press), 9–17.

CrossRef Full Text | Google Scholar

Chu, L., Liu, X., He, Z., Han, X., Yan, M., Qu, X., et al. (2020). Articular cartilage degradation and aberrant subchondral bone remodeling in patients with osteoarthritis and osteoporosis. J. Bone Mineral Res. 35 (3), 505–515. doi:10.1002/jbmr.3909

CrossRef Full Text | Google Scholar

Cooper, M. A., Fehniger, T. A., Turner, S. C., Chen, K. S., Ghaheri, B. A., Ghayur, T., et al. (2001). Human natural killer cells: a unique innate immunoregulatory role for the CD56(bright) subset. Blood 97 (10), 3146–3151. doi:10.1182/blood.v97.10.3146

CrossRef Full Text | Google Scholar

Dabrazhynetskaya, A., Ma, J., Guerreiro-Cacais, A. O., Arany, Z., Rudd, E., Henter, J. I., et al. (2012). Syntaxin 11 marks a distinct intracellular compartment recruited to the immunological synapse of NK cells to colocalize with cytotoxic granules. J. Cell Mol. Med. 16 (1), 129–141. doi:10.1111/j.1582-4934.2011.01280.x

CrossRef Full Text | Google Scholar

D'Orlando, O., Zhao, F., Kasper, B., Orinska, Z., Müller, J., Hermans-Borgmeyer, I., et al. (2013). Syntaxin 11 is required for NK and CD8⁺ T-cell cytotoxicity and neutrophil degranulation. Eur. J. Immunol. 43 (1), 194–208. doi:10.1002/eji.201142343

CrossRef Full Text | Google Scholar

Duruöz, M. T., Öz, N., Gürsoy, D. E., and Gezer, H. H. (2023). Clinical aspects and outcomes in osteoarthritis. Best Pract. and Res. Clin. Rheumatology 37, 101855. doi:10.1016/j.berh.2023.101855

CrossRef Full Text | Google Scholar

Falini, B., and Dillon, R. (2024). Criteria for diagnosis and molecular monitoring of NPM1-Mutated AML. Blood Cancer Discov. 5 (1), 8–20. doi:10.1158/2643-3230.BCD-23-0144

CrossRef Full Text | Google Scholar

Fang, Z., Li, J., Cao, F., and Li, F. (2022). Integration of scRNA-Seq and bulk RNA-seq reveals molecular characterization of the immune microenvironment in acute pancreatitis. Biomolecules 13 (1), 78. doi:10.3390/biom13010078

CrossRef Full Text | Google Scholar

Getz, G. S., Vanderlaan, P. A., and Reardon, C. A. (2007). The immune system and murine atherosclerosis. Curr. Drug Targets 8 (12), 1297–1306. doi:10.2174/138945007783220669

CrossRef Full Text | Google Scholar

Gu, Z. (2022). Complex heatmap visualization. Imeta 1 (3), e43. doi:10.1002/imt2.43

CrossRef Full Text | Google Scholar

Guan, M., Yu, Q., Zhou, G., Wang, Y., Yu, J., Yang, W., et al. (2024). Mechanisms of chondrocyte cell death in osteoarthritis: implications for disease progression and treatment. J. Orthop. Surg. Res. 19 (1), 550. doi:10.1186/s13018-024-05055-6

CrossRef Full Text | Google Scholar

Guan, Z., Jin, X., Guan, Z., Liu, S., Tao, K., and Luo, L. (2023). The gut microbiota metabolite capsiate regulate SLC2A1 expression by targeting HIF-1α to inhibit knee osteoarthritis-induced ferroptosis. Aging Cell 22 (6), e13807. doi:10.1111/acel.13807

CrossRef Full Text | Google Scholar

Hadzic, E., and Beier, F. (2023). Emerging therapeutic targets for osteoarthritis. Expert Opin. Ther. Targets 27 (2), 111–120. doi:10.1080/14728222.2023.2185133

CrossRef Full Text | Google Scholar

Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-Seq data. BMC Bioinforma. 14, 7. doi:10.1186/1471-2105-14-7

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Hu, Y., Gui, Z., Zhou, Y., Xia, L., Lin, K., and Xu, Y. (2019). Quercetin alleviates rat osteoarthritis by inhibiting inflammation and apoptosis of chondrocytes, modulating synovial macrophages polarization to M2 macrophages. Free Radic. Biol. Med. 145, 146–160. doi:10.1016/j.freeradbiomed.2019.09.024

CrossRef Full Text | Google Scholar

Jin, Z., Chang, B., Wei, Y., Yang, Y., Zhang, H., Liu, J., et al. (2022). Curcumin exerts chondroprotective effects against osteoarthritis by promoting AMPK/PINK1/Parkin-mediated mitophagy. Biomed. Pharmacother. 151, 113092. doi:10.1016/j.biopha.2022.113092

CrossRef Full Text | Google Scholar

Kong, R., Ji, L., Pang, Y., Zhao, D., and Gao, J. (2023). Exosomes from osteoarthritic fibroblast-like synoviocytes promote cartilage ferroptosis and damage via delivering microRNA-19b-3p to target SLC7A11 in osteoarthritis. Front. Immunol. 14, 1181156. doi:10.3389/fimmu.2023.1181156

CrossRef Full Text | Google Scholar

Kui, B., Pintér, J., Molontay, R., Nagy, M., Farkas, N., Gede, N., et al. (2022). EASY-APP: an artificial intelligence model and application for early and easy prediction of severity in acute pancreatitis. Clin. Transl. Med. 12 (6), e842. doi:10.1002/ctm2.842

CrossRef Full Text | Google Scholar

Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinforma. 9, 559. doi:10.1186/1471-2105-9-559

CrossRef Full Text | Google Scholar

Li, Z., Feng, J., Zhong, J., Lu, M., Gao, X., and Zhang, Y. (2022). Screening of the key genes and signalling pathways for diabetic nephropathy using bioinformatics analysis. Front. Endocrinol. (Lausanne) 13, 864407. doi:10.3389/fendo.2022.864407

CrossRef Full Text | Google Scholar

Liu, Z., Wang, L., Xing, Q., Liu, X., Hu, Y., Li, W., et al. (2022). Identification of GLS as a cuproptosis-related diagnosis gene in acute myocardial infarction. Front. Cardiovasc Med. 9, 1016081. doi:10.3389/fcvm.2022.1016081

CrossRef Full Text | Google Scholar

Njoto, I., Kalim, H., Soeatmadji, D. W., Handono, K., and Fatchiyah, F. (2019). Effect of hyperglycemia to the mRNA level and protein expression of perlecan at rat model of osteoarthritis with diabetes mellitus type 1. Med. Arch. 73 (3), 144–148. doi:10.5455/medarh.2019.73.144-148

CrossRef Full Text | Google Scholar

Offenhäuser, C., Lei, N., Roy, S., Collins, B. M., Stow, J. L., and Murray, R. Z. (2011). Syntaxin 11 binds Vti1b and regulates late endosome to lysosome fusion in macrophages. Traffic 12 (6), 762–773. doi:10.1111/j.1600-0854.2011.01189.x

CrossRef Full Text | Google Scholar

Pan, X., Jin, X., Wang, J., Hu, Q., and Dai, B. (2021). Placenta inflammation is closely associated with gestational diabetes mellitus. Am. J. Transl. Res. 13 (5), 4068–4079.

Google Scholar

Qiu, S., Munir, A., Malik, S. I., Khan, S., and Hassan, A. (2021). Identification of differentially expressed genes and pathways crosstalk analysis in rheumatoid and osteoarthritis using next-generation sequencing and protein-protein networks. Saudi J. Biol. Sci. 28 (8), 4656–4663. doi:10.1016/j.sjbs.2021.04.076

CrossRef Full Text | Google Scholar

Ramos, Y. F., den Hollander, W., Bovée, J. V., Bomer, N., van der Breggen, R., Lakenberg, N., et al. (2014). Genes involved in the osteoarthritis process identified through genome wide expression analysis in articular cartilage; the RAAK study. PloS one 9 (7), e103056. doi:10.1371/journal.pone.0103056

CrossRef Full Text | Google Scholar

Rios-Arce, N. D., Murugesh, D. K., Hum, N. R., Sebastian, A., Jbeily, E. H., Christiansen, B. A., et al. (2022). Preexisting type 1 diabetes mellitus blunts the development of posttraumatic osteoarthritis. JBMR Plus 6 (5), e10625. doi:10.1002/jbm4.10625

CrossRef Full Text | Google Scholar

Robin, X., Turck, N., Hainard, A., Tiberti, N., Lisacek, F., Sanchez, J. C., et al. (2011). pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinforma. 12, 77. doi:10.1186/1471-2105-12-77

CrossRef Full Text | Google Scholar

Shi, H., Yuan, X., Liu, G., and Fan, W. (2023). Identifying and validating GSTM5 as an immunogenic gene in diabetic foot ulcer using bioinformatics and machine learning. J. Inflamm. Res. 16, 6241–6256. doi:10.2147/JIR.S442388

CrossRef Full Text | Google Scholar

Steinmetz, J. D., Culbreth, G. T., Haile, L. M., Rafferty, Q., Lo, J., Fukutaki, K. G., et al. (2023). Global, regional, and national burden of osteoarthritis, 1990–2020 and projections to 2050: a systematic analysis for the global burden of disease study 2021. Lancet Rheumatology 5 (9), e508–e522. doi:10.1016/S2665-9913(23)00163-7

CrossRef Full Text | Google Scholar

Tan, Z., Chen, X., Zuo, J., Fu, S., Wang, H., and Wang, J. (2023). Comprehensive analysis of scRNA-Seq and bulk RNA-seq reveals dynamic changes in the tumor immune microenvironment of bladder cancer and establishes a prognostic model. J. Transl. Med. 21 (1), 223. doi:10.1186/s12967-023-04056-z

CrossRef Full Text | Google Scholar

Tibor, L. M., and Ganz, R. (2022). “Hip osteoarthritis: definition and etiology,” in Hip arthroscopy and hip joint preservation surgery (Springer), 161–174.

Google Scholar

Wang, H., Yan, Y., Pathak, J. L., Hong, W., Zeng, J., Qian, D., et al. (2023). Quercetin prevents osteoarthritis progression possibly via regulation of local and systemic inflammatory cascades. J. Cell Mol. Med. 27 (4), 515–528. doi:10.1111/jcmm.17672

CrossRef Full Text | Google Scholar

Yang, Y. Y., Gao, Z. X., Mao, Z. H., Liu, D. W., Liu, Z. S., and Wu, P. (2022). Identification of ULK1 as a novel mitophagy-related gene in diabetic nephropathy. Front. Endocrinol. (Lausanne) 13, 1079465. doi:10.3389/fendo.2022.1079465

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Yue, T., Chen, S., Zhu, J., Guo, S., Huang, Z., Wang, P., et al. (2021). The aging-related risk signature in colorectal cancer. Aging (Albany NY) 13 (5), 7330–7349. doi:10.18632/aging.202589

CrossRef Full Text | Google Scholar

Zhai, Y., Li, G., Li, R., Chang, Y., Feng, Y., Wang, D., et al. (2020). Single-cell RNA-sequencing shift in the interaction pattern between glioma stem cells and immune cells during tumorigenesis. Front. Immunol. 11, 581209. doi:10.3389/fimmu.2020.581209

CrossRef Full Text | Google Scholar

Zhang, D., Lu, W., Cui, S., Mei, H., Wu, X., and Zhuo, Z. (2022). Establishment of an ovarian cancer omentum metastasis-related prognostic model by integrated analysis of scRNA-seq and bulk RNA-Seq. J. Ovarian Res. 15 (1), 123. doi:10.1186/s13048-022-01059-0

CrossRef Full Text | Google Scholar

Zhang, L., Zhang, X., Liu, H., Yang, C., Yu, J., Zhao, W., et al. (2024). MTFR2-dependent mitochondrial fission promotes HCC progression. J. Transl. Med. 22 (1), 73. doi:10.1186/s12967-023-04845-6

CrossRef Full Text | Google Scholar

Zhu, R., Wang, Y., Ouyang, Z., Hao, W., Zhou, F., Lin, Y., et al. (2023). Targeting regulated chondrocyte death in osteoarthritis therapy. Biochem. Pharmacol. 215, 115707. doi:10.1016/j.bcp.2023.115707

CrossRef Full Text | Google Scholar