A total of 43 patients who met the diagnostic criteria for OP were enrolled in this study and randomly divided into two groups between June 2020 and December 2021 at the Department of Orthopedics, The Third Affiliated Hospital of Zhejiang Chinese Medical University. The study was approved by the ethics committee of The Third Affiliated Hospital of Zhejiang Chinese Medical University (approval No. ZSLL-KY-2021-017-01). Moreover, the study has been submitted for registration in the Chinese clinical trial registry using the following registration number: ChiCTR2200056265. All participants provided written informed consent before inclusion in the study.

The participants met the diagnostic criteria for a diagnosis of primary OP provided by the guidelines of the National Osteoporosis Foundation (Cosman et al., 2014). We included patients with a vertebral bone mineral density T-score ≤ -2.5 measured by dual-energy X-ray absorptiometry (DXA). Exclusion criteria were as follows: 1) patients with secondary OP; 2) patients with comorbid severe chronic functional disease; 3) patients who had been treated with bisphosphonates, denosumab, and calcitonin within 3 months; and 4) patients with digestive system or liver or kidney function diseases.

A total of 43 patients were included in the study, and their ages ranged from 71-87 years. Basic information such as the age, height, and weight of all participants was collected before treatment, and there was no significant difference in general characteristics between the two groups (P >0.05). Both groups were given 1 tablet of α-Calcitol (0.5 μg/capsule) orally per day for 3 months, and the treatment group was additionally given YGD orally for 3 months. Bone mineral density (BMD) measures were collected before starting treatment and 3 months after treatment, and fecal samples were collected at 3 months after treatment. All patients had fecal samples taken while in the fasted state and at similar time points in the morning.

Fecal samples were collected in sterile plastic cups and stored at -80°C before further processing. Total genomic DNA from samples was extracted using the SDS method. DNA concentration and purity were evaluated on 1% agarose gels. According to the concentration, DNA was diluted to 1 ng/µL using sterile water. Genomic DNA extraction and PCR amplification were performed using fecal samples, and PCR products were pooled and purified. Finally, sequencing libraries were generated using the TruSeq^® DNA PCR-Free Sample Preparation Kit (Illumina, USA) following the manufacturer’s recommendations, and index codes were added. The library quality was assessed on a Qubit@ 2.0 Fluorometer (Thermo Scientific) and Agilent Bioanalyzer 2100 system. Finally, the library was sequenced on an Illumina NovaSeq platform, and 250 bp paired-end reads were generated. The sequencing data were processed to obtain the effective tags (ET), and the uparse algorithm was used to cluster all the effective tags obtained from all samples, followed by operational taxonomic unit (OTU) clustering and species annotation; finally, alpha diversity and beta diversity were calculated, and functional annotation was performed.

Samples stored in a -80°C freezer were thawed on ice. A 400 μL solution (methanol: water = 7:3, V/V) containing an internal standard was added to a 20 mg sample, and the mixture was vortexed for 3 min. Each sample was sonicated in an ice bath for 10 min, vortexed for 1 min, and then placed at -20°C for 30 min. Each sample was then centrifuged at 12000 rpm for 10 min (4°C). The sediment was removed, and the supernatant was centrifuged at 12000 rpm for 3 min (4°C). A 200 μL aliquot of supernatant was obtained for Liquid chromatography–mass spectrometry (LC MS) analysis. All samples were analyzed by the LC MS system following the manufacturer’s instructions. The analytical conditions were as follows: UPLC: column, Waters ACQUITY UPLC HSS T3 C18 (1.8 µm, 2.1 mm*100 mm); column temperature, 40°C; flow rate, 0.4 mL/min; injection volume, 2 μL; solvent system, water (0.1% formic acid): acetonitrile (0.1% formic acid); gradient program, 95:5 V/V at 0 min, 10:90 V/V at 11.0 min, 10:90 V/V at 12.0 min, 95:5 V/V at 12.1 min, 95:5 V/V at 14.0 min.

Unsupervised PCA (principal component analysis) was performed by the statistics function prcomp using R software (www.r-project.org ). The data were unit variance scaled before unsupervised PCA was performed. The HCA (hierarchical cluster analysis) results obtained during the analysis of the samples and metabolites were presented as heatmaps with dendrograms, while Pearson correlation coefficients (PCCs) between the levels of molecules in the samples were calculated by the cor-function in R and presented as only heatmaps. Both HCA and PCC calculation were carried out by the R package Complex Heatmap. Identified metabolites were annotated using the Kyoto Encyclopedia of Genes and Genomes (KEGG) Compound database (http://www.kegg.jp/kegg/compound/ ), and annotated metabolites were then mapped to the KEGG Pathway database (http://www.kegg.jp/kegg/pathway.html ). Significantly enriched pathways were identified using a hypergeometric test p value in the analysis of a specific list of metabolites.

Regarding demographic and clinical parameters, two-tailed t tests were used for paired analyses (before and after intervention) by using the Statistical Package for the Social Sciences (SPSS) version 25.0 (SPSS Inc., Chicago, IL, United States). For the GM data, Wilcoxon’s signed-rank test was used for comparisons of paired samples, and the Mann Whitney U test was used to evaluate independent samples.

A total of 43 individuals were included in the study. They were randomly divided into two groups with 21 in the treatment group and 22 in the control group. The general characteristics of the two groups did not differ significantly and were comparable. After 3 months of treatment, patients in the treatment group showed significantly more improvements in BMD than those in the control group ( Table 1 ).

Paired-end sequencing of the GM samples obtained from both groups of patients was performed using the Illumina NovaSeq sequencing platform. Tiling and quality control were performed on the resulting data to obtain clean tags, which were then filtered for chimerism to yield effective tags that could be used for subsequent analysis. Clustering of OTUs was performed on the effective tags with 97% concordance (identity), followed by species annotation of the OTUs. A total of 2646 OTUs were identified, of which 755 were common to both groups, 1172 were unique to the treatment group and 719 were unique to the control group ( Figure 1A ). Based on the annotation results, the top 10 species ranked by the maximum abundance at the genus level for each sample were selected, and a column-shaped additive plot of the relative abundance of each species and their ratio in each sample at different taxonomic levels was generated ( Figure 1B ). The top 10 ranked species by maximum abundance at the genus level in each group (Escherichia-Shigella, Bacteroides, Faecalibacterium, Klebsiella, Megamonas, Akkermansia, Blautia, Bifidobacterium, Romboutsia, Holdemanella) are shown in Figure 1C . Based on the genus level data on species abundance, the top 35 genera were selected, the abundances of the microbes in each sample were used to cluster the microbes at the species and sample levels, and a heatmap of the species abundance was generated ( Figure 1D ). Of these genera, 15 had significantly more species clustered in the treatment group, and 20 had significantly more species clustered in the control group ( Figure 1E ). To further investigate the phylogenetic relationship of species at the genus level, the representative sequences of the top 100 genera were obtained by multiple sequence alignment ( Figure 1F ). Among them, the top five were Proteobacteria, Bacteroidetes, Firmicutes, Verrucomicrobiota, and Actinobacteria.

The dilution curve can directly reflect the sufficiency of the amount of sequencing data and indirectly reflect the abundance of species in the sample, and when the curve is flat, the amount of sequencing data is reasonable. The diversity of the microbes within the two patient groups was substantial and homogeneous ( Figure 2A ). The α-diversity index Chao1 was used to determine the ecological diversity within the microbial community. Analyses of the intergroup differences in the alpha diversity index indicated that the two groups of patients had significantly different flora species ( Figure 2B ). We performed principal coordinate analysis (PCoA, principal coordinates analysis) based on the weighted and unweighted UniFrac distances and selected the combinations of principal coordinates with the largest contribution for mapping ( Figure 2C ). The more similar species structure within the two patient groups is indicated in the figure. Box plots of the beta diversity data show significant differences in species between the two patient groups ( Figure 2D ). The multidimensionality confirmed the accuracy of the above results. To investigate the similarity among different samples, clustering trees of the samples can be generated by clustering the samples. The unweighted pair-group method with arithmetic mean (Lü et al., 2019) is a more commonly used method for cluster analysis. UPGMA clustering analysis was performed using the weighted and unweighted UniFrac distance matrices, and the clustering results were displayed by integrating them with the species relative abundances of each sample at the genus level (see Figures 2E , F ). Statistical analysis was used to identify the species with significantly different abundances among the subgroups and evaluate the enrichment of the differential species among the different subgroups ( Figure 2G ). The 193 species with significant differences in abundance at the genus level between the two groups are presented in the figure along with the top 47 differential species. Simper (similarity percentage) is a decomposition of the Bray Curtis dissimilarity index that is used to quantify the contribution of each species to the dissimilarity between two groups (Hamidi et al., 2019). The top 10 ranked species with contributions to the differences between the two groups and their abundances are presented ( Figure 2H ). The top ten species that contributed to the differences belonged to the genera Bacteroides, Escherichia-Shigella, Faecalibacterium, Megamonas, Blautia, Klebsiella, Romboutsia, Akkermansia, and Prevotella_9. Using LEfSe (LDA effect size) to compare the statistical significance and biological correlation of species differences between two groups, we identified biomarkers with significant differences between groups and identified features with different abundances and associated categories. The statistical results of LEfSe were visualized using three approaches: a histogram reflecting the distribution of the LDA value ( Figure 2I ), an evolutionary branch diagram (phylogenetic distribution) and figure illustrating the comparison of the abundance of biomarkers with significant differences between groups in different groups ( Figure 2J ).

Network analysis was performed after the data were filtered to yield valid data by performing a correlation index calculation across all samples ( Figure 3A ). Through the species co-occurrence network analysis, we can readily identify the species that dominate the interaction and the species groups that interact closely. These dominant species and species groups may play a unique and important role in maintaining the stability of the microbial community structure and function in the environment. The functional annotation of OTUs obtained from clustering was performed using Tax4Fun. According to the database annotation results, the top 10 functions of each sample or group with the largest abundance at each annotation level were selected, and a columnar stacking chart of the relative abundance of functions was generated to visually display the relative abundance of each microbe at different annotation level functions and their proportions ( Figure 3B ). The top 10 functional information points were analyzed and tested, and the functional information with significant differences was identified ( Figure 3C ). Based on the functional annotation and abundance information obtained from the samples using the database, the top 35 ranked features by abundance and their abundance information in each sample were selected to draw a heatmap and clustered at the level of functional differences ( Figure 3D ). The results showed that there were 4 functional clusters in the treatment group, and the enrichment of 3 functional clusters was significantly decreased.

To further investigate the differences in the sample metabolites between the treatment and control groups, after quantitatively identifying the statistically significant differential metabolites between different groups, we next aimed to investigate the link between the identified metabolites and biological processes or biological states by untargeted metabolomics. Then, the results were analyzed ( Figures 4A, B ) and indicated that the best model interpretation was obtained. Based on the OPLS-DA results, after the different varieties or differential metabolites between tissues were initially identified, the p value obtained through univariate analysis was combined to further identify the differential metabolites and draw a volcano plot of the results ( Figure 4C ). Fifty-five of the differential metabolites had higher levels and 19 had lower levels in the treatment group than in the control group, and 623 metabolites had no significant difference in levels. After qualitative and quantitative analysis of the detected metabolites, the fold changes in metabolite levels between groups were compared in addition to the grouping of specific samples ( Figure 4D ). The top 20 differential metabolites are presented in the figure. KEGG pathway enrichment was performed using the differential metabolites ( Figure 4E ).

To further investigate the connection between the flora and metabolites, we performed a combined analysis of the differences in the patients’ gut microbes and differential metabolites. As shown in Figure 5A , the differential metabolites were first subjected to PCA. PCA of the GM microbes was then performed ( Figure 5B ). Differential microbes and differential metabolites were subjected to correlation analysis, and Spearman correlation coefficients were calculated using the levels of microbes and metabolites; significant correlations were defined as those with a correlation |r| > = 0.8 and a p - value in the correlation coefficient significance test < 0.05. The results of these analyses of the differential microbes and metabolites are shown in Supplementary Table 2 . The relevant microorganisms and metabolites identified above were then plotted into a chord plot ( Figure 5C ). This correlation analysis resulted in the identification of significantly correlated microbes and metabolites at the genus level ( Supplementary Table 3 ). A chord diagram was drawn ( Figure 5D ). Correlation network plots were generated separately to demonstrate the correlations between microbes and metabolites based on two correlation analysis modalities ( Figures 5E , F ).