This study utilized three publicly available datasets from the Gene Expression Omnibus (GEO) database¹ to investigate molecular signatures associated with bone mineral density (BMD) and postmenopausal osteoporosis (PMO). GSE56815 includes circulating monocyte samples from 80 individuals, of which 40 postmenopausal women (20 with high BMD and 20 with low BMD) were analyzed to focus on the population most affected by PMO. GSE62402 comprises peripheral blood monocyte samples from 10 female participants (5 with high BMD and 5 with low BMD), although menopausal status is not specified, representing a limitation in its applicability to PMO-specific analyses. GSE230665 features femur tissue samples collected during hip arthroplasty from 12 postmenopausal osteoporosis patients and 3 healthy postmenopausal controls, providing direct insights into bone-specific molecular changes. These datasets, derived from Homo sapiens, were selected based on their relevance to BMD phenotypes, biological sample types pertinent to bone remodeling, and stratified data availability. Detailed sample characteristics are provided in Supplementary Table S1.

To investigate molecular differences related to osteoporosis, differential expression analysis was conducted on selected subsets from the GEO datasets using the limma package (6) in R software. A subset of 40 postmenopausal women from GSE56815, categorized into high and low BMD groups, was analyzed to identify differentially expressed genes (DEGs). Similarly, the GSE62402 dataset, containing peripheral blood monocyte samples from individuals in high and low BMD groups, was examined. The thresholds for DEG selection were set at p < 0.05 and |logFC| > 0.2 to balance the detection of biologically meaningful changes and the inclusion of a broader gene set for exploratory analyses. This choice reflects the exploratory nature of the study, aiming to capture subtle yet potentially significant changes in gene expression, which are particularly relevant in chronic metabolic diseases like osteoporosis. Using a less stringent threshold (p < 0.05) without adjustment for multiple comparisons increases the risk of false positives. However, applying adjusted p-values significantly reduced the number of DEGs, potentially excluding biologically relevant genes. To enhance the robustness of the findings, we performed an overlap analysis to identify genes consistently upregulated or downregulated across GSE56815 and GSE62402, highlighting key molecular signatures associated with BMD variation. These intersecting DEGs were further validated using the independent GSE230665 dataset, focusing on femur tissue samples, to ensure consistency across sample types and datasets.

To investigate the biological pathways associated with osteoporosis-related gene expression changes, we conducted Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis using Gene Set Enrichment Analysis (GSEA). This approach allows for the evaluation of predefined gene sets across the entire ranked list of genes, avoiding the need for arbitrary thresholds in defining differentially expressed genes (DEGs). Separate analyses were performed for genes upregulated and downregulated in the GSE56815 and GSE62402 datasets. The enrichment analysis was implemented using the gseKEGG function from the clusterProfiler package (7) in R, which evaluates pathway significance through permutation testing and adjusts p-values using the Benjamini-Hochberg (BH) method to control the false discovery rate (FDR). Pathways with FDR-adjusted p-values < 0.05 were considered significantly enriched. This approach ensures robust identification of key pathways underlying osteoporosis pathophysiology while minimizing biases associated with DEG selection.

The Comparative Toxicogenomics Database (CTD²) was leveraged to obtain comprehensive data on genes and environmental pollutants associated with the progression of postmenopausal osteoporosis (PMO). The CTD integrates manually curated information on chemical-gene-disease interactions, providing a valuable resource for exploring molecular mechanisms influenced by environmental exposures (8). From the CTD, we selected genes with an inference score > 25 to identify key targets and chemicals potentially associated with PMO. By intersecting these CTD-derived genes with the DEGs identified in our study, we pinpointed a subset of overlapping genes termed key environment-responsive genes (KERGs). These KERGs represent critical molecular targets influenced by environmental pollutants, providing a robust framework for understanding their role in PMO pathogenesis. To further explore the environmental factors influencing these KERGs, we searched for associated chemicals in CTD. We applied an additional selection criterion of reference count >10 to ensure the reliability of chemical-gene interactions. Only chemicals meeting this threshold were considered for network construction. Subsequently, we developed a disease-gene-chemical interaction network to visualize the relationships between PMO, identified genes, and environmental pollutants.

To further explore the interaction between environmental pollutants and the identified KERGs, molecular docking and molecular dynamics simulations were performed. The molecular docking was conducted using AutoDock Vina (version 1.5.6) to predict the binding affinity between selected chemicals and the KERGs. The three-dimensional (3D) structures of the chemicals were obtained from the PubChem database,³ while the protein structure of KERGs was retrieved from the Protein Data Bank (PDB)⁴ (9). Docking results were assessed based on binding energy, with values below −5.0 kcal/mol considered indicative of favorable interactions. Molecular dynamics (MD) simulations were carried out for the most promising chemical-KERG complexes using GROMACS software (version 2020.6-MODIFIED) with the Amber99sb-ildn force field. The simulations employed a TIP3P water model to accurately represent the aqueous environment, and periodic boundary conditions were applied to mimic physiological conditions. The prepared system underwent energy minimization to eliminate unfavorable contacts and optimize the initial structure. This was followed by equilibration in two stages: first under an NVT ensemble to stabilize the temperature, and subsequently under an NPT ensemble to ensure pressure stability. A 100 ns production run was then performed to observe the stability and dynamic behavior of the complexes. To analyze the stability of the chemical-KERG interactions, we analyzed three key parameters: root mean square deviation (RMSD) to monitor overall structural stability throughout the simulation, root mean square fluctuation (RMSF) to assess residue-level flexibility and identify regions with greater structural fluctuations, and radius of gyration (Rg) to determine the compactness of the protein-ligand complex, offering insights into both structural stability and conformational changes. In addition to these stability metrics, we generated a 3D free energy landscape (FEL) to visualize the conformational states and free energy minima of the complexes during the simulation.

Bisphenol A (BPA; CAS: 80–05-7, molecular formula: C₁₅H₁₆O₂, HY-18260, purity: 99.95%) was purchased from Med Chem Express (New Jersey, United States). A 1 mM stock solution was prepared in DMSO and subsequently diluted in culture medium to achieve the desired working concentrations. Bone marrow-derived macrophages (BMMs) were isolated from the tibias and femurs of 4-to 5-week-old C57BL/6 J mice. Bone marrow cells were flushed out and cultured in α-MEM complete medium supplemented with 10% fetal bovine serum (FBS), 1% penicillin–streptomycin, and 30 ng/mL macrophage colony-stimulating factor (M-CSF). Osteoclast differentiation was induced by treating BMMs with α-MEM complete medium containing 50 ng/mL receptor activator of nuclear factor kappa-B ligand (RANKL) and 30 ng/mL M-CSF, with RANKL-containing medium replaced every 2 days. Cytotoxicity was assessed by culturing BMMs in 96-well plates with complete medium containing various concentrations of BPA (0, 5, 10, 15, 20, and 25 μM) for 24 or 48 h. After the treatment period, CCK-8 reagent was added to each well, followed by incubation for 60 min. Absorbance at 450 nm was measured using a microplate reader to evaluate cell viability.

The concentration of BPA (20 μM) used in the intervention group was determined based on the cytotoxicity assay results, which confirmed that this dose was non-cytotoxic to bone marrow-derived macrophages (BMMs) while allowing for the evaluation of its biological effects on osteoclastogenesis. Osteoclast differentiation was analyzed using tartrate-resistant acid phosphatase (TRAcP) staining. Following fixation in 4% paraformaldehyde for 10 min, cells were stained for TRAcP and visualized using the EVOS FL Auto 2 imaging system (Thermo Fisher Scientific, United States). Quantification of TRAcP-positive multinucleated osteoclasts, including their number and size, was performed using ImageJ (v1.8.0) to assess the extent of osteoclastogenesis. To investigate the molecular effects of BPA on AKT1 expression, RNA was extracted from BMMs treated with BPA (20 μM) during the osteoclast differentiation process. Total RNA was isolated using Trizol reagent, and complementary DNA (cDNA) was synthesized according to the manufacturer’s instructions. RT-qPCR was conducted using the LightCycler 480 system (Roche, Switzerland) to measure the expression levels of AKT1, with GAPDH serving as an internal control. Relative expression levels were calculated using the ΔΔCt method. The primer sequences were: AKT1, 5′-TGTATGAGAAGAAGCTGAGCCC-3′ (Forward), 5′-AGTAGGAGAACTTGATCAGGCG-3′ (Reverse); GAPDH, 5′-CTCAGGAGAGTGTTTCCTCGTC-3′ (Forward), 5′-CCTTGACTGTGCCGTTGAATTT-3′ (Reverse).

After washing the cells with PBS, they were lysed in RIPA buffer containing 1% PMSF and 1% protease/phosphatase inhibitors. The lysates were then centrifuged at low temperature to extract total protein. Next, the protein samples were separated on a 10% SDS-PAGE gel at a constant voltage of 120 V for 60 min and transferred onto a PVDF membrane. The membrane was blocked with 5% skim milk for 90 min. Subsequently, the membrane was incubated at 4°C for 14 h with primary antibodies specific for p-AKT (Proteintech Group, 66444-1-Ig), AKT (Proteintech Group, 10176-2-AP), and GAPDH (Proteintech Group, 60004-1-Ig). This was followed by a 60-min incubation at room temperature with a fluorescently labeled secondary antibody (Goat anti-Rabbit IgG (H + L), Thermo Fisher Scientific, SA5-35571). Finally, the membrane was visualized using an infrared fluorescence imaging system (LI-COR Biosciences, United States), and band intensities were quantified with ImageJ (v1.8.0).

The comprehensive data screening strategy is outlined in the flowchart presented in Figure 1. From the GSE56815 dataset, which includes 40 postmenopausal women categorized into high and low BMD groups, a total of 425 DEGs were identified. Among these DEGs, 113 genes were found to be upregulated, while 312 genes were downregulated in the high BMD group compared to the low BMD group (Figure 2A). Similarly, the analysis of the GSE62402 dataset, consisting of peripheral blood monocyte samples from individuals with high and low BMD, identified 300 DEGs, with 29 genes upregulated and 271 genes downregulated in the high BMD group (Figure 2B). A total of four genes—AKT1, ATF1, PRNP, and SAR1A—were consistently downregulated across both datasets, indicating robust and reproducible molecular signatures related to BMD regulation (Figure 2C). Notably, no overlapping genes were identified in the upregulated category between the two datasets (Figure 2D). The expression levels of these four downregulated genes were further validated using the independent GSE230665 dataset, which contains clinical samples from PMO patients and healthy controls. As shown in the Figure 2E, the results confirmed the downregulation of AKT1, ATF1, PRNP, and SAR1A in osteoporosis samples compared to controls. Specifically, AKT1 and ATF1 showed statistically significant reductions in expression (p < 0.001 and p < 0.01, respectively), while PRNP exhibited no significant difference. SAR1A showed a moderate but statistically significant reduction in expression (p < 0.05).

The KEGG pathway enrichment analysis, performed using GSEA, provided insights into key biological processes associated with BMD variation. DEGs were ranked by logFC to generate gene lists, with enrichment direction reflecting the relative activity in the high BMD group compared to the low BMD (control) group. Both datasets revealed several unique pathways, as well as overlapping pathways, that were consistently activated or suppressed in the high BMD group. Specifically, in GSE56815, pathways such as cytokine-cytokine receptor interaction (NES = 1.54, adjusted p-value = 0.01) and alanine, aspartate, and glutamate metabolism (NES = 1.8, adjusted p-value = 0.03) were activated, indicating enhanced immune signaling and amino acid metabolism (Figure 3A). Conversely, oxidative phosphorylation (NES = −1.68, adjusted p-value = 0.02) and apoptosis (NES = −1.61, adjusted p-value = 0.03) were suppressed, reflecting decreased mitochondrial activity and apoptotic signaling (Figure 3B). In GSE62402, pathways including cytoskeleton in muscle cells (NES = 1.92, adjusted p-value < 0.001) and ECM-receptor interaction (NES = 2.13, adjusted p-value < 0.001) showed activation in the high BMD group, highlighting increased cellular structural reorganization (Figure 3C). Immune-related pathways such as Toll-like receptor signaling (NES = −2.04, adjusted p-value < 0.001) and neutrophil extracellular trap formation (NES = −2.18, adjusted p-value < 0.001) were suppressed, indicating a dampened inflammatory response (Figure 3D). Crucially, two pathways were consistently regulated across both datasets. The cytokine−cytokine receptor interaction pathway was activated in the high BMD group in both datasets, highlighting its essential role in immune signaling and bone remodeling. Conversely, the efferocytosis pathway was suppressed across both datasets, indicating a potential reduction in the clearance of apoptotic cells, which may contribute to changes in the bone microenvironment.

From the CTD, we retrieved a total of 16,758 data points linking genes and environmental pollutants to PMO (Supplementary Table S2). Using an inference score > 25 as the threshold, 150 genes associated with environmental chemicals influencing PMO were identified for further analysis. To pinpoint key molecular targets, these 150 genes were intersected with the DEGs identified from the GSE56815 and GSE62402 datasets. This intersection yielded a single overlapping gene, AKT1, which we identified as the KERG (Figures 2C,D). The identification of AKT1 suggests that it may serve as a critical molecular link between environmental pollutants and BMD regulation. To investigate the potential environmental modulators of AKT1, we filtered chemicals from the CTD using a reference count threshold of >10. The interactions were categorized into two primary groups: chemicals that increase AKT1 phosphorylation and those that decrease AKT1 phosphorylation. Several chemicals, such as bisphenol A (BPA) and lipopolysaccharides (LPS), were identified as increasing AKT1 phosphorylation, with BPA having a notable reference count of 24. Conversely, compounds like resveratrol and doxorubicin were linked to decreased AKT1 phosphorylation, suggesting an inhibitory influence on AKT1 activity. Interestingly, two chemicals, hydrogen peroxide and BPA, exhibited dual actions, capable of both increasing and decreasing AKT1 phosphorylation depending on the experimental context or biological conditions (Supplementary Figure S1; Supplementary Table S3). This duality highlights the complexity of AKT1 regulation and the context-dependent nature of environmental chemical interactions.

From the interactions identified between KERGs and environmental chemicals, we selected compounds that were categorized as increasing AKT1 phosphorylation for molecular docking studies. Due to the unavailability of 3D structures in the PubChem database, LPS, cadmium chloride, and particulate matter were excluded from the analysis. As a result, four key compounds—BPA, hydrogen peroxide, doxorubicin, and estradiol—were used for detailed molecular docking. The docking results (Table 1) revealed varying binding energies, with doxorubicin showing the strongest binding affinity to AKT1 (−7.6 kcal/mol) (Figure 4A), followed by estradiol (−6.2 kcal/mol) (Figure 4B), BPA (−5.9 kcal/mol) (Figure 4C), and hydrogen peroxide (−2.4 kcal/mol). Given that hydrogen peroxide exhibited a relatively weak interaction, it was excluded from further consideration for dynamic studies. For molecular dynamics simulations, both BPA and doxorubicin were selected. BPA was included due to its significant environmental relevance as a widespread contaminant and endocrine disruptor, as well as its novel potential impact on AKT1 and bone health, areas with limited prior investigation. Doxorubicin, while primarily a chemotherapeutic agent, was included due to its strong binding affinity to AKT1 and its potential to provide insights into AKT1 modulation mechanisms. The inclusion of doxorubicin adds depth to the study by allowing a comparative analysis between an environmental pollutant and a pharmacological compound, thereby broadening the understanding of AKT1 interactions. Molecular dynamics simulations were conducted to analyze the interactions of BPA and doxorubicin with AKT1, providing insights into their stability and dynamic behavior. For BPA, the root mean square deviation (RMSD) analysis revealed stable binding throughout the simulation, with minimal fluctuations in the AKT1-BPA complex compared to the free AKT1 structure. The root mean square fluctuation (RMSF) showed localized flexibility, particularly at residues interacting with BPA, indicating structural adaptability to accommodate the ligand. The radius of gyration (Rg) remained relatively constant, suggesting compactness and stability of the protein-ligand complex. The Gibbs free energy landscape demonstrated well-defined energy minima, reflecting stable conformational states of the AKT1-BPA complex (Figure 5A). For doxorubicin, the RMSD showed slightly higher fluctuations compared to BPA, indicating dynamic interactions within the AKT1-doxorubicin complex. The RMSF analysis highlighted regions of flexibility at binding residues, which may contribute to the observed higher fluctuations. The Rg values suggested sustained structural integrity of the AKT1-doxorubicin complex over the simulation period. Notably, the Gibbs free energy landscape displayed deeper and more distinct energy minima compared to BPA, indicating stronger binding affinity and conformational stability of the AKT1-doxorubicin complex (Figure 5B).

Concentrations of BPA up to 20 μM maintained acceptable levels of cell viability, while 25 μM significantly reduced viability, indicating potential cytotoxicity at higher doses (Figure 6A). Based on these results, 20 μM was selected as the working concentration for subsequent experiments. Osteoclastogenesis assays revealed that BPA significantly enhanced the formation of TRAcP-positive multinucleated osteoclasts compared to the RANKL-only group. The size of osteoclasts was markedly increased in the BPA-treated group, as visualized by TRAcP staining (Figure 6B) and quantified as a higher percentage of osteoclast area (Figure 6C). At the molecular level, real-time PCR analysis demonstrated that BPA treatment significantly upregulated AKT1 mRNA expression in BMMs during osteoclast differentiation (Figure 6D). Western blot analysis further showed that, following RANKL induction, BPA notably promoted AKT phosphorylation (Figures 6E,F). These findings suggest that BPA promotes osteoclastogenesis and enhances AKT1 expression, indicating its potential role in driving bone resorption and contributing to osteoporosis pathogenesis.