Utilizing the R Studio software, we conducted an intricate instrumental variable selection process. Initially, we select SNPs that showed a significant association with metformin treatment (P < 5 × 10−8) within the GWAS dataset. Subsequently, we established linkage disequilibrium parameters (r2 = 0.001, kb = 10000), culminating in the selection of 46 distinct SNPs. Next, we extracted data for the above-selected SNPs from the summary statistics of outcome trait (OP, LS-BMD, FN-BMD, and FA-BMD). Allele alignment was performed, and during this harmonization process, SNPs with inconsistent alleles and those with ambiguous palindromic SNPs that could not be corrected were removed. Subsequently, we employed the Phenoscanner database to interrogate phenotypes associated with the remaining SNPs. SNPs (specifically, rs34872471, rs780093, and rs849142) displaying associations with potential confounding factors such as pioglitazone, systemic lupus erythematosus, and high-density lipoprotein cholesterol were expunged, with a significance threshold set at P < 5E − 8. When exploring the causal effect of metformin treatment on forearm bone density, it was determined that SNPs linked to arm muscle mass and hand grip strength (rs780093, rs8756, rs4715207, rs1421085, rs76895963, rs947791, rs10195252, rs1483988, rs459193, and rs12146652) were excluded as possible confounding factors in addition to the three SNPs previously mentioned. The outliers identified by MR-PRESSO have been eliminated. Additionally, we computed the F-statistic for every SNP, and all of them exceeded the noteworthy threshold of 10. This observation attested to the substantial instrumental strength of these SNPs, thereby ensure the robustness of our ensuing MR analysis.

Ultimately, we employed 34, 34, 34, 24, 24, and 34 SNPs, respectively, as instrumental variables for the evaluation of metformin treatment's impact on outcomes, encompassing OP (Supplementary Table 1), LS-BMD (Supplementary Table 2), FN-BMD (Supplementary Table 3), FA-BMD (Supplementary Table 4), eBMD (Supplementary Table 5), and fracture (Supplementary Table 6).

Based on the inverse variance weighting (IVW) analysis outcomes, our findings unveiled a causal association between metformin treatment and a diminished risk of osteoporosis (OR: 0.859, 95% CI: 0.774–0.953, P = 0.004). This observation was corroborated by the weighted median (WM) method (OR: 0.830, 95% CI: 0.724–0.951, P = 0.007), Mendelian randomization-Egger (MR-Egger) analysis (OR: 0.646, 95% CI: 0.462–0.902, P = 0.015), penalized weighted Median (PWM) analysis (OR: 0.828, 95% CI: 0.722–0.950, P = 0.007), and maximum likelihood (ML) analysis (OR: 0.858, 95% CI: 0.781–0.943, P = 0.001). These diverse analytical approaches collectively yielded consistent outcomes (Fig. 1). Supplementary Fig. 1 shows the scatter plot of the above five methods for the association of metformin treatment with risk of OP.

Causal effects for metformin treatment on OP. MR-Egger, weighted median, inverse-variance weighted, penalized weighted median, and maximum likelihood estimates of Mendelian randomization are summarized. CI, confidence interval; nSNP, number of single nucleotide polymorphism; OR, odds ratio

Our comprehensive analysis of heterogeneity indicates the lack of significant heterogeneity in the IVW analysis (Cochran’s Q = 41.53, P = 0.146) and the MR-Egger analysis (Cochran’s Q = 37.89, P = 0.219) (Table 1). Furthermore, our analysis revealed no significant evidence of horizontal pleiotropy, as indicated by both the MR-Egger intercept test (P = 0.09), with P values exceeding the 0.05 threshold (Table 1). MR-PRESSO did not identify any outliers.

Furthermore, our study conducted leave-one-out tests and found that the causal influence of metformin treatment on osteoporosis remained constant even when individual SNPs were excluded (Supplementary Fig. 2). This confluence of results indicates the stability and reliability of our analysis regarding the causal relationship between metformin treatment and osteoporosis.

In our quest to delve deeper into the potential causal impact of metformin treatment on bone density at distinct anatomical sites, we conducted comprehensive MR analyses. The outcomes of our investigation revealed a positive correlation between metformin treatment and LS-BMD (IVW: OR: 1.063, 95% CI: 1.023–1.105, P = 0.002; MR-Egger: OR: 1.197, 95% CI: 1.047–1.367, P = 0.013; ML: OR: 1.066, 95% CI: 1.032–1.102, P = 0.017; Fig. 2), FN-BMD (IVW: OR: 1.034, 95% CI: 1.000–1.069, P = 0.049; MR-Egger: OR: 1.145, 95% CI: 1.020–1.284, P = 0.028; ML: OR: 1.035, 95% CI: 1.006–1.065, P = 0.016; Fig. 2) and eBMD (IVW: OR: 1.035, 95%CI: 1.023–1.047, P ≤ 0.001; MR-Egger: OR: 1.036, 95%CI: 1.010–1.062, P = 0.012; ML: OR: 1.036, 95%CI: 1.028–1.043, P ≤ 0.001; WM: OR: 1.032, 95%CI: 1.021–1.042, P ≤ 0.001; PWM: OR: 1.032, 95%CI: 1.021–1.042, P ≤ 0.001). The WM and PWM results for LS-BMD and FN-BMD analyses aligned with the IVW outcomes, albeit with somewhat diminished significance (LS-BMD: WM: OR: 1.043, 95% CI: 0.994–1.093, P = 0.084; PWM: OR: 1.033, 95% CI: 0.984–1.085, P = 0.186; Fig. 2. FN-BMD: WM: OR: 1.022, 95% CI: 0.980–1.066, P = 0.307; PWM: OR: 1.004, 95% CI: 0.963–1.046, P = 0.853; Fig. 2). However, no causal relationship was observed between metformin treatment and FA-BMD (IVW: OR: 1.050, 95% CI: 0.969–1.138, P = 0.237; Fig. 2). As shown in Fig. 2, WM, MR-Egger, PWM, and ML analyses consistently yielded analogous non-significant results (P < 0.05). It is worth noting that due to the detection of heterogeneity, we employed the random-effects IVW method to assess the causal associations between metformin treatment and LS-BMD, FN-BMD, FA-BMD, and eBMD. Supplementary Fig. 3 shows the scatter plot of the above five methods for the association of metformin treatment with LS-BMD, FN-BMD, FA-BMD, and eBMD.

Causal effects for metformin treatment on LS-BMD, FN-BMD, and FA-BMD and eBMD. MR-Egger, weighted median, inverse-variance weighted, penalized weighted median, and maximum likelihood estimates of Mendelian randomization are summarized. LS-BMD, lumbar spine bone mineral density; FN-BMD, femoral neck bone mineral density; FA-BMD, forearm bone mineral density; eBMD, estimated heel bone mineral density; CI, confidence interval; nSNP, number of single nucleotide polymorphism; OR, odds ratio

In the evaluation of heterogeneity, noteworthy findings emerged from the heterogeneity tests concerning LS-BMD (IVW: Cochran’s Q = 48.81, P = 0.038), FN-BMD (IVW: Cochran’s Q = 47.70, P = 0.047), FA-BMD (MR-Egger: Cochran’s Q = 34.45, P = 0.044), and eBMD (MR-Egger: Cochran’s Q = 62.578, P ≤ 0.001; IVW: Cochran’s Q = 62.593, P ≤ 0.001) (Table 2). However, no significant evidence of horizontal pleiotropy was observed based on the MR-Egger intercept tests (LS-BMD: P = 0.081; FN-BMD: P = 0.082; FA-BMD: P = 0.789; eBMD: P = 0.940) (P > 0.05) (Table 2). The MR-PRESSO analysis identified outliers in eBMD, leading to their exclusion from the dataset (rs10965246, rs11257655, rs11708067, rs1421085, rs2215383, rs459193, rs7177055, rs7615045, rs76895963, rs849142, rs8756, and rs947791). No outliers were detected in LS-BMD, FN-BMD, and FA-BMD. Furthermore, in the leave-one-out tests, no individual SNP was identified as exerting an influential effect on the causal relationship between metformin treatment and bone mineral density(Supplementary Fig. 4).

Our investigation has revealed that the metformin treatment is causally associated with a reduced risk of fractures (IVW: OR: 0.958, 95% CI: 0.928–0.989, P = 0.008; WM: OR: 0.957, 95% CI: 0.9223–0.993, P = 0.021; ML: OR: 0.957, 95% CI: 0.934–0.981, P ≤ 0.001; PWM: OR: 0.957, 95% CI: 0.922–0.994, P = 0.023; Fig. 3). The MR-Egger analyses concurred with the IVW results, albeit the significance was somewhat weakened (OR: 0.920, 95% CI: 0.824–1.029, P = 0.154; Fig. 3). Supplementary Fig. 5 shows the scatter plot of the above five methods for the association of metformin treatment with risk of fracture.

Causal effects for metformin treatment on fracture. MR-Egger, weighted median, inverse-variance weighted, penalized weighted median, and maximum likelihood estimates of Mendelian randomization are summarized. CI, confidence interval; nSNP, number of single nucleotide polymorphism; OR, odds ratio

There is a significant degree of heterogeneity present in both the IVW analysis (Cochran’s Q = 58.378, P = 0.003) and the MR-Egger analysis (Cochran’s Q = 59.358, P = 0.003) (Table 3). Nevertheless, the MR-Egger intercept tests showed no significant evidence of horizontal pleiotropy(P = 0.473) (Table 3). The MR-PRESSO analysis detected outliers and we subsequently excluded them from the dataset (rs1483988 and rs34872471).

Additionally, the leave-one-out tests did not identify any individual SNP that significantly affected the causal relationship between receiving metformin treatment and the risk of fractures (Supplementary Fig. 6).