Our prospective plan was to carry out a variety of two-sample univariable MR (UVMR) analyzes to examine the associations of each of the UKB biochemical biomarkers with overall, ER-positive, and ER-negative breast cancer liability. After our UVMR analyzes showed significant associations, we performed further multivariable MR (MVMR) analyzes to adjust for known risk factors and bidirectional analyses. We finally ranked our nominally significant biomarkers by importance using a multivariable Bayesian approach . Our analysis follows the guidelines for performing MR investigations  and our reporting follows the guidelines for strengthening the reporting of Mendelian randomization studies (STROBE-MR) (Additional file 2: Checklist S1) . We did not pre-register the study protocol.
Our study used summary-level exposure data from the UKB study  and summary-level outcome data from the Breast Cancer Association Consortium (BCAC) . The BCAC includes ~6000 samples from the UK , which amounts to, at most, a ~1.4% sample overlap between the exposure and outcome samples. Our data only includes women of European descent to reduce bias from population stratification.
We obtained publicly available summary-level genome-wide association study (GWAS) statistics on 34 serum, urine, and red blood cell biomarker levels; body mass index (BMI); and alcohol intake frequency from unrelated female participants of white-British ancestry (n = 194,174) in the UKB cohort study from Neale et al. . The genotypes and 34 biomarker levels were collected by the UKB study at baseline between 2006 and 2010 using various laboratory techniques and instruments by different suppliers [7, 10]. The GWASes were performed using age, age^2, and the first 20 principal components (PCs) as covariates . Inverse-rank normalized GWAS data was used because many of the quantitative biomarker traits were non-normally distributed. Most women (at least 59%) in the UKB cohort were post-menopausal . More information about the panel of UKB biomarkers and the original UKB study can be found elsewhere [3, 7].
Publicly available GWAS summary statistics on overall breast cancer cases (n = 122.977) and controls (n = 105,974) of European ancestry were obtained from the BCAC . Of the breast cancer cases, 69,501 were ER-positive, 21,468 were ER-negative, and the majority developed post-menopause. More details about the original studies are described elsewhere [8, 14, 15].
Selection of instrumental variables
For each exposure, we selected associated single-nucleotide polymorphisms (SNPs) at genome-wide significance (P < 5 × 10−8) and ensured their independence by removing those in linkage disequilibrium using the PLINK method (r
2 < 0.001, clumping distance = 10,000kb). We then harmonized the directions of the effect alleles between exposures and outcomes.
In all our MR analyses, SNPs must satisfy three assumptions to be considered valid IVs. Genetic variants must (1) strongly associate with the exposure (the relevance assumption), (2) be independent of confounders (the independence assumption), and (3) affect the outcome only through their effect on the exposure (the exclusion restriction assumption) .
The main univariable analysis consisted of inverse-variance weighted (IVW) MR between each exposure and each outcome. The IVW method first estimates the Wald ratio for each SNP by dividing the SNP-outcome association by the SNP-exposure association and then combines these ratios in a fixed effect meta-analysis where each ratio is weighted by the inverse of the variance of the SNP -outcome association . We used P < 0.05 as the nominal significance threshold. We also derived false discovery rate (FDR)-corrected P-values with the Benjamini-Hochberg (BH) method and used P < 0.05 as the FDR-corrected significance threshold. For exposures for which only 1 IV could be identified, we estimated the Wald ratio . Our results are reported as odds ratios (OR) per standard deviation (SD) change in the genetically predicted biomarker concentration.
A common violation of the exclusion restriction IV condition is caused by horizontal pleiotropy, where a genetic variant has an effect on the outcome that does not occur through the exposure . Therefore, we employed several additional univariable approaches with different underlying assumptions about the structure of the pleiotropy for all exposures, including the MR-Egger weighted median and MR Pleiotropy RESidual Sum and Outlier (MR-PRESSO) . The MR-Egger allows for some directional pleiotropy in its estimate of the causal effect by making the additional Instrument Strength Independent of Direct Effect (InSIDE) assumption, which states that across all instruments, the magnitude of the pleiotropic effect is independent of the strength of the genetic variant-exposure association . The weighted median allows for sparse or balanced pleiotropy by down-weighting outliers . The MR-PRESSO method allows for some directional pleiotropy by identifying and adjusting for outliers .
We tested the robustness of our univariable findings by performing MVMR [22, 23] and bidirectional MRI. MVMR was used to adjust for previously reported risk factors, while bidirectional MR was employed to rule out potential reverse causation.
We performed two-sample MVMR analyzes for all seven biomarkers that were nominally significantly associated with overall breast cancer in IVW MR. We searched for associations that P < 10−8 of all variants used as IVs in Phenoscanner [24, 25] (Additional file 3: T1-T7), a database providing summarized GWASes, and adjusted for traits that could be considered reasons for horizontal pleiotropy. MVMR assumes that pleiotropic pathways operate through the risk factors included in the model . For all MVMR analyses, we included SNPs that were genome-wide significantly associated (P < 5 × 10−8) with any exposure or risk factor that was taken into consideration in an MVMR model and not in linkage disequilibrium (r
2 < 0.001, clumping distance = 10,000kb).
As lipids are correlated we included HDL cholesterol, low-density lipoprotein (LDL) cholesterol, triglycerides, and lipoprotein A in MVMR models to observe the direct associations of each lipid with each outcome.
As BMI  and alcohol intake  are associated with breast cancer risk, we included BMI and alcohol intake frequency in MVMR models for each of the seven biomarkers that we found to be nominally significantly associated with overall breast cancer in IVW MR.
As oestrogen decreases alkaline phosphatase (ALP) expression and activity in breast cancer cells  and we could not obtain enough genetic variants for oestradiol, we adjusted for testosterone and SHBG in an MVMR model with ALP.
After adjusting for BMI in MVMR, significant associations between SHBG and breast cancer risk have been found so  we included BMI and SHBG in MVMR models.
Due to the low prior probability of association between ALP and breast cancer, we performed a bidirectional univariable MR analysis of genetically predicted overall, ER-positive, and ER-negative breast cancer liability and ALP levels.
We used MR Bayesian model averaging (MR-BMA) to agnostically rank the causal importance of the seven biomarkers found to be nominally significantly associated with overall breast cancer in IVW MR while accounting for potential pleiotropy . Empirical P-values were calculated using a permutation approach  and adjusted for multiple testing using the BH method with P < 0.05 as the significance threshold. All independent (r
2 < 0.001) genetic variants associated with any of the biomarkers at genome-wide significance were included in the analysis (n = 460).
We used MR-BMA to consider each combination of biomarkers (all single biomarkers, all pairs of biomarkers, all triplets, and so on) as a candidate model in an MVMR analysis using weighted regression. Each candidate model was assigned a posterior probability (PP) that expresses the likelihood that the candidate model contains the true set of causal biomarkers using the regression’s goodness-of-fit measure.
Then, we used MR-BMA to perform model-averaging to assign each biomarker a marginal inclusion probability (MIP) and report each biomarker’s model-averaged causal effect (MACE) on overall breast cancer. The MIP represents the probability that the biomarker is a causal determinant of breast cancer risk, and the MACE represents the biomarker’s weighted average direct causal effect on risk across all candidate models. The MIP was calculated by summing up the posterior probabilities of all candidate models where the biomarker is present. The MACE underestimates the true causal effect of a biomarker on overall breast cancer and should not be interpreted in absolute terms, but as an indication of the effect direction and to compare the relative causal effects among biomarkers.
We used 0.5 as the prior probability for inclusion in the main analysis, which reflected an a priori belief that half of the candidate models or that half of the nominally significant biomarkers are causal, and priors of 0.25 and 0.75 as sensitivity analyses.
We employed the TwoSampleMR MendelianRandomization MRPRESSO and ieugwasr  R packages, as well as the GitHub repository https://github.com/verena-zuber/ for MR-BMA for our analyzes using R (version 4.0.5). We searched for secondary trait associations using Phenoscanner [24, 25].