Your new experience awaits. Try the new design now and help us make it even better

ORIGINAL RESEARCH article

Front. Toxicol., 03 September 2025

Sec. Environmental Toxicology

Volume 7 - 2025 | https://doi.org/10.3389/ftox.2025.1617663

Exploring plasticisers-osteoporosis links and mechanisms: a cohort and network toxicology study

Xingyao YangXingyao Yang1Xin WangXin Wang2Shifu BaoShifu Bao1Zhengjiang LiZhengjiang Li1Shuxing XingShuxing Xing1Zhangzhen Du,
Zhangzhen Du1,3*
  • 1Department of Orthopaedics, Chengdu Fifth People’s Hospital, Chengdu, China
  • 2Grade 2021 of School of Medical Imaging, North Sichuan Medical College, Nanchong, China
  • 3Beijing Institute of Basic Medical Science, Beijing, China

Background: Plasticisers, widely present in daily life, have been linked to osteoporosis (OP), though the precise mechanisms remain unclear.

Methods: This study examined the association between mono (2-ethylhexyl) phthalate (MEHP) and OP using multivariate logistic regression based on NHANES data. Network toxicology identified key targets and pathways involved in MEHP-induced OP. Molecular docking and dynamics simulations validated the stability of MEHP-target interactions. The effects of MEHP on osteogenic differentiation were further assessed in mouse bone marrow stromal cells (BMSCs).

Results: All logistic regression models confirmed a significant positive correlation between MEHP levels and OP. Network toxicology analysis identified CTSD, SOAT1, and VCP as key targets and the apoptosis pathway as a key mechanism in MEHP-induced OP. Molecular simulations demonstrated stable MEHP binding to these targets. Cellular experiments revealed that MEHP significantly inhibited BMSC osteogenesis by downregulating CTSD and VCP, while SOAT1 showed a weaker correlation.

Conclusion: MEHP exposure is positively associated with OP risk, with CTSD, VCP, and the apoptosis pathway potentially playing key roles in impairing BMSC osteogenesis.

1 Introduction

Plasticisers are a class of chemicals widely used in industry and daily life. In some market financial reports, the use of phthalates accounts for more than 70% of the total phthalate plasticiser market. Among phthalate plasticiser, di (2-ethylhexyl) phthalate (DEHP) once dominated the market, accounting for over 70%. However, due to concerns regarding health and environmental impacts, its usage has gradually decreased in recent years, and it now accounts for about 50% of the global market share. Phthalates are primarily used to enhance the flexibility and processing properties of plastic products. These chemicals are most commonly found in plastic products (such as PVC pipes, food packaging, and children’s toys), as well as construction materials (such as flooring, wallpaper), medical devices (such as infusion tubes and blood bags), and cosmetics. Due to their widespread application in both production and consumption, humans inevitably become exposed to phthalates through skin contact, food intake, or air inhalation (Mondal et al., 2022). In the human body, DEHP is metabolized to mono (2-ethylhexyl) phthalate (MEHP) (Yuan et al., 2022).

In recent years, studies have found that phthalates and some alternative plasticisers may have a range of adverse effects on human health, particularly endocrine disruption and developmental effects. These effects are closely associated with various diseases. Phthalates have been shown to be linked to reproductive system disorders (such as infertility), metabolic diseases (such as obesity and diabetes), cardiovascular diseases, and certain cancers (Mariana et al., 2023; Hlisníková et al., 2020; Guo et al., 2023; Basso et al., 2022). However, there are few studies on the relationship between phthalates and bone health, especially OP, and existing research lacks systematic exploration.

Given that OP is a public health issue closely related to an aging population, characterized by reduced bone density and increased fracture risk, it is particularly important to study the potential impact of environmental exposure factors (such as MEHP) on OP. This study aims to investigate whether MEHP may affect bone metabolism through endocrine disruption, metabolic disorders, or other mechanisms, providing new scientific evidence for the prevention and treatment of OP. The flow chart of this study is shown in Figure 1.

Figure 1
Flowchart illustrating a study on the association between Mono(2-ethylhexyl) phthalate (MEHP) and osteoporosis (OP). It includes steps like cohort study, chemical and disease target mining, enrichment analysis, and identification of key targets. Key elements include GO and KEGG analyses, CytoHubba and machine learning algorithms, and molecular docking. Experimental validation involves CCK8 assay and ALP staining, with results on MEHP concentration and its effect on protein levels of CTSD and VCP. Diagrams and charts detail ligand-protein interactions, nomograms, and ROC curves.

Figure 1. experimental procedure.

2 Methods

2.1 The relationship between MEHP exposure and osteoporosis risk: insights from NHANES

2.1.1 Clinical sample screening

This cross-sectional study leveraged data from the U.S. National Health and Nutrition Examination Survey (NHANES; http://www.cdc.gov/nchs/nhanes.htm), a program providing nationally representative insights into the health and nutritional status of the U.S. population. The study protocol received approval from the National Center for Health Statistics (NCHS) Ethics Review Board, with the latest approval obtained in August 2022, and all participants provided written informed consent. We analyzed data from the 2013–2018 NHANES cycles, covering three survey waves. Out of an initial 29,400 participants, exclusions were applied for individuals under 20 years of age (n = 12,343), missing MEHP measurements (n = 11,853), and missing OP data (n = 2,358). The resulting analytical cohort included 2,846 individuals aged 20 years and older.

2.1.2 Measurement of variables, criteria, screening

Urine specimens were analyzed for phthalate content using high-performance liquid chromatography coupled with electrospray ionization tandem mass spectrometry (HPLC-ESI-MS/MS). The primary DEHP metabolite measured in our study was MEHP, a sensitive and representative biomarker of DEHP exposure (Wang et al., 2019). For phthalate concentrations below the lower limit of detection (LLOD), values were imputed as LLOD divided by the square root of 2, with the LLOD for MEHP established at 0.8 ng/mL. Detailed laboratory protocols are available in the Centers for Disease Control and Prevention (CDC) Laboratory Procedure Manual: Phthalates and Phthalate Alternative Metabolites.

Bone mineral density (BMD), reported in grams/cm2, was assessed using dual-energy X-ray absorptiometry (DXA) with Hologic QDR-4500A fan-beam densitometers at NHANES mobile examination centers. Scans focused on the proximal femur of the left hip, measuring BMD at the total femur, femoral neck, and trochanter; the right hip was scanned if left-side replacements or metal implants were reported. Exclusion criteria aligned with NHANES guidelines, omitting participants who were pregnant, weighed over 300 pounds, or had bilateral hip conditions such as fractures, replacements, or metal pins.

OP was defined in accordance with World Health Organization (WHO) guidelines as bone mineral density (BMD) values that are 2.5 standard deviations or more below the mean BMD of a young adult reference population.

The covariates in this study included age (year), sex (male/female), race (American/White/Black/other), marital tatus (married/living with a partner/widowed/divorced/separated/never married), education attainment (Less than 9th grade/high school graduate or equivalent/some college or AA degree/9–11th grade and college graduate or above), poverty income ratio (PIR), smoking status (never, former and now), drinking status (never, former, mild, moderate, heavy), hypertension (yes/no), diabetes (yes/no), Cardiovascular diseases (CVD, yes/no) and body mass index (BMI,kg/m2). Details of each variable were publicly available at www.cdc.gov/nchs/nhanes/.

2.1.3 Statistical analysis

Statistical analyses followed Centers for Disease Control and Prevention (CDC) guidelines, applying NHANES sampling weights recalculated according to National Center for Health Statistics (NCHS) recommendations to accommodate the survey’s complex design. Weighted chi-square tests were used for categorical variables, and one-way analysis of variance (ANOVA) was applied for continuous variables. The relationship between MEHP and OP was explored through multivariable logistic regression across three models. Model 1 was unadjusted; Model 2 adjusted for sex, age, and race; Model 3 included additional adjustments for marital status, education, PIR, smoking status, drinking status, hypertension, diabetes, CVD and BMI. The investigation of potential dose–response associations between the MEHP and OP was conducted using a restricted cubic spline (RCS) model. For two-sided tests, a p-value less than <0.05 was deemed to denote statistical significance. Data processing and statistical analyses were performed in R version 4.1.3. For cell experiments, data are presented as mean ± standard deviation (SD). For Western blot and ALP staining experiments, each assay was independently repeated three times. For qPCR experiments, a sample size of n = 3 was used. Statistical analysis was performed using one-way or two-way ANOVA, followed by Tukey’s post hoc test. A p-value <0.05 was considered statistically significant. All statistical analyses were conducted using GraphPad Prism software (version 7.0).

2.2 Advanced exploration of the potential targets and pathways for MEHP in OP

The dataset “GSE56815” was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The dataset contains 40 control samples and 40 disease samples. The data were normalized, and DEGs between the OP group and the control group were identified using the “Limma” package based on the following criteria: p-value <0.05 and |log2FC| >0.1. A volcano plot was generated using the “ggplot2” R package to visualize the screening results. The “SMILES” structure of MEHP was obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/), and its potential targets were predicted using the SwissTargetPrediction (http://swisstargetprediction.ch/), SEA (https://sea.bkslab.org/), and SuperPred (https://prediction.charite.de/subpages/target_prediction.php) databases. The targets were then standardized using the String database (https://cn.string-db.org/). A Venn diagram of the intersecting targets between the DEGs and MEHP targets was drawn using the “VennDiagram” package. Finally, the obtained targets were subjected to GO (Gene Ontology) analysis and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway prediction, and their potential involvement in metabolic or signaling pathways was mapped.

2.3 Integrating protein-protein interaction networks with machine learning for key targets screening

Based on the intersecting targets, a PPI network was constructed in the String database (https://cn.string-db.org/) with a minimum required interaction skey of 0.15. The results of the PPI network were analyzed using the cytoHubba plugin in Cytoscape software. The top 10 genes in the network were identified using four algorithms: “Betweenness,” “Degree,” “MCC,” and “Stress” (Du et al., 2024). Finally, the hub targets obtained from the four algorithms were visualized. Based on hub targets, we utilized three machine learning algorithms—Boruta, LASSO, and SVM—to screen for key targets involved in MEHP-induced osteoporosis (Jiang et al., 2024). The Boruta algorithm is a supervised classification feature selection method that identifies all relevant features for the classification task (Lyu et al., 2023). Least Absolute Shrinkage and Selection Operator (LASSO) logistic regression is used to identify and retain the most important genes, thereby improving the predictive accuracy and interpretability of the model. This algorithm was implemented using the “glmnet” package. To reduce the risk of overfitting, we used Support Vector Machine Recursive Feature Elimination (SVM-RFE) to carefully select the optimal genes from the metadata collection. A Venn diagram of the intersecting genes identified by the three machine learning methods was drawn using the “VennDiagram” package as key targets.

2.4 Construction of a nomogram model for key targets

To further evaluate the diagnostic predictive ability of key targets, we constructed a nomogram using the rms software package (Ma et al., 2025). Subsequently, calibration curves and Receiver Operating Characteristic (ROC) curves were utilized to assess the accuracy of the nomogram’s predictions. Additionally, clinical impact curve (CIC) was employed to evaluate the predictive value of the model.

2.5 Molecular docking for MEHP and key target of OP

Molecular docking was employed to investigate the binding affinity between MEHP and the key target proteins (Yang et al., 2025). The three-dimensional structures of the target proteins were retrieved from the UniProt (https://www.uniprot.org/) and Protein Data Bank (PDB, https://www.rcsb.org/) databases. The chemical structure of MEHP in SDF format was obtained from the PubChem database and subsequently converted to mol2 format using Open Babel GUI software. Both protein receptors and ligand molecules were preprocessed using AutoDock Tools, including the removal of water molecules, the addition of hydrogen atoms, and the assignment of partial atomic charges. Blind docking was performed using AutoDock Vina, with the grid box defined large enough to encompass the entire protein surface, allowing for unbiased exploration of potential binding sites. Grid center coordinates and dimensions were determined based on the overall structure of the protein. For each docking simulation, five binding conformations were generated, and the conformation with the lowest binding energy was selected for further analysis. All docking results were visualized using PyMOL.

2.6 Molecular dynamics simulation

To minimize the discrepancies between the conformations obtained from protein-ligand docking and the actual complexes, we conducted molecular dynamics simulations (MD) on key targets-MEHP complexes using GROMACS 2020.6 software (Tuo et al., 2024). We adopted the AMBER99SB-ILDN/GAFF force field for each simulation system, implemented by Sobtop. The initial system was constructed in a dodecahedron box with a 1.2 nm layer between the protein surface and the edge of the box, populated by the TIP3P water model. Each system was neutralized by adding appropriate amounts of Na+ and Cl counter ions. Prior to the MD simulation, energy minimization was executed with the steepest descent algorithm. Then, the canonical (NVT) and isothermal-isobaric (NPT) ensembles were implemented to equilibrate the system for 100 ps. The state-balanced system was configured to maintain a constant temperature of 310 K and a standard pressure of 1.0 bar, along with periodic boundary conditions.

Finally, the system underwent a 100 ns MD simulation to assess complex stability. RMSD, radius of gyration (Rg), and solvent accessible surface area (SASA) were calculated from the trajectory data. The binding free energy of the complexes was estimated using the gmx_MMPBSA tool. Additionally, free energy landscape (FEL) plots were generated to visualize conformational transitions and energy states of the protein-ligand complexes. The FEL visualization script was sourced from Charles Hahn’s open-source repository (https://github.com/CharlesHahn/DuIvy/blob/master/sources/other/PCA_FEL/xpm2png.py).

2.7 Isolation and culture of BMSCs

All animal procedures followed the Animal Welfare Act and were approved by the Ethics Committee of the Beijing Institute of Basic Medical Sciences. Two-day-old neonatal mice were euthanized with an overdose of isoflurane. After 5 min of immersion in 75% ethanol, bones were isolated and cleaned of surrounding tissues. The tibiae and femora were digested with 0.1% type I collagenase to remove soft tissue, then minced and placed in T25 flasks with DMEM (high-glucose) containing 20% FBS and 1% PS. After 2 days, cells growing around the bone fragments were detached with 0.25% trypsin and transferred to new T25 flasks. Cells were cultured until the P3 generation for further experiments.

2.8 Quantitative real-time PCR (qRT-PCR) and alkaline phosphatase (ALP) staining

Cells were cultured in six-well plates, washed twice with 2 mL of cold PBS, and then lysed with 1 mL of TRIzol reagent. After a 5-min incubation at room temperature, the cells were collected into Eppendorf tubes. For micromass or tissue samples, these were quickly frozen in liquid nitrogen, followed by disruption using a multibead shocker for RNA extraction. For the qRT-PCR, 500 ng of RNA was reverse-transcribed into cDNA using ReverTra Ace (Toyobo). PCR amplification was performed in a total reaction volume of 20 μL, including 1 μL of cDNA and 4 μL of SYBR FAST qPCR Master Mix (Tsingke). The primer sequences used for qRT-PCR are provided below. Gene expression was quantified using the comparative Ct method (2−ΔΔCT), with data normalized to β-actin. Each sample was analyzed in triplicate to ensure reliability. Primer sequences are shown in Table 1.

Table 1
www.frontiersin.org

Table 1. Primer sequences for the tested genes.

BMSCs were fixed with 4% paraformaldehyde for 30 min. ALP staining was then performed using a kit (SCR004, Sigma-Aldrich, United States) according to the manufacturer’s instructions.

2.9 Western blot analysis (WB), siRNA transfection and osteogenic differentiation

Cells were lysed using RIPA buffer 1. The extracted proteins were boiled at 100°C for 10 min, separated by 10% PAGE gel, and subsequently transferred onto a PVDF membrane. The membrane was then blocked with 5% skim milk at room temperature for 1 h. Following blocking, it was incubated overnight at 4°C with primary antibodies against BMP2 (bs-0514R, Bioss, China), RUNX2 (BSM-52672R, Bioss, China), GAPDH (10494-1-AP, Proteintech, China), COL1 (14695-1-AP, Proteintech, China), VCP (82463-1-RR, Proteintech, China), and CTSD (21327-1-AP, Proteintech, China). After primary antibody incubation, the membrane was treated with goat-derived HRP-conjugated secondary antibodies at room temperature for 1 h. The membrane was washed five times with TBST between steps. Finally, the protein bands were visualized using the Super ECL Plus Chemiluminescent HRP Substrate (Solarbio, China). siRNA was purchased from Sangon (Shanghai, China) and transfected into BMSCs (50 nM) using Lipofectamine 3000 (Invitrogen, United States). BMSCs were seeded in six-well plates (2 × 10^5 cells/well) 1 day prior. The siRNA and transfection reagent were mixed following the manufacturer’s instructions, incubated at room temperature for 15–20 min to form complexes, and added to the cells. Transfection efficiency was assessed via Western blot. After 36 h, the transfected BMSCs were induced for Osteogenic Differentiation. siRNA sequences are listed below (Table 2).

Table 2
www.frontiersin.org

Table 2. siRNA sequences.

3 Results

3.1 The results of the cohort study

3.1.1 Baseline characteristics of participants

The baseline characteristics of 2,846 participants were summarized according to MEHP quartiles (Supplementary Table S1). The group was evenly divided between men and women, with an average age of 39.046 ± 0.325 years. Significant differences were observed across MEHP quartiles in several demographic and health-related variables, including sex, race, PIR, BMI, education and drinking status (all P < 0.05). Participants in the higher MEHP quartiles were more likely to be male, smokers, and had higher BMI, lower PIR and OP.

3.1.2 Association between MEHP and OP

Table 3 Presents the multivariable regression analysis examining the association between MEHP and OP. The results indicate that higher MEHP is associated with the prevalence of OP. In every Model, participants in Quartile 4 had significantly higher prevalence of osteoporosis compared to those in Quartile 1 (Model1:OR = 3.651, 95% CI: 2.018 to 6.603, p < 0.001, Model2:OR = 4.069, 95% CI: 2.261 to 7.324, p < 0.001, Model3:OR = 4.228, 95% CI: 2.112 to 8.461, p < 0.001).These results indicate a strong positive association between MEHP levels and OP, underscoring the potential role of MEHP in the increased prevalence of OP.

Table 3
www.frontiersin.org

Table 3. Multivariate logistic analysis of the association between mono (2-ethylhexyl) phthalate (MEHP) and osteoporosis (OP). Data are presented as β coefcient (95% CI). Model 1 was unadjusted; Model 2 adjusted for gender, age, and race; Model 3 included additional adjustments for marital status, education, poverty income ratio, body mass index, smoking status, drinking status, hypertension, diabetes, cardiovascular. OR, odd ratio.

3.1.3 Analysis of restricted cubic spline regression

After accounting for multiple covariates, we identified a significant non-linear relationship between MEHP and OP in RCS regression (Figure 2). MEHP were nonlinearly and positively correlated with the prevalence of OP (p for nonlinearity = 0.0089).

Figure 2
A line graph shows the relationship between MEHP levels and osteoporosis, depicting a nonlinear upward trend. The red line curves initially then rises steadily. P-value for nonlinearity is 0.0089, indicating statistical significance. Horizontal axis represents MEHP, and the vertical axis represents osteoporosis levels.

Figure 2. Analysis of restricted cubic spline regression. Restricted cubic spline analysis of the association between mono (2-ethylhexyl) phthalate (MEHP) and osteoporosis.

3.2 Advanced exploration into the identification of potential targets and pathways for MEHP in OP

The volcano plot shows that 440 upregulated genes and 482 downregulated genes were identified through differential analysis (Figure 3A). After obtaining MEHP-related targets from multiple databases, a total of 204 MEHP-associated targets were identified. The Venn diagram shows that there are 19 intersecting targets between the DEGs and MEHP-related targets (Figure 3B) (Supplementary Table S2). The GO functional annotation results of the targets indicate their involvement in 456 biological processes (BP), including the collagen catabolic process (GO:0030574). Within the cellular component (CC) category, significant enrichment was observed in the secretory granule lumen (GO:0034774). In the molecular function (MF) category, serine-type peptidase activity (GO:0008236) was predominantly enriched (Figure 3C). KEGG pathway enrichment analysis identified 17 significantly enriched pathways (P < 0.05), among which the apoptosis pathway may play a crucial role in MEHP-induced OP (Figure 3D). Figure 3E illustrates the pathway map of the intersecting genes involved in apoptosis.

Figure 3
A composite image consists of five panels:A) A volcano plot showing gene expression changes with dots representing genes, separated into downregulated (blue) and upregulated (red) categories.B) A Venn diagram illustrating the overlap between differentially expressed genes (DEGs) and MEHP, with numbers indicating gene counts.C) Three bar charts showing gene ontology categories: biological process (BP), cellular component (CC), and molecular function (MF), with specific terms listed.D) A bubble chart displaying pathways affected, with color indicating category and size for gene count.E) A pathway map with various cellular processes and highlighted areas in red boxes.

Figure 3. Identification of Related Target Genes (A) Volcano plot of DEGs in “GSE56815”. (B) Identification of 19 related target genes. (C) Gene Ontology analysis. BP (biological process), CC (cell composition) and MF (molecular function). (D) Kyoto Encyclopedia of Genes and Genomes pathway prediction. (E) Apoptotic Pathway Diagram with Red-Marked Intersection Targets.

3.3 Identification of key targets

The protein-protein interaction (PPI) network comprises 17 nodes and 57 edges. (Figure 4A). Four algorithms of the cytoHubba plugin were used: “Betweenness,” “Degree,” “MCC,” and “Stress” to identify the top 10 genes in each (Figure 4B). As shown in the Figure, the intersection of the four algorithms resulted in eight hub targets, including MMP8, LDHA, NFKB1, VCP, CTSK, CTSD, CTNNB1, and SOAT1. (Figure 4C). LASSO, SVM-RFE, and Boruta algorithms were used to identify key targets. When applying LASSO based on 10-fold cross-validation, the minimum error value corresponded to 8 key targets, including MMP8, LDHA, NFKB1, VCP, CTSK, CTSD, CTNNB1, and SOAT1 (Figure 4D). The Boruta algorithm identified 6 key targets: MMP8, VCP, CTSK, CTSD, CTNNB1, and SOAT1 (Figure 4E). The SVM-RFE algorithm was also validated using 10-fold cross-validation. The most accurate algorithm with the lowest estimated error identified 5 key targets: CTSD, NFKB1, VCP, SOAT1, and LDHA (Figure 4F). Therefore, based on the results of the above machine learning models, 3 key targets (CTSD, VCP, and SOAT1) were identified (Figure 4G).

Figure 4
A composite image featuring several scientific visualizations: (A) a network graph with connected nodes. (B) Four network diagrams labeled Betweenness, Degree, MCC, and Stress, with nodes and connecting lines. (C) A Venn diagram showing overlaps among four sets, labeled Betweenness, Degree, MCC, and Stress. (D) Two plots: the left shows log lambda values against binomial deviance; the right is a coefficient versus L1 norm plot. (E) A box plot showing attribute importance. (F) Two line graphs: the left shows accuracy versus the number of features; the right shows error versus the number of features. (G) A Venn diagram comparing LASSO, Boruta, and SVM-REF methods.

Figure 4. Identification of Key Targets. (A) PPI network diagram of intersecting target. (B) “Betweenness,” “Degree,” “MCC,” and “Stress” to identify the top 10 targets. Darker Colours indicate larger values for the relevant algorithms. (C) Venn diagram, four algorithms for screening targets. (D) Lasso, plots of ten-fold cross-validation curves and trajectories of their coefficients for the LASSO model, with vertical dashed lines denoting optimal lambda values. (E) Boruta, boxplot of importance distribution of each target in Boruta algorithm. Green represents the screened key targets. (F) SVM-RFE, trend plots of accuracy and cross-validation error for SVM models for optimizing hyperparameter selection. (G) Venn diagram for determining key targets.

3.4 Constructing a nomogram and validating its performance

A nomogram was developed based on three key targets (Figure 5A). The calibration curve results revealed a strong agreement between the observed incidence and the nomogram-predicted probabilities (Figure 5B). The CIC demonstrated that the model provides a favorable overall net benefit in clinical practice (Figure 5C). In the final ROC analysis, an Area Under Curve (AUC) of 0.868 indicated excellent model accuracy (Figure 5D).

Figure 5
Panel A shows a nomogram with variables SOAT1, CTSD, and VCP contributing to total points, predicting a probability of 0.127. Panel B presents a calibration plot contrasting predicted and observed probabilities. Panel C is a decision curve illustrating the number deemed high risk at various thresholds and cost-benefit ratios. Panel D displays a ROC curve with an AUC of 0.868, indicating predictive performance.

Figure 5. Construct a nomogram and validate its performance (A) nomogram of key targets. (B) The calibration curve of nomogram. (C) CIC of nomogram. (D) Receiver Operating Characteristic (ROC) curves of nomogram.

3.5 Molecular docking

Molecular docking helps to understand the binding of proteins expressed by key targets and the ligand MEHP. The binding energies of MEHP with CTSD, SOAT1, and VCP were −5.7, −6.7, and −6.4 (kcal/mol), respectively (Figures 6A–C). These binding energies are all less than −5 (kcal/mol), indicating good binding between the key targets and MEHP (Tuo et al., 2024).

Figure 6
Molecular docking diagrams showing interactions between MEHP and enzymes. A: CTSD-MEHP complex with binding energy (BE) of -5.7 kcal/mol, highlighting interactions at GLY-79. B: SOAT1-MEHP complex with BE of -6.7 kcal/mol, highlighting interactions at SER-456 and HIS-425. C: VCP-MEHP complex with BE of -6.4 kcal/mol, highlighting interactions at GLY-521, GLY-523, SER-524, and THR-525. Each panel includes a 3D structure, detailed interaction view, and a 2D interaction diagram with labeled types: van der Waals, hydrogen bonds, and pi interactions.

Figure 6. Molecular Docking Analysis (A) CTSD-MEHP. (B) SOAT1-MEHP. (C) VCP-MEHP. Purple - proteins, green - MEHP, yellow - hydrogen bonds, orange - hydrogen bond linked amino acids.

3.6 Molecular dynamics simulations

To further investigate the stability of protein–ligand interactions, molecular dynamics simulations were performed on three protein–ligand complexes: CTSD-MEHP, SOAT1-MEHP, and VCP-MEHP. The root mean square deviation (RMSD) was used to assess whether the simulation systems had reached equilibrium, with RMSD values within 1 nm indicating the relative stability of protein–ligand interactions under physiological conditions. As shown in Figure 7A, the RMSD values of the three complexes rapidly stabilized at 0.217 ± 0.036 nm, 0.604 ± 0.085 nm, and 0.317 ± 0.063 nm, respectively, suggesting that all three complexes maintained structural stability throughout the simulation. The radius of gyration (Rg) was analyzed to evaluate the compactness of receptor–ligand binding. As depicted in Figure 7B, the Rg values of the three complexes remained stable throughout the simulation, with final values of 2.078 ± 0.022 nm, 3.337 ± 0.044 nm, and 2.236 ± 0.039 nm, respectively. Additionally, the solvent-accessible surface area (SASA), a crucial parameter reflecting protein folding and stability, was assessed. As shown in Figure 7C, the SASA values remained stable, with average values of 168.925 ± 4.921 nm2, 420.744 ± 14.031 nm2, and 158.723 ± 3.147 nm2, respectively, indicating that the overall structural integrity of the protein–ligand complexes was preserved during the simulation. Furthermore, the binding free energy (∆G_bind) of the three protein–ligand complexes was calculated using the MM/GBSA method over the final 40 ns of the simulation. Lower ∆G_bind values indicate stronger receptor–ligand binding affinity. As shown in Figure 7D, the ∆G_bind values for the CTSD-MEHP, SOAT1-MEHP, and VCP-MEHP complexes were −33.56 kcal/mol, −51.58 kcal/mol, and −29.79 kcal/mol, respectively, suggesting strong binding affinities for all three complexes. In addition, free energy landscape (FEL) analysis revealed the presence of multiple low-energy states during the last 30 ns of the simulation, further confirming the stability of the protein–ligand interactions (Figure 7E).

Figure 7
Graphs and visualizations illustrate molecular dynamics analysis. Panels A and B show line graphs for RSMD and Rg over time for CTSD, SOAT1, and VCP. Panel C presents a graph of SASA versus time for the same proteins. Panel D displays a bar chart of binding free energy components for these proteins, indicating ΔGGAS, ΔGSOLV, and ΔTOTAL values. Panel E includes three 3D surface plots showing free energy landscapes for CTSD-MEHP, SOAT1-MEHP, and VCP-MEHP, with a color gradient indicating energy intensity.

Figure 7. Molecular Dynamics Simulations (A) RMSD values of the three complexes. (B) Radius of gyration (Rg) values of the three complexes. (C) Solvent-accessible surface area (SASA) values of the three complexes. (D) Binding free energies (∆Gbind) of the four complexes. (E) Free Energy Landscape (FEL).

3.7 The effect of MEHP on osteogenic differentiation of BMSCs

The osteogenic differentiation ability of bone marrow mesenchymal stem cells (BMSCs) plays a crucial role in OP. To assess the effect of MEHP on OP, we first performed a CCK8 assay to investigate the inhibitory effect of MEHP on BMSC proliferation. The results showed that 500 ng/mL of MEHP significantly inhibited BMSC proliferation (Figure 8A), with the inhibitory effect increasing as the concentration increased. To further confirm these results, we treated BMSCs with 500, 600, and 700 ng/mL MEHP, and after 1, 3, 5, and 7 days of treatment, a marked decrease in cell proliferation was observed compared to the control group (Figure 8B). Given that excessive inhibition of cell viability could potentially lead to excessive cell death during subsequent osteogenic differentiation, we selected 500 ng/mL of MEHP to investigate its effect on BMSC osteogenesis. After 7 days of osteogenic induction, alkaline phosphatase (ALP) staining, qRT-PCR, and Western blot (WB) analyses revealed that MEHP significantly suppressed osteogenesis in BMSCs, as evidenced by the decreased expression of ALP, type I collagen (COL1), RUNX2, and BMP2 (Figures 8C–G). Molecular docking results indicated that MEHP may influence osteoporosis by interacting with three proteins: CTSD, VCP, and SOAT. To validate this, we examined the correlation between CTSD, VCP, SOAT, and the osteogenic process in BMSCs. qRT-PCR results showed that during BMSC osteogenesis, the expression of COL1, BMP2, and RUNX2 genes significantly increased (Figures 8H–J). We also observed a strong correlation between the gene expression of CTSD and VCP with BMSC osteogenic differentiation, but not with SOAT (Figures 8K–M). Therefore, we focused on the effect of MEHP on CTSD and VCP. WB analysis demonstrated that while MEHP inhibited osteogenesis, it also reduced the protein levels of CTSD and VCP (Figure Figure 8G). To further confirm whether MEHP-mediated suppression of CTSD and VCP proteins correlates with osteogenesis in BMSCs, we used siRNA to knock down CTSD and VCP. The results showed that siRNA treatment significantly impaired BMSC osteogenesis, as indicated by the decreased protein levels of BMP2, COL1, and RUNX2 after 7 days of osteogenic induction (Figures 8N,O), as well as a reduction in ALP expression (Figure 8P). In conclusion, MEHP significantly inhibits BMSC osteogenic differentiation by downregulating CTSD and VCP expression, and the reduction in CTSD and VCP levels is inversely correlated with the osteogenic potential of BMSCs.

Figure 8
A scientific figure with multiple panels showing experimental data. Panel A: Bar graph of CCK-8 assay results at various MEHP concentrations, showing decreased OD values with higher concentrations. Panel B: Line graph of CCK-8 assay results over seven days, showing OD values at different MEHP concentrations. Panel C: Images of ALP staining for NC, OST, and OST+MEHP conditions. Panels D-F: Bar graphs of mRNA expression levels for COL1, RUNX2, and BMP2 across NC, OST, and OST+MEHP groups.Panels G and N: Western blot images showing protein levels of COL1, CTSD, RUNX2, VCP, BMP2, and GAPDH under different conditions.Panels H-M: Bar graphs showing mRNA expression of various genes over different time periods.Panel O: Western blot of proteins under NC, OST, OST+siVCP conditions.Panel P: ALP staining images for NC, OST, OST+siVCP, and OST+siCTSD.

Figure 8. Validation of Key Targets’ Osteogenic Relevance (A,B) CCK8 assay to determine MEHP concentration. (C) Alkaline phosphatase (ALP) staining of BMSCs after MEHP intervention. (D–F) Decreased expression of type I collagen (COL1), RUNX2, and BMP2 after MEHP intervention. (G) MEHP reduced the protein levels of CTSD and VCP, (H–M) Significant increase in the expression of COL1, BMP2, and RUNX2 genes during osteogenesis of BMSCs. (N, O) MEHP-mediated inhibition of CTSD and VCP proteins impaired the osteogenic capacity of BMSCs. (P) ALP staining after siRNA treatment.

4 Discussion

In previous years, reports have indicated that sports equipment in elementary and middle schools has been found to contain phthalates far exceeding safe levels, raising concerns about the impact of phthalates on children’s growth and development. However, such chemicals are inevitably encountered in modern life. Regardless of age, the question remains whether phthalates affect bone growth and development, and how they exert this effect needs further investigation.

In our study, we first selected 2,846 OP patients from the NHANES database to observe the relationship between the human metabolite MEHP and OP. In three predictive models adjusted for different covariates, we found a positive correlation between MEHP levels and OP. To further explore the mechanism by which MEHP contributes to OP, we downloaded the OP dataset GSE56815, which includes multiple samples, and identified DEGs. We then intersected these DEGs with MEHP target genes, identifying 19 intersecting genes. The primary functions of these 19 genes are involvement in the catabolic processes of collagen and osteoclast differentiation, as well as participation in pathways related to apoptosis. Oxidative stress-induced DNA damage, apoptosis, and cellular senescence are potential contributors to osteoporosis. These mechanisms, particularly apoptosis, offer promising new avenues for osteoporosis treatment (Chandra and Rajawat, 2021). In previous studies, scholars have identified two pathways. There was KEGG apoptosis and activation of REACTOME caspase activation via the extrinsic apoptotic signaling pathway (Huang HL. et al., 2023). The “Betweenness,” “Degree,” “MCC,” and “Stress” algorithms were used to determine the top 10 key genes, and finally, Boruta, LASSO, and SVM-RFE algorithms identified three key targets. Additionally, a nomogram was developed based on the key targets associated with MEHP-induced OP. The model exhibited strong predictive performance, facilitating the translation from basic research to clinical practice, and further underscored the critical roles of the three key targets in disease pathogenesis and progression. Molecular docking and dynamics analysis further confirmed that the MEHP ligand could easily dock with the proteins expressed by CTSD, SOAT1, and VCP, and that the binding interactions were stable.

Currently, CTSD has attracted significant attention in cancer research due to its overexpression in various cancers, where it promotes tumor cell proliferation, invasion, and metastasis. As a potential tumor biomarker and therapeutic target, it has been widely studied (Wang X. et al., 2023). In a study by Zheng et al., it was found that CTSD regulates the radiosensitivity of glioblastoma by affecting the fusion of autophagosomes and lysosomes (Zheng et al., 2020). Additionally, CTSD plays a crucial role in neurodegenerative diseases, and alterations in its expression may contribute to Parkinson’s disease, Alzheimer’s disease, and Huntington’s disease (Drobny et al., 2022). In some bioinformatics analyses, CTSD expression has been shown to differ significantly in patients with OP (Wang X. et al., 2022; Huang Y. et al., 2023). Furthermore, experiments involving rat BMSCs have shown that CTSD promotes osteogenesis (Sun et al., 2024). Our experiments also confirmed that CTSD and VCP have some osteogenic promoting abilities.

SOAT1 is currently primarily studied for its role in maintaining cholesterol homeostasis and mediating cholesterol esterification metabolism. It has been reported in diseases such as liver cancer, atherosclerosis, and Alzheimer’s disease (Song et al., 2021). SOAT1 expression is significantly elevated in most cancers and is strongly correlated with prognosis. By increasing cholesterol esterification, it promotes epithelial-mesenchymal transition (EMT) in cancer and contributes to liver cancer development (Fu et al., 2024). SOAT1 is also upregulated in colon cancer, enhancing the migration and invasion abilities of colon cancer cells to promote its progression (Wang XC. et al., 2022). In a study by Sun et al., it was shown that inhibiting SOAT1 could enhance the efficacy of glioma radiotherapy both in vivo and in vitro (Sun et al., 2023). SOAT1 plays a critical role not only in metabolic diseases and cancer but also potentially in growth development and bone health. In our experimental validation, we found no significant association between SOAT1 and osteogenic differentiation. We speculate that SOAT1 may indirectly affect the function of osteoblasts and osteoclasts by regulating intracellular cholesterol levels and lipid metabolism (Kim et al., 2021; Akhmetshina et al., 2023), thereby participating in the formation and remodeling of bones. Previous studies have reported that conditions such as hyperglycemia, hyperlipidemia, and obesity are high-risk factors for OP (Kruse, 2018; Hussain et al., 2023; Guimarães et al., 2024), possibly due to adipose tissue competitively inhibiting bone growth (Deng et al., 2024).

Mutations in the VCP gene are associated with multisystem proteinopathies, a disease spectrum that includes myopathies, bone diseases, neurodegenerative diseases, and motor neuron diseases (Columbres et al., 2023). VCP myopathy often presents with uncommon clinical manifestations, such as ocular muscle paralysis, ptosis, and dysphagia (Guo et al., 2019). As one of the most common metabolic bone diseases, Paget’s disease of bone has a high incidence among Caucasians. The overactivity or dysfunction of osteoblasts may lead to reduced bone density or skeletal deformities, and mutations in the VCP gene have been confirmed (Chung and Van Hul, 2012). VCP dysfunction can cause protein accumulation, which in turn affects osteoblast differentiation and mineralization. Normal osteoblast differentiation relies on the autophagic process (Wang J. et al., 2023; Yoshida et al., 2022), and VCP plays a key role in the lysosome-mediated autophagy pathway (Arhzaouy et al., 2019). With the combined effect of these two factors, mutated VCP severely impacts osteoblast function. Inhibition of VCP, while regulating the NF-κB signaling pathway, also suppresses osteoclast differentiation (Wei et al., 2023).

Although many organizations have recognized the widespread use of phthalates and its potential impact on human health, numerous measures have been taken to limit its usage or find harmless alternatives. In developed countries, phthalates containing DEHP are classified as “Substances of Very High Concern” (SVHC). The EU’s REACH regulation restricts the use of certain types of phthalates (such as DEHP) in toys and children’s care products. However, in developing countries, due to cost factors and delayed regulations, the use of phthalates still dominates. Although there has been a clear global trend toward restrictions and alternatives, and efforts are actively being made to find substitutes, alternative plasticisers (such as citrates and epoxides) are considered safer in certain applications, their long-term exposure effects still require further evaluation.

This study focused on the effects of MEHP on bone health, specifically examining its impact on bone marrow stromal cells (BMSCs). BMSCs were selected due to their essential role as progenitor cells in bone formation and remodeling. They not only directly contribute to osteogenesis but also regulate osteoblast and osteoclast activity, making them central to bone homeostasis. Given their sensitivity to endocrine disruptors like MEHP, which can interfere with key pathways such as PPARγ and Wnt/β-catenin (Zhang et al., 2024; Huo et al., 2024), BMSCs provide an ideal model for understanding the early cellular effects of MEHP on bone metabolism. While osteoblasts and osteoclasts play crucial roles in bone resorption and formation (McNamara, 2021), our focus on BMSCs allows for a deeper understanding of the upstream disruptions in bone homeostasis. Stem cells are known to be particularly sensitive to toxic exposures (Gao et al., 2025; Lin et al., 2024), making BMSCs a suitable target for investigating the initial impacts of MEHP. Although the effects of MEHP on osteoblasts and osteoclasts are important, research in this area remains limited. Therefore, our study centered on BMSCs, with plans to extend future research to include osteoblast mineralization and osteoclast differentiation to provide a more comprehensive understanding of MEHP’s role in bone metabolism.

Ultimately, disease pathogenesis is a complex process. Within the intricate human body, the impact of environmental pollutants and endocrine disruptors warrants further attention from various perspectives. For example, air pollutants contain numerous toxins, which, upon entering the human body, activate immune and inflammatory responses to varying degrees, potentially leading to osteoporosis (Allen et al., 2024). Moreover, air pollutants and phthalates likely share overlapping genes and pathways in the development of OP (Yu et al., 2023). We believe that with the growing awareness of health and environmental protection, as well as the promotion of alternative plasticizers, the usage of phthalates is expected to gradually decrease in the future.

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

Ethical approval was not required for the studies involving humans because NHANES data is publicly available and is not required to obtain ethical approval. The studies were conducted in accordance with the local legislation and institutional requirements. The human samples used in this study were acquired from The NHANES data sample is derived from the U.S. non-institutionalised (civilian) population, and a complex multi-stage probability sampling design ensures that the sample is nationally representative. Written informed consent to participate in this study was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and the institutional requirements. The animal study was approved by the Beijing Institute of Basic Medical Sciences. The study was conducted in accordance with the local legislation and institutional requirements.

Author contributions

XY: Methodology, Visualization, Conceptualization, Software, Validation, Writing – original draft, Writing – review and editing, Data curation. XW: Methodology, Validation, Conceptualization, Software, Visualization, Writing – original draft. SB: Methodology, Software, Writing – original draft, Visualization, Validation. ZL: Supervision, Software, Writing – review and editing, Visualization, Validation, Data curation. SX: Resources, Funding acquisition, Project administration, Supervision, Methodology, Writing – review and editing. ZD: Visualization, Software, Writing – original draft, Formal Analysis, Resources, Data curation, Funding acquisition, Conceptualization, Methodology, Writing – review and editing, Validation, Supervision.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. This work was funded by the Sichuan Provincial Department of Science and Technology (24NSFSC1606) for Zhangzhen Du and the Chengdu Science and Technology Bureau Key Research and Development Project (2023ZYFS0235) for Shuxing Xing.

Conflict of interest

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

Generative AI statement

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

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/ftox.2025.1617663/full#supplementary-material

References

Akhmetshina, A., Kratky, D., and Rendina-Ruedy, E. (2023). Influence of cholesterol on the regulation of osteoblast function. Metabolites 13 (4), 578. doi:10.3390/metabo13040578

PubMed Abstract | CrossRef Full Text | Google Scholar

Allen, O., Knight, M. M., and Verbruggen, S. W. (2024). Air pollution and osteoporosis. Curr. Osteoporos. Rep. 22 (6), 590–598. doi:10.1007/s11914-024-00889-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Arhzaouy, K., Papadopoulos, C., Schulze, N., Pittman, S. K., Meyer, H., and Weihl, C. C. (2019). VCP maintains lysosomal homeostasis and TFEB activity in differentiated skeletal muscle. Autophagy 15 (6), 1082–1099. doi:10.1080/15548627.2019.1569933

PubMed Abstract | CrossRef Full Text | Google Scholar

Basso, C. G., de Araújo-Ramos, A. T., and Martino-Andrade, A. J. (2022). Exposure to phthalates and female reproductive health: a literature review. Reprod. Toxicol. 109, 61–79. doi:10.1016/j.reprotox.2022.02.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Chandra, A., and Rajawat, J. (2021). Skeletal aging and osteoporosis: mechanisms and therapeutics. Int. J. Mol. Sci. 22 (7), 3553. doi:10.3390/ijms22073553

PubMed Abstract | CrossRef Full Text | Google Scholar

Chung, P. Y., and Van Hul, W. (2012). Paget's disease of bone: evidence for complex pathogenetic interactions. Semin. Arthritis Rheum. 41 (5), 619–641. doi:10.1016/j.semarthrit.2011.07.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Columbres, R. C. A., Chin, Y., Pratti, S., Quinn, C., Gonzalez-Cuyar, L. F., Weiss, M., et al. (2023). Novel variants in the VCP gene causing multisystem proteinopathy 1. Genes. (Basel) 14 (3), 676. doi:10.3390/genes14030676

PubMed Abstract | CrossRef Full Text | Google Scholar

Deng, X., Wu, X., Sun, Z., Liu, Q., and Yuan, G. (2024). Associations between new obesity indices and abnormal bone density in type 2 diabetes mellitus patients. Osteoporos. Int. 35 (10), 1807–1815. doi:10.1007/s00198-024-07163-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Drobny, A., Prieto, H. S., Dobert, J., Kluge, A., Bunk, J., Schlothauer, T., et al. (2022). The role of lysosomal cathepsins in neurodegeneration: mechanistic insights, diagnostic potential and therapeutic approaches. Biochim. Biophys. Acta Mol. Cell. Res. 1869 (7), 119243. doi:10.1016/j.bbamcr.2022.119243

PubMed Abstract | CrossRef Full Text | Google Scholar

Du, J., Zhou, T., Zhang, W., and Peng, W. (2024). Developing the new diagnostic model by integrating bioinformatics and machine learning for osteoarthritis. J. Orthop. Surg. Res. 19 (1), 832. doi:10.1186/s13018-024-05340-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, R., Xue, W., Liang, J., Li, X., Zheng, J., Wang, L., et al. (2024). SOAT1 regulates cholesterol metabolism to induce EMT in hepatocellular carcinoma. Cell. Death Dis. 15 (5), 325. doi:10.1038/s41419-024-06711-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, H., Yang, X., Pan, P., Liu, X., Ma, Y., Chen, Y., et al. (2025). Pubertal low dose exposure to benzophenone-3 (BP-3) alters murine mammary stem cell functions. Ecotoxicol. Environ. Saf. 292, 117982. doi:10.1016/j.ecoenv.2025.117982

PubMed Abstract | CrossRef Full Text | Google Scholar

Guimarães, G. C., Coelho, J. B. C., Silva, J. G. O., de Sant'Ana, A. C. C., de Sá, C. A. C., Moreno, J. M., et al. (2024). Obesity, diabetes and risk of bone fragility: how BMAT behavior is affected by metabolic disturbances and its influence on bone health. Osteoporos. Int. 35 (4), 575–588. doi:10.1007/s00198-023-06991-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, X., Zhao, Z., Shen, H., Qi, B., Li, N., and Hu, J. (2019). VCP myopathy: a family with unusual clinical manifestations. Muscle Nerve 59 (3), 365–369. doi:10.1002/mus.26389

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, T., Meng, X., Liu, X., Wang, J., Yan, S., Zhang, X., et al. (2023). Associations of phthalates with prostate cancer among the US population. Reprod. Toxicol. 116, 108337. doi:10.1016/j.reprotox.2023.108337

PubMed Abstract | CrossRef Full Text | Google Scholar

Hlisníková, H., Petrovičová, I., Kolena, B., Šidlovská, M., and Sirotkin, A. (2020). Effects and mechanisms of phthalates' action on reproductive processes and reproductive health: a literature review. Int. J. Environ. Res. Public Health 17 (18), 6811. doi:10.3390/ijerph17186811

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, H. L., Wu, C. K., Wu, D. J., Liu, W. H., Lee, Y. S., and Wu, C. L. (2023a). Apoptosis pathways and osteoporosis: an approach to genomic analysis. J. Gene Med. 25 (11), e3555. doi:10.1002/jgm.3555

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, Y., Zhou, X., Li, X., Huang, D., Fang, Z., and Ding, R. (2023b). A pan-cancer analysis identifies SOAT1 as an immunological and prognostic biomarker. Oncol. Res. 31 (2), 193–205. doi:10.32604/or.2023.027112

PubMed Abstract | CrossRef Full Text | Google Scholar

Huo, S., Tang, X., Chen, W., Gan, D., Guo, H., Yao, Q., et al. (2024). Epigenetic regulations of cellular senescence in osteoporosis. Ageing Res. Rev. 99, 102235. doi:10.1016/j.arr.2024.102235

PubMed Abstract | CrossRef Full Text | Google Scholar

Hussain, S. M., Ebeling, P. R., Barker, A. L., Beilin, L. J., Tonkin, A. M., and McNeil, J. J. (2023). Association of plasma high-density lipoprotein cholesterol level with risk of fractures in healthy older adults. JAMA Cardiol. 8 (3), 268–272. doi:10.1001/jamacardio.2022.5124

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, J., Chen, Y., Su, Y., Zhang, L., Qian, H., Song, X., et al. (2024). Identification and experimental validation of diagnostic and prognostic genes CX3CR1, PID1 and PTGDS in sepsis and ARDS using bulk and single-cell transcriptomic analysis and machine learning. Front. Immunol. 15, 1480542. doi:10.3389/fimmu.2024.1480542

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, H., Oh, B., and Park-Min, K.-H. (2021). Regulation of osteoclast differentiation and activity by lipid metabolism. Cells 10 (1), 89. doi:10.3390/cells10010089

PubMed Abstract | CrossRef Full Text | Google Scholar

Kruse, C. (2018). The new possibilities from “Big Data” to overlooked associations between diabetes, biochemical parameters, glucose control, and osteoporosis. Curr. Osteoporos. Rep. 16 (3), 320–324. doi:10.1007/s11914-018-0445-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, Y. C., Wuputra, K., Kato, K., Ku, C. C., Saito, S., Noguchi, M., et al. (2024). Di-n-butyl phthalate promotes the neural differentiation of mouse embryonic stem cells through neurogenic differentiation 1. Environ. Pollut. 347, 123722. doi:10.1016/j.envpol.2024.123722

PubMed Abstract | CrossRef Full Text | Google Scholar

Lyu, F., Wang, L., Jia, Y., Wang, Y., Qi, H., Dai, Z., et al. (2023). Analysis of zinc and stromal immunity in disuse osteoporosis: mendelian randomization and transcriptomic analysis. Orthop. Surg. 15 (11), 2947–2959. doi:10.1111/os.13840

PubMed Abstract | CrossRef Full Text | Google Scholar

Ma, B., Wang, X., Xu, C., Xu, Z., Zhang, F., and Cheng, W. (2025). Identification of MEG3 and MAPK3 as potential therapeutic targets for osteoarthritis through multiomics integration and machine learning. Sci. Rep. 15 (1), 23240. doi:10.1038/s41598-025-06175-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Mariana, M., Castelo-Branco, M., Soares, A. M., and Cairrao, E. (2023). Phthalates' exposure leads to an increasing concern on cardiovascular health. J. Hazard Mater 457, 131680. doi:10.1016/j.jhazmat.2023.131680

PubMed Abstract | CrossRef Full Text | Google Scholar

McNamara, L. M. (2021). Osteocytes and Estrogen deficiency. Curr. Osteoporos. Rep. 19 (6), 592–603. doi:10.1007/s11914-021-00702-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Mondal, T., Mondal, S., Ghosh, S. K., Pal, P., Soren, T., Pandey, S., et al. (2022). Phthalates - a family of plasticizers, their health risks, phytotoxic effects, and microbial bioaugmentation approaches. Environ. Res. 214 (Pt 3), 114059. doi:10.1016/j.envres.2022.114059

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Y., Liu, J., Zhao, K., Gao, L., and Zhao, J. (2021). Cholesterol-induced toxicity: an integrated view of the role of cholesterol in multiple diseases. Cell. Metab. 33 (10), 1911–1925. doi:10.1016/j.cmet.2021.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, S., Qi, G., Chen, H., He, D., Ma, D., Bie, Y., et al. (2023). Ferroptosis sensitization in glioma: exploring the regulatory mechanism of SOAT1 and its therapeutic implications. Cell. Death Dis. 14 (11), 754. doi:10.1038/s41419-023-06282-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, C., He, W., Wang, L., Hao, T., Yang, X., Feng, W., et al. (2024). Studies on the role of MAP4K2, SPI1, and CTSD in osteoporosis. Cell. Biochem. Biophys. 83, 2115–2126. doi:10.1007/s12013-024-01621-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Tuo, Y., Lu, X., Tao, F., Tukhvatshin, M., Xiang, F., Wang, X., et al. (2024). The potential mechanisms of catechins in tea for Anti-Hypertension: an integration of network pharmacology, molecular docking, and molecular dynamics simulation. Foods 13 (17), 2685. doi:10.3390/foods13172685

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Y., Zhu, H., and Kannan, K. (2019). A review of biomonitoring of phthalate exposures. Toxics 7 (2), 21. doi:10.3390/toxics7020021

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Pei, Z., Hao, T., Ariben, J., Li, S., He, W., et al. (2022a). Prognostic analysis and validation of diagnostic marker genes in patients with osteoporosis. Front. Immunol. 13, 987937. doi:10.3389/fimmu.2022.987937

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X. C., Luo, L. M., Huang, T. S., and Feng, L. F. (2022b). SOAT1 is a new prognostic factor of colorectal cancer. Ir. J. Med. Sci. 191 (4), 1549–1554. doi:10.1007/s11845-021-02746-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, X., Guo, S., Zhou, H., Sun, Y., Gan, J., Zhang, Y., et al. (2023a). Pan-cancer transcriptomic analysis identified six classes of immunosenescence genes revealed molecular links between aging, immune system and cancer. Genes. Immun. 24 (2), 81–91. doi:10.1038/s41435-023-00197-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Zhang, Y., Cao, J., Wang, Y., Anwar, N., Zhang, Z., et al. (2023b). The role of autophagy in bone metabolism and clinical significance. Autophagy 19 (9), 2409–2427. doi:10.1080/15548627.2023.2186112

PubMed Abstract | CrossRef Full Text | Google Scholar

Wei, R., Cao, Y., Wu, H., Liu, X., Jiang, M., Luo, X., et al. (2023). Inhibition of VCP modulates NF-κB signaling pathway to suppress multiple myeloma cell proliferation and osteoclast differentiation. Aging (Albany NY) 15 (16), 8220–8236. doi:10.18632/aging.204965

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, X., Liang, W., Feng, Z., Li, G., Chen, X., and Zhang, J. (2025). Molecular mechanisms of polychlorinated biphenyls in breast cancer: insights from network toxicology and molecular docking approaches. Front. Pharmacol. 16, 1604993. doi:10.3389/fphar.2025.1604993

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshida, G., Kawabata, T., Takamatsu, H., Saita, S., Nakamura, S., Nishikawa, K., et al. (2022). Degradation of the NOTCH intracellular domain by elevated autophagy in osteoblasts promotes osteoblast differentiation and alleviates osteoporosis. Autophagy 18 (10), 2323–2332. doi:10.1080/15548627.2021.2017587

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, X. H., Cao, H. W., Bo, L., Lei, S. F., and Deng, F. Y. (2023). Air pollution, genetic factors and the risk of osteoporosis: a prospective study in the UK biobank. Front. Public Health 11, 1119774. doi:10.3389/fpubh.2023.1119774

PubMed Abstract | CrossRef Full Text | Google Scholar

Yuan, Y. Z., Ye, C., Sun, J. H., Hu, M. Y., Huo, S. J., Zhu, Y. T., et al. (2022). Toxicokinetics of mono-(2-ethylhexyl) phthalate with low-dose exposure applying fluorescence tracing technique. Toxicol. Appl. Pharmacol. 434, 115814. doi:10.1016/j.taap.2021.115814

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Zheng, L., Cheng, D., Lei, C., Li, H., Zhou, J., et al. (2024). Chronic di(2-ethylhexyl) phthalate exposure at environmental-relevant doses induces osteoporosis by disturbing the differentiation of bone marrow mesenchymal stem cells. Sci. Total Environ. 914, 169918. doi:10.1016/j.scitotenv.2024.169918

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, W., Chen, Q., Wang, C., Yao, D., Zhu, L., Pan, Y., et al. (2020). Inhibition of Cathepsin D (CTSD) enhances radiosensitivity of glioblastoma cells by attenuating autophagy. Mol. Carcinog. 59 (6), 651–660. doi:10.1002/mc.23194

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: plasticiser, osteoporosis, NHANES, molecular docking, molecular dynamics simulation

Citation: Yang X, Wang X, Bao S, Li Z, Xing S and Du Z (2025) Exploring plasticisers-osteoporosis links and mechanisms: a cohort and network toxicology study. Front. Toxicol. 7:1617663. doi: 10.3389/ftox.2025.1617663

Received: 24 April 2025; Accepted: 04 August 2025;
Published: 03 September 2025.

Edited by:

Youji Wang, Shanghai Ocean University, China

Reviewed by:

Stefaan Verbruggen, Queen Mary University of London, United Kingdom
Johana Marquez, Corporacion Universitaria Rafael Nunez Facultad de Ciencias de la Salud, Colombia

Copyright © 2025 Yang, Wang, Bao, Li, Xing and Du. 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: Zhangzhen Du, ZHV6aGFuZ3poZW4xOTg5QGNkdXRjbS5lZHUuY24=

Disclaimer: 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.