This research employed a cross-sectional research design. From August to October 2024, community-dwelling adults were recruited from Xinglinyuan Community in Haikou City, Hainan Province, China, through a combination of convenience sampling and volunteer recruitment.

Inclusion criteria were as follows: (1) aged ≥ 50 years; (2) residence in the target community for more than 6 months; (3) voluntary participation with a clear understanding of the study objectives; (4) ability to communicate effectively and independently complete all required assessments. Exclusion criteria included: (1) a history of lumbar spine surgery that could interfere with accurate BMD measurement; (2) use of antibiotics within the previous 3 months; (3) severe physical disabilities or terminal illnesses that limited participation; (4) serious mental disorders preventing completion of the questionnaire; (5) refusal to provide informed consent or withdrawal during data collection. The study protocol was approved by the Ethics Committee of Hainan Medical University (Approval No. HYLL-2023-098). All participants provided written informed consent prior to enrollment, and the study was conducted in accordance with the principles of the Declaration of Helsinki. The study finally obtained 39 older adult participants with 19 and 20 in the OP group and control group, respectively.

Fasting venous blood samples were collected from all participants in the morning under standardized conditions. Approximately 5 mL of whole blood was drawn into vacuum tubes without anticoagulants and centrifuged at 3000 rpm for 10 min to separate the serum.

Serum levels of S-equol were measured using a commercial ELISA kit (Shanghai Enzyme-Linked Biotechnology Co., Ltd., China). Inflammatory cytokines, including interleukin-6 (IL-6), IL-1β, tumor necrosis factor-α (TNF-α), monocyte chemoattractant protein-1 (MCP-1), and high-sensitivity C-reactive protein (hs-CRP), as well as malondialdehyde (MDA) and total antioxidant status (TAS), were quantified using ELISA kits purchased from Nanjing Boyan Biotechnology Co., Ltd., China, in accordance with the manufacturer’s protocols. All assays were performed in duplicate, and the mean values were used for subsequent statistical analyses. According to the manufacturer’s instructions, the intra-assay and inter-assay coefficients of variation (CVs) for the serum S-equol assay were <8 and <10%, respectively. The oxidative stress index (OSI) was calculated as the ratio of MDA to TAS (25).

BMD of the lumbar spine (L1–L4 vertebrae) was measured for all participants using dual-energy X-ray absorptiometry (DXA) scanner (GE Lunar Prodigy, USA). Prior to scanning, participants were instructed to remove any metal-containing items that might interfere with image quality, and they were positioned in a standard supine posture on the DXA examination table to align the lumbar spine with the scanner’s detection range. Based on the World Health Organization (WHO) diagnostic criteria, participants were classified into two groups according to their T-scores: the OP group with T-scores ≤ − 2.5 (26), and the control group with T-scores > − 2.5. The measurements were performed at the Department of Clinical Laboratory, the First Affiliated Hospital of Hainan Medical University, and all scans were performed on the same device to minimize inter-instrument variability and ensure data consistency.

Fresh fecal samples were collected from each participant using sterile collection tubes and immediately transported to the laboratory in a dry ice box, where they were stored at – 80 °C until further processing. Total microbial DNA was extracted from frozen fecal samples using the QIAamp PowerFecal Pro DNA Kit (QIAGEN, China), following the manufacturer’s instructions. DNA concentration and purity were assessed using a NanoDrop One spectrophotometer (Thermo Fisher Scientific, USA).

Total DNA extracted from fecal samples was used for whole-metagenome shotgun sequencing. Sequencing libraries were prepared following standard protocols and sequenced on the Illumina NovaSeq XPlus platform (Illumina Inc., USA), generating 150 bp paired-end reads. Raw sequencing data were subjected to quality control using FastQC and Trimmomatic, with removal of low-quality reads and host DNA contamination. Clean reads were then assembled using MEGAHIT and annotated using the MetaPhlAn3 pipeline for taxonomic profiling at the species level. Functional annotation was performed based on the Kyoto Encyclopedia of Genes and Genomes (KEGG) and EggNOG databases using HUMAnN3. Relative abundances of microbial taxa and functional pathways were calculated and compared between the control and OP groups. Alpha-diversity and beta-diversity metrics were computed using QIIME2. Differentially abundant species and functions were identified using LEfSe and STAMP, with a Linear Discriminant Analysis (LDA) score threshold set at 2.0 and p < 0.05 considered statistically significant.

Untargeted metabolomic profiling of fecal samples was conducted by LC–MS/MS using an Ultra-High-Performance Liquid Chromatography (UHPLC) system (Thermo Dionex Ultimate 3,000) coupled with a Q Exactive HF-X Orbitrap Mass Spectrometer (Thermo Fisher Scientific, USA). Metabolite extraction was performed using a methanol-based protocol, followed by centrifugation and supernatant collection. Samples were analyzed in both positive and negative electrospray ionization (ESI) modes.

Raw MS data were processed using Compound Discoverer (v3.1, Thermo, USA) for peak detection, alignment, and normalization. Metabolite identification was based on accurate mass-to-charge ratio (m/z), retention time, and MS/MS fragmentation patterns by searching against HMDB, KEGG, and mzCloud databases. Multivariate statistical analyses including principal component analysis (PCA) and orthogonal partial least squares discriminant analysis (OPLS-DA) were performed using SIMCA (v14.1) to visualize group separation and detect differential metabolites.

Differential metabolites between OP and control groups were identified based on VIP > 1, fold change ≥ 1.2, and adjusted p < 0.05. KEGG pathway enrichment analysis of differential metabolites was conducted to explore their potential biological relevance, particularly focusing on metabolic pathways associated with bone metabolism, inflammation, and microbial-host interactions.

To explore the relationship between the gut microbiota and fecal metabolites, multi-omics correlation analysis was conducted using the Spearman correlation method. First, the relative abundance of differential gut microbial species obtained from metagenomic sequencing and the normalized peak intensities of differential fecal metabolites identified through untargeted metabolomics were selected. Spearman’s rank correlation coefficients were calculated between microbial species and metabolites across all participants.

Significant correlations were defined based on both correlation coefficient thresholds (|rho| > 0.5) and statistical significance (p < 0.05). The resulting correlation matrix was visualized as a heatmap using the “pheatmap” package in R (version 4.2). Furthermore, microbial-metabolite interaction networks were constructed and visualized using Cytoscape software (version 3.9.1) to better elucidate the potential biological interactions between dominant gut taxa and their associated metabolites.

Descriptive statistics were applied to summarize demographic variables and clinical characteristics. Continuous variables were tested for normality using the Shapiro–Wilk test. Normally distributed data were presented as means ± standard deviations (SD), whereas non-normally distributed data were presented as medians with interquartile ranges (IQR). Between-group comparisons were performed using independent sample t-tests or Mann–Whitney U tests, as appropriate. Categorical variables were reported as counts and percentages and compared using Fisher’s exact tests. All statistical analyses were performed using SPSS version 20.0 (IBM Corp., Armonk, NY, USA), and R software version 4.2. p < 0.05 was considered statistically significant unless otherwise specified.

The demographic characteristics and baseline clinical parameters of the OP and control groups are summarized in Table 1. No significant differences were observed between the two groups in terms of age, sex, smoking status, alcohol consumption, or other sociodemographic variables (all p > 0.05), showing that the groups were well matched and comparable.

Of particular interest, S-equol was significantly lower levels in the OP group compared to controls (3,561 ± 304 vs. 3,855 ± 469 pg/mL, p = 0.026). In addition, IL-1β levels were significantly higher in the OP group than in the control group (45.0 ± 16.5 vs. 35.8 ± 10.3, p = 0.043), and indicated higher inflammation in individuals with OP.

Additionally, participants with OP showed significantly lower lumbar spine BMD values across all vertebrae (L1–L4) (all p < 0.001). While the average T-score in the OP group was markedly lower (−2.88 ± 0.31) than that in the control group (−1.35 ± 0.08, p < 0.001), consistent with the diagnostic classification of osteoporotic bone loss.

To comprehensively characterize the gut microbial composition in all participants, whole-metagenome shotgun sequencing was performed on fecal samples. As shown in Figure 1A, there was no statistically significant differences in the Chao1 index between the OP and control groups, although the control group exhibited a slightly higher median value. In addition, Shannon and Simpson indices were also calculated to assess α-diversity, and no significant differences were observed between groups (Supplementary Figure S1). These findings suggest that the overall richness and evenness of the gut microbial communities were comparable between osteoporotic and non-osteoporotic individuals.

To assess inter-individual variability in gut microbial composition, beta diversity was analyzed using principal coordinates analysis (PCoA) based on Bray-Curtis distance metrics. The PCoA plot (Figure 1B) showed a partial separation trend between the OP and control groups along the primary coordinate axis (PCoA1, explaining 21.88% of the variance), although the difference was not statistically significant (p = 0.122). The 95% confidence ellipses illustrate considerable overlap between the two groups, showing that overall community structure was broadly similar, with some degree of compositional divergence. These findings imply that while α-diversity remained comparable, subtle differences in microbial community composition may exist between osteoporotic and non-osteoporotic individuals, potentially contributing to functional disparities in gut metabolism or host-microbiota interactions.

The Firmicutes-to-Bacteroidetes (F/B) ratio, a commonly used indicator of gut microbiota stability and metabolic state, has been associated with host energy metabolism, inflammation, and immune regulation. Therefore, we evaluated the F/B ratio as shown in Figure 1C, but no statistically significant difference was observed between the OP and control groups. To further characterize the gut microbial profiles of the study participants, genus-level taxonomic composition was assessed and compared between the OP and control groups. As shown in Figure 1D, both groups mainly comprised genera from the phyla Bacteroidota and Bacillota, with notable differences in the relative abundance of several key taxa. Streptococcus was relatively enriched in the OP group, while Ruminococcus exhibited slightly higher abundance in the control group. Other SCFA-producing genera such as Faecalibacterium, Roseburia, and Blautia were also more abundant in the control group.

To identify key microbial taxa associated with OP, we employed both LEfSe analysis and species-level differential abundance testing. LEfSe analysis (Figure 2A) showed significant enrichment of taxa such as Klebsiella pneumoniae, Gammaproteobacteria, and Akkermansia in the OP group, while Faecalibacterium prausnitzii, Roseburia, and Bacteroides uniformis were more abundant in the control group. In parallel, we conducted unpaired t-tests at the species level to compare relative abundances between groups. As shown in Figure 2B, the top 20 species with the most significant differences displayed group-specific enrichment trends. Among these, species including Phocaeicola salanitronis, Butyricicoccus intestinisimiae, and Bacteroides sp. ET71 exhibited a higher abundance in the control group. Notably, Akkermansia sp. Marseille-P6666 and Oscillibacter sp. CU971, known to participate in lipid metabolism, anti-inflammatory functions and gut barrier integrity, were present at significantly lower relative abundances in individuals with OP compared with controls.

To further explore the ecological interactions among the differential species, we constructed group-specific microbial co-occurrence networks based on Spearman correlations (Figure 2C: Control; Figure 2D: OP). In these networks, red-colored nodes indicate species with higher degree centrality, such as Clostridium sp. CAG433 and Paenibacillus oryzisoli, showing that these taxa frequently co-occur with multiple other species and may serve as ecological hub species within the gut microbiome. Comparison of the two network structures revealed clear differences between the OP and control groups. Notably, species such as Neisseria sp. HMSC065C04 and Providencia sp. PROV120 appear in both networks; however, their connection strengths and interaction partners vary between the groups. Although these species are present across conditions, their roles within the microbial ecosystem may differ, possibly contributing distinctively to the development or modulation of OP.

To elucidate the functional implications of gut microbiota alterations in OP, we first performed microbial metabolic pathway analysis based on predicted KEGG orthologs (Figure 3A). The OP group exhibited significant upregulation in pathways related to ABC transporters, glutathione metabolism, and antimicrobial resistance, including bacterial secretion systems and biofilm formation. In contrast, the control group exhibited enrichment in core biosynthetic pathways such as ribosome biogenesis, RNA polymerase, and methanogenesis, showing distinct microbial functional associated with host bone health status.

Subsequently, a spearman correlation heatmap (Figure 3B) showed associations between host characteristics (age, sex, BMI), the phytoestrogen metabolite S-equol, serum biomarkers (e.g., IL-6, IL-1β, TNF-α, hs-CRP, MCP-1, MDA, TAS), and gut microbial species. Notably, S-equol exhibited positive correlations with several microbial taxa, including Paenibacillus oryzisoli, Thiopseudomonas, and Chloroflexi bacterium. Conversely, negative correlations were observed between S-equol and taxa such as Citrobacter sp. Lgbk14 and Klebsiella sp. HSTU Sny5. Additionally, inflammatory markers and oxidative stress indicators showed distinct associations with specific bacterial species, indicating associations between microbial composition and host redox and inflammatory status. These correlations suggest potential microbial contributors to both metabolic phenotypes and S-equol bioavailability in the context of OP.

Based on these findings, we applied a Random Forest algorithm to identify microbial features most strongly associated with OP status (Figure 3C). The top-ranked species included Phocaeicola salanitronis, Cytobacillus gottheilii, Prevotella sp. CAG255, and Pediococcus acidilactici, which contributed notably to the model’s classification accuracy. Of these, several species showed clear enrichment in the OP group, while others were more abundant in controls, reflecting disease-specific microbial signatures. To evaluate the diagnostic utility of the top discriminatory taxa, receiver operating characteristic (ROC) analysis was performed (Figure 3D). Cytobacillus gottheilii achieved the highest area under the curve (AUC = 0.8575), followed by Pediococcus acidilactici (AUC = 0.8077) and Prevotella sp. CAG255 (AUC = 0.7817), showing their potential as non-invasive biomarkers for identifying OP.

To assess the global variation in fecal metabolite composition between osteoporotic individuals and controls, a partial least squares discriminant analysis (PLS-DA) was conducted. As shown in Figure 4A, the PLS-DA plot demonstrated a clear separation between the OP and control groups along the first two components (PC1 and PC2), which explained 4.34% and 7.48% of the total variance respectively, and it indicated that OP was accompanied by notable alterations in gut microbial metabolism, warranting further identification of differential fecal metabolites and associated biochemical pathways.

To further identify specific fecal metabolites that differed significantly between osteoporotic individuals and controls, differential analysis was conducted and visualized using a volcano plot. As shown in Figure 4B, a total of 83 metabolites showed significant differences between the OP and control groups. 43 upregulated metabolites (red dots) were predominantly enriched in the OP group, while 40 downregulated metabolites (green dots) were more abundant in control individuals. The majority of metabolites (gray dots, n = 1,119) did not exhibit statistically significant differences between the two groups. The identified differentially abundant metabolites provide a foundation for subsequent pathway enrichment and correlation analyses to explore their potential involvement in bone metabolism or S-equol production. To further visualize the expression patterns of the most significantly different metabolites, a hierarchical clustering heatmap was constructed based on the top-ranked differential features (Figure 4C). Each column represents an individual sample, and each row corresponds to a significantly dysregulated metabolite. The heatmap revealed distinct clustering patterns between the OP and control groups, supporting the separation observed in the PLS-DA and volcano plots. Several metabolites, including D-tryptophan, 3-indoleacrylic acid, decanedioic acid, and guanidinoacetate, exhibited markedly higher intensities in the control individuals. In contrast, several amino acid metabolites such as Gln-Val-Ile-Asp, Leu-Ser-Asn, acetamide, and avenanthramide D, were enriched in the OP individuals. The clear metabolic differentiation and specific fecal metabolites may be served as potential biomarkers or mechanistic indicators of gut metabolic dysbiosis in OP.

To explore potential functional implications, differential fecal metabolites between osteoporotic and control groups were mapped to the KEGG database. The KEGG bar plot categorizes the enriched pathways into six major functional classes, as showed in Figure 4D. Among these, metabolism-related pathways were the most dominant, particularly those involved in lipid and amino acid metabolism, underscoring their relevance to membrane lipid remodeling and nitrogen metabolism, which could contribute to osteoporotic pathophysiology. Complementarily, the KEGG scatter plot (Figure 4E) offers an overview of enrichment significance and the number of differential metabolites involved. Pathways with the highest enrichment scores and lowest p-values included glycerophospholipid metabolism, glycine/serine/threonine metabolism, and retrograde endocannabinoid signaling, showing a functional link to host-microbiota interactions and skeletal health. Additional pathways such as bile acid biosynthesis, gonadotropin-releasing hormone (GnRH) signaling, cAMP signaling, and sphingolipid signaling were also enriched, implicating alterations in hormonal regulation and immune signaling cascades.

To further assess the group-specific differences in fecal metabolites, boxplot analysis was performed on 15 representative discriminatory metabolites between the OP and control groups (Figure 5A). Notably, several metabolites, such as D-tryptophan, 3-indoleacrylic acid, and 2-hydroxyindoleacetic acid, were present at lower levels in the OP group and were associated with differences in microbial tryptophan and amino acid metabolic profiles. To further identify key metabolites that best distinguish osteoporotic individuals from controls, a random forest classification model was applied. As shown in Figure 5A, the top 15 discriminatory metabolites were ranked according to their importance, measured by the mean decrease in classification accuracy. Among these, 5-oxooctanoic acid, atenolol, and Gln-Val-Ile-Asp, were identified as the most informative features contributing to the model’s predictive performance. Color-coding of relative abundance patterns indicated clear group-specific differences, with several metabolites markedly enriched in the OP group (red), while others showed higher levels in controls (blue). These metabolites are involved in diverse biological processes including lipid metabolism, bile acid conjugation, tryptophan degradation, and amino acid turnover, pathways that may be relevant to gut microbial activity and bone metabolism. To evaluate the diagnostic potential of key metabolites identified through random forest analysis, ROC curve analyses were conducted using the top three features ranked by MeanDecreaseGini. As shown in Figure 5C, atenolol exhibited the highest discriminative ability with an AUC of 0.7904, followed closely by Gln-Val-Ile-Asp (AUC = 0.7787) and 5-oxooctanoic acid (AUC = 0.6774).

To explore the potential clinical relevance of fecal metabolic alterations, we performed a correlation analysis between significantly different fecal metabolites and various host parameters, including serum S-equol, inflammatory cytokines, and oxidative stress markers (Figure 5D). Notably, S-equol showed negative correlations with Asn-Trp and acetamide, and a strong positive correlation with lysylproline, showing a potential involvement in modulating amino acid metabolism. Furthermore, inflammatory and oxidative stress markers exhibited significant correlations with some distinct metabolites, highlighting the complex interplay among host immune responses, oxidative stress, and microbial metabolic outputs.

To investigate the potential interplay between different gut microbial taxa and fecal metabolites, we performed Spearman correlation analyses and constructed a metabolite-microbe interaction network. As shown in Figures 6A,B multiple gut microbial species exhibited significant correlations with distinct metabolites including amino acid derivatives, plant secondary metabolites, and phospholipids. These findings highlight intricate connections between specific microbial taxa and host fecal metabolites, and may contribute to host metabolic phenotypes relevant to bone metabolism and warrant further mechanistic validation.