ACADL and ADH1B signify ketone body metabolic reprogramming in osteoarthritic synovium: insights from bioinformatics and animal model studies

TL;DR

This study identifies ACADL and ADH1B as key genes in synovial ketone body metabolism linked to osteoarthritis, offering potential for early diagnosis and targeted therapies.

Contribution

The study discovers ACADL and ADH1B as novel OA biomarkers and proposes progesterone and fomepizole as potential targeted therapies.

Findings

ACADL and ADH1B are associated with synovial inflammation and lipid oxidation in OA.

Progesterone and fomepizole show strong binding affinities to ACADL and ADH1B, respectively.

Mouse models confirmed reduced protein levels in OA synovium.

Abstract

Osteoarthritis (OA) is characterized by articular degeneration and chronic joint pain, partly resulting from synovial inflammation. Accumulating evidence suggests that alterations in the synovial ketone body metabolism (KBM) are closely associated with OA pathogenesis. This study aimed to investigate the metabolic changes in synovial tissues to optimize the treatment of clinical OA.

Analysis of OA and normal control synovial transcriptomic datasets extracted from the Gene Expression Omnibus (GEO) identified 808 differentially expressed genes (DEGs). These DEGs were integrated with KBM-related genes from the metabolic databases, yielding 50 candidates related to OA progression. Following enrichment analysis, protein-protein interaction network construction via STRING, and machine learning with expression analysis, two genes were identified as OA biomarkers: ACADL, encoding long-chain acyl-CoA dehydrogenase and ADH1B, alcohol dehydrogenase 1 B.

The nomogram based on these data revealed high accuracy in the training and validation sets. Functional analysis revealed that these genes function in lipid oxidation, a process critical for synovial cell energy metabolism, as well as in the redox balance that prevents oxidative stress from worsening OA inflammation. Immune infiltration analysis revealed that their expression significantly correlated with 21 immune subtypes, including pro-inflammatory M1 macrophages and Th17 cells, which drive synovial inflammation. Molecular docking analysis identified progesterone and fomepizole as potential agents with satisfactory affinities for ACADL and ADH1B, respectively. Assessment of mouse models of OA confirmed a significant reduction in the synovium at the protein level.

ACADL and ADH1B link KBM abnormalities to immune dysregulation in the OA synovium. The nomogram enables the precise early diagnosis of OA, and progesterone and fomepizole are promising targeted therapies. These findings deepen the current understanding of OA pathogenesis and support the advancement of personalized treatments for clinical translation.

Figures

Candidate gene acquisition of ketone body metabolism in osteoarthritic synovium. (A) Heatmap of top 50 differentially expressed genes in the transcriptomic data from the osteoarthritic and control synovium. The color of the upper density heatmap represents the gene expression density of each sample, with redder colors indicating higher density; in the lower heatmap, the vertical axis represents genes, where purple denotes high expression and yellow denotes low expression. (B) Volcano plot of DEGs between OA and control groups (adjust p < 0.05 and | log2FC| > 1). (C) Venn diagram of the intersection between KBM-related genes and DEGs.

Function and interaction analysis of the candidate genes. (A) The GO enrichment analysis results indicated that the candidate genes are primarily involved in the regulation of systemic metabolism and lipid homeostasis (p < 0.05). (B) KEGG enrichment of the candidate genes (p < 0.05). (C) PPI network of the candidate genes (confidence score > 0.4).

Two machine learning algorithms identified osteoarthritis biomarkers. (A) Based on the 31 feature genes screened by the Boruta algorithm, importance ranking was performed. (B) Based on the five feature genes selected by the XGBoost algorithm, the horizontal axis represents the contribution of each gene to the model’s gain, with larger values indicating a greater contribution to improving model performance. The vertical axis lists the names of the respective genes. (C) The intersection of the two algorithms identified five consensus genes. (D,E) Validations of the consensus genes in two independent cohorts (Wilcoxon rank-sum tests). Data were presented as Median with Q1, Q3 (D,E). NS, not significant; *p < 0.05, p < 0.01, *p < 0.001, p < 0.0001.

The evaluation of clinical predicting value, functional signaling pathways, and molecular crosstalk of ACADL and ADH1B. (A) The nomogram constructed based on the two biomarkers. “Points” indicates the individual scores corresponding to the values of each variable; “Total Points” represent the sum of all individual scores, which corresponds to the “Odds” axis at the bottom and can be used to estimate the predicted probability of the target event (e.g., disease occurrence). (B,C) The predictive performance of the nomogram model was jointly validated by ROC analysis (B) and calibration curves (C), with results demonstrating its strong discriminative ability and calibration accuracy. (D,E) The primary pathways identified in GSEA of ACADL (D) and ADH1B (E). (F,G) GeneMANIA-based interaction networks of ACADL (F) and ADH1B (G).

Immune microenvironment changes in OA synovium. (A) Heatmap of immune cell infiltration. (B) There exists a statistically significant difference in the infiltration levels of 21 immune cell types between the OA group and the control group (Wilcoxon rank-sum tests). (C) Correlation analysis among the 21 differential immune cell types (Spearman correlation analysis). (D) Correlation analysis between the biomarkers and differentially infiltrated immune cells (Spearman correlation analysis). Data were presented as Median with Q1, Q3 (B). NS, not significant; *p < 0.05, p < 0.01, *p < 0.001, p < 0.0001.

Pharmacological prediction and molecular docking of ACADL and ADH1B. (A) Drug-biomarker regulatory network. Orange nodes represent the two biomarkers, while purple nodes represent the predicted drugs. (B,C) Molecular docking simulations of ACADL-progesterone (B) and ADH1B-fomepizole (C).

ACADL and ADH1B expression in murine synovium. (A–D) The representative images of immunohistochemistry staining of ACADL (A) and ADH1B (B) in the synovium of mouse model and statistical analysis. (E) Western blot analysis of ACADL and ADH1B of mouse model. CTRL, control; OA, osteoarthritis. Data were presented as Mean ± SD (B,D,E). *p < 0.05, *p < 0.001, p < 0.0001, non-paired Student t-test. Scale bar = 200μm (left) or 50μm (right).

Tables

TABLE 1. Molecular docking information of ACADL and ADH1B.

Protein	PDB ID	Compound	Compound CID	Docking score
ACADL	8W0T	Progesterone	5,994	−8.8 kcal/mol
ADH1B	1DEH	Fomepizole	3,406	−9.3 kcal/mol

Full text

Introduction

1

Osteoarthritis (OA) is a chronic joint disease that occurs among individuals worldwide. The latest data from the Global Burden of Disease indicate that OA affected 7.6% of the global population in 2020, with its prevalence variably projected to increase by 132.2% over the next three decades and by 60es a% by 2050 (1). OA is one of the most common causes of chronic pain and disability among the elderly (2). Pathologically, OA is characterized by cartilage degradation and synovial inflammation, ultimately leading to joint dysfunction (3, 4). Clinically, management of early-to mid-stage OA includes weight loss, physical therapy, analgesic medications, and chondroprotective agents to alleviate joint symptoms. Conversely, in advanced OA, surgical interventions, such as osteotomy and joint replacement are required (5). Insufficient understanding of the molecular mechanisms underlying OA pathogenesis has rendered disease progression irreversible, with current strategies limited to symptomatic treatments (6). Therefore, exploring the aberrant molecular signatures in OA tissues may facilitate the identification of therapeutic targets that influence disease progression, consequently promoting the development of etiological interventions.

Arthroscopic and imaging examinations of affected patients have revealed that approximately 50% of patients with OA with knee pain presented with synovitis, with its severity positively correlated with the degree of pain (1, 7, 8). Infiltrated macrophages and T lymphocytes can release interleukin-1β, tumor necrosis factor-α, prostaglandin E2, and nerve growth factor. These factors induce peripheral and central sensitization, ultimately contributing to pain (9). Reciprocally, they can interrupt local synovial cells during lipid metabolism, which is involved in the alternative energy supply and inflammatory cytokine synthesis (10, 11), creating a vicious cycle driving disease progression.

Specifically, ketone bodies generated from fatty acid β-oxidation under stress conditions play core roles in both metabolic reprogramming and cellular signaling. Studies have demonstrated that elevated blood ketone body levels ameliorate inflammation and cartilage loss in rat models of OA (12). A recent study further confirmed that β-hydroxybutyrate enhances chondrocyte mitophagy via the HCAR2/AMPK/PINK1/Parkin pathway, thereby reducing chondrocyte senescence and inflammatory cytokine secretion (13). Consequently, we hypothesized that aberrant ketone body metabolism (KBM) in the synovial cells induces inflammation and contributes to joint pain.

This study aimed to integrate transcriptomic datasets, employing bioinformatics approaches to comprehensively identify KBM-associated genes in the synovium of patients with OA, elucidate their biological implications, and validate the findings in animal models. Ultimately, this study provides a novel theoretical foundation for deciphering the KBM-related mechanisms in OA and formulating targeted therapeutic strategies.

Materials and methods

2

Data collection

2.1

Gene expression profiles of normal and OA synovial tissues were obtained from the Gene Expression Omnibus (GEO) database.^1^ The GSE55235 dataset (platform: GPL96) was selected as the training cohort, which comprised transcriptomic data from 10 OA samples and 10 healthy controls. For validation purposes, the GSE82107 dataset (platform: GPL570) containing 10 OA and seven healthy synovial tissue transcriptomes was retrieved from the same repository. KBM-related genes were obtained from the Molecular Signatures Database (MSigDB)^2^ using the keywords “Ketone Body” and “homo sapiens.” After removing duplicates, 485 unique KBM-related genes were retained for subsequent analysis (Supplementary Table 1).

Differential expression analysis

2.2

Identification of DEGs between OA and control samples from the GSE55235 training set was achieved using the “limma” R package (v 3.54.0) (14). Genes with an adjusted p < 0.05 and an absolute log_2_Fold Change (| log_2_FC|) > 1 were considered statistically significant. The Benjamini-Hochberg method was employed for p-value adjustment. A chart displaying the top 10 genes with the highest dysregulation (arranged by decreasing | log_2_FC| values) was generated with the R package “ggplot2” (v 3.4.1) to visualize expression patterns (15). Expression profiles of the top 25 up-regulated and down-regulated DEGs were mapped using the R package “pheatmap” (v 1.0.12) (16).

Candidate genes acquisition and functional enrichment analysis

2.3

For the purpose of identifying genes related to ketone body metabolism, a Venn diagram analysis was performed utilizing the R package “ggvenn” (v 0.1.9) (17) to determine overlapping genes between the two gene sets, and the intersection genes of the two partial genes were defined as candidate genes for subsequent analysis.

To delve into the biological processes and signaling mechanisms that contribute to the advancement of OA disease, enrichment analyses for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were carried out on the candidate genes set employing the R package “clusterProfiler” (v 4.7.1.3) (18). Statistical significance was set at p < 0.05. The top 10 terms with the lowest p-values were depicted using the R package “enrichplot” (v 1.18.3) (19).

Protein-protein interaction network

2.4

To understand the interaction between proteins encoded by candidate genes, PPI network was constructed for candidate genes, candidate genes were inputted into the STRING database^3^ with a minimum interaction confidence score set at 0.4. The PPI network was then constructed using standard parameters and subsequently visualized using Cytoscape software (v 3.8.2) (20).

Feature gene selection using machine learning

2.5

To further refine candidate genes, two independent machine learning algorithms-Boruta analysis and eXtreme Gradient Boosting (XGBoost)-were applied to screen feature genes in the training set (OA vs. Control samples). The Boruta algorithm, a random forest-based global feature selection method, was implemented using the R package “Boruta” (v 8.0.0) (21). This method evaluates the relative importance of all candidate genes by comparing their significance with randomly permuted “shadow features.” Genes with importance scores statistically higher than shadow features (p < 0.05) were retained as significant features. The XGBoost algorithm, an optimized gradient-boosting framework, was employed through the R package “XGBoost” (v 2.1.1.1) (22). Gene importance was ranked based on the “gain” metric, which quantifies the contribution of each gene to model accuracy improvement. The Venn diagram was produced with the aid of the R package “Venn Diagram” (v 1.7.3) (23) to identify overlapping genes derived from Boruta and XGBoost. Genes common to both methods were defined as final feature genes for subsequent validation and functional studies.

Identification of biomarker

2.6

To further pinpoint and confirm the expression of characteristic genes, a Wilcoxon rank test was conducted on disease and control samples in the training set (GSE55235) and the validation set (GSE82107). A two-tailed p < 0.05 was considered statistically significant. The expression levels of the final feature genes (derived from the intersection of Boruta and XGBoost) were compared between OA and control groups. Biomarkers were characterized as genes that display statistically significant differential expression (p < 0.05) and concordant trends across both the training and validation datasets.

Construction and validation of the clinical nomogram

2.7

In order to elevate clinical applicability, a nomogram was developed employing the R package “rms” (v 6.7.0) (24) based on the biomarkers identified in previous analyses. In this model, “Points” represented individual scores assigned to each gene, while “Total Points” denoted the cumulative score derived from all selected genes. The total score was utilized to predict OA prevalence. Calibration curve: the nomogram’s predictive accuracy was evaluated employing the R package “regplot” (v 1.1) (25). A calibration plot comparing predicted probabilities with observed outcomes was generated, where a slope close to 1 indicated high model reliability. Receiver Operating Characteristic curve analysis: the nomogram’s discriminative power was confirmed through the utilization of the R package “pROC” (v 1.18.0) (26). Area under the curve > 0.7 was considered indicative of robust predictive performance.

Functional analysis of biomarkers

2.8

Gene Set Enrichment Analysis (GSEA) was executed to delve into the biological functions of the identified biomarkers in the development of OA. The steps taken in the analysis included: Performing Spearman correlation analysis between each biomarker and all other genes in the OA training set (GSE55235) with the aid of the R package “psych.” (v 2.1.6) (27). Genes were ranked based on correlation coefficients (from highest to lowest). The gene list was analyzed through GSEA employing the R package “clusterProfiler” (version 4.6.2) (28), in conjunction with the KEGG pathway gene set (c2.cp.kegg.v7.0.entrez.gmt) from the MSigDB repository.^4^ Enrichment was deemed significant when it met the criteria of a nominal p < 0.05, a false discovery rate (FDR) < 0.25, and an absolute normalized enrichment score (| NES|) > 1. The five pathways with the highest enrichment scores, were graphically represented using the R package “enrichplot” (v 1.18.3) (19).

In order to further evaluate the interactions between biomarkers and their functionally similar genes and related functions, the GeneMANIA database^5^ was utilized. The analysis followed these steps: biomarkers identified in previous analyses were submitted to GeneMANIA. GeneMANIA’s built-in algorithms assigned weights to interaction types based on their relevance to the input genes, and functionally similar genes were retained for network visualization.

Immune infiltration analysis

2.9

The xCell algorithm was used to assess the infiltration levels of 34 immune cell types in the OA training set (GSE55235) employing the R package “xCell” (v 1.1.0) (29). Heatmaps and boxplots were generated to visualize differential immune infiltration patterns between OA patients and normal controls. Wilcoxon rank-sum tests were conducted employing the R package “stats” (v 4.2.2) (30) to contrast immune cell infiltration scores between the OA and control groups. Immune cell types with significantly altered infiltration levels (p < 0.05) were designated as differentially infiltrated immune cells. Boxplots illustrating these differences were created using the R package “ggplot2” (v 3.4.1) (15). Spearman correlation analysis between differentially infiltrated immune cells was conducted across all training set samples using the R package “psych” (v 2.1.6) (27), the correlation coefficients (| cor| > 0.3) and (p < 0.05) were considered significant. Additionally, the package was utilized to analyze the Spearman correlation between biomarkers and differentially infiltrating immune cells (| cor| > 0.3, p < 0.05), and the correlation results were visualized using a heatmap.

Drug prediction and molecular docking analysis

2.10

To identify potential therapeutic agents targeting the biomarkers for OA treatment, the Drug-Gene Interaction Database^6^ was queried to retrieve drug-gene interactions. Drug-gene interaction networks were constructed and visualized using Cytoscape software (v 3.9.1) (20). Subsequently, the Comparative Toxicogenomics Database^7^ was utilized to obtained drugs associated with OA. Drugs overlapping between Drug-Gene Interaction Database and Comparative Toxicogenomics Database predictions were prioritized. OA-relevant drugs targeting biomarkers were selected for further molecular docking analysis. Their chemical structures were downloaded from the PubChem database.^8^ Molecular docking was performed using CB-Dock2.^9^

Animal model validation

2.11

We adopted surgery-induced OA mice model in this study, which was established by destabilization of the medial meniscus (DMM). 18 mice (C57BL/6, aged 8-week male, weighed ∼20 g) were purchased, acclimatized for 1 week under standard housing conditions (specific pathogen free, 22°C with a 12 h light/dark cycle and free access to standard rodent chow and water), and then randomly divided into two groups. For the OA group, we anesthetized the mice and made a 5 mm longitudinal incision on the medial side of the knee joint, then exposed the joint capsule. Next, we gently everted the patella to visualize the medial meniscus and its tibial attachment. Then we precisely transected the medial meniscal-tibial ligament, confirmed by slight meniscal displacement. Then we closed and disinfected the wound. For the control group, we performed the same procedures except meniscus destabilization. The success of the surgically induced model was confirmed 8 weeks post-surgery using behavioral and histopathological assessment (data not shown) (31).

For immunohistochemistry staining, we harvested 12 knee joints from the OA and control group, respectively. They were decalcified in 10% ethylenediaminetetraacetic acid (pH 7.4) at 37°C for 72 h with daily solution change prior to paraffin embedding. 5μm-thick slides were baked at 60°C for 2 h, deparaffinized in xylene, rehydrated through graded ethanol, and rinsed in distilled water. Antigen retrieval was conducted by microwave heating in 0.01 M citrate buffer (pH 6.0) (boiling for 5 min, simmering for 15 min), followed by cooling and phosphate-buffered saline washes. Endogenous peroxidase was blocked with 3% H^2^O^2^ (15 min), and non-specific binding with 5% bovine serum albumin (30 min). Then the slides were incubated with ACADL polyclonal antibody (Proteintech, Wuhan, China, 17442-1-AP, 1:200) or ADH1B polyclonal antibody (Affinity, Liyang, China, DF12809, 1:50) overnight at 4°C, followed by Goat anti-rabbit IgG H&L (HRP) secondary antibody (Abcam, Cambridge, UK, ab6721, 1:1,000) (60 min at room temperature). Signals were visualized with 3,3’-diaminobenzidine tetrahydrochloride (Beyotime, Shanghai, China, P0203), nuclei counterstained with hematoxylin, and sections dehydrated, cleared, and mounted. Positive cell proportion was calculated using ImageJ software.

For Western blot, we harvested 6 knee joints from the other mice of the two groups. Frozen synovial tissue was first weighed, minced, and homogenized on ice in radio-immunoprecipitation assay buffer lysis buffer containing 1% protease and phosphatase inhibitors; the homogenate was incubated on ice for 30 min, then centrifuged at 12,000 g for 15 min at 4°C to collect the supernatant. Protein concentration was quantified by bicinchoninic acid assay, and equal amounts of protein (30μg) were denatured by boiling with 5 × sodium dodecyl sulfate—polyacrylamide gel electrophoresis loading buffer at 95°C for 10 min. Samples were then subjected to above electrophoresis, followed by transfer to polyvinylidene fluoride membranes. The membranes were blocked with 5% skim milk in Tris-buffered saline with Tween-20 (TBST) for 60 min at room temperature, incubated with ACADL polyclonal antibody (Proteintech, Wuhan, China, 17442-1-AP, 1:1,000) or ADH1B polyclonal antibody (Affinity, Liyang, China, DF12809, 1:2,000) or Vinculin monoclonal antibody (Proteintech, Wuhan, China, 66305-2-Ig, 1:10,000) overnight at 4°C, washed with TBST, then incubated with Goat anti-rabbit IgG H&L (HRP) secondary antibody (Abcam, Cambridge, United Kingdom, ab6721, 1:10,000) or Rabbit anti-mouse IgG H&L (HRP) secondary antibody (Abcam, Cambridge, United Kingdom, ab6728, 1:20,000) for 60 min at room temperature. After final TBST washes, Signals were visualized with 3,3’-diaminobenzidine tetrahydrochloride (Beyotime, Shanghai, China, P0203) to be evaluated.

Statistical analysis

2.12

Bioinformatics analyses were conducted with the aid of the R programming language (v 4.2.2). The Wilcoxon rank sum test or non-paired Student’s t-test were applied to contrast the differences between two groups. Fisher’s exact test was used to compute the significance p-value. p < 0.05 was deemed statistically significant.

Results

3

Acquisition of 50 candidate genes for OA

3.1

Transcriptomic analysis revealed 808 DEGs in OA samples in contrast to controls (| log_2_FC| > 1, adjusted p < 0.05), with 496 genes upregulated and 312 downregulated (Figure 1A). Subsequent intersection analysis with known KBM-related genes revealed 50 differentially expressed KBM-related genes (Figures 1B,C). These overlapping genes demonstrated significant dysregulation patterns, suggesting their potential involvement in OA pathogenesis.

Candidate gene acquisition of ketone body metabolism in osteoarthritic synovium. (A) Heatmap of top 50 differentially expressed genes in the transcriptomic data from the osteoarthritic and control synovium. The color of the upper density heatmap represents the gene expression density of each sample, with redder colors indicating higher density; in the lower heatmap, the vertical axis represents genes, where purple denotes high expression and yellow denotes low expression. (B) Volcano plot of DEGs between OA and control groups (adjust p < 0.05 and | log2FC| > 1). (C) Venn diagram of the intersection between KBM-related genes and DEGs.

Exploration of the functions and pathways associated with the candidate genes

3.2

To elucidate the biological roles of 50 candidate genes in OA pathogenesis, comprehensive functional annotation was performed. GO analysis revealed significant enrichment in 1,078 biological processes, 24 cellular components, and 104 molecular functions (Supplementary Tables 2–4). The top five biological processes terms were cellular ketone metabolic process, regulation of cellular ketone metabolic process, regulation of small molecule metabolic process, regulation of cellular ketone metabolic process, regulation of lipid metabolic process, and fatty acid metabolic process. These findings indicated the candidate genes’ predominant involvement in systemic metabolic regulation and lipid homeostasis.

Cellular components analysis highlighted enrichment in endoplasmic reticulum lumen, membrane raft, membrane microdomain, clathrin-coated vesicle membrane and coated vesicle membrane, suggesting their roles in vesicle-mediated transport and subcellular metabolic compartmentalization. Molecular functions analysis identified DNA-binding transcription factor binding, flavin adenine dinucleotide binding and oxidoreductase activities as dominant features and encoding NAD (P)-dependent enzymes played a critical role in redox balance maintenance during cartilage degeneration (Figure 2A).

Function and interaction analysis of the candidate genes. (A) The GO enrichment analysis results indicated that the candidate genes are primarily involved in the regulation of systemic metabolism and lipid homeostasis (p < 0.05). (B) KEGG enrichment of the candidate genes (p < 0.05). (C) PPI network of the candidate genes (confidence score > 0.4).

KEGG pathway analysis further identified 62 significantly enriched pathways (Supplementary Table 5). The top five pathways, non-alcoholic fatty liver disease, drug metabolism-cytochrome P450, foxO signaling pathway, alcoholic liver disease, and retinol metabolism-highlighted crosstalk between local joint metabolism and systemic energy dysregulation (Figure 2B).

The PPI network comprised 38 nodes (proteins) and 123 edges (interactions) with a confidence threshold of 0.4 (Figure 2C). Notably, Interleukin-6 was found to interact with multiple proteins, including ACADL, EGR1, PPARGC1A, and FOXO3, highlighting its potential central role in the network. The dense interconnectivity indicated robust coordination among these genes. This structural organization suggested a tightly regulated metabolic cascade potentially disrupted in OA pathophysiology.

Identification of ACADL and ADH1B as biomarkers for OA

3.3

Two machine learning algorithms were systematically employed to refine the candidate gene pool. The Boruta algorithm identified 31 feature genes (Figure 3A), while XGBoost identified 5 genes (Figure 3B). Intersection analysis revealed 5 consensus genes (ADM, ADH1B, CBR3, ACADL, DHRS9) shared by both methods (Figure 3C). Subsequent validation in two independent cohorts (GSE55235 and GSE82107) demonstrated significant downregulation of ACADL (p < 0.05) and ADH1B (p < 0.05) in OA samples compared to controls (Figures 3D,E). These findings position ACADL and ADH1B as critical molecular players in OA pathogenesis.

*Two machine learning algorithms identified osteoarthritis biomarkers. (A) Based on the 31 feature genes screened by the Boruta algorithm, importance ranking was performed. (B) Based on the five feature genes selected by the XGBoost algorithm, the horizontal axis represents the contribution of each gene to the model’s gain, with larger values indicating a greater contribution to improving model performance. The vertical axis lists the names of the respective genes. (C) The intersection of the two algorithms identified five consensus genes. (D,E) Validations of the consensus genes in two independent cohorts (Wilcoxon rank-sum tests). Data were presented as Median with Q1, Q3 (D,E). NS, not significant; *p < 0.05, **p < 0.01, ***p < 0.001, ***p < 0.0001.

Clinical value of ACADL and ADH1B

3.4

A dynamic nomogram was developed to quantify OA progression risk, incorporating two biomarkers (ACADL and ADH1B). The model assigned differential point values to distinct expression ranges of each biomarker, enabling cumulative risk calculation through variable summation (Figure 4A). At a total score (total points = 167), the nomogram predicted an increase in OA progression odds (odds = 8,290). Receiver Operating Characteristic curve analysis achieved exceptional discriminative performance (Area under the curve = 0.99), maintaining diagnostic accuracy across validation cohorts (Figure 4B). Calibration curves revealed close alignment between predicted and observed OA probabilities (Hosmer-Lemeshow test, p = 0.229), confirming high prediction accuracy (Figure 4C).

The evaluation of clinical predicting value, functional signaling pathways, and molecular crosstalk of ACADL and ADH1B. (A) The nomogram constructed based on the two biomarkers. “Points” indicates the individual scores corresponding to the values of each variable; “Total Points” represent the sum of all individual scores, which corresponds to the “Odds” axis at the bottom and can be used to estimate the predicted probability of the target event (e.g., disease occurrence). (B,C) The predictive performance of the nomogram model was jointly validated by ROC analysis (B) and calibration curves (C), with results demonstrating its strong discriminative ability and calibration accuracy. (D,E) The primary pathways identified in GSEA of ACADL (D) and ADH1B (E). (F,G) GeneMANIA-based interaction networks of ACADL (F) and ADH1B (G).

Exploration of the functions and pathways associated with ACADL and ADH1B

3.5

GSEA was conducted to clarify the functional roles of ACADL and ADH1B in OA progression. Spearman correlation-based ranking revealed distinct pathway activation patterns for each biomarker. ACADL displayed a notable enrichment in 108 biological processes (p < 0.05, FDR < 0.25, | NES| > 1), as outlined in Supplementary Table 6. The primary pathways identified were organic anion transport, organic acid transport, monocarboxylic acid transport, basement membrane organization, and collagen fibril organization, as depicted in Figure 4D. These findings suggest ACADL’s dual involvement in metabolic transport and extracellular matrix remodeling, processes critical for cartilage homeostasis.

ADH1B demonstrated enrichment in 63 pathways (Supplementary Table 7), predominantly associated with ion homeostasis and neurodevelopmental processes. The most significant pathways included monoatomic ion transport, neuron development, cell morphogenesis, monoatomic cation transport, and ion transmembrane transport (Figure 4E). The enrichment of neurodevelopmental pathways implies potential crosstalk between OA-related metabolic dysregulation and sensory nerve remodeling in joint tissues.

To elucidate the functional cooperativity of ACADL and ADH1B, GeneMANIA-based interaction networks were constructed for both biomarkers (Figures 4F,G). The ACADL-centric network comprised 20 functionally associated genes, with physical interactions dominating the connectivity. Co-expression and predicted interactions constituted secondary connection modes, while co-localization, genetic interactions, pathway sharing, and domain homology contributed minimally. The ADH1B network similarly contained 20 interactors with identical edge-type distribution. Both networks displayed characteristic hub-and-spoke architecture, suggesting these biomarkers serve as critical nodes in their respective biological systems. The high prevalence of physical Interactions implied direct molecular binding between the biomarkers and their partners. These findings collectively highlight the central regulatory roles of ACADL and ADH1B in their corresponding metabolic pathways through extensive molecular interactions.

Correlation between ACADL/ADH1B and immune microenvironment in OA synovium

3.6

Comparative analysis of immune cell composition revealed distinct microenvironment alterations between OA and control groups. Stacked bar plots demonstrated differential distribution patterns of 34 immune cell subtypes across cohorts (Figure 5A). Immune infiltration analysis identified 21 immune cell subtypes with significant differences in infiltration levels between the OA and control groups (p < 0.05) (Figure 5B). Macrophages and B cells exhibited the most pronounced proportional increase in OA samples (p < 0.0001).

*Immune microenvironment changes in OA synovium. (A) Heatmap of immune cell infiltration. (B) There exists a statistically significant difference in the infiltration levels of 21 immune cell types between the OA group and the control group (Wilcoxon rank-sum tests). (C) Correlation analysis among the 21 differential immune cell types (Spearman correlation analysis). (D) Correlation analysis between the biomarkers and differentially infiltrated immune cells (Spearman correlation analysis). Data were presented as Median with Q1, Q3 (B). NS, not significant; *p < 0.05, **p < 0.01, ***p < 0.001, ***p < 0.0001.

To investigate the coordinated regulation of immune components in disease progression, Spearman correlation analysis was conducted to evaluate intercellular interactions among immune cell subtypes. Cross-cell-type correlation analysis revealed a robust positive association (cor = 0.86, p < 0.0001) between monocytes and cDC transcriptional profiles (Figure 5C and Supplementary Table 8), indicating biologically meaningful co-regulation. This strong correlation pattern suggested synchronized activation states and functional synergy between the two cell populations.

Further analysis revealed distinct correlation patterns between key biomarkers (ACADL and ADH1B) and immune cell infiltration. Both genes exhibited predominant negative correlations with most immune cell types (p < 0.01), particularly with macrophages (ACADL: cor = −0.76; ADH1B: cor = −0.77) and CD4^+^ Tem (ACADL: cor = −0.74; ADH1B: cor = −0.75). It is worth noting that paradoxical positive correlations were observed in certain immune cell subsets, such as regulatory T cells (Tregs; ACADL: cor = 0.60; ADH1B: cor = 0.46), basophils (ACADL: cor = 0.48; ADH1B: cor = 0.56) and CD8^+^ naive T cells (ACADL: cor = 0.65; ADH1B: cor = 0.64), with all correlations reaching statistical significance (p < 0.05) (Figure 5D and Supplementary Table 9).

These results collectively suggest that ACADL and ADH1B may function through cell type-specific interaction mechanisms, influencing disease pathogenesis via immune-microenvironment modulation.

Pharmacological prediction and molecular docking analysis

3.7

A total of 192 drugs targeting two biomarkers (ACADL and ADH1B) were identified. Notably, elenium and nickel exhibited dual-targeting effects on both biomarkers. For ADH1B-specific targeting, fomepizole as prioritized, while progesterone was identified as the primary candidate for ADH1B and ACADL modulation (Figure 6A). Molecular docking simulations were performed to evaluate binding affinities between drugs and their targets. Progesterone demonstrated robust binding affinity with ACADL (docking score = −8.8 kcal/mol), while fomepizole exhibited strong interaction with ADH1B (docking score = −9.3 kcal/mol) (Figures 6B,C and Table 1).

Pharmacological prediction and molecular docking of ACADL and ADH1B. (A) Drug-biomarker regulatory network. Orange nodes represent the two biomarkers, while purple nodes represent the predicted drugs. (B,C) Molecular docking simulations of ACADL-progesterone (B) and ADH1B-fomepizole (C).

Animal model validation for Acadl and Adh1b expression

3.8

To confirm the expressional level of Acadl and Adh1b in osteoarthritic synovium, we established post-traumatic OA model in mice by DMM and performed immunohistochemistry staining of their knee joint. The results showed that the positive cell rate of ACADL and ADH1B were significantly lower in OA synovium, respectively (Figures 7A–D). Consistently, ACADL and ADH1B protein levels were lower in OA synovium revealed by Western blot analysis (Figure 7E).

*ACADL and ADH1B expression in murine synovium. (A–D) The representative images of immunohistochemistry staining of ACADL (A) and ADH1B (B) in the synovium of mouse model and statistical analysis. (E) Western blot analysis of ACADL and ADH1B of mouse model. CTRL, control; OA, osteoarthritis. Data were presented as Mean ± SD (B,D,E). *p < 0.05, ***p < 0.001, ***p < 0.0001, non-paired Student t-test. Scale bar = 200μm (left) or 50μm (right).

Discussion

4

OA is known to be closely associated with metabolic disorders. Previous studies have demonstrated that an abnormal KBM in chondrocytes promotes OA progression by disrupting energy homeostasis (32). In this study, for the first time, we integrated multi-omics data and machine learning algorithms to systematically identify the KBM-related genes ACADL and ADH1B in synovial tissues as key biomarkers for OA diagnosis, further validating their decreased expression in OA animal models. The clinical nomogram model constructed using these two genes exhibited excellent predictive efficacy (area under the curve = 0.99), highlighting the potential of these genes as therapeutic targets for OA.

ACADL encodes the rate-limiting enzyme of mitochondrial fatty acid β-oxidation, governing ATP generation to sustain the primary energy supply under stress (33, 34). In our study, both bioinformatic analyses and animal models demonstrated significant downregulation of ACADL expression in the OA synovium, concurrent with lipid droplet accumulation. This is similar to the lipid accumulation observed in chondrocytes (35, 36). PPI network analysis revealed that ACADL formed a metabolic module with HADH, suggesting that ACADL may impair ketone body synthesis under hypoxic conditions by reducing acetyl-CoA production. This could result in a compromised energy supply to synoviocytes, as well as an elevation of lipid byproducts, including inflammatory factors. Consequently, decreased ACADL may contribute to the pathogenesis of pain and joint swelling in patients by enhancing neutrophil chemotaxis and macrophage activation in the synovium, accompanied by the release of inflammatory factors. In the literature, Chung et al. confirmed that ACADL deficiency in type II alveolar epithelial cells inhibits the production of the neutrophil chemokine C-X-C motif chemokine ligand 2 by impairing mitochondrial long-chain fatty acid β-oxidation, consequently attenuating alveolar inflammation in acute lung injury (37). Further, Zhang et al. verified that the marine cyclopeptide Stylissatin A promotes fatty acid accumulation in macrophages by inhibiting ACADL-mediated fatty acid β-oxidation, thereby blocking NF-κB signaling and inducible nitric oxide synthase expression to exert anti-inflammatory effects (38). These contradictory tissue-specific roles of ACADL require further in-depth investigation.

ADH1B, a member of the alcohol dehydrogenase family, is highly expressed in the adipose tissue and exerts antioxidant effects in OA (39, 40). Consistent with this knowledge, we identified ADH1B downregulation as a signature of OA, which may impair the NAD^+^/NADH balance, leading to the accumulation of reactive oxygen species. The accumulation of ADH1B in synoviocytes may activate the NF-κB and MAPK pathways, inducing massive secretion of proinflammatory cytokines and chemokines. This drives the recruitment of immune cells, including neutrophils and macrophages, into the joint cavity, thereby forming a vicious cycle of inflammation. Reactive oxygen species can also directly or indirectly activate matrix metalloproteinases and aggrecanases, which excessively degrade cartilage collagen and proteoglycans, ultimately resulting in joint destruction. Furthermore, a significant negative correlation was identified between ADH1B expression and macrophage infiltration. Importantly, prior research has shown that inhibition of M1 macrophage polarization in the synovium can delay OA progression, possibly due to the regulation of the secretion of inflammatory factors and remodeling of the local immune microenvironment of the joint (41–43). These results indicate that modulation of the local immune microenvironment to restore articular homeostasis is viable.

Acetaminophen is a first-line analgesic for OA. However, its clinical utility is limited by its dose-dependent hepatotoxicity (44). Treatment with fomepizole can mitigate this risk by inhibiting CYP2E1-mediated NAPQI formation and reducing JNK-associated oxidative stress (42). When combined with N-acetylcysteine, this represents a promising strategy for managing acetaminophen overdoses (42). On the other hand, progesterone acts as a potential disease-modifying agent for OA, particularly in postmenopausal women, by promoting chondrocyte survival and exerting anti-inflammatory effects (45). However, their clinical translations require further validation through more rigorous preclinical studies and systematic clinical trials to fully elucidate the associated therapeutic potential and safety profiles.

Collectively, this study proposes that the novel biomarkers ACADL and ADH1B could serve as objective and efficient diagnostic indicators of OA. Candidate compounds, such as progesterone and fomepizole, have laid a theoretical foundation for the development of OA-modifying drugs. Regarding the mechanism, this study further clarified how KBM affects joint homeostasis by clarifying the metabolism-immunity axis of the OA synovium. Future research is needed to explore the role of the abnormal synovium in driving OA development and to evaluate the potential of targeted interventions to alleviate OA synovial metabolic disorders, thereby more effectively relieving patient symptoms and improving prognosis.

Conclusion

5

Bioinformatics analysis and animal experiments confirmed that ACADL and ADH1B are KBM-related genes downregulated in OA synovial tissues. They exhibit excellent clinical predictive value and are significantly associated with the infiltration of various immune cells. Progesterone and fomepizole, as the predicted drugs targeting ACADL and ADH1B respectively, hold potential for OA treatment.