We searched for published genome-wide association studies (GWASs) evaluating individuals of European ancestry on the GWAS Catalog and PubMed (the last search was performed in September 2023). Summary statistics of Ca and P were obtained from the Neale Laboratory with a sample size of 315,153 UK Biobank participants (16). The median, decile 1, and decile 9 values for Ca were 2.376 mmol/L, 2.268 mmol/L, and 2.497 mmol/L, respectively, whereas for P, the corresponding values were 1.165 mmol/L, 0.954 mmol/L, and 1.363 mmol/L. Summary statistics of 25(OH)D traits with a sample size of 417,580, were obtained from the UK Biobank, and the median, mean, and interquartile range of 47.9 mmol/L, 49.6 mmol/L, and 33.5–63.2 nmol/L, respectively (17). Summary statistics of PTH from the INTERVAL study with a sample size of less than 5,000 were used, and the range of PTH was not publicly available (18).

Summary statistics for kidney function traits, including CKD, CKD_Rapid3, CKDi25, eGFRcrea, eGFRcys (estimated glomerular filtration rate using serum cystatin C), BUN (blood urea nitrogen), UACR (urinary albumin–creatinine ratio), and microalbuminuria (MA), were obtained from several published GWAS meta-analyses (19–22). CKD was defined as eGFRcrea of<60 mL min⁻¹ per 1.73 m². CKD_Rapid3 was characterized by an annual eGFRcrea decline of 3 mL min⁻¹ per 1.73 m² or more and CKDi25 represented a decline of 25% or more in eGFRcrea with a final eGFRcrea below 60 mL min⁻¹ per 1.73 m², among individuals with an initial eGFRcrea exceeding 60 mL min⁻¹ per 1.73 m². MA was defined as UACR > 25 mg/g in women and > 17 mg/g in men. The majority of individual studies included in these meta-analyses were population-based, with less than 10% originating from diabetic populations, representing less than 5% of the total population in each meta-analysis. The GWASs for CKD, CKD_Rapid3, CKDi25, eGFRcrea, eGFRcys, BUN, UACR, and MA, included a subset of participants, representing 0.16%, 5.8%, 7.7%, 0.16%, 0.16%, 0.16%, 0%, and 0% respectively, from the UK Biobank. These subsets of participants overlapped with those included in GWAS analyses for other exposure factors. Detailed information regarding these studies can be found in Supplementary Table S1 or by referring to the original studies.

As all these GWASs were approved by their respective local research ethics committees and institutional review boards, no additional ethical approval was necessary for this study.

LDSC regression regressed SNV GWAS χ2 statistics for one phenotype to infer SNV-based heritability (h2) or χ2 statistics cross products for a pair of phenotypes to infer SNV-based coheritability on LDSC scores (23). To ensure sufficient statistical power, only phenotypes with a heritability z-score > 4, indicating a substantial proportion of phenotype variance explained by heredity, were included in subsequent analyses. The presence of an intercept close to 1, indicating polygenicity rather than population stratification or cryptic relatedness among samples, accounted for the majority of the observed increase in mean χ2 statistics. The LDSC v1.0.1 software (https://github.com/bulik/ldsc) was employed for this purpose.

The first step in univariable MR analysis was to construct effective IVs. As some exposure factors had a limited number of candidate IVs (P-value < 5E−08) or more than 65% of those (P-value < 5E−08) were not available in the GWAS results of the specific outcome, a relatively relaxed significance threshold was set at 5E−06 to improve statistical power. Hence, the first assumption of MR was satisfied. Subsequently, these IVs were grouped based on the 1000 Genomes Project LD structure and independent SNPs (R2 < 0.001 with any other SNP within 10,000 kb) with the most significant P-value. Moreover, to satisfy the latter two assumptions as much as possible, IVs associated with known confounders, such as hypertension, diabetes mellitus, lipid traits, body mass index, and gout, available at http://www.phenoscanner.medschl.cam.ac.uk/, were excluded. Following the removal of palindromic SNVs with intermediate allele frequencies, the IVs related to kidney traits (P-value [kidney traits] < P-value [exposure factors]), and weak IVs with an F value (beta2/se2) of <10, the initially selected IVs were constructed.

The next step was to perform an MR analysis with the initially selected IVs. In the absence of heterogeneity and pleiotropy, the IVW (inverse variance weighted) method combined information on all uncorrelated IVs into a single causal estimate with larger weights given to more precise IVs based on their inverse variances. In particular, the IVW radial as the main approach improved the performance of the IVW estimate and Cochran’s Q statistics by using modified second-order weights (28). Heterogeneity was checked with the Cochran Q test and pleiotropy was tested with MR-Egger. The median-based and median-based approaches were conducted as a supplement given significant heterogeneity, whereas the Robust Adjusted Profile Score (RAPS) (29), which addressed idiosyncratic pleiotropy by robustifying the adjusted profile score, and MR-Egger were conducted as a supplement given significant pleiotropy. Causal directionality was evaluated with the Steiger test (30). The symmetrical funnel plot was generated to identify bias. Additionally, sensitivity analyses were conducted, such as the leave-one-out analysis with the IVW method to check whether the overall estimate was driven by a single SNP. In instances of substantial heterogeneity, the Mendelian Randomization Pleiotropy Residual Sum and Outlier (MR-PRESSO) method was employed for sensitivity analysis. This approach assessed heterogeneity by iteratively excluding each variant to identify horizontal pleiotropic outliers using the Residual Sum of Squares (RSS) (global test). Moreover, the method addressed horizontal pleiotropy by eliminating outliers and examined for significant differences in causal estimates before and after outlier removal (distortion test). When the global test results were significant, the corrected results of MR-PRESSO were considered more reliable than the raw results of that (31).

To further mitigate the risk of false positives, global heterogeneity, and reverse causality, pleiotropic SNPs identified from the pleiotropic analyses, outliers detected by IVW radial (28), and IVs potentially associated with kidney traits (with a P-value [kidney traits] < 0.05 if the IV threshold was 5E−08) were excluded to establish the final set of IVs. Subsequently, a new MR analysis was conducted using these final IVs, and the results were compared with those obtained using the initially selected IVs.

The multivariable MR analysis included sets of exposure factors that exhibited genetic correlations among them. IVs were constructed by using multiple sets of overlapping SNPs from the GWASs for relevant exposure factors, which met the univariable MR inclusion criteria described previously. In the absence of heterogeneity and pleiotropy, the multivariable IVW method was the primary choice. If pleiotropy was detected, the multivariable Egger method was used, whereas the multivariable median method was recommended for addressing heterogeneity. Similar to univariable MR, we conducted multivariable MR analyses before and after removing pleiotropic SNPs.

The MR analyses mentioned were conducted using R software v4.1.0 and the following R packages: “TwoSampleMR” v0.5.5, “RadialMR” v2.2.1, and “MendelianRandomization” v0.4.4. The significance threshold for MR analyses mentioned was set at P-value = 0.05/12 = 4.2E−03 with Bonferroni correction applied due to the multiple testing (3 exposures × 4 outcomes = 12 tests).

Phenotypes with z-scores < 4, including PTH, CKD_Rapid3, CKDi25, eGFRcys, and MA, were excluded from subsequent analyses ( Table 1 ). Genetic correlations were observed between serum 25(OH)D and eGFRcrea (rg −0.094; P-value 1.4E−05), as well as BUN (rg 0.127; P-value 1.7E−06). Serum Ca showed a significant genetic correlation with UACR (rg 0.202, P-value 5.0E−04). However, there was no statistically significant genetic correlation between serum P and CKD (P-value 0.03). Among the exposure factors, only serum 25(OH)D was correlated with serum P (rg 0.076; P-value 0.002) ( Table 2 ).

The characteristics of the annotated genes and related LD-independent loci are described in Supplementary Tables S2 and S3 . Figure 2 illustrates the proportion of horizontally pleiotropic genes for each exposure among all the genes associated with that particular exposure, as well as the percentage of horizontally pleiotropic LD-independent loci for each exposure among all the loci associated with that exposure. The potential pleiotropic SNPs located within the pleiotropic genes or loci are presented in Supplementary Table S4 .

Due to the unavailability of more than 65% of candidate IVs with a P-value less than 5E−08 in the GWAS results of UACR, a relatively relaxed significance threshold of 5E−06 was adopted when investigating UACR as the outcome.

The characteristics of the initially selected effective IVs after removing confounding IVs ( Supplementary Table S4 ), and the final set of effective IVs chosen after additional elimination of pleiotropic IVs ( Supplementary Table S5 ), outliers ( Supplementary Table S6 ), and IVs associated with the outcome, were detailed in Supplementary Tables S5–S30 . Over 85% of the variants overlapped between outliers and IVs associated with the outcome.

There was no evidence of a violation of the causal directionality assumption in any of the univariable MR analyses ( Supplementary Tables S31, S32 ). Initially selected IVs exhibited widespread heterogeneity but no evidence of horizontal pleiotropy. Upon using the finally selected IVs, the heterogeneity disappeared. However, horizontal pleiotropy was only evident for the relationship between 25(OH)D and BUN (P-value 0.04) ( Supplementary Tables S33–S36 ). Although the global test of MR-PRESSO identified horizontal pleiotropic outliers only with initially selected IVs, except for the relationships involving Ca or 25(OH)D and UACR, the distortion test did not reveal significant differences in causal estimates before and after the removal of outliers ( Supplementary Tables S37, S38 ). All symmetrical funnel plots indicated no disproportionate effects from selected IVs ( Supplementary Figures S1–S4 ).

None of the MR analysis methods found any exposure factors associated with higher odds of CKD after the Bonferroni correction (P-value>4.2E-03) ( Supplementary Tables S39, 40 and Supplementary Figures S5, S6 ).

In the analysis conducted with the initially selected IVs ( Supplementary Table S41 ), the primary method, weighted median, supported both the associations between P and eGFRcrea (beta −0.01; P-value 2.0E−04), as well as between Ca and eGFRcrea (beta 0.009; P-value 8.0E−04). This finding was further confirmed by the main sensitivity analysis, MR-PRESSO ( Figure 3 ) ( Supplementary Table S37 ), rather than the leave-one-out analysis ( Supplementary Figure S7 ). In the analysis utilizing the finally selected IVs ( Supplementary Table S42 ), the IVW radial method consistently affirmed the association between P and eGFRcrea (beta −0.005; P-value 2.0E−03), supported by sensitivity analyses with MR-PRESSO ( Figure 3 ) and leave-one-out analysis ( Supplementary Figure S8 ).

Causal estimates for BUN ( Supplementary Tables S43, 44 ) indicated that with the initially selected IVs, only higher serum P was linked to a higher BUN, as evidenced by the weighted median (beta 0.02; P-value 3.0E−03), as well as sensitivity analyses using the MR-PRESSO ( Figure 3 ) and the leave-one-out analysis ( Supplementary Figure S9 ). Furthermore, in the analysis with the final IV selection, the primary analysis method, IVW radial approach, also validated the causal relationship between P and BUN (beta 0.02; P-value 7.5E−07), supported by sensitivity analyses using MR-PRESSO ( Figure 3 ) and leave-one-out analysis ( Supplementary Figure S10 ).

Regarding UACR ( Supplementary Tables S45, S46 ), in the analysis using either the initially selected or the finally selected IVs, it was found that only genetically predicted higher Ca was related to the increase in UACR, as confirmed by the IVW radial method (beta 0.11, P-value 2.0E−03; beta 0.13, P-value 2.0E−04) and sensitivity analyses using the MR-PRESSO ( Figure 3 ) and the leave-one-out analysis ( Supplementary Figures S11, S12 ).

The concurrent utilization of serum P and 25(OH)D as exposure variables in a multivariable MR analysis was driven by the identified genetic correlation between the two. The characteristics of the selected SNPs after excluding confounding IVs ( Supplementary Table S47 ) and those finally selected after further removing pleiotropic IVs ( Supplementary Table S48 ) are detailed in Supplementary Tables S49–S54 . Both the initial and final IV sets showed extensive heterogeneity but no evidence of horizontal pleiotropy ( Supplementary Tables S55, S56 ). The multivariable median method considered as the main method only showed that the associations between serum P and eGFRcrea or BUN were still significant (P-value < 4.2E-03) ( Figure 4 ), regardless of the exclusion of pleiotropic IVs, consistent with the results in univariable MR analysis.