The murine pre-osteoblastic MC3T3-E1 cell line (passages 5–15) was obtained from Wuhan PriCells Biotechnology Co., Ltd. The cells were routinely cultured in α-MEM complete medium (PriCells, Wuhan) containing 10% premium-grade fetal bovine serum (FBS, PriCells, Wuhan) and 1% penicillin-streptomycin solution (100 U/mL, PriCells, Wuhan). Cultures were maintained at 37°C in a humidified atmosphere with 5% CO₂. TPhP Treatment: Triphenyl phosphate (TPhP) (Sigma-Aldrich) was dissolved in DMSO (Sigma-Aldrich) to create stock solutions and diluted in α-MEM complete medium to the required concentrations. MC3T3-E1 cells were treated with TPhP at 0 (vehicle control), 10, 20, and 40 μM for subsequent functional assays. For pathway validation experiments, cells were treated with TPhP (40 μM) alone or in combination with the MEK/ERK activator C16-PAF (MCE, at 1 μM) or the MEK inhibitor MEK-IN-6 (MCE 10 μM). Cells were co-incubated for 48 hours before collection for Transwell, CCK-8, or Western blot analysis.

The standardized chemical structure and SMILES formula of TPhP were retrieved from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). A ligand-based multi-target prediction strategy was employed, utilizing the SWISS Target Prediction platform (http://www.swisstargetprediction.ch/) and the SEA-DOCK reverse docking system for high-throughput screening of potential targets. The species parameter was set to “Homo sapiens,” with a prediction probability threshold >0.6. The identified targets were imported into the UniProt database (https://www.uniprot.org/) for standard annotation and functional classification, establishing a dataset of potential TPhP targets.

The DisGeNET disease-gene association database (https://www.disgenet.org/) and the TTD therapeutic target database (http://db.idrblab.net/ttd/) were integrated for systematic searches using keywords “osteoporosis,” “bone metabolism,” and “bone remodeling.” Highly relevant genes involved in bone metabolism regulation were screened based on evidence-level scoring (relevance score ≥0.4) to establish a bone metabolism-specific target library. Overlapping genes between TPhP targets and bone metabolism-specific targets were identified using Venn diagram analysis on the BioVenn online platform, forming the candidate key targets for TPhP-mediated regulation of bone metabolism.

The overlapping genes identified were imported into the STRING database (https://string-db.org/), with a combined confidence threshold >0.4 to construct the protein-protein interaction (PPI) network. Cytoscape software (v3.10.1) was used for network visualization and topological analysis. Functional modules were identified using the MCODE clustering algorithm (degree cutoff = 2, node score cutoff = 0.2, k-core = 2). Core hub genes were determined using 12 ranking algorithms provided by the CytoHubba plugin, including Degree, Betweenness, and Closeness centrality measures.

Functional enrichment analysis of TPhP’s potential targets was performed using the Metascape platform (https://metascape.org/) and the DAVID annotation system (https://david.ncifcrf.gov/). GO analysis covered three aspects: Biological Processes (BP), Cellular Components (CC), and Molecular Functions (MF). Fisher’s exact test was used to evaluate significance (P < 0.05) with false discovery rate (FDR) correction. KEGG pathway enrichment focused on signaling pathways relevant to bone metabolism regulation, such as the Wnt, BMP, and MAPK pathways.

Molecular docking was performed to investigate the binding characteristics between TPHP and six core target proteins (IGF1R, NR3C1, MAP3K1, BRAF, WNK4, and CNR2) using AutoDock software. The three-dimensional structures of these proteins were obtained from the Protein Data Bank (PDB). Prior to docking, all protein structures were preprocessed by removing water molecules and adding hydrogen atoms. The TPHP structure was energy-minimized before docking simulation. The binding sites were defined based on the surface cavities of each protein. The docking parameters were set to achieve optimal conformational sampling, and the most stable binding conformations were selected based on the lowest binding energy. The binding energy threshold was set at -5 kcal/mol to indicate excellent binding affinity. PyMOL software was utilized to visualize and analyze the three-dimensional configurations of the most stable binding conformations between TPHP and each target protein. The molecular interactions and binding modes were thoroughly examined to evaluate the stability of the protein-ligand complexes.

Core osteoporosis-related genes were identified using the LASSO (Least Absolute Shrinkage and Selection Operator) regression algorithm. The LASSO regression model was constructed using R software (version 4.1.2) with the “glmnet” package. The penalty parameter λ was determined through 10-fold cross-validation to minimize mean squared error. The optimal λ value was selected when the model achieved the best balance between bias and variance. To validate the screening results, receiver operating characteristic (ROC) curve analysis was performed for both the identified candidate genes and the established osteoporosis model. The area under the curve (AUC) was calculated to evaluate the predictive accuracy of the model, with an AUC threshold of 0.7 considered to indicate significant discriminatory power. The validation process was conducted using the “pROC” package in R.

A wound healing assay was conducted to assess the effect of TPhP on MC3T3-E1 cell migration. Cells were cultured to 90–95% confluence, and a uniform scratch was created on the monolayer using a sterile 200 μL pipette tip. Detached cells were gently removed by washing twice with PBS. The cells were then treated with TPhP at 0 (vehicle control), 10, 25, 50, and 100 μM. Images of wound closure at 0 h and 48 h were captured using an inverted microscope equipped with a CCD camera (Nikon Ti2-U, Japan). Image-Pro Plus 6.0 software was used to measure wound areas, and cell migration rates.

The CCK-8 assay was used to assess the effects of TPhP on cell proliferation. MC3T3-E1 cells in the logarithmic growth phase were seeded into 96-well plates at a density of 1×10³ cells/well. After adherence, the cells were treated with TPhP at concentrations ranging from 0 to 100 μM. At 24, 48, 72, 96, and 120 hours, 10 μL of CCK-8 reagent (Beyotime, Wuhan, China) was added to each well and incubated at 37°C for 2 hours in the dark. Absorbance was measured at 450 nm using a microplate reader. Each concentration was tested in triplicate, and the experiments were independently repeated three times.

A Transwell chamber assay was used to evaluate TPhP’s effects on cell invasion and migration. For the invasion assay, Matrigel (BD, 1:8 dilution) was precoated onto polycarbonate membranes (8 μm pores, Corning) and incubated overnight at 4°C. Cells pre-treated with serum starvation for 12 h were seeded into the upper chamber at 2×10⁴ cells/well in 200 μL serum-free medium containing different concentrations of TPhP. The lower chamber was filled with 600 μL complete medium containing 10% FBS as a chemoattractant. After 24 h, cells on the membrane’s upper surface were gently removed, while invaded cells on the lower surface were fixed with methanol, stained with crystal violet, and counted in five randomly selected fields under a microscope. Migration assays followed the same procedure without Matrigel coating.

Cells treated with TPhP were lysed in RIPA buffer containing protease and phosphatase inhibitors (Beyotime, Wuhan) on ice for 30 min. Lysates were centrifuged at 12,000g for 15 min to collect the supernatant. Protein concentration was quantified using a BCA assay, and 30 μg of denatured protein per sample was resolved via SDS-PAGE and transferred onto PVDF membranes (0.45 μm, Millipore) using a semi-dry transfer system. After blocking with 5% BSA in TBST for 2 h, membranes were incubated overnight at 4°C with specific primary antibodies. After washing, HRP-conjugated secondary antibodies (1:5000) were applied for 1 h at room temperature. Signals were developed using enhanced chemiluminescence substrate (Thermo Scientific™ SuperSignal™) and visualized on a ChemiDoc™ MP Imaging System (Bio-Rad). Densitometric analysis was performed using ImageJ software, with β-actin as the loading control. All primary antibodies were obtained from Cell Signaling Technology (CST): p-ERK1/2 (CST, 1:1000); ERK1/2 (CST,1:1000); p-MEK (CST,1:1000); MEK (CST,1:1000); p-P38 (CST,1:1000);P38 (CST,1:1000); p-JNK (CST,1:1000); JNK (CST,1:1000; E-Cadherin (CST,1:1000); N-Cadherin (CST,1:1000); GAPDH (CST,1:1000)

Total RNA was extracted using the TRIzol method (Yeasen, Shanghai), and 1 μg of RNA was reverse-transcribed into cDNA using the PrimeScript™ RT Master Mix (Yeasen, Shanghai). Gene expression was quantified using SYBR Green dye on an ABI QuantStudio™ 3 system. GAPDH was used as an internal control, and relative expression was calculated using the 2^−ΔΔCt method. Each sample was tested in triplicate, with experiments independently repeated three times. The primers for all targets are listed in Table 1 .

Quantitative data were analyzed and plotted using GraphPad Prism 10.0 software. Results were expressed as mean ± standard deviation (Mean ± SD). Statistical comparisons between two groups were performed using independent sample tt-tests, while one-way analysis of variance (ANOVA) followed by Tukey’s post hoc test was used for multiple group comparisons. A two-tailed P-value < 0.05 was considered statistically significant. Significance levels were denoted as follows: *P < 0.05, **P < 0.01, and ***P < 0.001. All experiments were conducted at least three times independently.

Through comprehensive research, we conducted systematic screening, initially predicting and evaluating TPHP toxicity in PubChem, which revealed active responses to estrogen receptor α (ER) ( Supplementary Figure 1 ). It is well established that endogenous estrogen disruption is closely linked to osteoporosis, and previous studies have verified TPHP’s impact on human hormone levels (21, 22). These findings indicate a strong correlation between TPHP and osteoporosis. We initially identified 1,659 potential TPHP-related targets through CHEMBL, STTCH, and Swiss Target Prediction databases ( Figure 1A ). Subsequently, we screened 970 osteoporosis-related genes from GeneCard, OMIM, and DisGeNET databases ( Figure 1B ). After thorough analysis and deduplication, we identified 78 overlapping target molecules, which were considered crucial candidates for TPHP-regulated bone metabolism. The distribution of TPHP and osteoporosis target molecules was visualized using a Venn diagram ( Figure 1C ).

Using the STRING database, we constructed a PPI network of osteoporosis-related proteins and visualized their topological features using Cytoscape software. The network diagram clearly illustrated the correlation between relevant proteins, TPHP, and osteoporosis ( Figure 2A ). To further identify key proteins in TPHP-induced osteoporosis, we first conducted an in-depth analysis of osteoporosis-related gene sets in the GEO database to identify differentially expressed genes ( Figure 2B ). By intersecting these differential genes with the aforementioned osteoporosis-related genes, we ultimately identified 47 target molecules significantly influenced by TPHP in osteoporosis ( Figure 2C ).

To investigate the biological functions and signaling pathways of the 47 differentially expressed target proteins ( Figure 3A ), we performed detailed Gene Ontology (GO) functional analysis using the DAVID database, focusing on human-specific aspects. The analysis revealed three major categories: biological processes, cellular components, and molecular functions. We ranked these functions by significance level and visualized the top 10 most significant entries from each category ( Figures 3B, C ). The results indicated that these molecules primarily participate in biological processes including enzyme activity regulation, organ development, signal transduction, and protein modification; in terms of cellular structure, they are mainly distributed in specialized membrane structures, enzyme complexes, and cell junctions; regarding molecular functions, they primarily involve various kinase activities and growth factor binding. Additionally, KEGG pathway enrichment analysis visualized the signaling pathways involving these molecules, which were primarily enriched in MAPK signaling pathway, cellular signal transduction, protein modification, and metabolic regulation pathways ( Figures 3D, E ).

We employed LASSO regression machine learning algorithm to screen core osteoporosis-related genes. The LASSO regression algorithm established an osteoporosis-related model, identifying six disease-related candidate core target genes (IGF1R, NR3C1, MAP3K1, BRAF, WNK4, CNR2) ( Figures 4A, B ). To verify the accuracy of the selected candidate genes, we conducted ROC tests on both the candidate genes and the osteoporosis model. The ROC curve results showed AUC values greater than 0.7, confirming the scientific validity of the osteoporosis model construction ( Figures 4C, D ).

To explore the interaction characteristics between TPHP and the six key target proteins (IGF1R, NR3C1, MAP3K1, BRAF, WNK4, and CNR2), we conducted simulation analysis using AutoDock software. It is well known that binding energy below 0 kcal/mol indicates spontaneous binding between receptor and ligand without external energy, while binding energy below -5 kcal/mol indicates excellent binding. The molecular docking results showed that TPHP could tightly bind to the surface cavities of IGF1R, NR3C1, MAP3K1, BRAF, WNK4, and CNR2, with binding energies all below -5 kcal/mol (specifically -7.5, -6.2, -6.1, -8.9, -8.9, -6.4, -9.4 kcal/mol, respectively), indicating high binding stability. Using PyMOL software’s three-dimensional visualization technology, we clearly demonstrated the most stable binding conformations ( Figure 5 ), further confirming their tight interactions.

To investigate the direct effects of TPHP on bone cells, we selected MC3T3-E1 mouse osteoblasts as our research model. Cells were exposed to TPHP at concentrations of 0, 10, 20, and 40 μM for 48 hours, and the effects were evaluated through scratch assays and CCK-8 detection. Results showed that TPHP significantly inhibited MC3T3-E1 cell migration in a concentration-dependent manner ( Figures 6A, B ). Cell viability assays further confirmed that TPHP exhibited significant inhibitory effects on cell proliferation, particularly pronounced at higher concentrations ( Figure 6C ). Migration and invasion assays verified TPHP’s inhibitory effect on bone cell migration and invasion capabilities, while Western blot results indicated that TPHP could regulate EMT-related protein expression ( Figure 6D ). Based on previous pathway enrichment analysis, we further investigated whether TPHP affects bone cell function through the MAPK signaling pathway. Next, we explored its molecular mechanisms. Western blot analysis ( Figure 6F ) showed that TPhP treatment led to a significant decrease in the phosphorylation levels of p-MEK and p-ERK1/2 within the MAPK pathway, while p-P38 and p-JNK levels remained unchanged, suggesting that TPhP might specifically inhibit the MEK/ERK signaling axis. Concurrently, TPhP treatment upregulated the expression of the epithelial marker E-Cadherin and downregulated the mesenchymal marker N-Cadherin, indicating a potential impact on the cellular Epithelial-Mesenchymal Transition (EMT) process.

To validate the key role of the MAPK/ERK pathway, we employed a MEK/ERK activator (C16-PAF) and an inhibitor (MEK-IN-6). Results from Transwell ( Figure 6G ) and CCK-8 ( Figure 6H ) assays demonstrated that C16-PAF significantly reversed the TPhP-induced inhibition of cell invasion/migration and proliferation, indicating that reactivating the MEK/ERK pathway could “rescue” the negative effects of TPhP. In contrast, the MEK inhibitor MEK-IN-6 failed to reverse TPhP’s effects. Correspondingly, Western blot results ( Figure 6I ) confirmed that C16-PAF restored the TPhP-suppressed phosphorylation levels of p-ERK1/2 and p-MEK1/2 and reversed the changes in E-Cadherin and N-Cadherin expression. Finally, we used qPCR to examine the effect of TPhP on the expression of the six previously identified key target genes. The results showed that, compared to the control group, 40 μM TPhP treatment significantly downregulated the mRNA levels of NR3C1, IGF1R, MAP3K1, BRAF, WNK4, and CNR2 ( Figure 6J ).

In summary, the results presented in Figure 6 visually demonstrate that TPhP impairs the proliferation, migration, and invasion functions of MC3T3-E1 osteoblasts by inhibiting the MAPK/ERK signaling pathway and downregulating associated gene expression.