Ziwei Shen
Suming Li
Teng Li1
Lining Wang- 1Department of Orthopedics and Traumatology, Nanjing Hospital of Traditional Chinese Medicine, Nanjing, Jiangsu, China
- 2Department of Orthopedics and Traumatology, School of Integrated Traditional Chinese and Western Medicine, Nanjing University of Chinese Medicine, Nanjing, Jiangsu, China
Background: Bitong Mixture (BM) has shown efficacy in alleviating pain in knee osteoarthritis (OA) in clinical practice; however, the molecular mechanisms underlying its therapeutic effects remain to be fully elucidated. This study aimed to identified BM-related OA biomarkers and explore their functional implications.
Methods: An integrative strategy combining bioinformatics prediction and experimental validation was used. Biomarkers were screened from public OA transcriptomic data using differential expression analysis, network pharmacology, and machine learning. Their functions were explored via enrichment and immune infiltration analyses. Molecular docking predicted interactions between herbal compounds and targets. Single-cell analysis characterized biomarker expression in chondrocyte subsets. A rat OA model and reverse transcription quantitative polymerase chain reaction (RT-qPCR) were employed for in vivo validation.
Results: Bioinformatic prediction identified three potential biomarkers: MMP9, MMP2, and SPP1. They demonstrated certain diagnostic performance for OA and were implicated in pathways related to extracellular matrix organization and immune regulation. Immune analysis revealed significant correlations, notably between MMP2 and activated dendritic cells (cor = 0.66) and between SPP1 and CD4+ central memory T cells (cor = -0.75). Molecular docking suggested strong binding affinity between luteolin (a BM component) and MMP9. Single-cell analysis indicated high expression of these potential biomarkers in hypertrophic chondrocytes, inflammatory chondrocytes, and fibrochondrocytes. In vivo validation confirmed that BM alleviated OA symptoms and histopathological damage in rats. RT-qPCR results showed that BM treatment alleviated the OA-induced upregulation of MMP9, MMP2, and SPP1 expression.
Conclusion: MMP9, MMP2, and SPP1 are potential therapeutic biomarkers for BM in OA. The efficacy of BM may be attributed to its regulation of extracellular matrix remodeling and immune responses, which provides a possible mechanistic explanation for its clinical use.
1 Introduction
Osteoarthritis (OA) is a chronic degenerative joint disorder characterized by progressive cartilage destruction, synovial inflammation, and subchondral bone sclerosis (1). This multifaceted disease represents one of the most prevalent musculoskeletal conditions worldwide, significantly impacting quality of life and imposing substantial socioeconomic burdens on healthcare systems globally (2). The pathogenesis of OA involves complex multifactorial interactions, including aging, mechanical overload, metabolic dysfunction (e.g., obesity), genetic susceptibility, and inflammatory cascades that perpetuate tissue damage (3–5). These interconnected pathways contribute to the heterogeneous nature of OA, making it challenging to develop universally effective therapeutic interventions (6). Epidemiologically, OA affects over 300 million people globally, with a disproportionate impact on postmenopausal women and the elderly population, demonstrating rising prevalence due to aging demographics and increasing obesity rates (7). The disease burden is projected to escalate further as life expectancy increases and sedentary lifestyles become more prevalent. Current therapeutic approaches, including nonsteroidal anti-inflammatory drugs (NSAIDs), intra-articular corticosteroids, hyaluronic acid injections, and joint replacement surgery, primarily focus on symptom management rather than addressing underlying disease mechanisms, thus failing to halt or reverse disease progression (8). Notably, chronic NSAID administration increases risks of peptic ulcers, cardiovascular complications, and renal impairment, while surgical options are limited by prosthesis longevity, postoperative complications, and patient-specific factors (9, 10). These therapeutic limitations underscore the urgent need to identify novel prognostic biomarkers and innovative therapeutic targets to better elucidate OA pathogenesis and provide more effective approaches for clinical diagnosis, prognosis assessment, and targeted treatment strategies.
Traditional Chinese medicine (TCM), a comprehensive holistic therapeutic system encompassing herbal compounds, acupuncture, manual therapies, and dietary interventions, has garnered increasing scientific attention for its multi-target efficacy and remarkably low toxicity profile in managing chronic diseases, including OA (11). This ancient medical paradigm operates on the principle of restoring balance within the body’s energy systems while addressing root causes rather than merely symptomatic relief. TCM formulations, systematically classified into standardized patent medicines and individualized customized decoctions, demonstrate remarkable therapeutic versatility by modulating complex pathophysiological processes, including inflammation, oxidative stress, cartilage metabolism, and immune dysfunction, through sophisticated synergistic interactions of multiple bioactive components (12–14). Recent pharmacological investigations have highlighted the significant therapeutic potential of specific TCM herbs in OA management, with Curcuma longa exhibiting potent anti-inflammatory properties through NF-κB pathway inhibition, and Eucommia ulmoides demonstrating cartilage-protective effects via enhanced collagen synthesis and chondrocyte proliferation (15). These findings validate the traditional empirical knowledge with modern scientific evidence. The Bitong mixture (BM) is a hospital-prepared medicine of Nanjing Hospital of Traditional Chinese Medicine (Approval No.: Su Yao Zhi Zi Z04000864). It has been clinically used for the treatment of osteoarthritis for decades and achieved favorable therapeutic outcomes., It is a classical TCM formula comprising Aconitum carmichaelii, Paeonia lactiflora, and multiple complementary herbs, has demonstrated significant analgesic and anti-arthritic effects in preliminary clinical observations and experimental studies (16–18). Despite these promising therapeutic outcomes, the specific bioactive compounds, precise molecular targets, and underlying mechanisms responsible for BM’s OA improvement remain poorly characterized and inadequately understood. Therefore, systematic investigation of BM’s comprehensive pharmacological basis is imperative to bridge the gap between traditional empirical knowledge and contemporary evidence-based clinical application, facilitating its integration into modern therapeutic protocols.
In this study, We adopted a strategy integrating bioinformatics, network pharmacology, and experimental verification to screen candidate biomarkers related to BM treatment of OA, predict their potential molecular targets and pathways, and conduct preliminary validation through animal models and molecular experiments, thereby providing clues and hypotheses for the subsequent in-depth mechanism study and clinical application of BM.
2 Materials and methods
2.1 Study design and workflow
To identify potential biomarkers and explore the therapeutic mechanisms of BM in OA, we pinpointed candidate biomarkers through data mining, experimentally validated their relevance and BM’s in vivo efficacy, and conducted mechanistic and cellular-scale analysis to further dissect the underlying contexts (Supplementary Figure 1).
2.2 Establishment of osteoarthritis rat model
In this experiment, a total of 15 8-week-old male Sprague-Dawley (SD) rats weighing 250 g were obtained from SPF(Beijing)BIOTECHNOLOGY Co., Ltd. These rats were first subjected to one week of adaptive feeding. The sample size was determined with reference to similar pharmacological studies of OA (19) and based on a sample size calculation using G*Power software, which indicated that a minimum of 4 animals per group was required to achieve statistical significance (effect size d = 2.0, α = 0.05, β = 0.8). To ensure robust statistical power, five rats per group were ultimately adopted in this study. Additionally, the monosodium iodoacetate (MIA) model was selected for its ability to rapidly and reproducibly induce key OA pathological features including cartilage degradation, synovitis, and pain behavior, which is suitable for short-term intervention studies and has been widely used in OA drug screening and mechanistic research (20). Subsequently, 15 rats were selected and randomly divided into 3 groups: the Control group, the OA group, and the OA+Arthralgia Compound group. Except for the rats in the Control group, the rats in other groups were anesthetized with 2.5-3% isoflurane. An 8% sodium iodoacetate solution was prepared (2 mg of sodium iodoacetate was dissolved in 50 μl of normal saline). With the rat’s knee joint slightly flexed, a 1-ml insulin syringe was inserted slightly medial to the patella and punctured obliquely upward and medially. When a distinct sense of giving way was felt, indicating that the needle tip had pierced into the joint cavity, the corresponding volume of sodium iodoacetate solution was injected into the bilateral knee joint cavities of each rat to establish the model. Intervention began 14 days after model establishment. For the rats in the OA + BM group, they were intragastrically administered with the Bitong mixture (BM) (containing various traditional Chinese medicine components such as Caulis Sargentodoxae 1.35g/kg, Radix Saposhnikoviae 0.9g/kg, Processed Aconitum carmichaeli Debx. 0.9g/kg, Processed Aconitum kusnezoffii Reichb. 0.9g/kg, Radix Angelicae Pubescentis 0.9g/kg, Radix Dipsaci 1.08g/kg, Radix Achyranthis Bidentatae 1.35g/kg, Rhizoma Cibotii 0.9g/kg, Rhizoma Polygoni Cuspidati 1.35g/kg, Radix Glycyrrhizae 0.9g/kg, Stir-fried Rhizoma Atractylodis Macrocephalae 1.08g/kg, Semen Coicis 1.35g/kg, Radix Clematidis 0.9g/kg, Pericarpium Citri Reticulatae 0.9g/kg, Cortex Periplocae 0.9g/kg, Eupolyphaga Seu Steleophaga 0.54g/kg, Radix Cynanchi Paniculati 0.9g/kg, Poria 1.35g/kg, Ramulus Cinnamomi 0.9g/kg) at a dose of 19.35 g/kg, which was derived from the clinical human dose (assuming a body weight of 70 kg) using a human-to-rat equivalent dose conversion factor of 6.3 based on body surface area (21), once a day for 14 consecutive days. BM is a standardized hospital preparation produced by Nanjing Hospital of Traditional Chinese Medicine with the approval number [Su Yao Zhi Zi Z04000864]. Its production strictly complies with the quality standards outlined in the Pharmacopoeia of the People’s Republic of China (22) and the Good Preparation Practice for Pharmaceutical Preparations in Medical Institutions (https://www.nmpa.gov.cn/yaopin/ypfgwj/ypfgbmgzh/20010313010101514.html), ensuring regulatory oversight and clinical applicability. To guarantee batch-to-batch consistency, the manufacturing process follows established standardized operating procedures. Key steps include soaking multiple Chinese herbal pieces in water, decoction twice, combining filtrates, allowing sedimentation, concentrating, adding flavoring and preservative agents, filtering, filling, and sterilizing, thereby ensuring both chemical stability and pharmacological consistency between batches. The weight and Lequesne MG score were collected to assess the treatment effect of BM in OA. The articular cartilage tissue samples and serum samples of 15 experimental rats were collected for subsequent experiments. All experimental procedures involving animals were approved by the Institutional Animal Care and Use Committee (IACUC).
2.3 Histological staining
After the articular cartilage tissue of 15 rats was removed, it was decalcified using an ethylenediaminetetraacetic acid (EDTA) decalcification solution. After the decalcification was completed, the tissue was taken out, washed thoroughly with phosphate-buffered saline (PBS), and then immersed in a 4% paraformaldehyde solution for fixation. Subsequently, the tissue was dehydrated successively using 70% ethanol, 80% ethanol, 95% ethanol, and 100% ethanol. After the dehydration steps were finished, the tissue was soaked in a mixture of xylene and absolute ethanol until it became transparent. Finally, the transparent tissue was placed in paraffin for soaking and infiltration.
In this experiment, hematoxylin and eosin (H&E) staining was performed on tissue samples from three groups. The 5-μm sliced tissue was placed in alcohol for flattening. After it was flattened, the sliced tissue on the slide was put into an oven for baking. Once the paraffin melted, the slide was taken out and immersed in xylene for dewaxing. Subsequently, the slide was hydrated successively with 100% alcohol, 95% alcohol, 80% alcohol, and 70% alcohol, and then rinsed with PBS solution. Next, the slide was stained with hematoxylin and then rinsed with distilled water. It was then put into an alcohol-hydrochloric acid solution for differentiation and placed in tap water to make it turn blue again. After that, the slide was stained with eosin and rinsed in distilled water. Subsequently, it was dehydrated successively with 95% alcohol and 100% alcohol, and then soaked in xylene until it became transparent. Finally, the slide was sealed with neutral balsam.
Similar to the H&E staining, the sliced tissue with a 3-μm sliced tissue was stained with safranin. Then it was stained with Weigert staining solution and differentiated with an acidic differentiating solution. Subsequently, it was immersed in fast green staining solution. After the immersion, the slice was washed with a weak acidic solution, and then it was immersed in safranin staining solution again. After that, it was successively soaked in alcohol for dehydration and in xylene until it reached a transparent state. Finally, it was sealed with optical resin.
The histological scoring of all tissue sections was performed under blinded conditions. Specifically, all tissue sections were assigned anonymous numerical codes prior to evaluation. The researcher responsible for the histological assessment was unaware of the group assignments (Control, OA, or OA+BM) of each sample throughout the entire scoring process. The codes were only revealed and linked to the experimental groups after the analysis was completed.
2.4 Data acquisition
In this study, three OA-related human datasets were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds). Specifically, GSE114007 (GPL11154) included 18 normal (control group) and 20 OA (case group) human knee cartilage tissue samples, and GSE169077 (GPL96) comprised 5 normal and 6 OA knee cartilage tissue samples. Additionally, the single-cell dataset GSE255460 (GPL24676) consisted of 16 OA and 3 human knee cartilage tissue samples. Among them, GSE114007 and GSE169077 were used as the training set and validation set, respectively. The classification criteria were relatively consistent across datasets, and sample classification in all datasets adhered to the original studies’ diagnostic criteria. Batch correction was not performed since GSE114007 and GSE169077 were analyzed separately as independent training and validation sets, and GSE255460 was processed independently as a single-cell dataset without merging into a unified expression matrix. The detailed characteristics of these datasets, including species, tissue source, disease stage, and platform, are provided in Supplementary Table 1.
2.5 Differential expression and enrichment analysis
In the human OA transcriptomic dataset GSE114007, differential expression analysis was employed to screen for differentially expressed genes (DEGs) between case and control groups by DEseq2 package (v 1.34.0) (19), the results were presented in volcano maps and heat maps by ggplot2 (v 3.3.6) (20) and pheatmap3 (v 1.1.9) (21) packages, respectively. Briefly, raw read counts were preprocessed by removing genes with an average expression below 0.5. A pseudocount of 1 was added to all counts for subsequent log transformation. Differential expression analysis was then conducted using DESeq2 based on a negative binomial generalized linear model. The design matrix was constructed with the sample condition (OA vs. Control) as the variable. Gene-wise dispersion estimates were calculated and a model was fitted using the DESeq() function. The Benjamini-Hochberg procedure was applied to adjust P-values for multiple testing control. DEGs were identified based on thresholds of |log2 fold change (FC)| ≥ 1 and an adjusted p-value (false discovery rate, FDR) < 0.05 (23).
Meanwhile, drug targets corresponding to active compounds in BM (oral bioavailability (OB) ≧ 30% and drug likeness (DL) ≧ 0.18) (24, 25) were predicted in the Traditional Chinese Medicine Systems Pharmacology (TCMSP) database (https://tcmsp-e.com/tcmsp.php). After that, the drug targets obtained above were intersected with the human OA-related DEGs to further identify candidate genes. Subsequently, the biological functions and pathways associated with candidate genes were explored through the clusterProfiler package (v 4.8.2) (22), which included Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses (adj. p < 0.05) (26).
2.6 Construction of protein-protein interaction network
To further explore the interactions among the human OA-derived candidate genes at the protein level, their PPI network was constructed using the STRING database (http://string.embl.de/) (confidence score > 0.4) (27). Subsequently, the top 1 connectivity gene cluster was selected using the MCODE plugin. Next, hub genes in the top 1 gene cluster were further identified using 3 algorithms (Degree, Betweenness, Closeness) in the cytoHubba plugin. The intersection of hub genes identified by these 3 algorithms was employed to mine candidate key genes.
2.7 Machine learning analysis
Firstly, the human OA-derived candidate key genes were subjected to least absolute shrinkage and selection operator (LASSO) regression analysis according to the glmnet package (v 4.1-2) (23), where 10-fold cross validation was used to identify feature genes 1. Meanwhile, the Boruta algorithm was also used to screen the feature genes 2. The key feature genes were obtained by taking the intersection of feature genes from these two machine learning algorithms. Subsequently, the expression levels of key feature genes were investigated in the human OA dataset GSE114007 and GSE169077 using the Wilcoxon test (p < 0.05) (28). The genes showing significant differences in case and control groups and consistent expression trends across both datasets were considered potential biomarkers. Following this, receiver operating characteristic (ROC) curves for these potential biomarkers were plotted in both datasets using the pROC package (v 1.18.0) (24) to further evaluate their diagnostic performance for OA (area under curve (AUC) values > 0.7) (29).
2.8 Construction of a nomogram
In order to further explore the predictive performance of potential biomarkers for OA patients, a potential biomarker nomogram was constructed, where the occurrence of OA was predicted according to the total points corresponding to each potential biomarker. After that, the accuracy of the nomogram model prediction was further evaluated by calibration curve, decision curve analysis (DCA), and ROC curve.
2.9 Functional analysis of potential biomarkers
In the human OA dataset GSE114007, potential biomarkers related to signaling pathways were explored. Specifically, the correlation of potential biomarkers with other genes was calculated separately and sequenced. After that, the human C2: KEGG gene set (Homo sapiens) was downloaded through the msigdbr package (v 7.5.1) (30) as the background set. Afterwards, the sequenced genes were subjected to gene set enrichment analysis (GSEA) by the GSEA function (adj. p < 0.05). In addition, in order to explore the roles played by potential biomarkers in vivo as well as the associated interacting proteins, GeneMANIA (http://genemania.org) was employed to construct the gene-gene interaction (GGI) network of potential biomarkers, where the top 20 genes and the top 7 signaling pathways were presented.
2.10 Immune-related analyses
To understand whether there was a significant difference in the enrichment of immune cells during the progression of OA compared with normal people, xcell (v 1.1.0) (31) was used to perform immune cell infiltration analysis in the human OA dataset GSE114007. Besides, the infiltration scores of 64 immune cells were obtained, samples that were significantly enriched in immune cells were screened, and all immune cells that were not enriched in the samples were removed. Thereafter, the Wilcoxon test was used to analyze the difference in immune cells enrichment scores between the case and control groups (p < 0.05). Eventually, the relationship between potential biomarkers and 64 immune cells was revealed by the Spearman method.
2.11 Regulatory network analysis
In order to deeply explore the expression regulation mode of potential biomarkers, their regulatory networks were constructed. Firstly, miRNAs that may regulate potential biomarker expression were predicted using the miRDB database (http://www.mirbase.org). Among them, the top 10 miRNAs were used for lncRNA prediction in the Starbase database (http://starbase.sysu.edu.cn)(pancancerNum≧10). After that, the top 5 lncRNAs sorted by pancancerNum associated with each miRNA were used for lncRNA-miRNA-mRNA network construction via Cytoscape (v 3.9.1) (32). Also, the transcription factors (TFs) of the potential biomarkers were predicted in the CistromeDB database (Http://cistrome.org/db).
2.12 Molecular docking
In this study, the herbal-active compounds-gene - gene network of potential biomarkers was constructed via Cytoscape software (v 3.9.1) (32). Besides, the common targets of the compounds were selected for molecular docking. In detail, mRNA-encoded proteins in PDB format files were downloaded from Uniprot Protein Database (https://www.uniprot.org), while SDF format files of components were obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). After that, the drug active compounds and gene active compounds results were pretreated with the CB-DOCK2 tool (33). Finally, the molecular docking of the active compounds of the drug with the mRNA-encoded proteins was performed via AutoDock (v 1.1.2) (34), and the results were visualized via PyMol (v2.4.1) (35).
2.13 Single-cell analysis
After completing the molecular docking analysis, the GSE255460 dataset (human OA single-cell RNA-seq) was further analyzed using the Seurat package for single-cell RNA sequencing data. First, low-quality cells were filtered out based on quality control criteria: only cells with nFeature_RNA between 200 and 4,000, nCount_RNA below 20,000, and percent.mt under 15% were retained for further analysis. The retained data was then normalized using the “LogNormalize” method. After log normalization, the “FindVariableFeatures” function with the “vst” method was used to identify the top 2,000 highly variable genes to optimize computational efficiency. Next, principal component analysis (PCA) was conducted on individual samples using the “RunPCA” function, and sample integration was achieved using the “CCAIntegration” method within the “IntegrateLayers” function. After PCA dimensionality reduction, unsupervised clustering was performed using “FindNeighbors” and “FindClusters” functions in “Seurat” (v 5.0.1) (36) (resolution = 0.5) (37). Subsequently, cell annotation was performed based on the TSNE clustering results and marker genes provided in the literature (38). The distribution of potential biomarkers across cell clusters was visualized. Finally, the “CellChat” (v 1.6.1) (39) was employed to analyze cell-to-cell communication at the molecular level. The number of receptor-ligand pairs was counted, and chord diagrams were employed to visualize the communication and regulatory networks between cells.
2.14 Reverse transcription-quantitative PCR
Because articular cartilage is the primary site of degeneration in osteoarthritis, it was selected for RT-qPCR analysis to directly evaluate gene expression changes at the disease epicenter. Articular cartilage tissues were harvested from the rat models described in Section 2.1. RNAs of 5 control, 5 OA, and 5 treatment (using BM to treat OA) group articular cartilage tissue samples were isolated via TRIzol reagent (Ambion, America). After being detected by electrophoresis and analyzed for purity on NanoPhotometer N50, RNAs were subjected to cDNA synthesis using the SureScript-First-strand-cDNA-synthesis kit (Servicebio, China). Then RT-qPCR was conducted by Universal Blue SYBR Green qPCR Master Mix (Servicebio, China) on a CFX connect RT-qPCR detection system (BIO-RAD, America). Primers of potential biomarkers and internal reference gene (GAPDH) were shown in Supplementary Table S1. Besides, the reaction system and program for reverse transcription and RT-qPCR were shown in Supplementary Tables S2-S5. After RT-qPCR reaction, the 2-ΔΔCт method was employed to calculate relative expression. Eventually, t-test (P < 0.05) (40) was conducted to assess intergroup differences, and Graphpad Prism 5 software (v 8.0) (41) was applied for visualization.
2.15 Statistical analysis
Bioinformatic analyses were conducted in the R program (v 4.2.2). Additionally, employing the Wilcoxon rank sum test or t-test to determine group differences, with a significance threshold defined at P < 0.05.
3 Results
3.1 BM alleviates OA symptoms and joint damage (rat model)
To initially validate the therapeutic efficacy of BM, we assessed its effects on pain behavior and joint histopathology in an MIA-induced rat OA model. In the animal model, the knee joints of the rat samples in the OA group exhibited obvious hyperplasia and inflammatory phenotypes when compared with those in the control group and the OA + BM group (Figure 1A). It is speculated that this was caused by the excessive infiltration of inflammatory cells. In addition, the Lequesne MG scores of the three groups of animal models showed that the Lequesne MG scores reached the highest level from 3 to 7 days after the injection of sodium iodoacetate. This indicated that the damage caused by sodium iodoacetate to the rats was most significant during this 3- to 7-day period. However, the scoring results on the 28th day showed that after the intervention of the OA+BM group, the scores decreased significantly (Figure 1B). These results suggested that BM can effectively alleviate the damage caused by OA. Furthermore, the results of H&E staining and safranin staining showed that, compared with the control group and the OA+BM group, the articular cartilage tissue of the samples in the OA group exhibited hyperplasia, an increase in the number of cell layers, thickening of the fibrous tissue, indicating a large number of infiltrating inflammatory cells, and severe cartilage loss (Figures 1C, D). These results further indicated at the histopathological level that BM could alleviate the symptoms of OA and have a protective effect on joint tissue.
Figure 1. (A) Three groups of rat knee joint photos, showing significant hyperplasia and inflammatory phenotype in the OA group knees. (B) In the line chart of Lequesne MG score of the three groups over time, the score of the 3–7 days OA group was the highest, and the score of the 28 days OA + Bitong mixture group was significantly decreased. (C) H&E staining images of articular cartilage tissues in three groups of rats. Compared with the control group and the OA+BM group, in the OA group, the cartilage tissue proliferates, the number of cell layers increases, and the fibrous tissue thickens. ((1:control, 2:OA+BM, 3:OA). (D) Safranin staining images of articular cartilage tissues in three groups of rats, Compared with the control group and the OA+BM group, Cartilage loss was severe in the OA group(1:control, 2:OA+BM, 3:OA).
3.2 Identification and functional enrichment of 64 candidate genes (human: GSE114007)
Through the identification of candidate genes targeted by BM and functional enrichment analysis, we obtained preliminary insights into their biological relevance to OA. In GSE114007, 2,476 DEGs were identified, comprising 1,329 upregulated and 1,147 downregulated genes (Figures 2A, B), which overlapped with 261 drug targets of BM, further identifying 64 candidate genes (Figure 2C). These candidates were collectively enriched into 1,347 GO terms (55 cellular components (CCs), 112 molecular functions (MFs), and 1,180 biological processes (BPs)), primarily involving functions such as response to oxygen levels, response to decreased oxygen levels, among others (Figure 2D). Additionally, 116 signaling pathways were excavated for candidate genes, such as the TNF signaling pathway, the AGE−RAGE signaling pathway in diabetic complications, and so on (Figure 2E).
Figure 2. (A) Volcano plot of differentially expressed genes in the GSE114007 dataset, showing the log-fold change and adjusted P-values of genes; (B) Heatmap of differentially expressed genes, intuitively reflecting the expression levels of genes in different samples; (C) Venn diagram showing the intersection of differentially expressed genes and the drug targets of Bi Pain Mixture obtained 64 candidate genes. (D) The GO enrichment analysis results for candidate genes, displaying the enriched cellular components, molecular functions, and biological processes related entries. (E) The KEGG enrichment analysis results of candidate genes, displaying the relevant signaling pathways.
3.3 Protein interaction network identifies 7 candidate key genes (human: GSE114007)
Building upon the identified candidate genes, we performed a PPI network analysis to elucidate their interconnectivity and prioritize central regulators within the OA pathological network. In the PPI network, 60 candidate genes formed 640 regulatory relationships. For example, top2 interacted at the protein level with multiple genes such as BIRC, NUF2, and RASSF1 (Figure 3A). Subsequently, the top 1 connectivity gene cluster identified by the MCODE plugin (Figure 3B) was further analyzed using Degree, Betweenness, and Closeness algorithms (Figures 3C–E). The intersection of genes identified by these algorithms yielded 7 candidate key genes, namely MMP9, MMP2, SPP1, PTGS2, SERPINE1, TNF, and FN1 (Figure 3F).
Figure 3. (A) Schematic diagram of the interaction between some genes in the PPI network, showing the interaction between top2 and genes such as BIRC, NUF2, and RASSF1. (B) Gene clusters identified by the MCODE plugin with the highest connectivity. (C–E) Degree, Betweenness, Closeness Algorithm Analysis Results. (F) The intersection of results from 3 algorithms yields 7 candidate key genes.
3.4 Machine learning identifies MMP9, MMP2, and SPP1 as potential biomarkers (human: GSE114007, GSE169077)
To distill the diagnostic relevance of the candidate key genes, we conducted machine learning-based screening for potential biomarkers. In LASSO regression analysis, 6 feature genes 1 (MMP9, MMP2, SPP1, PTGS2, SERPINE1, and FN1) were identified when the error was minimized (lambda.min = 0) (Figure 4A). Simultaneously, 5 feature genes 2 (MMP9, MMP2, SPP1, PTGS2, and FN1) were also mined in Boruta analysis (Figure 4B). The intersection of these feature genes was taken to identify 5 key feature genes (MMP9, MMP2, SPP1, PTGS2, and FN1) (Figure 4C). Among them, MMP9, MMP2, and SPP1 were considered potential biomarkers for OA as they were significantly upregulated in case samples of both GSE114007 and GSE169077 (Figure 4D). Subsequently, these potential biomarkers exhibited certain diagnostic performance for OA, with their area under curve (AUC) values for ROC curves exceeding 0.7 in both datasets (Figure 4E).
Figure 4. (A) LASSO regression analysis results show that the 6 feature genes identified when the error is minimal. (B) Boruta analysis results, displaying the 5 extracted feature genes. (C) Venn diagram, presenting the intersection of LASSO regression and Boruta analysis results yields 5 key feature genes. (D) MMP9, MMP2, and SPP1 expression levels comparison between case and control groups in the GSE114007 and GSE169077 datasets. (E) MMP9, MMP2, and SPP1 in the two datasets, showing their diagnostic performance.
3.5 Diagnostic nomogram based on potential biomarkers (human: GSE114007)
To preliminarily assess the integrated diagnostic potential of MMP9, MMP2, and SPP1 for OA, we constructed an exploratory nomogram based on the GSE114007 dataset (Figure 5A). Afterwards, its predictive performance was further evaluated. The calibration curve showed C-index of 0.861, which exceeded 0.7 (Figure 5B). In DCA analysis, the net benefit of the nomogram was greater than that of individual potential biomarkers (Figure 5C). Moreover, the AUC value of the ROC curve was 0.861, also exceeding 0.7 (Figure 5D). These results consistently indicated that the nomogram of potential biomarkers had certain predictive performance for OA. However, as this nomogram is an exploratory tool constructed solely on transcriptomic data, its clinical applicability requires further validation in future prospective cohorts that incorporate clinical variables.
Figure 5. (A) Potential biomarker nomogram, showing the score axes corresponding to the three potential biomarkers MMP9, MMP2 and SPP1, as well as the relationship between the total score and the probability of osteoarthritis occurrence; (B) Calibration curve, showing that the C index is 0.861; (C) The DCA analysis results indicated that the net benefit of the nomogram model was greater than that of a single biomarker; (D) The ROC curve shows that the AUC value is 0.861.
3.6 Signaling pathways and gene networks associated with potential biomarkers (human: GSE114007)
Elucidating the broader functional context of these potential biomarkers, we performed pathway and network enrichment analyses. GSEA was employed to further provide insight into the function of potential biomarkers. The results indicated that all biomarkers were collectively involved in systemic lupus erythematosus. Additionally, MMP9 and MMP2 were enriched together in the ribosome pathway. MMP9 and SPP1 were significantly enriched in the spliceosome and neuroactive ligand receptor interaction pathways (Figures 6A–C). Subsequently, the top 20 genes associated with the potential biomarkers (e.g., TIMP1, ITGA5) and seven signaling pathways were identified in the GGI network. These pathways included extracellular matrix organization, collagen metabolic process, and others (Figure 6D).
Figure 6. (A–C) GSEA analysis showed that potential biomarkers were jointly involved in the signaling pathways of systemic lupus erythematosus; MMP9 and MMP2 are enriched in the ribosomal signaling pathway; MMP9 and SPP1 are enriched in the spliceosome and neuroactive ligand-receptor interaction signaling pathways. (D) The top 20 genes and 7 signaling pathways related to potential biomarkers in the GGI network.
3.7 Altered immune infiltration in OA and correlation with potential biomarkers (human: GSE114007)
To elucidate the immunological context in which these potential biomarkers operate, we assessed immune cell infiltration patterns in OA. In this study, the infiltration of 64 immune cells was explored according to xcell (Figure 7A), of which, 22 immune cells, such as adipocytes, astrocytes, had significantly higher infiltration scores in case samples, whereas the remaining 8 immune cells (e.g., basophils, mast cells), were the opposite (Figure 7B). Furthermore, in the correlation analysis of potential biomarkers with 64 immune cells, MMP2 showed the largest positive correlation with aDC (cor=0.66, p ≤ 0.0001), while SPP1 showed the largest negative correlation with CD4 Tcm (cor = -0.75, p ≤ 0.0001) (Figure 7C).
Figure 7. (A) The results of xcell analysis on the infiltration of 64 types of immune cells are presented; (B) Comparison of 22 immune cells with higher infiltration scores and 8 immune cells with lower infiltration scores; (C) In the correlation analysis of potential biomarkers and 64 types of immune cells, the correlations between MMP2 and aDC, and SPP1 and CD4_Tcm were displayed.
3.8 Regulatory networks and molecular docking predict interactions between BM compounds and potential biomarkers (database prediction)
Building upon these potential biomarkers, we explored their upstream regulatory networks and performed molecular docking to predict their mechanistic roles in OA. In the lncRNA-miRNA-mRNA network of potential biomarkers, 25 miRNAs, 24 lncRNAs, and 3 potential biomarkers were included. They formed regulatory relationships included PVT1 - has-miR-181c-5p - SPP1, LINC02863 - has-miR-149-5Pp - MMP9 and so on (Figure 8A). After that, TFs corresponding to potential biomarkers were also predicted, and the results showed that FOS, EP300, and TRIM28 were the most critical TFs for MMP9, MMP2, and SPP1, respectively (Figures 8B–D). Furthermore, 4 herbal (Jixueteng, Chenpi, Gancao, and Huzhang), and 3 active ingredients (luteolin, quercetin, and nobiletin) were used to construct the herbal-active compounds-gene network with potential biomarkers. The relationship pairs formed included Jixueteng-luteolin-MMP9, Huzhang-quercetin-SPP1, and so on (Figure 8E). Finally, these active compounds were applied for molecular docking with potential biomarkers to predict the binding affinity, in which luteolin demonstrated the most favorable docking effect with MMP9 in the simulation, and its predicted docking binding energy was -10.1 kcal/mol (Supplementary Figure 2).
Figure 8. (A) The lncRNA-miRNA-mRNA network shows the regulatory relationships among 25 mirnas, 24 lncrnas and 3 potential biomarkers; (B–D) The prediction results of the key transcription factors of MMP9, MMP2 and SPP1; (E) The figure shows the herb - active compound - gene network, presenting the pairs of relationships among four herbs, three active components and potential biomarkers.
3.9 Single-cell profiling of OA cartilage identifies key cell types (human: GSE255460)
We next used scRNA-seq to resolve the cellular expression patterns of these potential biomarkers and the associated changes in the joint communication network during OA. Moving on to the single-cell RNA sequencing analysis, the distribution of nFeature RNA (the number of genes detected per cell), nCount RNA (the total gene expression per cell), and percent.mt (the proportion of mitochondrial genes per cell) in the quality-controlled samples was visualized using violin plots (Figure 9A). In addition, highly variable genes within the GSE255460 dataset were identified and displayed (Figure 9B). To explore the heterogeneity among cells, PCA was conducted (Figure 9C), while the integrated data indicated that cells were evenly mixed without significant outliers. Subsequently, the top 30 principal components were selected for further analysis (Figure 9D). Before cell annotation, TSNE clustering analysis was performed, identifying a total of 14 cell clusters (Figure 10A). The subsequent cell annotation identified 10 cell types: regulatory chondrocytes (RegC), proliferative chondrocytes (ProC), homeostatic chondrocytes (HomC), prehypertrophic chondrocytes (preHTC), effector chondrocytes (EC), fibrochondrocytes (FC), hypertrophic chondrocytes (HTC), pre-fibrochondrocytes (preFC), pre-inflammatory chondrocytes (preInfC), and inflammatory chondrocytes (InfC) (Figure 10B). Annotation results revealed that the SPP1 gene is highly expressed in HTC and InfC cells, MMP2 is abundantly expressed in FC cells, and MMP9 is relatively enriched in InfC cells. Thus, HTC, InfC, and FC were selected as key cell types, suggesting that these cells may play significant roles in the drug effects on OA (Figure 10C). Moreover, the analysis of cell-cell interactions revealed that in the control group, HTC cells played a key role in communication with other cells. However, in osteoarthritis, their interactions weakened, except with preHTC and FC cells, where communication remained stable. Meanwhile, IfnC cells showed no significant changes between the groups. In contrast, FC cells had stronger communication with preHTC cells and weaker interactions with EC and RegC cells in osteoarthritis, suggesting that preHTC cells may be crucial in the progression of the disease (Figure 10D). In contrast, the interaction patterns of InfC cells did not change notably between the control and disease groups. Additionally, the interactions between FC and preHTC cells were strengthened in the disease group, whereas their interactions with EC and RegC were weakened. These findings suggest that pre-HTC cells may also play a key role in the progression of osteoarthritis.
Figure 9. (A) Violin plot, showing the distribution of nFeature RNA, nCount RNA and percent.mt in the samples after quality control; (B) Display diagram of highly variable genes; (C) Principal Component analysis (PCA) plot, showing the heterogeneity between cells; (D) The analysis graph of the integrated data shows that the cell distribution is uniform without obvious outliers.
Figure 10. (A) The T-distribution random neighborhood embedding (tSNE) clustering analysis diagram identified 14 cell clusters; (B) The results of cell annotation determined 10 types of cells; (C) The distribution of potential biomarkers in different cell types; (D) The results of intercellular communication analysis show the communication changes among different cell types.
3.10 BM downregulates potential biomarker expression (rat model)
RT-qPCR analysis was subsequently performed in the rat model to assess the expression of these biomarkers under BM intervention. Compared to the control group, the expression of MMP9, MMP2, and SPP1 in OA group was significantly upregulated (P < 0.01). Their expression trends were consistent with those in GSE114007 and GSE169077 datasets, revealing the reliability of aforementioned results. Besides, compared to the OA group, MMP9, MMP2, and SPP1 were significantly downregulated in the treatment group (P < 0.01). Notably, the expression level in the treatment group was closer to control group, implying that BM might have a corrective or normalizing effect on the pathological processes associated with OA (Figure 11).
Figure 11. It is a bar chart of the relative expression levels of MMP2, MMP9, and SPP1 in the control group, OA group, and treatment group. The expression in the OA group is significantly higher than that in the control group, and the expression in the treatment group is significantly lower than that in the OA group and close to the control group.**: p < 0.01, ***: < 0.001.
4 Discussion
OA is a degenerative joint disease characterized by cartilage degradation, synovial inflammation, osteophyte formation, and pain. Current therapies focus on symptom management, but targeted interventions remain limited (42). The Bitong mixture (BM), a traditional herbal formula, has shown empirical efficacy in alleviating OA symptoms, though its molecular mechanisms are poorly understood (16–18). This study systematically explored the potential mechanisms and biological markers of BM in treating OA by integrating bioinformatics, animal experiments, and computational simulations. First, MMP9, MMP2, and SPP1 were screened as BM-related candidate markers based on human transcriptome data; subsequently, a rat OA model was constructed to preliminarily verify that BM could alleviate joint pathology and downregulate the expression of biological markers, transitioning from correlation analysis based on human data to validation through animal experiments. Network pharmacology and molecular docking suggested that the active components of BM (such as luteolin) may directly act on these targets. Single-cell analysis further revealed that these markers are primarily enriched in hypertrophic and inflammatory chondrocyte subpopulations, indicating that BM may influence matrix metabolism and the immune microenvironment by regulating these cell functions. In summary, this study suggests that BM may alleviate the disease progression of OA by regulating MMP9, MMP2, and SPP1 through a synergistic effect of multiple components, thereby affecting extracellular matrix metabolism and the immune microenvironment. These results provide preliminary evidence and direction for further elucidating its mechanisms of action and exploring potential clinical applications.
4.1 The role and mechanism of potential biomarkers in osteoarthritis
Matrix metalloproteinase 9 (MMP9) belongs to the matrix metalloproteinase family. Its main function is to degrade extracellular matrix components such as collagen and gelatin (43).Beyond ECM remodeling, MMP9 participates in inflammatory cascades by facilitating leukocyte infiltration and amplifying pro-inflammatory cytokine release (e.g., IL-1β, TNF-α), which perpetuates synovial inflammation and joint damage (44). Under inflammatory stimulation, the expression of MMP9 will be induced to increase, forming a vicious cycle and exacerbating the development of osteoarthritis (45). In addition, MMP9 can regulate osteogenesis/osteolysis by inhibiting bone resorption and promoting bone formation (46). In osteoarthritis, the expression of MMP9 is significantly upregulated (45), and its expression level is positively correlated with the severity of the disease. This indicates that MMP9 plays a key role in the pathogenesis of osteoarthritis and can be used as an important indicator for evaluating the condition of osteoarthritis (47, 48). Its potential mechanisms may involve degrading articular cartilage, promoting inflammatory responses and the release of inflammatory mediators, while playing a significant role in osteophyte formation.
Matrix metalloproteinase 2 (MMP2), a crucial member of the MMP family, degrades extracellular matrix components including type IV collagen and gelatin (49, 50). In OA, MMP2 expression is markedly upregulated from early disease stages, with levels escalating in correlation with radiographic and histological severity scores (e.g., OARSI grading). Studies demonstrate its capacity to degrade cartilage and synovial ECM, disrupting joint integrity and function. This elevation from early OA stages suggests MMP2’s promotive role in disease progression, highlighting its potential as a potential diagnostic biomarker and therapeutic target (51–53). Moreover, MMP2 regulates cellular proliferation, migration, and differentiation, potentially mediating OA pathology by modulating chondrocyte and synoviocyte biological behaviors.
Secreted phosphoprotein 1 (SPP1/osteopontin), a multifunctional extracellular matrix phosphoglycoprotein involved in cell adhesion, migration, and signaling, shows significant upregulation in OA (54). Research reveals SPP1 enhances inflammatory cell adhesion/migration, amplifies inflammatory responses, and modulates chondrocyte proliferation/differentiation while suppressing cartilage matrix synthesis. Elevated SPP1 levels in OA patients’ articular cartilage and synovial tissue correlate strongly with inflammatory severity and cartilage degeneration, indicating its pivotal role in OA pathogenesis through inflammation regulation and chondrocyte functional modulation (55). Additionally, SPP1’s involvement in bone metabolism suggests its potential contribution to osteophyte formation and joint deformity via bone remodeling alterations in OA.
Notably, the upregulation of these potential biomarkers in OA patients and their normalization following BM intervention suggest that BM may inhibit MMP activity and osteopontin-mediated pathological calcification, thereby preserving cartilage integrity (56). These potential biomarkers show preliminary diagnostic efficacy (AUC > 0.7) in different data sets, suggesting they have good diagnostic potential, and future efforts can focus on further validating their clinical application value. The discovery of these potential biomarkers not only serves as an early assessment indicator with potential value in OA diagnosis but also provides a theoretical basis for developing new treatment strategies based on BM.
4.2 Pathway analysis reveals key BM-regulated mechanisms
The GSEA results demonstrated that all three potential biomarkers are collectively involved in systemic lupus erythematosus (SLE), suggesting autoimmune involvement in OA pathogenesis. In OA patients, synovial immune dysregulation resembling SLE may trigger inflammatory cytokine overproduction (e.g., IL-6, TNF-α) and cartilage breakdown. Autoimmune components may promote synovial immune cell infiltration, accelerating joint destruction through cytokine-mediated pathways (57, 58), highlighting immunomodulation as a potential OA therapeutic target.
MMP9 and MMP2 co-enrich in the ribosome pathway, a central protein synthesis site crucial for cellular homeostasis. This enrichment implies dual regulatory mechanisms in osteoarthritis: 1) MMP9/MMP2 may modulate ribosomal biogenesis and function, potentially influencing ribosome assembly or mRNA translation efficiency; 2) Their expression could conversely depend on ribosomal signaling. In OA chondrocytes, ribosomal dysfunction may impair synthesis of cartilage matrix proteins like collagen, reducing tissue elasticity and accelerating degradation through compromised structural integrity (59, 60).
MMP9 and SPP1 exhibited dual-pathway enrichment in spliceosome function and neuroactive ligand-receptor interactions. The spliceosome pathway enrichment suggests their potential regulatory roles in mRNA splicing, where aberrant processes in OA could lead to chondrocyte dysfunction through impaired metabolic/proliferative/differentiation pathways (61, 62). Concurrently, their neuroactive ligand-receptor association implies possible involvement in pain modulation, particularly through regulating OA-elevated nociceptive mediators like substance P and CGRP (63, 64). These findings propose that MMP9/SPP1 may mediate OA progression through mRNA processing accuracy and pain signal transduction, offering potential therapeutic targets for both cartilage degeneration and pain management in OA (65).
In summary, the signaling pathways associated with these potential biomarkers encompass immune regulation, protein synthesis, gene expression regulation, and neural signal transduction at multiple levels. This not only provides a new perspective for understanding the pathogenic mechanism of OA but also offers potential clues for the mechanistic research of BM in treating OA.
4.3 Association analysis of immune cells and potential biomarkers
Using xCell, we conducted immune cell infiltration analysis on the GSE114007 dataset to systematically profile 64 immune cell types in OA pathogenesis. Significant infiltration differences emerged between OA cases and controls, with 22 cell types—including adipocytes and astrocytes—exhibiting higher infiltration scores in OA samples. Enhanced adipocyte infiltration may modify inflammatory microenvironments through fat-derived mediators like leptin (pro-inflammatory) and adiponectin, which regulate immune responses (66). For astrocytes, their increased presence suggests neuro-immune crosstalk in OA, potentially mediated by neurotransmitters and cytokines influencing immune cell activity and inflammation progression (67).
The infiltration scores of eight immune cells, including basophils and mast cells, were significantly lower in OA samples than in controls. Basophils and mast cells regulate allergic responses and inflammation by releasing histamine, leukotrienes, and other mediators (68). Their reduced infiltration in OA may impair inflammation regulation, exacerbating uncontrolled inflammatory responses and disease progression.
To explore potential relationships between potential biomarkers and immune cells, Spearman correlation analysis was conducted. Results revealed the strongest positive correlation between MMP2 and activated dendritic cells (aDCs) (cor = 0.66, P ≤ 0.0001). As crucial antigen-presenting cells, aDCs process antigens and activate T cells to initiate immune responses. The MMP2-aDC correlation suggests MMP2 may interact with aDCs to modulate immune regulation and inflammatory responses in osteoarthritis pathogenesis. Functionally, MMP2 degrades extracellular matrix components, potentially releasing cartilage breakdown products that activate aDCs. Conversely, cytokines secreted by activated aDCs might upregulate MMP2 expression, forming a positive feedback loop that exacerbates inflammation and joint damage (69–71).
SPP1 exhibited the strongest negative correlation with CD4_Tcm (central memory CD4+ T cells) (cor = -0.75, P ≤ 0.0001). CD4_Tcm plays crucial roles in maintaining immune memory and regulating immune responses (72). The strong negative correlation suggests that SPP1 might impair immune memory maintenance and immune balance through suppressing CD4_Tcm functions, potentially driving OA pathogenesis. The encoded osteopontin (OPN) interacts with multiple cell surface receptors to regulate cell adhesion, migration, and signaling (73, 74). In summary, through immune infiltration analysis, we speculate that BM may regulate potential biomarkers to inhibit the excessive activation of dendritic cells and reduce the functional suppression of CD4+ central memory T cells, alleviating the vicious cycle of inflammation and tissue destruction in osteoarthritis, which helps restore local immune homeostasis in the joints and ultimately exerts a therapeutic effect in alleviating the progression of OA. However, the specific related effects still need further verification through subsequent experiments (75, 76).
4.4 Analysis of molecular docking results and mechanisms
In molecular docking analysis, we constructed an herb-active compound-gene network linking potential biomarkers and herbal active compounds to investigate their interaction mechanisms. The network identified four herbs (Jixueteng, Chenpi, Gancao, and Huzhang) and three active compounds (luteolin, quercetin, nobiletin) associated with the potential biomarkers.
Among them, luteolin demonstrated the strongest binding affinity with MMP9 at −10.1 kcal/mol. Binding energy (negative values indicate stability) revealed their high-affinity interaction and capability to form a stable complex. The excellent binding affinity between luteolin and MMP9 underscores luteolin’s potential as a direct MMP9 inhibitor. Studies have shown that luteolin exhibits multiple biological activities including anti-inflammatory, antioxidant, and anti-tumor effects (77–79). In osteoarthritis treatment, luteolin may competitively block MMP9’s catalytic site, thereby reducing collagenase activity and cartilage degradation (80). This aligns with BM’s observed efficacy in reversing OA potential biomarker expression toward normal levels.
Additionally, luteolin may exert therapeutic effects on osteoarthritis by modulating other signaling pathways. Studies indicate that luteolin inhibits nuclear factor-κB (NF-κB) pathway activation, reducing the expression and release of inflammatory factors (81–83). In osteoarthritis, inflammatory response significantly contributes to joint damage. NF-κB pathway activation promotes production of inflammatory cytokines like iIL-1β and TNF-α, exacerbating joint inflammation and cartilage degradation (84). By suppressing the NF-κB pathway, luteolin likely alleviates inflammatory responses, thereby demonstrating therapeutic potential for osteoarthritis.
Studies indicate that luteolin promotes chondrocyte proliferation and differentiation while inhibiting apoptosis. In osteoarthritis, chondrocyte damage and apoptosis reduce cartilage matrix synthesis. By enhancing chondrocyte proliferation and differentiation, luteolin may increase cartilage matrix production, aiding the repair of damaged articular cartilage (85, 86).
In summary, molecular docking simulation shows that there is a potential high-affinity binding pattern between luteolin and MMP9, suggesting that it may play a potential role in alleviating extracellular matrix degradation, regulating inflammatory response, and promoting cartilage homeostasis by modulating MMP9 activity and downstream related signaling pathways. This provides a preliminary theoretical basis for the hypothesis of “BM regulating the improvement of OA pathological process through potential biomarkers” in this study, and future exploration of the mechanism of action of BM and the development of new drugs can be based on this. However, since molecular docking analysis is only a predictive study, its results need to be further validated through functional experiments.
4.5 Conclusion and limitations
This study first systematically explored the possible mechanisms of BM in OA by integrating multi-omics analysis and experimental validation. The results showed that BM may affect the OA process and produce therapeutic effects by regulating extracellular matrix remodeling and immune microenvironment modulation mediated by MMP9, MMP2, and SPP1, providing references and new insights for understanding the potential molecular mechanisms of traditional Chinese medicine compound treatment for OA. Although this study integrates bioinformatics and animal experiments, there are still the following limitations: (1) The sample size, heterogeneity, and late-stage OA single-cell data of public datasets limit the understanding of the disease as a whole, and there is a risk of model overfitting; (2) The MIA model has the defect of not being able to fully simulate the complex disease course of human OA, and currently lacks solvent and positive drug controls; (3) The findings of this study are mainly based on computational predictions and mRNA levels, lacking direct evidence, protein levels, and functional experimental validation in OA patients. Based on this, future research should deepen in the following directions: (1) Conduct larger-scale prospective clinical studies to systematically collect a series of samples from OA patients receiving BM treatment, and validate the function of markers at the species level; (2) Supplement solvent and positive drug controls, and validate the role of markers and key interaction mechanisms at the protein and functional levels through experiments such as Western blotting, immunohistochemistry, and gene knockdown/overexpression; (3) Utilize technologies such as organoids and spatial transcriptomics to elucidate the regulatory network of BM on chondrocyte subpopulations; (4) Integrate computational and experimental approaches to clarify the compound material basis and optimize its design, ultimately achieving the transformation towards evidence-based precision treatment.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics statement
The animal study was approved by The Ethics Review Committee of Nanjing Municipal Hospital of Traditional Chinese Medicine. The study was conducted in accordance with the local legislation and institutional requirements.
Author contributions
ZS: Writing – review & editing, Writing – original draft. SL: Supervision, Conceptualization, Writing – review & editing, Methodology, Funding acquisition, Resources. TL: Validation, Visualization, Writing – review & editing, Investigation, Software. LW: Validation, Formal analysis, Data curation, Writing – review & editing.
Funding
The author(s) declared that financial support was received for this work and/or its publication. This work was supported by the Science and Technology Development Fund of Nanjing Hospital of Traditional Chinese Medicine (No. 0030402).
Acknowledgments
We would like to express our sincere gratitude to all individuals and organizations who supported and assisted us throughout this research. Special thanks to the following authors: Teng Li, Lining Wang. In conclusion, we extend our thanks to everyone who has supported and assisted us along the way. Without your support, this research would not have been possible.
Conflict of interest
The author(s) declared that this work 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) declared that generative AI was not used in the creation of this manuscript.
Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.
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/fimmu.2026.1739355/full#supplementary-material
Supplementary Figure 1 | Schematic overview illustrating the integrated study design for investigating BM in OA.
Supplementary Figure 2 | The figure shows the visual display of the docking result between luteolin and MMP9 molecules, indicating that the docking binding energy is -10.1 kcal/mol. It should be noted that this analysis is predictive and has not been functionally validated; therefore, these findings are hypothesis-generating and intended only for guidance in future mechanistic research.
Supplementary Table 1 | Summary of the GEO datasets used for bioinformatics analysis, including dataset accession numbers, sample labels, disease status, species, tissue source, and technology platform.
References
1. Xu Y, Wang JY, Meng T, Ma XW, Li H, and Li K. Role of hydrogels in osteoarthritis: A comprehensive review. Int J Rheum Dis. (2023) 26:2390–401. doi: 10.1111/1756-185X.14968
2. GBD 2021 Osteoarthritis Collaborators. 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 Rheumatol. (2023) 5:e508–22. doi: 10.1016/S2665-9913(23)00163-7
3. Xi Y, Wang Z, Wei Y, NI QX, Li D, Ting Z, et al. Gut microbiota and osteoarthritis: from pathogenesis to novel therapeutic opportunities. Am J Chin Med. (2025) 53:43–66. doi: 10.1142/S0192415X2550003X
4. Sampath SJP, Venkatesan V, Ghosh S, and Kotikalapudi N. Obesity, metabolic syndrome, and osteoarthritis-an updated review. Curr Obes Rep. (2023) 12:308–31. doi: 10.1007/s13679-023-00520-5
5. Nedunchezhiyan U, Varughese I, Sun AR, Wu X, Crawford R, and Prasadam I. Obesity, inflammation, and immune system in osteoarthritis. Front Immunol. 13:907750. doi: 10.3389/fimmu.2022.907750
6. Hu J and Wang J. Macrophages: pivotal orchestrators and emerging therapeutic targets in osteoarthritis pathogenesis. J Transl Med. (2025) 23:1250. doi: 10.1186/s12967-025-07340-2
7. Foster NE, Eriksson L, Deveza L, and Hall M. Osteoarthritis year in review 2022: Epidemiology & therapy. Osteoarthritis Cartilage. (2023) 31:876–83. doi: 10.1016/j.joca.2023.03.008
8. Motaung A, Rants’o TA, Walvekar P, Luliński P, and Choonara YE. Recent advances in intra-articular bioactive delivery for the treatment of osteoarthritis. Int J Pharm. 687:126401. doi: 10.1016/j.ijpharm.2025.126401
9. Popelka S, Sosna A, Vavřík P, Jahoda D, Barták V, and Landor I. Eleven-year experience with total ankle arthroplasty. Acta Chir Orthop Traumatol Cech. (2016) 83:74–83. doi: 10.55095/achot2016/012
10. Zzhang ZH, Qi YS, Wei BG, Bao HRC, and Xu YS. Application strategy of finite element analysis in artificial knee arthroplasty. Front Bioeng Biotechnol. (2023) 17:1127289. doi: 10.3389/fbioe.2023.1127289
11. Valsamidou E, Gioxari A, Amerikanou C, Zoumpoulakis P, Skarpas G, and Kaliora AC. Dietary interventions with polyphenols in osteoarthritis: A systematic review directed from the preclinical data to randomized clinical studies. Nutrients. (2021) 13. doi: 10.3390/nu13051420
12. Chen JX, Liu AF, Zhou QX, Yu WJ, Guo TC, Jia YZ, et al. Acupuncture for the treatment of knee osteoarthritis: an overview of systematic reviews. Int J Gen Med. (2021) 19:8481–94. doi: 10.2147/IJGM.S342435
13. Wu BY, Yang L, Chen LY, Ma L, and Guo YT. Traditional Chinese medicine therapies for patients with knee osteoarthritis: A protocol for systematic review and network meta-analysis. Med (Baltimore). (2022) 101:e29404. doi: 10.1097/MD.0000000000029404
14. Zhou X, Xiang KM, Yuan XY, Wang ZP, and Li KL. Chinese herbal medicine Wutou decoction for knee osteoarthritis: A protocol for systematic review and meta-analysis. Med (Baltimore). (2020) 99:e22767. doi: 10.1097/MD.0000000000022767
15. Wang PE, Xu JB, Sun Q, Ge QW, Qiu M, Zou KA, et al. Chondroprotective Mechanism of Eucommia ulmoides Oliv.- Glycyrrhiza uralensis Fisch. Couplet Medicines in Knee Osteoarthritis via Experimental Study and Network Pharmacology Analysis. Drug Des Devel Ther. (2023) 17:633–46. doi: 10.2147/DDDT.S397185
16. Tan N, Jian G, Peng J, Tian X, and Chen B. Chishao - Fuzi herbal pair restore the macrophage M1/M2 balance in acute-on-chronic liver failure. J Ethnopharmacol. 328:118010. doi: 10.1016/j.jep.2024.118010
17. Tan X, Zheng DH, Lin Q, Wang LL, Zhu ZS, Huang YF, et al. Confirmation of pain-related neuromodulation mechanism of Bushen Zhuangjin Decoction on knee osteoarthritis. J Ethnopharmacol. (2022) 324:117772. doi: 10.1016/j.jep.2024.117772
18. Peng MZ, Yao ZL, Zhang JL, Lin Y, Xu L, Zhang Q, et al. Discovery and validation of anti-arthritic ingredients and mechanisms of Qingfu Juanbi Tang, a Chinese herbal formulation, on rheumatoid arthritis. J Ethnopharmacol. (2024) 329:118140. doi: 10.1016/j.jep.2024.118140
19. Li S, Yan Z, Zhi X, Wei HZ, Zi QZ, Zhe ND, et al. Extracellular vesicles from IPFP-MSCs trigger osteoarthritis by transferring mtDNA. Bioact Mater. (2025) 58:252–73. doi: 10.1016/j.bioactmat.2025.11.046
20. AbuBakr N, Fares AE, Mostafa A, and Farag DBE. Mesenchymal stem cells-derived microvesicles versus platelet-rich plasma in the treatment of monoiodoacetate-induced temporomandibular joint osteoarthritis in Albino rats. Heliyon. (2022) 8:e10857. doi: 10.1016/j.heliyon.2022.e10857
21. Xu SY, Bian RL, and Chen X. Pharmacological Experimental Methodology. Beijing: People’s Medical Publishing House (2002).
22. National Pharmacopoeia Commission. Pharmacopoeia of the People’s Republic of China. 2020 Edition. Beijing: China Medical Science Press (2020).
23. Sebastiano M, Chastel O, Eens M, and Costantini D. Gene expression provides mechanistic insights into a viral disease in seabirds. Sci Total Environ. (2024) 957:177478. doi: 10.1016/j.scitotenv.2024.177478
24. Zhang B, Fu D, Wang X, and Hu X. Network-based analysis and experimental validation of identified natural compounds from Yinchen Wuling San for acute myeloid leukemia. Front Pharmacol. (2025) 16:1591164. doi: 10.3389/fphar.2025.1591164
25. Huang R, Lu J, Yang X, Sheng G, Qin F, and Yang X. Molecular mechanism of the effect of BixiezelanYin on knee osteoarthritis based on network pharmacology and molecular docking. Med (Baltimore). (2025) 104:e41459. doi: 10.1097/MD.0000000000041459
26. Yu T, Fang Z, Cheng Y, Zhou Y, Yan CJ, Chang L, et al. A novel palmitoylation-based molecular signature reveals COX6A1 as a key regulator in metabolic dysfunction-associated steatotic liver disease. J Transl Med. (2025) 23:1212. doi: 10.1186/s12967-025-07253-0
27. Zhao H, Wang D, Fu D, and Xue L. Predicting the potential ankylosing spondylitis-related genes utilizing bioinformatics approaches. Rheumatol Int. (2015) 35:973–9. doi: 10.1007/s00296-014-3178-9
28. Deng Y, Wang D, Wang C, Kai YG, Ming L, Yan ZH, et al. Serum CHI3L1 levels correlate with disease activity in rheumatoid arthritis and reveal potential molecular mechanisms. Front Immunol. (2025) 16:1729989. doi: 10.3389/fimmu.2025.1729989
29. Mao S, Wang Y, Gu M, Kai L, Ji BM, Jun M, et al. Study on differentially expressed genes and pattern recognition receptors in osteoporosis based on bioinformatics analysis. Sci Rep. (2025) 15:31287. doi: 10.1038/s41598-025-16891-9
30. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, and Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell Syst. (2015) 1:417–25. doi: 10.1016/j.cels.2015.12.004
31. Aran D, Hu ZC, and Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. (2017) 18:220. doi: 10.1186/s13059-017-1349-1
32. Shannon P, Markiel A, Ozier Q, Baliga N, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. (2003) 13:2498–504. doi: 10.1101/gr.1239303
33. Liu Y, Grimm M, Dai WT, Hou MC, Xiao ZX, and Cao Y. CB-Dock: a web server for cavity detection-guided protein-ligand blind docking. Acta Pharmacol Sin. (2020) 41:138–44. doi: 10.1038/s41401-019-0228-6
34. Trott O and Olson AJ. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem. (2010) 31:455–61. doi: 10.1002/jcc.21334
35. Seeliger D and Groot BLD. Ligand docking and binding site analysis with PyMOL and Autodock/Vina. J Comput Aided Mol Des. (2010) 24:417–22. doi: 10.1007/s10822-010-9352-6
36. Satija R, Farrell JA, Gennert D, Schier AF, and Regev A. Spatial reconstruction of single-cell gene expression data. Nat Biotechnol. (2015) 33:495–502. doi: 10.1038/nbt.3192
37. Stuart T, Srivastava A, Madad S, Lareau CA, and Satija R. Single-cell chromatin state analysis with Signac. Nat Methods. (2021) 18:1333–41. doi: 10.1038/s41592-021-01282-5
38. Fan Y, Bian XH, Meng XG, Li L, Fu LY, Zhang YN, et al. Unveiling inflammatory and prehypertrophic cell populations as key contributors to knee cartilage degeneration in osteoarthritis using multi-omics data integration. Ann Rheum Dis. (2024) 83:926–44. doi: 10.1136/ard-2023-224420
39. Jin SQ, Guerrero-Juarez CF, Zhang LH, Chang IV, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. (2021) 12:1088. doi: 10.1038/s41467-021-21246-9
40. Zhang SL, Zhang KS, Wang JF, Yi HW, Hui Z, Chong T, et al. Corresponding changes of autophagy-related genes and proteins in different stages of knee osteoarthritis: an animal model study. Orthop Surg. (2022) 14:595–604. doi: 10.1111/os.13057
41. Chang JJ, Wu H, Wu J, Liu M, Zhang WT, Hu YF, et al. Constructing a novel mitochondrial-related gene signature for evaluating the tumor immune microenvironment and predicting survival in stomach adenocarcinoma. J Transl Med. (2023) 21:191. doi: 10.1186/s12967-023-04033-6
42. Zhang Y and Ji Q. Macrophage polarization in osteoarthritis progression: a promising therapeutic target. Front Cell Dev Biol Front Cell Dev Biol. (2023) 11:1269724. doi: 10.3389/fcell.2023.1269724
43. Zhang J, Xu L, Zhang JJ, Liu Y, Li X, Ren T, et al. Pan-cancer analysis of the prognostic and immunological role of matrix metalloproteinase 9. Med (Baltimore). (2023) 102:e34499. doi: 10.1097/MD.0000000000034499
44. Arab HH, Abd El-Aal SA, Ashour AM, EI-Sheikh AK, and AI Khabbaz HU. Targeting inflammation and redox perturbations by lisinopril mitigates Freund’s adjuvant-induced arthritis in rats: role of JAK-2/STAT-3/RANKL axis, MMPs, and VEGF. Inflammopharmacology. (2022) 30:1909–26. doi: 10.1007/s10787-022-00998-w
45. Qi WZ, Jin L, Huang SQ, Aikebaier A, Xue S, and Wang QY. Modulating synovial macrophage pyroptosis and mitophagy interactions to mitigate osteoarthritis progression using functionalized nanoparticles. Acta Biomater. (2024) 181:425–39. doi: 10.1016/j.actbio.2024.05.014
46. Mao W, Zheng Y, Zhang W, Yin JR, Liu ZY, He PL, et al. Enocyanin promotes osteogenesis and bone regeneration by inhibiting MMP9. Int J Mol Med. (2025) 55:9. doi: 10.3892/ijmm.2024.5450
47. Jarecki J, Małecka-Masalska T, Kosior-Jarecka E, Widuchowski W, Krasowski P, and Gutbier M. Concentration of selected metalloproteinases and osteocalcin in the serum and synovial fluid of obese women with advanced knee osteoarthritis. Int J Environ Res Public Health. (2022) 19:3530. doi: 10.3390/ijerph19063530
48. Luo S, Li W, Wu W, and Shi Q. Elevated expression of MMP8 and MMP9 contributes to diabetic osteoarthritis progression in a rat model. J Orthop Surg Res. (2021) 16:64. doi: 10.1186/s13018-021-02208-9
49. Wen D, Chen Z, Zhang Z, and Jia Q. The expression, purification, and substrate analysis of matrix metalloproteinases in Drosophila melanogaster. Protein Expr Purif. (2020) 171:105629. doi: 10.1016/j.pep.2020.105629
50. Chen YJ, Jeng JH, Chang HH, Huang MY, Tsai FF, and Yao CC. Differential regulation of collagen, lysyl oxidase and MMP-2 in human periodontal ligament cells by low- and high-level mechanical stretching. J Periodontal Res. (2013) 48:466–74. doi: 10.1111/jre.12028
51. Zheng LL, Chen WS, Xian GY, Pan BQ, Ye YY, and Gu MH. Identification of abnormally methylated-differentially expressed genes and pathways in osteoarthritis: a comprehensive bioinformatic study. Clin Rheumatol. (2021) 40:3247–56. doi: 10.1007/s10067-020-05539-w
52. Guo SL, Han CT, Jung JL, Chen WJ, Mei JJF, Lee HC, et al. Cystatin C in cerebrospinal fluid is upregulated in elderly patients with chronic osteoarthritis pain and modulated through matrix metalloproteinase 9-specific pathway. Clin J Pain. (2014) 30:331–9. doi: 10.1097/AJP.0b013e31829ca60b
53. Conte R, Finicelli M, Borrone A, Margarucci S, Peluso G, Calarco A, et al. MMP-2 Silencing through siRNA Loaded Positively-Charged Nanoparticles (AcPEI-NPs) Counteracts Chondrocyte De-Differentiation. Polymers (Basel). (2023) 15:1172. doi: 10.3390/polym15051172
54. Ermakov S, Leonov A, Trofimov S, Malkin I, and Livshits G. Quantitative genetic study of the circulating osteopontin in community-selected families. Osteoporos Int. (2011) 22:2261–71. doi: 10.1007/s00198-010-1451-7
55. Kang X, Zhang K, Wang Y, Zhao Y, and Lu Y. Single-cell RNA sequencing analysis of human chondrocytes reveals cell-cell communication alterations mediated by interactive signaling pathways in osteoarthritis. Front Cell Dev Biol. (2023) 11:1099287. doi: 10.3389/fcell.2023.1099287
56. Li M, Peng Y, Xiao Q, Li MC, Peng Y, and Xiao QQ. Pharmacological inhibition of PGK1 by genipin reduces inflammation and oxidative damage in osteoarthritis. Phytomedicine. (2026) 150:157729. doi: 10.1016/j.phymed.2025.157729
57. Brune Z, Rice MR, and Barnes BJ. Potential T cell-intrinsic regulatory roles for IRF5 via cytokine modulation in T helper subset differentiation and function. Front Immunol. (2020) 11:1143. doi: 10.3389/fimmu.2020.01143
58. Milenkovic I and Novoa EM. Ribosomal protein paralogues in ribosome specialization. Philos Trans R Soc Lond B Biol Sci. (2025) 380:20230387. doi: 10.1098/rstb.2023.0387
59. Zhou X, Liao WJ, Liao JM, Liao P, and Lu H. Ribosomal proteins: functions beyond the ribosome. J Mol Cell Biol. (2015) 7:92–104. doi: 10.1093/jmcb/mjv014
60. Chabronova A, van den Akker GGH, Housmans BAC, Caron MMJ, Cremers A, Surtel DAM, et al. Ribosomal RNA-based epitranscriptomic regulation of chondrocyte translation and proteome in osteoarthritis. Osteoarthritis Cartilage. (2023) 31:374–85. doi: 10.1016/j.joca.2022.12.010
61. He M, Tay WT, Song N, Liu T, Vincent Ren SX, Zhu CS, et al. Cardioprotective role of RBFox1 in myocardial infarction-induced heart failure. Cardiovasc Res. (2025) 121:2534–48. doi: 10.1093/cvr/cvaf206
62. Ismaeel A, Peck BD, Montgomery MM, Burke BI, Goh J, Franco AB, et al. microRNA-1 regulates metabolic flexibility by programming adult skeletal muscle pyruvate metabolism. Mol Metab. (2026) 98:102182. doi: 10.1016/j.molmet.2025.102182
63. Fu HY, Jie LS, Gong ZJ, Huang ZL, Zhu ZS, Liu JY, et al. Mechanism study on liquiritin alleviating nociceptive sensitization in knee osteoarthritis via promoting M2 macrophage polarization through regulation of the Rap1/PI3K/Akt signaling pathway. Phytother Res. (2026) 40:1–12. doi: 10.1002/ptr.70174
64. Liebmann K, Castillo M, Jergova S, Rahimi B, Kaplan LD, Best TM, et al. Functional assessment of genetically modified infrapatellar fat pad mesenchymal stem/stromal cell-derived extracellular vesicles (EVs): potential implications for inflammation/pain reversal in osteoarthritis. Cells. (2025) 14:1952. doi: 10.3390/cells14241952
65. Komori T. Runx2, an inducer of osteoblast and chondrocyte differentiation. Histochem Cell Biol. (2018) 149:313–23. doi: 10.1007/s00418-018-1640-6
66. Kredel L, Batra A, and Siegmund B. Role of fat and adipokines in intestinal inflammation. Curr Opin Gastroenterol. (2014) 30:559–65. doi: 10.1097/MOG.0000000000000116
67. MacLaren R, Cui W, and Cianflone K. Adipokines and the immune system: an adipocentric view. Adv Exp Med Biol. (2008) 632:1–21. doi: 10.1007/978-0-387-78952-1_1
68. Greene D, Moore Fried J, and Wang J. IgE in allergic diseases. Immunol Rev. (2025) 334:e70057. doi: 10.1111/imr.70057
69. Huang R, Wang WJ, Chen ZY, Chai J, Qi Q, Zheng HD, et al. Identifying immune cell infiltration and effective diagnostic biomarkers in Crohn’s disease by bioinformatics analysis. Front Immunol. (2023) 14:1162473. doi: 10.3389/fimmu.2023.1162473
70. Xia H, Zhang LL, Dai J, Liu XM, Zhang X, Zeng Z, et al. Effect of selenium and peroxynitrite on immune function of immature dendritic cells in humans. Med Sci Monit. (2021) 27:e929004. doi: 10.12659/MSM.929004
71. Baek JH, Birchmeier C, Zenke M, and Hieronymus T. The HGF receptor/Met tyrosine kinase is a key regulator of dendritic cell migration in skin immunity. J Immunol. (2012) 189:1699–707. doi: 10.4049/jimmunol.1200729
72. Ge T, He G, Cui Q, Wang SC, Wang ZK, Xie YY, et al. Decoupled dynamics of absolute and relative lymphocyte counts and age-polarized CD4+/CD8+ ratio in infants versus older adults. Front Immunol. (2025) 16:1599515. doi: 10.3389/fimmu.2025.1599515
73. Du J, Wang G, Zhao Y, Li M, Yang S, and Liu L. Osteopontin induces airway smooth muscle cells proliferation and migration by modulating FAK/Src/YAP axis. Int Immunopharmacol. (2025) 169:115886. doi: 10.1016/j.intimp.2025.115886
74. Johnson GA, Burghardt RC, and Bazer FW. Osteopontin: a leading candidate adhesion molecule for implantation in pigs and sheep. J Anim Sci Biotechnol. (2014) 5:56. doi: 10.1186/2049-1891-5-56
75. Liu Y, Jiang W, Huang J, and Zhong L. Bioinformatic analysis combined with immune infiltration to explore osteoarthritis biomarkers and drug prediction. Med (Baltimore). (2024) 103:e38430. doi: 10.1097/MD.0000000000038430
76. Han X, Tan S, Du B, Liu J, and Qu L. Discovery and experimentally mice model validation of CFI, a natural killer T cell-related gene, as a key biomarker in osteoarthritis. Clin Exp Immunol. (2025) 219:uxaf054. doi: 10.1093/cei/uxaf054
77. Ren FJ, Li Y, Luo HY, Gao S, Jiang SS, Yang J, et al. Extraction, detection, bioactivity, and product development of luteolin: A review. Heliyon. (2024) 10:e41068. doi: 10.1016/j.heliyon.2024.e41068
78. Gao XY, Liu ZX, Chen JX, Zhu DS, Liu H, Li JL, et al. Encapsulation of luteolin by self-assembled Rha/SSPS/SPI nano complexes: Characterization, stability, and gastrointestinal digestion in vitro. Food Res Int. (2024) 188:114532. doi: 10.1016/j.foodres.2024.114532
79. Wang R, Li X, Xu YH, Li YY, Zhang WS, Guo RQ, et al. Progress, pharmacokinetics and future perspectives of luteolin modulating signaling pathways to exert anticancer effects: A review. Med (Baltimore). (2024) 103:e39398. doi: 10.1097/MD.0000000000039398
80. Wang X, Zhang YH, Chang X, Wen XD, Tian F, Yu HJ, et al. The inhibitory effect of the active ingredients in the Bushen Huoxue formula on the IL-17A signaling pathway and its alleviating effect on osteoarthritis. J Inflammation Res. (2025) 18:6505–27. doi: 10.2147/JIR.S506716
81. Guo YF, Xu NN, Sun W, Zhao Y, Li CY, and Guo MY. Luteolin reduces inflammation in Staphylococcus aureus-induced mastitis by inhibiting NF-kB activation and MMPs expression. Oncotarget. (2017) 8:28481–93. doi: 10.18632/oncotarget.16092
82. Hsuan CF, Hsu HF, Tseng WK, Lee TL, Wei YF, Hsu KL, et al. Glossogyne tenuifolia Extract Inhibits TNF-α-Induced Expression of Adhesion Molecules in Human Umbilical Vein Endothelial Cells via Blocking the NF-kB Signaling Pathway. Molecules. (2015) 20:16908–23. doi: 10.3390/molecules200916908
83. Shehnaz SI, Roy A, Vijayaraghavan R, Sivanesan S, and Pazhanivel N. Modulation of PPAR-γ, SREBP-1c and inflammatory mediators by luteolin ameliorates β-cell dysfunction and renal damage in a rat model of type-2 diabetes mellitus. Mol Biol Rep. (2023) 50:9129–42. doi: 10.1007/s11033-023-08804-8
84. Li Y, Fu T, Yu W, Wen HY, Wang ZQ, Lyu ZX, et al. Mesenchymal stem cells and extracellular vesicles for knee osteoarthritis: clinical application, mechanism exploration and prospect. Stem Cell Res Ther. (2025) 16:688. doi: 10.1186/s13287-025-04783-8
85. Shivnath N, Siddiqui S, Rawat V, Khan MS, and Arshad M. Solanum xanthocarpum fruit extract promotes chondrocyte proliferation in vitro and protects cartilage damage in collagenase induced osteoarthritic rats (article reference number: JEP 114028). J Ethnopharmacol. (2021) 274:114028. doi: 10.1016/j.jep.2021.114028
Keywords: osteoarthritis, Bitong mixture, biomarkers, single-cell analysis, bioinformatics, network pharmacology, molecular docking, RT-qPCR
Citation: Shen Z, Li S, Li T and Wang L (2026) Identification of potential biomarkers related to the Bitong mixture in osteoarthritis based on bioinformatics and network pharmacology, and exploration of the mechanism involved. Front. Immunol. 17:1739355. doi: 10.3389/fimmu.2026.1739355
Received: 04 November 2025; Accepted: 27 January 2026; Revised: 14 January 2026;
Published: 13 February 2026.
Edited by:
Balachandran Ravindran, Institute of Life Sciences (ILS), IndiaReviewed by:
Hadida Yasmin, Cooch Behar Panchanan Barma University, IndiaMadhu Baghel, Henry Ford Health System, United States
Copyright © 2026 Shen, Li, Li and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Suming Li, b3J0aG9kb2N0b3JAMTYzLmNvbQ==