This study systematically investigates the causal relationship between the plasma metabolome (exposures) and PCa (outcome) using MR. The methods of this study were in compliance with the STROBE-MR checklist (Skrivankova et al., 2021). Three foundational assumptions necessary for MR analysis were satisfied: (1) genetically determined IVs are strongly associated with the plasma metabolome; (2) the genetic IVs are not associated with any confounding factors related to the plasma metabolome and PCa; and (3) the genetic IVs influence PCa only through their effect on the plasma metabolome. This study utilized publicly available datasets which had already received ethical approval and informed consent; thus, no additional ethical declaration or consent was required. An overview of this study is depicted in Figure 1.

To mitigate non-independence arising from the overlap of exposure and outcome samples, datasets for exposure and outcome were sourced from different databases. In selecting exposure and outcome GWAS summary statistics, criteria encompassed the European origin of the population, the larger size of population samples and SNPs, and the time of data publication. Recently, Chen et al. (Chen et al., 2023) reported a highly comprehensive exploration of the genetic influences on the plasma metabolome, with a focus on metabolites implicated in human diseases. In their study, a total of 8,299 unrelated individuals of European ancestries were recruited from the Canadian Longitudinal Study on Aging (CLSA) cohort (Raina et al., 2019), all of whom had genome-wide genotyping and measurements of circulating plasma metabolites. Following rigorous quality control and data curation, 1,091 metabolite levels were quantified in plasma samples by Metabolon, Inc. (https://www.metabolon.com/) using the Ultrahigh Performance Liquid Chromatography-Tandem Mass Spectroscopy (UPLC-MS/MS) platform. Of these, 850 known metabolites were categorized into eight super pathways: lipid, amino acid, nucleotide, xenobiotics, cofactor and vitamins, carbohydrate, peptide, and energy. Furthermore, to elucidate biological processes not apparent from the study of individual plasma metabolites alone, 309 genetic determinants of metabolite ratios were calculated, including substrate to product ratios and metabolite pairs sharing an enzyme or transporter. Notably, the most recent GWAS summary statistics for PCa, with the largest samples to date (122,188 cases; 604,640 controls; 25,146,978 SNPs) (Wang et al., 2023), was selected for the primary MR analysis. To validate the conclusions drawn, GWAS data were obtained from the Prostate Cancer Association Group to Investigate Cancer Associated Alterations in the Genome (PRACTICAL) Consortium, comprising 79,148 cases, 61,106 controls, and 19,716,640 SNPs (Schumacher et al., 2018). An overview of the GWAS summary statistics utilized in this study is available in the GWAS Catalog (https://www.ebi.ac.uk/gwas/) (Table 1).

To identify the IVs for metabolite levels and ratios, several procedures were undertaken to verify the validity of the first assumption. Initially, genetic variants were identified through association thresholds of p < 1 × 10⁻⁵, a standard practice in MR analysis to encompass a broader range of variation when SNPs available for exposure are limited (Liu et al., 2024). Subsequently, SNPs showing independent inheritance and minimal linkage disequilibrium (LD) (r² < 0.001, kb = 10,000) within the 1,000 Genomes European subset (Genomes Project et al., 2015) were identified through a clumping procedure in R software. SNPs that were ambiguous or palindromic were also excluded. To mitigate the risk of weak instrument bias, the F statistic (F > 10) was employed to ascertain if the identified IVs possessed sufficient power to represent metabolites or metabolite ratios (Burgess et al., 2011). We also utilized the NHGRI-EBI GWAS Catalog (www.ebi.ac.uk/gwas), which adheres to the findable, accessible, interoperable, and reusable (FAIR) data principles, to determine whether selected SNPs (or those in LD with them) were previously associated with other traits or diseases (P < 1 × 10⁻⁵) (Sollis et al., 2023). Using the LDtrait tool, SNPs associated with potential confounders, such as type 2 diabetes mellitus, non-PCa cancer, smoking, alcohol use, stroke, autoimmune diseases, coronary artery disease, and hyperlipidemia, were excluded to minimize confounding between significantly associated metabolites and PCa. To augment the reliability of our findings, an online power calculator (mRnd) was utilized to calculate the statistical power for MR (https://cnsgenomics.shinyapps.io/mRnd/).

In this study, five complementary MR methods were utilized to estimate the causal relationship between exposure and outcome, namely MR-Egger (Bowden et al., 2015), weighted median (Bowden et al., 2016), inverse variance weighted (IVW), simple mode, and weighted mode. Among these, the IVW method, which provides the most accurate causal estimation between exposure and outcome, served as the primary approach (Burgess et al., 2013). Furthermore, several supplementary analyses were performed to validate the robustness of our MR results. The MR-Egger method not only identifies violations of the IVs assumption but also calculates intercepts to assess potential horizontal pleiotropy and bias from ineffective IVs, providing effect estimates unaffected by these violations (Bowden et al., 2015). Cochran’s Q statistic was applied to estimate heterogeneity among the variables (SNPs) (Cohen et al., 2015). Moreover, the MR-PRESSO global test, an additional MR method designed to identify and correct horizontal pleiotropic outliers, was utilized to assess the existence of horizontal pleiotropy (Verbanck et al., 2018). A “leave-one-out (LOO)” sensitivity analysis was conducted to determine whether results were affected by a single SNP (Burgess and Thompson, 2017). Finally, scatter plots and funnel plots were used to visually display the relationships and interplay between each genetic instrument.

Notably, all analyses utilized the “TwoSampleMR” (version 0.4.22) and MR-PRESSO packages in R software (version 3.6.0). A P-value of <0.05 was deemed statistically significant. Moreover, in light of multiple testing concerns, the false discovery rate (FDR) and Bonferroni correction were applied to identify statistically significant results in multiple comparisons. According to prior studies, FDR < 0.05 and Bonferroni correction are regarded as indicative of significant causal relationships, with Bonferroni correction being notably more stringent.

To investigate whether the researched outcomes had an impact on the plasma metabolome, a reverse MR analysis was performed. In this reverse analysis, SNPs (P < 5 × 10⁻⁸, r² < 0.001, kb = 10,000) selected from GWAS data on PCa as IVs, with blood plasma metabolite levels as the outcome, were used to assess the bidirectionality of the previously determined causal relationship.

Metabolic pathway analysis was conducted using the online tool Metaconflict 6.0 (https://www.metaboanalyst.ca/) (Chong et al., 2018). Functional enrichment analyses and the pathway analyses module were utilized to identify metabolite groups or pathways potentially relevant to the biological processes underlying PCa. Importantly, this study focused solely on metabolites that surpassed the recommended association threshold as determined by IVW (P _IVW-mre < 0.05).

Gene expression data were obtained from the publicly available microarray dataset GSE46602, hosted in the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). This dataset consists of 14 benign prostate gland samples and 36 PCa samples (Mortensen et al., 2015). Subsequently, Gene Set Enrichment Analysis (GSEA) was conducted using the R package “clusterProfiler” to identify significantly enriched biological pathways. Given the focus of our study, we specifically examined gene sets related to glutathione metabolism. The analysis was performed using default parameters, and the resulting p-values were adjusted using the FDR method to control for multiple hypothesis testing.

After a rigorous selection process, MR analysis was conducted to evaluate the causality of 1,091 metabolite levels and 309 metabolite ratios on PCa using independent GWAS summary datasets. The IVs generated for the plasma metabolome ranged from 11 to 93 SNPs, with a minimum F statistic of 19.50, indicating that all IVs were sufficiently robust for MR analysis. Supplementary Tables S1, S2 displays the characteristics of SNPs and their genetic associations with the 1,400 plasma metabolome and PCa.

The IVW method of MR was utilized to establish causality between the 1,400 plasma metabolome components and PCa, using independent GWAS summary data. A total of 152 suggestive and unique causal associations (P _IVW-mre < 0.05) were identified, comprising 117 plasma metabolite levels and 35 metabolite ratios, as detailed in Supplementary Tables S3, S4. Of the 117 plasma metabolites, 91 were identified as known metabolites (including 40 lipids, 33 amino acids, 7 xenobiotics, 4 nucleotides, 4 cofactor and vitamins, 2 carbohydrates, and 1 peptide), while 26 were unknown or partially characterized, as classified by the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa and Goto, 2000). Importantly, all overlapping components of the plasma metabolome exhibited consistent causal effects on PCa across various GWAS dataset pairs, as shown in Figure 2. Additionally, the outcomes of heterogeneity analysis employing Cochran’s Q statistic (applied to both IVW and MR-Egger methods) and pleiotropy analysis using MR-Egger intercepts and the MR-PRESSO global test for these suggestive causal associations are detailed in Supplementary Tables S5, S6.

To confirm the robustness of our findings, further screening was performed on the 152 suggestive associations. Associations within the plasma metabolome consistently identified across all five complementary MR methods (IVW, simple mode, weighted mode, weighted median, and MR-Egger) were included. Figure 3 demonstrates that 20 causal associations were consistent across two distinct GWAS datasets of PCa. All three causative associations among the overlapped metabolite ratios were negatively correlated with PCa risk, specifically: arginine to glutamate ratio, phosphate to uridine ratio, and glycerol to mannitol/sorbitol ratio, with corresponding odds ratios (OR) and confidence intervals detailed in the text. Additionally, no evidence of pleiotropy or heterogeneity was observed among these associations. However, after applying FDR correction to these P-values, no significant confirmative associations were observed, categorizing the results as suggestive causal associations (Figure 3).

Among the 17 overlapping metabolite levels, 14 causative associations passed all sensitivity analyses without evidence of heterogeneity or horizontal pleiotropy, with three showing significant association with PCa risk (FDR < 0.05). Significant associations were observed with cysteinylglycine disulfide levels negatively correlated with PCa risk, and uridine and N-acetyl-L-glutamine (NAG) levels positively correlated, as detailed with corresponding OR and confidence intervals in Figure 3. Bonferroni correction (P < 0.05/1,400; 3.571 × 10⁻⁵) was utilized to identify the most significant causal relationships among the 1,400 plasma metabolome components and PCa. Results indicated that NAG levels are significantly associated with an increased risk of PCa. The statistical power for the association of NAG levels with PCa was 1.0, demonstrating robustness against a type I error rate of 0.05. Fixed-effect IVW estimates confirmed the association of NAG levels with an increased risk of PCa across two datasets. Specifically, additional methods including the weighted median, weighted mode, simple mode, and MR-Egger, indicated consistent results, as shown in Figure 4 and Supplementary Figure S1, S2. Furthermore, based on the “leave-one-out” analysis, none of the associations were driven solely by a single SNP, as evidenced in Supplementary Figures S3–S6. The MR Steiger directionality test results confirmed the accuracy of our causal direction estimate (P < 0.001). Reverse MR analysis showed no causal relationship from PCa to NAG levels, indicating directionality in the observed association, as detailed in Supplementary Table S7. Consequently, it was determined that NAG levels might be causally associated with PCa, underscoring the reliability of the results.

Metabolic pathway analysis identified five pathways significantly associated with the development of PCa (P < 0.05): Glutathione metabolism, Biosynthesis of unsaturated fatty acids, Arginine and proline metabolism, Tyrosine metabolism, and Pyrimidine metabolism. Notably, the Glutathione metabolism pathway emerged as a recurrent feature across two independent PCa datasets, underscoring its potential role in PCa pathophysiology (Table 2). To further delineate the molecular underpinnings of this finding, we performed GSEA focused on glutathione-related gene sets. This approach enabled a more integrated assessment of dysregulated metabolic signatures in PCa, transcending the metabolite-level observations. The GSEA results revealed a consistent downregulation of multiple glutathione metabolism-associated gene sets in malignant prostate tissues compared to benign controls. Among these, GOMF_Glutathione_Transferase_Activity and WP_Glutathione_Metabolism showed statistically significant enrichment patterns (FDR q < 0.05) (Supplementary Figure S7).