Mechanistic Insights into Active Components of Rosa Roxburghii Juice Against Fluoride-Induced Osteoarthritis

TL;DR

This study explores how Rosa roxburghii juice protects against fluoride-induced osteoarthritis by identifying key compounds and molecular mechanisms.

Contribution

The study systematically identifies active components and molecular mechanisms of Rosa roxburghii juice in preventing fluoride-induced osteoarthritis.

Findings

Quercetin, Epicatechin, Emodin, and Ellagic acid are key bioactive components in Rosa roxburghii juice that protect against F-OA.

These compounds suppress p53 mRNA expression and restore proteasome-mediated p53 degradation to mitigate fluoride-induced damage.

LiP–MS analysis revealed significant protein conformational changes, particularly in TP53, following treatment with the active compounds.

Abstract

Fluoride-induced osteoarthritis (F-OA) is a debilitating manifestation of endemic fluorosis, with limited preventive or therapeutic strategies. Rosa roxburghii juice (RRJ), a traditional medicinal/edible product, has shown protective effects against skeletal fluorosis, yet its active constituents and molecular mechanisms are not fully understood. In this study, an integrated strategy combining bioinformatics analysis, network pharmacology, molecular docking and dynamics simulations, limited proteolysis–mass spectrometry (LiP–MS), and in vitro experiments was employed to systematically elucidate the protective mechanisms of RRJ against F-OA. Forty-four core F-OA-associated genes were identified, with TP53 and the p53 signaling pathway emerging as central regulatory hubs. Quercetin, Epicatechin, Emodin, and Ellagic acid were screened as key bioactive components of RRJ and demonstrated strong binding affinity toward core targets, including TP53. Cellular experiments showed that these compounds significantly attenuated sodium fluoride-induced cellular injury. LiP–MS analysis further revealed widespread protein conformational remodeling following treatment, with TP53 exhibiting pronounced structural sensitivity. Mechanistically, these active compounds mitigated fluoride-induced pathological changes by suppressing p53 mRNA expression and restoring proteasome-mediated p53 degradation. This study provides systematic pharmacological evidence supporting Rosa roxburghii fruit as a promising functional food for the prevention and management of skeletal fluorosis and F-OA.

Figures

Full text

1. Introduction

Chronic systemic fluorosis is induced by prolonged and excessive fluoride intake, with skeletal fluorosis representing the most severe clinical subtype. Excessive fluoride accumulates in mineralized tissues and progressively disrupts bone homeostasis, resulting in osteosclerosis, joint pain, restricted mobility, skeletal deformity, and, in advanced stages, neurological compression-related paralysis [1]. Dental fluorosis is widely recognized as the earliest biomarker of fluoride overexposure; however, skeletal fluorosis reflects much deeper metabolic damage characterized by dysregulated osteoblast and osteoclast activity, oxidative stress, inflammatory injury, and impaired bone remodeling [2,3]. When fluoride toxicity involves articular tissues, secondary osteoarthritis frequently develops as a major clinical manifestation. Unlike primary osteoarthritis, which is usually driven by age-related degeneration, mechanical overload, or obesity, fluoride-induced osteoarthritis (F-OA) originates from chronic toxic insult. Recent mechanistic studies have demonstrated that fluoride disrupts cartilage extracellular matrix homeostasis, induces chondrocyte apoptosis, alters subchondral bone turnover, and upregulates catabolic mediators such as matrix metalloproteinase-13 (MMP-13) and inflammatory cytokine [4,5]. Due to its toxic etiology, conventional symptomatic treatments such as NSAIDs or intra-articular corticosteroids provide limited therapeutic benefit, highlighting the need for multi-target intervention strategies capable of attenuating fluoride-related oxidative, inflammatory, and metabolic disturbances.

Bioinformatics and network pharmacology demonstrate powerful complementary strengths for uncovering the deep toxicological mechanisms by which fluoride induces osteoarthritis. They can systematically integrate massive multi-omics data—from the genome and transcriptome to the metabolome—and, through efficient data mining and network analysis, reconstruct the full chain of “fluoride exposure–molecular targets–signaling pathways–pathological phenotypes.” This system-level mapping guides multi-target intervention strategies and thus enables a holistic dissection of the complex network mechanisms underlying fluorosis-associated osteoarthritic lesions [6]. Rosa roxburghii Tratt (RRT), a fruit belonging to the Rosaceae family, is indigenous to southwestern China, serves as a significant local agricultural product, and is traditionally used to aid digestion. Its fruit is rich in flavonoids, polyphenols, triterpenes, and vitamin C and exhibits marked antioxidant, anti-inflammatory, and cytoprotective activities. With both notable medicinal properties and edibility, it serves as a representative example of a medicine–food homology (MFH) substance [7]. Moreover, recent studies report that Rosa roxburghii juice (RRJ) has strong antioxidant effects and can improve indicators of liver and kidney function, protect collagenous tissues, and modulate the gut-microbiota–gut–bone axis; these findings together point to the therapeutic and development potential of RRJ against endemic fluorosis-related skeletal disease, including fluoride-induced osteoarthritis [8,9]. Accordingly, an integrated pharmacological approach was adopted to elucidate RRJ’s protective mechanisms against F-OA, alongside a preliminary in-depth exploration of its bioactive components and underlying actions. Through bioinformatic prediction, network toxicology/pharmacology analyses, molecular docking, liquid chromatography–mass spectrometry (Lip-MS), ligand fishing, and molecular dynamics simulations, we constructed an interaction network between the active constituents of RRJ and OA-related targets, followed by validation using Rat synovial fibroblast models to confirm the computational findings. This study is the first to systematically demonstrate the potential of RRJ as a multi-component, multi-target MFH intervention for F-OA, thereby laying the groundwork for subsequent translational studies and fostering the advancement of therapeutic or functional food products for skeletal fluorosis.

2. Materials and Methods

2.1. Materials and Reagents

Sodium fluoride standard solution (S885178) and ELISA kits for TNF-α, IL-1β, IL-6, and IL-10 were sourced from Macklin Biochemical Technology Co., Ltd. (Shanghai, China) and Shanghai Fanyin Biotechnology Co., Ltd. (Shanghai, China), respectively. The Reactive Oxygen Species (ROS) Assay Kit (S0033S) was obtained from Beyotime Biotechnology Co., Ltd. (Shanghai, China). The Cell Cycle Detection Kit (KGA9101) was purchased from KeyGEN BioTECH Co., Ltd. (Nanjing, China). The antibodies used in the experiment included the following: anti-P53 (10442-1-AP), anti-P21 (28248-1-AP), and anti-P-P53 (67826-1-lg) were all procured from Wuhan Sanying Biotechnology Co., Ltd. (Wuhan, China); anti-α-tubulin (ER130905) was acquired from Hangzhou Hua’an Biotechnology Co., Ltd. (Hangzhou, China); and Rat synovial fibroblast cells (iCell-0135a) and the cell culture medium (11095080) were obtained from iCell Bioscience Inc. (Shanghai, China) and Thermo Fisher Scientific Inc. (Waltham, MA, USA), respectively. The RNA-easy Isolation Reagent kit was purchased from Vazyme Biotech Co., Ltd. (Nanjing, China). The cDNA Synthesis SuperMix kit and PerfectStart Green qPCR SuperMix were procured from TransGen Biotech Co., Ltd. (Beijing, China). Emodin (C_15_H_10_O_5_, CAS 518-82-1, HPLC ≥ 98%), Quercetin (C_15_H_10_O_7_, CAS 117-39-5, HPLC ≥ 98%), Epicatechin (C_15_H_14_O_6_, CAS 490-46-0, HPLC ≥ 98%), and Ellagic acid (C_14_H_6_O_8_, CAS 476-66-4, HPLC ≥ 98%) were purchased from Shanghai Yuanye Bio-Technology Co., Ltd. (Shanghai, China).

2.2. Bioinformatic Data Download and Preprocessing

In this study, microarray datasets were analyzed to obtain gene expression profiles and clinical phenotypic information associated with the keywords “osteoarthritis,” “Homo sapiens,” and “gene expression.” All relevant gene expression datasets and platform annotation files were accessed from the publicly available Gene Expression Omnibus (GEO) repository (https://www.ncbi.nlm.nih.gov/geo/). As such, this study did not require separate ethical approval or informed consent. To ensure analytical consistency, two datasets (GSE82107 and GSE206848) sharing the GPL570-55999 platform were selected. GSE82107 comprised samples from 7 healthy controls and 10 OA patients, and GSE206848 included 7 controls alongside 9 OA patients.

To obtain accurate mRNA expression profiles, Perl scripts were used to match and align the transcriptomic data from GSE82107 and GSE206848. The sva package (version 3.50.0) in R (version 4.3.2) [10] was employed to identify and correct batch effects using its default parameter settings. Background correction was applied using the robust multi-array average (RMA) method [11] to ensure data consistency across platforms. The ComBat algorithm was subsequently employed to further eliminate batch effects between samples. The processed data were transformed using a natural logarithm to optimize data distribution and enhance analytical performance. Post-correction expression data were visualized using box plots and principal component analysis (PCA), yielding a final dataset of 14 control samples and 17 OA samples for downstream analysis.

Furthermore, to explore pathogenic genes associated with OA, we retrieved information on genes known to be involved in OA progression (as of October 2024) from the Comparative Toxicogenomics Database (CTD, http://ctdbase.org). Candidate genes related to sodium fluoride were predicted by integrating data from ChEMBL (https://www.ebi.ac.uk/chembl/), STITCH (http://stitch.embl.de/), PubChem (https://pubchem.ncbi.nlm.nih.gov/), and the CTD database.

2.3. Identification and Processing of Co-Expressed Genes Related to OA and Fluorosis

To identify differentially expressed genes (DEGs) between osteoarthritis (OA) samples and healthy controls, differential expression analysis was conducted on the integrated dataset using the R package limma (version 3.58.1) [12]. We defined statistical significance using a p-value cutoff of 0.05 and required an absolute log_2_ fold change (|log_2_FC|) of more than 0.585 for gene selection. Subsequently, high-confidence OA-associated genes with an Inference Score > 10 were obtained from the CTD. These genes were merged with sodium fluoride-associated genes predicted from STITCH, ChEMBL, PubChem, and CTD databases. After removing duplicates, a candidate gene set related to fluoride toxicity was compiled. The intersection among OA-related DEGs, key OA genes, and fluoride-related genes was then defined, yielding a set of candidate genes potentially linking OA to fluorosis.

The expression patterns of DEGs were visualized through heatmaps using the ggplot2 (version 4.0.1) and pheatmap (version 1.0.13) packages [13]. To further explore gene-gene correlations, the corrplot package was used to calculate Pearson correlation coefficients, helping to uncover potential regulatory modules and network structures. These overlapping target genes were then input into the STRING database [14] (https://cn.string-db.org/) for protein-protein interaction (PPI) analysis. The resulting network was imported into Cytoscape software (version 3.8.2) [15] for further topological analysis and hub gene identification.

2.4. Pathway and Functional Profiling of Fluorosis-Related OA Genes

To elucidate the biological roles and pathways of the candidate genes, we conducted Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses with the R package clusterProfiler (version 4.10.1) [16]. A p value of < 0.05 was considered statistically significant. This analysis sought to clarify the biological processes and signaling pathways by which sodium fluoride may contribute to the initiation and progression of osteoarthritis (OA), offering preliminary mechanistic insights.

2.5. Preparation of Rosa Roxburghii Juice (RRJ)

Fresh fruits of Rosa roxburghii Tratt. cultivar ‘Guinong 5’ were purchased from the Jinsha area of Bijie City, Guizhou Province, China. The plant material was authenticated as fresh fruits of the current season by the School of Pharmacy, Guizhou Medical University. A photograph of the representative fruit samples used in this study is provided in Supplementary Figure S1. After collection, the fruits were placed in sealed plastic ziplock bags, transported to the laboratory in a portable refrigerated container maintained at 4 °C, thoroughly washed with distilled water, and air-dried. The juice was extracted using a juice extractor, and the freshly obtained juice was directly used for subsequent analyses and is hereafter referred to as Rosa roxburghii juice (RRJ).

Sample preparation for LC–MS analysis: A 600 μL aliquot of the RRJ suspension was transferred into a 1.5 mL EP tube, followed by the addition of 400 μL of methanol. The mixture was vortexed for 10 s, after which 200 μL of the resulting solution was mixed with 200 μL of 40% (v/v) aqueous methanol and vortexed again for 10 s. For LC-MS analysis, samples were centrifuged at 16,000× g for 15 min at 4 °C, and the supernatant was collected for analysis. For HPLC analysis, the sample was prepared by diluting 1 mL of the RRJ suspension to 10 mL with 70% (v/v) methanol in a volumetric flask. After standing, the supernatant was similarly collected for HPLC analysis.

2.6. Establishment of LC–MS and HPLC Analytical Methods for RRJ

LC–MS analysis: Chemical profiling of RRJ was performed using an ultra-high-performance liquid chromatography (UHPLC) system coupled to a Q-Exactive HFX mass spectrometer (Thermo Fisher Scientific, Waltham, MA, USA). Chromatographic separation was achieved at a column temperature of 35 °C with a flow rate of 0.3 mL/min. The mobile phase consisted of 0.1% formic acid in water (A) and 0.1% formic acid in acetonitrile (B). Gradient elution was conducted according to the program shown in Table 1. Mass spectrometric detection was carried out using an electrospray ionization (ESI) source operated in both positive and negative ion modes. Data acquisition was performed in Full-MS/dd-MS^2^ mode with resolutions set to 60,000 (MS^1^) and 15,000 (MS^2^). The stepped normalized collision energies were 20, 40, and 60, and the MS^1^ scan range covered m/z 90 to 1300.

For quantitative analysis, an HPLC-DAD method was developed using a Dionex Ultimate 3000 system (Thermo Fisher Scientific, Waltham, MA, USA). Chromatographic separation was achieved on a Gemini NX-C18 reversed-phase column (Phenomenex) maintained at 30 °C, with an injection volume of 5 μL. The mobile phase consisted of solvent A (0.1% phosphoric acid in water) and solvent B (acetonitrile), applied under a linear gradient program over 30 min.

For ellagic acid analysis, the gradient was programmed from 80% A/20% B to 70% A/30% B, with detection at 254 nm. For epicatechin determination, the gradient was adjusted from 85% A/15% B to 80% A/20% B, and detection was performed at 280 nm. Reference standards of ellagic acid (23 mg/mL) and epicatechin (68.6 mg/mL) were dissolved in methanol to prepare stock solutions.

Method validation was conducted in accordance with the International Council for Harmonisation guideline ICH Q2(R1). Linearity was evaluated at six concentration levels for each analyte by constructing calibration curves of peak area versus concentration. The corresponding calibration curves are shown in Supplementary Figure S2A (epicatechin) and Supplementary Figure S2B (ellagic acid). The limits of detection (LODs) and quantification (LOQs) were determined based on signal-to-noise ratios of 3 and 10, respectively, and are illustrated in Supplementary Figure S2C–F.

Precision was assessed through intra-day and inter-day analyses at three concentration levels, while repeatability was evaluated using six independently prepared sample solutions (Supplementary Table S1). Accuracy was determined by recovery experiments using the standard addition method at spiking levels of 80%, 100%, and 120%, with results presented in Supplementary Table S2. Comprehensive validation parameters, including regression equations, correlation coefficients (R^2^), LODs, LOQs, precision (RSD%), and recovery values, are summarized in Supplementary Table S3.

2.7. Network Pharmacology Analysis

The active compounds present in RRJ were sourced based on a comprehensive review of previously published literature [17] and their molecular structure information was retrieved from the PubChem database. Potential targets of these compounds were predicted using the SwissTargetPrediction platform [18,19] (https://www.swisstargetprediction.ch/) and the Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform (TCMSP) (https://www.tcmsp-e.com/load_intro.php?id=43). These predicted targets were then intersected with the previously identified fluorosis-related osteoarthritis (OA) gene set to obtain shared potential regulatory targets.

The shared targets were submitted to the STRING database to construct a protein–protein interaction (PPI) network. A Sankey diagram of “drug-compound-target-disease” relationships was constructed using the Microbioinformatics platform. Additionally, GO functional annotation and KEGG pathway enrichment analyses of the common targets were conducted using the clusterProfiler package in R [16] to investigate the biological processes and signaling pathways potentially underlying the effects of RRJ on fluorosis-associated osteoarthritis (OA). Enrichment outcomes were visualized using bubble plots, with statistical significance defined as p < 0.05.

2.8. Molecular Docking

Compound structures were retrieved from PubChem (https://pubchem.ncbi.nlm.nih.gov/) and converted to PDBQT format using Open Babel 2.3.2. Crystal structures of the target proteins, including AKT1 (7MYX), BAX (8SRX), FGR (7JT9), FOS (1A02), HP (4 × 0 L), IL1A (5UC6), MMP12 (1RMZ), PTGS1 (6Y3C), and TP53 (3DCY), were obtained from the Protein Data Bank (PDB).

All receptor proteins underwent structural preparation in PyMOL 2.3.0, including removal of solvent molecules and existing ligands, after which hydrogen addition and charge assignment were performed with AutoDockTools. The finalized protein models were exported in PDBQT format for subsequent docking analysis. Receptor–ligand molecular docking was executed using AutoDock Vina 1.5.7 [20]. Binding interactions were characterized using the Protein–Ligand Interaction Profiler (PLIP) [21] (https://plip-tool.biotec.tu-dresden.de/plip-web/plip/index), while docking poses were rendered with PyMOL. Two-dimensional representations of protein–ligand interactions were generated using Discovery Studio 2019.

2.9. Molecular Dynamics (MD) Simulation and Binding Free Energy Calculation

Molecular dynamics (MD) simulations were conducted using GROMACS 2022.3 [22]. The CHARMM36 force field (charmm36-jul2022.ff) was applied to generate topology files for the TP53 protein. The simulation system was solvated with the TIP3P water model, and Na^+^ and Cl^−^ ions were introduced to ensure overall charge neutrality. Equilibration was performed under both NVT (constant volume and temperature) and NPT (constant pressure and temperature) ensembles for 100 ps at 300 K. This was followed by an 80 ns production MD simulation for each protein–ligand complex using a 2 fs integration time step.

To investigate the binding affinity between ligands (Quercetin, Ellagic Acid, Epicatechin, and Emodin) and TP53, the molecular mechanics/Poisson-Boltzmann surface area (MM/PBSA) method was applied. The binding free energy was calculated using the gmx_MMPBSA (version 1.9.1) [23] module (Python-based interface) from the GROMACS suite [24]. The total binding free energy was derived from the difference between the free energy of solvation and the molecular mechanical potential energy of the complex.

Furthermore, the binding energy contributions of individual amino acid residues at the protein-ligand interface were analyzed to identify key residues involved in stabilizing the complexes with TP53.

2.10. Sample Preparation and Extraction in Limited Proteolysis Mass Spectrometry (Lip-MS) Analysis

Rat knee synovial fibroblasts were homogenized and proteins were extracted using cold 1× PBS, with protein quantification performed by BCA assay. After activation, cells were randomly divided into untreated and drug-treated groups (treated with a mixture of 5 μM Quercetin, 20 μM Epicatechin, 5 μM Emodin, and 1 μM Ellagic acid). During drug treatment, the control group received the corresponding solvent.

The treated samples were digested with Proteinase K (PK), and 15 μg of protein was used for SDS-PAGE. Following PK treatment, denaturants (UA/DOC) were added, followed by DTT at a final concentration of 20 mM and incubation at 30 °C for 2 h. The samples were then cooled to room temperature, and IAA was added to a final concentration of 25 mM. The mixture was shaken at 600 rpm for 1 min and then incubated in the dark at room temperature for 30 min.

Ammonium bicarbonate buffer (50 mM) was added to reduce the UA/DOC concentration to below 1.5 M. Then, 2 μg of trypsin was added, and the samples were incubated at 37 °C for 16 h. After desalting and lyophilization, the samples were reconstituted in 0.1% formic acid. Peptide concentration was measured by OD280. Finally, 2 μg of peptides were mixed with an appropriate amount of iRT standard peptides for DIA mass spectrometry analysis.

2.11. Establishment of Lip-MS Analysis Method

Chromatographic separation was performed using a Vanquish Neo system (Thermo Fisher, Waltham, MA, USA) followed by Data-Independent Acquisition (DIA) on an Astral high-resolution mass spectrometer (Thermo Scientific, Waltham, MA, USA) operated in positive ion mode. The precursor ion range was set from 380 to 980 m/z. MS1 spectra were acquired at a 240,000 resolution (at 200 m/z) with a normalized AGC target of 500% and a maximum injection time of 5 ms. For DIA MS2 scans, a 2 m/z isolation window was applied, HCD collision energy was set to 25 eV, and a normalized AGC target of 500% with a 3 ms maximum injection time was used. DIA data were analyzed with Spectronaut software (https://biognosys.com/software/spectronaut/).

2.12. Cell Culture and Treatment

Rat synovial fibroblasts (RSFs) were cultured in complete DMEM medium (Gibco, China) supplemented with 5.5 mmol/L D-glucose, 10% fetal bovine serum (FBS, TransSerum, Australia), and 1% penicillin-streptomycin. Prior to experiments, cells were synchronized by serum starvation for 24 h under standard culture conditions (37 °C, 95% humidity, 5% CO_2_). The cells were then pretreated for 2 h with either sodium fluoride or the respective test compounds.

The experimental groups were designated as follows: Control, Sodium Fluoride (600 μM NaF), Quercetin (600 μM NaF + 5 μM Quercetin), Epicatechin (600 μM NaF + 20 μM Epicatechin), Emodin (600 μM NaF + 5 μM Emodin), and Ellagic Acid (600 μM NaF + 1 μM Ellagic acid).

Cell viability and proliferation were assessed using the MTT assay. Briefly, RSFs were seeded in 96-well plates at an initial density of 2.5–3 × 10^4^ cells per well, including blank control wells (without cells) and untreated control wells (cells with medium only). Each treatment group consisted of six replicate wells. Following 48 h of treatment with NaF (Guangrun Bio, Nanjing, China) at the specified concentrations, 10 μL of MTT reagent was added to each well and incubated for 4 h at 37 °C. Optical density at 490 nm was determined using a Varioskan™ LUX multimode microplate reader (Thermo Fisher Scientific, Waltham, MA, USA).

2.13. Real-Time Quantitative PCR (RT-qPCR) Analysis

Total RNA from RSFs was extracted using the RNeasy Total RNA Kit (Vazyme, Nanjing, China). Following genomic DNA removal, cDNA was generated using the EasyScript^®^ Genomic DNA Removal and cDNA Synthesis Kit (TransGen Biotech, Beijing, China). Quantitative real-time PCR was performed on the resulting cDNA with PerfectStart^®^ Green qPCR SuperMix (TransGen Biotech, Beijing, China) on a Bio-Rad CFX Manager 3.1 system (Bio-Rad Laboratories). Relative mRNA expression was calculated via the 2^(−ΔΔCt) method, using β-Actin as the reference gene. All primers were custom-synthesized by Sangon Biotech (Shanghai, China); the primer sequences are listed in Table 2.

2.14. Determination of Intracellular Reactive Oxygen Species (ROS) Levels

Intracellular reactive oxygen species (ROS) levels were measured using the fluorescent probe 2’,7’-dichlorofluorescin diacetate (DCFH-DA). After the indicated treatments, cells were incubated with 10 μM DCFH-DA in serum-free medium at 37 °C for 30 min in the dark. Following incubation, cells were washed three times with ice-cold phosphate-buffered saline (PBS) to remove excess probe. Fluorescence intensity was immediately analyzed using a NovoCyte flow cytometer (ACEA Biosciences, Hangzhou, China). Data acquisition and analysis were performed with NovoExpress software (https://www.agilent.com/product/novoexpress-software) (Agilent Technologies, Santa Clara, CA, USA). The mean fluorescence intensity (MFI) was calculated to quantify intracellular ROS levels, reflecting the degree of oxidative stress.

2.15. Detection of Inflammatory Cytokine Protein Expression (ELISA)

Cell culture supernatants were collected, and the secretion levels of inflammatory cytokines (TNF-α, IL-1β, IL-6) were quantitatively measured using commercial ELISA kits (R&D Systems, Fanyin, Shanghai, China) according to the manufacturer’s instructions. The changes in cytokine secretion across experimental groups were assessed to validate the corresponding mRNA expression results obtained from RT-qPCR.

2.16. Cell Cycle Analysis

Cell cycle distribution was analyzed by flow cytometry. After treatment, cells were harvested, washed with PBS, and fixed in 70% cold ethanol at 4 °C overnight. Fixed cells were washed twice with PBS and incubated with propidium iodide (PI) staining solution containing RNase A (Keygen Biotech, Jiangsu, China) at 37 °C for 30 min in the dark. Cell cycle profiles were determined using a NovoCyte flow cytometer (ACEA Biosciences, Hangzhou, China). Data were analyzed with NovoExpress software (Agilent Technologies, Santa Clara, CA, USA). The percentages of cells in the G0/G1, S, and G2/M phases were calculated to evaluate the effects of NaF and drug treatments on cell cycle progression and proliferation.

2.17. Immunofluorescence Staining and Analysis

RSFs were cultured on glass coverslips in 24-well plates at a density of 1 × 10^5^ cells per well for 24 h. After the designated treatments, cells were fixed with 4% paraformaldehyde (PFA) in PBS for 15 min at room temperature and washed three times with PBS. Cells were permeabilized using 0.2% Triton X-100 (Solarbio, Beijing, China) for 10 min and then blocked with 5% bovine serum albumin (BSA) in PBS for 1 h at room temperature to minimize non-specific binding.

Cells were treated with a 1:200 dilution of anti-p53 primary antibody in blocking solution and incubated overnight at 4 °C. Cells were washed three times with PBS and then incubated with a Cy3-conjugated secondary antibody for 1 h at room temperature in the dark. After washing, coverslips were mounted onto glass slides with an anti-fade medium containing 4′,6-diamidino-2-phenylindole (DAPI) to visualize nuclei.

Images were obtained using a Nikon Eclipse C1 fluorescence microscope (Nikon, Tokyo, Japan), and the fluorescence signal of p53 was quantified with ImageJ software (v1.52a, NIH, USA). For each coverslip, a minimum of five randomly selected fields were analyzed, with three independent experiments performed for statistical reliability.

2.18. Detection of Key Signaling Pathway Protein Expression by Western Blot

Total protein was isolated from cells using RIPA lysis buffer, and concentrations were assessed via BCA assay. Denaturation was performed by heating protein samples at 100 °C for 10 min in loading buffer. Equivalent amounts of protein were separated on 10% SDS-PAGE gels and subsequently transferred to 0.45 μm PVDF membranes for immunoblotting. Membranes were blocked with 5% skim milk at room temperature for 2 h and then incubated overnight at 4 °C with the following primary antibodies: anti-P53 (1:5000), anti-P-P53 (1:1000), anti-P21 (1:2000), and anti-α-tubulin (1:5000). The next day, secondary antibody incubation was carried out for 2 h at room temperature. Protein signals were detected using the NcmECL High Kit (NCM Biotech, Suzhou, China) and quantified via densitometric analysis with Image Lab software (Bio-Rad, Hercules, CA, USA, version 5.2).

2.19. Statistical Analysis

All algorithms and data analyses were performed using R software (version 4.3.2). Comparisons of continuous variables between two groups were conducted using Student’s t-test in GraphPad Prism (version 9.0), while comparisons among multiple groups were carried out using one-way ANOVA followed by Dunnett’s post hoc test. The significance levels were defined as follows: * p < 0.05, ** p < 0.01, *** p < 0.001 (indicating between-group differences); # p < 0.05, ## p < 0.01, ### p < 0.001 (denoting other specific comparisons).

3. Results

3.1. Integrated Bioinformatic Analysis Identifies Core Genes Associated with Fluoride-Related Osteoarthritis

We analyzed two microarray datasets, GSE82107 and GSE206848, to investigate the potential relationship between fluoride exposure and the development of OA. These datasets consisted of 14 healthy control samples and 19 OA samples. Initially, batch effects between the two datasets were identified and corrected using the “sva” (version 3.50.0) package in R (version 4.3.2). Principal Component Analysis (PCA) of the raw data revealed a clear separation between the two datasets before correction (Figure 1A). After the removal of batch effects, however, the samples from both groups were distributed homogeneously, indicating successful integration of the datasets for further analysis (Figure 1B). This confirmed that the combined dataset was reliable for subsequent analyses.

Comparative analysis between OA and control samples identified 1165 differentially expressed genes (DEGs), of which 675 were upregulated and 490 downregulated (Figure 1C). Further integration of these DEGs with gene sets from the GEO database, known OA-related genes, and NaF-associated genes allowed us to identify 44 overlapping candidate genes (Figure 1D). The heatmap demonstrated that OA and control groups exhibited markedly different expression patterns for these genes (Figure 1E), while correlation analysis revealed close functional interactions among these genes (Figure 1F). Consequently, a total of 44 core F-OA genes were ultimately identified, with their detailed information presented in Table 3.

3.2. Prediction of the Potential Mechanisms of RRJ Active Ingredients Against F-OA Based on Network Pharmacology Analysis

To systematically elucidate the molecular mechanism of RRJ active components in the treatment of F-OA, we first constructed a protein-protein interaction (PPI) network using the 44 previously identified core genes through the STRING database (confidence score > 0.4; Figure 2A). Subsequently, based on the previously established RRJ bioactive ingredient database [17], we collected 351 potential targets associated with the 44 ingredients from SwissTargetPrediction and TCMSP and performed an intersection analysis with the core genes. The results showed that these overlapping targets occupied central hub positions in the PPI network (Figure 2B), including nine core targets such as TP53, PTGS1, MMP12, and IL1A. Further analysis revealed that seven of the 44 RRJ ingredients, including Quercetin, Epicatechin, Emodin, and Ellagic acid, were interconnected with these core targets, as depicted in Figure 2C.

Pathway enrichment analysis of these target genes indicated their significant enrichment in oxidative stress, the p53 signaling pathway, and the TNF signaling pathway (Figure 2D). Gene Ontology (GO) analysis indicated that these genes were primarily involved in biological processes such as response to oxidative stress, regulation of cytochrome c release from mitochondria, and T cell apoptotic processes, with molecular functions including DNA binding, antioxidant activity, and protein folding chaperone binding (Figure 2E). Gene Set Enrichment Analysis (GSEA) further demonstrated that gene sets from the control group were predominantly enriched in pathways related to mitochondrial energy metabolism, while the OA group showed significant enrichment in inflammation and stress response pathways (Figure 2F,G). In summary, integrating evidence from core gene identification, multi-faceted enrichment analyses, and the existing literature, we found that TP53 and its associated p53 signaling pathway play pivotal roles in the aforementioned processes of oxidative stress, inflammation, immune response, and apoptosis. The bioactive components related to RRJ may exert a systemic protective effect against F-OA by modulating multiple pathways, including oxidative stress, apoptosis, and inflammation, through the TP53/p53 signaling axis.

3.3. Molecular Docking Analysis of Active Ingredients Related to RRJ

Molecular docking was used to assess the interaction mechanisms between the principal bioactive constituents of RRJ and core target proteins. The binding energy, a crucial thermodynamic parameter, indicates the stability of protein–ligand complexes, with lower values corresponding to more stable binding. Energies below 0 kcal/mol signify a spontaneous interaction, whereas values below −5.0 kcal/mol are indicative of strong binding affinity.

The molecular docking results revealed that the binding energies between the active components of RRJ and the core targets ranged from −5.0 to −9.0 kcal/mol (Figure 3A), reflecting potent and stable intermolecular interactions between the bioactive compounds and target proteins. Representative docking poses were also shown (Figure 3B–J), illustrating the binding modes and key interaction sites between the ligands and receptors. These results suggest that the active components of RRJ may exert their therapeutic effects by achieving high-affinity binding with disease-related target proteins, thereby modulating the pathological processes of F-OA. Finally, based on the average binding free energy from molecular docking, we further predicted four components—Quercetin, Epicatechin, Emodin, and Ellagic acid—to be the highly core target-associated bioactive ingredients in RRJ. These results suggest that these highly associated RRJ components may regulate the pathological process of fluoride-induced osteoarthritis through high-affinity binding to disease-relevant target proteins.

3.4. RRJ Highly Correlated Active Ingredients Can Restore the Viability and Function of NaF-Damaged Synovial Fibroblasts

To elucidate the protective effects of the four highly RRJ-associated active components—Quercetin, Epicatechin, Emodin, and Ellagic acid—in alleviating fluoride-induced osteoarthritic cell damage, we first evaluated the cytotoxicity of sodium fluoride (NaF) on rat synovial fibroblasts. RRJ NaF caused a clear dose-dependent decline in cell viability, with viability reduced to below 50% at 1000 and 1400 μM (Figure 4A). We next assessed the intrinsic cytotoxicity of Quercetin, Epicatechin, Emodin, and Ellagic acid (Figure 4B) and established a cellular injury model using 600 μM NaF. Under this condition, all four active components significantly rescued NaF-induced loss of cell viability (Figure 4C). Accordingly, the following treatment regimens were applied: Sodium Fluoride (600 μM NaF), Quercetin (600 μM NaF + 5 μM Quercetin), Epicatechin (600 μM NaF + 20 μM Epicatechin), Emodin (600 μM NaF + 5 μM Emodin), and Ellagic Acid (600 μM NaF + 1 μM Ellagic acid). Because NaF-mediated cytotoxicity was accompanied by alterations in cell-cycle progression, we further examined cell-cycle distribution. NaF exposure induced pronounced G1-phase arrest, reflected by an elevated G1 fraction and reduced S and G2/M fractions (Figure 4D), whereas each active component effectively reversed this arrest and restored normal cell-cycle dynamics (Figure 4E–G).

To further explore the anti-inflammatory and antioxidant mechanisms involved, we quantified the mRNA levels of inflammatory cytokines and ROS. qPCR analysis revealed that NaF markedly upregulated the mRNA expression levels of IL-1β, IL-6, and TNF-α, whereas all four active components of RRJ significantly counteracted this effect (Figure 5A–C). NaF also induced a substantial increase in ROS production, which was effectively scavenged by these active components (Figure 5D,E). At the protein level, NaF stimulated the secretion of IL-1β, IL-6, TNF-α and IL-10, and this induction was strongly inhibited by treatment with the active components (Figure 5F–I).

Collectively, these findings demonstrate that Quercetin, Epicatechin, Emodin, and Ellagic acid can counteract fluoride ion toxicity by enhancing cell viability, suppressing inflammatory responses, reducing oxidative stress, and restoring normal cell cycle regulation. These combined pharmacological effects suggest their potential in effectively alleviating fluoride poisoning-induced osteoarthritis.

3.5. Protein Conformational Remodeling and Functional Network Reprogramming Induced by Core Components of RRJ Identified Through LiP–MS

To determine whether Quercetin, Epicatechin, Emodin, and Ellagic acid can induce protein conformational remodeling in synovial fibroblasts, we employed Limited Proteolysis–Mass Spectrometry (LiP–MS) to compare conformation-specific peptides between control and treated samples under both native and denatured conditions (Figure 6A). Quantitative analysis revealed significant changes in peptide abundance between the groups treated with Quercetin, Epicatechin, Emodin, and Ellagic acid and the control group. A total of 8853 peptides were significantly upregulated and 3678 were downregulated, corresponding to 2396 upregulated and 1492 downregulated proteins, respectively (Figure 6B). These results indicate that the highly RRJ-associated active components induced extensive alterations in protein structural states. Unsupervised clustering demonstrated that the treatment and control groups were clearly segregated, with significantly increased expression of proteins in the treated samples (Figure 6C). This distinct clustering pattern highlights the substantial remodeling of both protein conformation and expression profiles induced by the highly RRJ-associated active components.

Functional enrichment analysis of the differentially expressed proteins showed significant enrichment in extracellular regions, vesicles, cytoplasmic components, and adhesion junctions (Figure 6D). KEGG pathway analysis demonstrated significant enrichment in pathways related to neurodegenerative diseases and oxidative stress, including oxidative phosphorylation, diabetic cardiomyopathy, and endoplasmic reticulum protein processing (Figure 6E). These results indicate that the highly RRJ-associated active components regulate processes related to cellular homeostasis, stress response, immunity, and inflammation through extensive remodeling of protein conformations and expression networks. Notably, targeted fishing for the TP53 protein revealed consistency between the LiP-MS results and the predicted core target TP53, with multiple peptides showing significant upregulation after treatment. This suggests that these active components may directly enhance the structural stability of TP53, thereby perturbing the P53 pathway and exerting key regulatory effects at the molecular level.

3.6. Comprehensive MD and MMPBSA Analyses Reveal Stable and Favorable Binding of Highly Relevant RRJ Active Compounds to TP53

To further elucidate the molecular mechanisms underlying the binding of Ellagic acid, Emodin, Epicatechin, and Quercetin to TP53, we performed 80 ns molecular dynamics (MD) simulations to systematically evaluate the conformational dynamics, structural stability, and interaction characteristics of the complexes. Global structural analysis showed that the backbone root mean square deviation (RMSD) of TP53 in all complexes stabilized between 0.15–0.30 nm, indicating overall structural stability, with the TP53–Quercetin complex exhibiting slightly higher conformational flexibility during the later stages of the simulation (Figure 7A). Ligand RMSD analysis revealed that Ellagic acid and epicatechin remained highly stable within the binding pocket (RMSD < 0.3 nm), whereas emodin and Quercetin displayed larger fluctuations (up to ~0.8 nm), suggesting potential conformational rearrangements or relatively weaker binding stability (Figure 7B).

Assessment of the radius of gyration (Rg) demonstrated that TP53 maintained a compact and stable structure upon ligand binding (Figure 7C). Solvent-accessible surface area (SASA) values exhibited a slight decreasing trend across all systems, while the TP53–Quercetin complex maintained relatively higher SASA, implying a more solvent-exposed binding conformation (Figure 7D). Per-residue root mean square fluctuation (RMSF) analysis indicated elevated flexibility in the N- and C-terminal regions of TP53. Notably, moderate fluctuations were observed near residues 150 and 250 in the TP53–Ellagic acid and TP53–Quercetin complexes, reflecting ligand-induced increases in local dynamic flexibility (Figure 7E).

To evaluate binding stability, we further analyzed the time-dependent evolution of hydrogen bonds, atomic contacts within 0.4 nm, and the minimum distance between TP53 and each ligand. The TP53–Ellagic acid complex maintained 1–4 hydrogen bonds throughout the simulation, with a marked increase after 40 ns, suggesting progressive stabilization (Figure 7F). The TP53–Emodin complex exhibited a consistently higher hydrogen-bond count (2–6), whereas the TP53–epicatechin complex fluctuated between 1–4 hydrogen bonds (Figure 7F). The TP53–Quercetin complex formed the most extensive hydrogen-bond network, reaching up to 7–8 interactions, highlighting its strong polar binding characteristics (Figure 7F).

Atomic contact analysis showed that the TP53–Quercetin and TP53–Emodin complexes displayed the highest numbers of close contacts (~250–350 atoms), indicative of compact binding interfaces, while TP53–Ellagic acid and TP53–Epicatechin formed fewer contacts (~150–250 atoms) (Figure 7G). The minimum distance between TP53 and all four ligands remained <0.3 nm, confirming tight binding throughout the simulation. Among these, the TP53–Ellagic acid complex showed larger fluctuations in minimum distance, implying possible binding-mode adjustments, whereas the TP53–Quercetin complex maintained a relatively stable minimum distance, further supporting its strong binding affinity (Figure 7H).

To assess thermodynamic stability and conformational diversity, Gibbs free-energy (GFE) landscapes were constructed using principal component analysis (PCA). Distinct low-energy basins were observed for all complexes (Figure 8A–D). The TP53–Ellagic acid complex exhibited the narrowest and deepest free-energy basin, reflecting high conformational stability. In contrast, the TP53–Emodin complex displayed a broader basin, suggesting the presence of multiple favorable low-energy states and greater conformational adaptability. The TP53–Epicatechin and TP53–Quercetin complexes showed basins of moderate depth and width, representing an intermediate balance between flexibility and stability.

We further quantified the binding energetics of TP53–ligand interactions using MMPBSA decomposition analysis (Figure 8E–H). All four ligands exhibited negative total binding free energies, indicating thermodynamically favorable binding: Emodin (ΔG_total = −20.44 kcal·mol^−1^), Ellagic acid (−11.76 kcal·mol^−1^), Epicatechin (−13.71 kcal·mol^−1^), and Quercetin (−11.64 kcal·mol^−1^). Gas-phase contributions—including van der Waals and electrostatic interactions—were the primary sources of favorable binding (ΔG_gas: Ellagic acid = −38.90, Emodin = −68.16, Epicatechin = −44.09, Quercetin = −56.49 kcal·mol^−1^). The polar solvation term (ΔEPB) provided an opposing, unfavorable contribution (30.08, 50.49, 33.69, and 47.91 kcal·mol^−1^, respectively), while nonpolar solvation energies contributed slightly favorable effects (≈−2.9 to −3.3 kcal·mol^−1^). Overall, the balance between strong gas-phase interactions and partially offsetting polar solvation energies led to consistently favorable total binding energies for all four complexes.

3.7. Highly Correlated Active Ingredients in RRJ Alleviate NaF-Induced p53 Activation by Inhibiting p53 Expression and Promoting Proteasomal Degradation

Finally, we evaluated the transcriptional and translational regulation of p53 following NaF exposure and examined the modulatory effects of the highly relevant active compounds in RRJ (Quercetin, Epicatechin, Emodin, and Ellagic acid). qPCR analysis showed that NaF markedly increased p53 mRNA expression, whereas all four compounds significantly suppressed this transcriptional upregulation (Figure 9A). Western blot results further confirmed that NaF substantially elevated p53 protein levels, while treatment with any of the four active components significantly attenuated this increase (Figure 9B,C). To determine whether NaF-induced p53 accumulation is transcriptionally driven, Pearson correlation analysis was performed between p53 mRNA and protein levels. A significant positive correlation was observed (r = 0.603, p = 0.0081; Figure 9D), indicating that NaF-induced p53 upregulation is partly transcriptionally mediated, with RRJ-derived compounds attenuating p53 activation through coordinated mRNA and protein suppression. Consistently, immunofluorescence showed pronounced nuclear p53 accumulation in NaF-treated cells, which was attenuated by antioxidant intervention (Figure 9E,F).

To further elucidate the mechanism underlying NaF-induced p53 accumulation, we performed CHX chase and CHX+MG132 inhibition assays to assess p53 protein stability and proteasome-dependent degradation. The CHX chase experiments demonstrated that NaF markedly extended the half-life of p53, indicating impaired degradation (Figure 10A–D), whereas Quercetin, Epicatechin, Emodin, and Ellagic acid restored the normal degradation rate and reversed NaF-induced p53 stabilization. Moreover, results from the CHX+MG132 assays showed that NaF inhibited proteasome-mediated p53 degradation, while the highly relevant active compounds in RRJ effectively reinstated this degradation process (Figure 10E–H).

Collectively, these findings demonstrate that NaF induces abnormal p53 accumulation by simultaneously enhancing its transcription and suppressing its proteasome-dependent degradation. In contrast, the highly correlated active ingredients in RRJ effectively counteract this pathological process through a dual regulatory mechanism, thereby protecting cells from fluoride-induced injury.

3.8. Qualitative and Quantitative Analysis of Four Active Ingredients in RRJ

Based on the integrated mechanistic analyses described above, four RRJ constituents showing the highest correlation with fluoride-responsive targets—Quercetin, Emodin, Epicatechin, and Ellagic acid—were subjected to qualitative and quantitative characterization. Authentic reference standards were employed for compound annotation. Quantitative HPLC analysis confirmed the presence of epicatechin and Ellagic acid in RRJ. Method validation results demonstrated excellent linearity for both epicatechin and ellagic acid, with correlation coefficients (R^2^) greater than 0.999. The LOD and LOQ values indicated high sensitivity of the method. Intra-day precision (RSD%) was 2.26% for epicate-chin and 8.31% for ellagic acid, while inter-day precision was below 10% for both compounds. The average recoveries ranged from 97% to 103%, indicating satisfactory accuracy. Detailed validation parameters are summarized in Supplementary Table S3. Finally, the concentrations of Epicatechin and Ellagic acid were measured to be 41.11 mg/mL and 0.63 mg/mL, respectively (Figure 11). In contrast, Quercetin and Emodin were not directly detected.

Furthermore, qualitative analysis was conducted on Epicatechin, Ellagic acid, as well as Quercetin and Emodin. RRJ was analyzed using ultra-high-performance liquid chromatography coupled with high-resolution mass spectrometry (UHPLC–HRMS) in both positive and negative ionization modes. A total of 130 chemical constituents were identified in RRJ. Epicatechin and Ellagic acid were consistently annotated in the UHPLC–HRMS dataset; qualitative analysis confirmed the presence of Epicatechin and Ellagic. In the qualitative detection, Quercetin and Emodin were not directly detected, and a derivative of Quercetin, namely, Quercetin 3,4′-diglucoside, was detected. Its glycosylated derivative is likely to be metabolized by liver enzymes to form free quercetin. It is highly likely that quercetin exists in the blood when the RRJ is absorbed by the human body, as detailed in Table A1 (Appendix A).

Collectively, these data establish epicatechin and Ellagic acid as the predominant quantifiable constituents among the highly correlated RRJ compounds, while Quercetin- and Emodin-associated species were present mainly in conjugated forms, providing a chemical basis for subsequent functional and mechanistic investigations.

4. Discussion

Our study demonstrates that polyphenolic constituents of Rosa roxburghii juice (RRJ) confer significant protection against fluoride-induced osteoarthritic damage through modulating oxidative stress, inflammatory pathways, and p53 proteostasis. These findings provide mechanistic insights into the therapeutic potential of medicine–food homology agents for managing fluoride-associated osteoarthritis (OA).

Integrated transcriptomic profiling identified TP53, FOS, BAX, and AKT1 as key fluoride-responsive hub genes, consistent with previous reports that implicate p53 activation in chondrocyte senescence, mitochondrial dysfunction, and cartilage degeneration in OA [25]. This convergence highlights p53 as a central stress response mediator in both fluoride toxicity and cartilage pathology.

Computational docking analyses further revealed that Quercetin, Epicatechin, Emodin, and Ellagic acid form stable interactions with the p53 protein, with Emodin exhibiting the strongest predicted binding affinity, corroborating structural studies indicating that natural polyphenols can directly modulate p53 conformation and stability [26]. In parallel, limited proteolysis mass spectrometry (LiP–MS) demonstrated widespread protein structural remodeling following RRJ treatment, particularly in pathways related to oxidative phosphorylation and proteostasis, consistent with the sensitivity of LiP–MS to capture ligand-induced conformational dynamics [27].

Functionally, sodium fluoride (NaF) exposure induced ROS accumulation, G1 phase arrest, and cytotoxicity in rat synovial fibroblasts, corroborating earlier evidence of fluoride-induced oxidative injury and mitochondrial damage [28,29,30,31].

Mechanistically, NaF exposure stabilized p53 by impairing its proteasomal degra-dation, a finding consistent with the literature, which suggests that fluoride disrupts proteostasis and cellular stress responses [32,33,34]

Several limitations of the present study should be acknowledged. First, the in vitro model employed does not fully recapitulate the biomechanical loading, immune interactions, and multicellular complexity characteristic of in vivo osteoarthritic tissues. Second, only a subset of LiP–MS-predicted targets was experimentally validated, and broader target confirmation is required in future investigations. Although in silico modeling suggested potential direct interactions between RRJ-derived polyphenols and p53, definitive structural validation, including biophysical binding assays and site-directed mutagenesis, remains necessary.

In addition, it should be noted that Quercetin and Emodin were not detected in the final quantitative analyses of RRJ, nor were they identified as free aglycones in UHPLC–HRMS profiling. This absence does not necessarily indicate that these constituents are precluded from their presence in RRJ, as differences in juice processing, sample preparation, or compound stability may influence their detectability. Accordingly, further targeted analyses and optimized extraction strategies are warranted to clarify the presence, forms, and bioavailability of Quercetin- and Emodin-related constituents in RRJ. Future studies should validate these findings using rodent models of fluoride-induced OA and employ structural proteomics to delineate specific polyphenol–protein interactions. Pharmacokinetic profiling will also be essential to determine which bioactive compounds most effectively accumulate in joint tissues. Collectively, our results support the promise of multi-component natural products such as RRJ as safe and multi-target therapeutic agents against fluoride-associated OA.

5. Conclusions

This study systematically elucidates the mechanistic basis by which RRJ mitigates F-OA, highlighting Epicatechin and Ellagic acid as key correlated constituents. Integrated transcriptomic, proteomic, and biochemical analyses identified TP53 as a potential core hub gene in F-OA pathogenesis. RRJ-associated constituents exhibited high predicted affinity toward TP53/p53 and induced protein conformational remodeling. Functionally, these compounds attenuated fluoride toxicity by suppressing p53 transcription and promoting proteasome-mediated p53 degradation, thereby limiting aberrant p53 activation.

Collectively, these findings provide a mechanistic rationale for the development of RRJ as a functional food for the prevention and intervention of fluoride-associated skeletal disorders and underscore the multi-target therapeutic potential of medicine–food homology substances in chronic toxicological diseases.