Exploiting the mediating role of the metabolome to unravel transcript-to-phenotype associations

Abstract
Editor's evaluation
Introduction
Results
Discussion
Methods
Data availability
References
Article and author information
Metrics

Abstract

Despite the success of genome-wide association studies (GWASs) in identifying genetic variants associated with complex traits, understanding the mechanisms behind these statistical associations remains challenging. Several methods that integrate methylation, gene expression, and protein quantitative trait loci (QTLs) with GWAS data to determine their causal role in the path from genotype to phenotype have been proposed. Here, we developed and applied a multi-omics Mendelian randomization (MR) framework to study how metabolites mediate the effect of gene expression on complex traits. We identified 216 transcript-metabolite-trait causal triplets involving 26 medically relevant phenotypes. Among these associations, 58% were missed by classical transcriptome-wide MR, which only uses gene expression and GWAS data. This allowed the identification of biologically relevant pathways, such as between ANKH and calcium levels mediated by citrate levels and SLC6A12 and serum creatinine through modulation of the levels of the renal osmolyte betaine. We show that the signals missed by transcriptome-wide MR are found, thanks to the increase in power conferred by integrating multiple omics layer. Simulation analyses show that with larger molecular QTL studies and in case of mediated effects, our multi-omics MR framework outperforms classical MR approaches designed to detect causal relationships between single molecular traits and complex phenotypes.

Editor's evaluation

The reviewers found that your article brings important new methods and insight for how to analyze large, complex, multi-omic datasets in order to highlight specific molecular hypotheses for follow-up validation. As realistic/affordable sample sizes for population studies with omics data have recently exploded in size, this has raised the clear need for new or adapted statistical methods for best exploiting the increased statistical power. This work is compelling and we believe it should be of general interest to biologists and biostatisticians, and of particular interest to those working on (or those who could work on) large cohorts.

https://doi.org/10.7554/eLife.81097.sa0

Introduction

Genome-wide association studies (GWASs) have identified thousands of single nucleotide polymorphisms (SNPs) associated with a wide range of complex traits (MacArthur et al., 2017; Visscher et al., 2017). However, the path from GWAS to biology is not straightforward as most SNPs implicated by GWASs reside in non-coding regions of the genome (MacArthur et al., 2017) and do not directly inform on the functional mechanism through which variants exert their effect on phenotypes.

GWASs have been performed on gene expression (Võsa et al., 2021), DNA methylation (Min et al., 2021), protein (Sun et al., 2018), and metabolite (Shin et al., 2014; Lotta et al., 2021) levels, identifying genetic variants influencing molecular traits, commonly referred to as molecular quantitative trait loci (molQTLs). The large overlap between complex and molecular trait-associated variants suggests that integrating these data can help interpreting GWAS loci (Vandiedonck, 2018; Taylor et al., 2019; Ongen et al., 2017). Advances in the field of transcriptomics make gene expression the best studied molecular phenotype, thanks to the presence of large expression QTL (eQTL) studies (e.g., eQTLGen Consortium [Võsa et al., 2021] N>30,000). Availability of these datasets fostered the development of summary statistic-based statistical approaches aiming at identifying associations between transcripts and complex traits (Zhu et al., 2016; Porcu et al., 2019; Hormozdiari et al., 2014; Gusev et al., 2016), prioritizing genes from known GWAS loci for functional follow-up, and inferring the directionality of these relations (Porcu et al., 2019; Porcu et al., 2021b). However, the cascade of events that mediates the effect of genetic variants on complex traits involves more than one molecular trait. Although approaches used for gene expression can be extended to other molecular data, investigating whether these molecular traits reside along the same causal pathway remains under-explored and only recently have studies applied colocalization and Mendelian randomization (MR) to methylation, gene expression, and protein levels data (Giambartolomei et al., 2018; Wu et al., 2018; Gleason et al., 2020; Sadler et al., 2022) and to a lesser extent to metabolic QTLs (mQTL) (Yin et al., 2022).

Metabolites are often the final products of cellular regulatory processes and the most proximal omic layer to complex phenotypes. Their levels could thus represent the ultimate response of biological systems to genetic and environmental changes. For instance, the metabolic status of organisms reflects disease progression, as metabolic disturbances can often be observed several years prior to the symptomatic phase (Shah et al., 2010; Wang et al., 2011; Sabatine et al., 2005). Therefore, using metabolomics to identify early-stage biomarkers of complex phenotypes, such as prediabetes and COVID-19 susceptibility, has gained increased interest (Wang-Sattler et al., 2012; Julkunen et al., 2021). While two-sample MR approaches using metabolites as single exposure have revealed biomarkers for several diseases (Qian et al., 2021; Lord et al., 2021; Porcu et al., 2021a), these analyses focused on the prediction of disease risk rather than on deciphering the mechanisms of discovered associations.

Integrating transcriptomics with metabolomics data can provide insights into how metabolites are regulated, elucidating targetable functional mechanisms. Here, we develop a framework based on established MR methodology that hypothesizes a mediating role of the metabolome in the transcript-to-phenotype axis, with the primary exposure being defined as an upstream omic layer, namely gene expression. Specifically, our integrative MR analysis combines summary-level multi-omics (i.e., GWAS, eQTL, and mQTL) data to compute the indirect effect of gene expression on complex traits mediated by metabolites in three steps (Figure 1). First, we map the transcriptome to the metabolome by identifying causal associations between transcripts and metabolites. Next, we screen metabolites for downstream causal effects on 28 complex phenotypes, resulting in the identification of gene expression → metabolite → phenotype cascades (Figure 1A). In parallel, we prioritize trait-associated genes by testing the association of transcripts with phenotypes (Figure 1B). Third, for transcripts found to causally influence either a metabolite (A) or a complex phenotype (B), we test whether the identified target genes exert their effect on the phenotype through the metabolite using multivariable MR (MVMR; Figure 1C). Finally, we carried out extensive power analyses to determine under which conditions the mediation analysis (Figure 1C) outperforms the conventional exposure-outcome MR framework (Figure 1B).

Figure 1 with 1 supplement see all

Download asset Open asset

Workflow overview.

(A) Estimation of the causal transcript-to-metabolite and metabolite-to-phenotype effects through univariable Mendelian randomization (MR). (B) Estimation of the causal transcript-to-phenotype effects through univariable transcriptome-wide MR (TWMR). (C) Estimation of the direct (i.e., not mediated by the metabolites) and mediated effect of transcripts on phenotypes through multivariable MR (MVMR) by accounting for mediation through the metabolome.

Results

Mapping the transcriptome onto the metabolome

We applied univariable MR to identify metabolites whose levels are causally influenced by transcript levels in whole blood (Figure 1A). Summary statistics for cis-eQTLs stem from the eQTLGen Consortium meta-analysis of 19,942 transcripts in 31,684 individuals (Võsa et al., 2021), while summary statistics for mQTLs originate from a meta-analysis of 453 metabolites in 7824 individuals from two independent European cohorts: TwinsUK (N = 6056) and KORA (N = 1768) (Shin et al., 2014). After selecting SNPs included in both the eQTL and mQTL studies, our analysis was restricted to 7884 transcripts with ≥3 instrumental variables (IVs) (see Methods, Figure 1—figure supplement 1A) and 242 metabolites with an identifier in the Human Metabolome Database (HMDB) (Wishart et al., 2022) (see Methods, Supplementary file 1a). By testing each gene for association with the 242 metabolites, we detected 96 genes whose transcript levels causally impacted 75 metabolites, resulting in 133 unique transcript-metabolite associations (FDR 5% considering all 1,907,690 instrumentable gene-metabolite pairs; Supplementary file 1b). Most involved genes (86%; 83/96) were causally influencing the level of a single metabolite, with some notable exceptions acting as mQTL hubs, such as TMEM258 and FADS2 both affecting the same 11 metabolites, followed by FADS1 affecting a subset of six metabolites. While only 5 (3.8%) of the 133 associations were reported in HMDB, an automated literature review (see Methods) identified a match for 22 (16.5%) of the identified transcript-metabolite pairs (Supplementary file 1b).

Mapping the metabolome onto complex phenotypes

Univariable metabolome-wide MR (MWMR; Figure 1A) was used to identify causal relationships between 48 metabolites with ≥3 IVs (Figure 1—figure supplement 1B) and 28 complex phenotypes. The latter include a wide range of anthropometric traits, cardiovascular assessments, and blood biomarkers, whose summary statistics originate from the UK Biobank (UKB) (Bycroft et al., 2018;Supplementary file 1c). Overall, 34 metabolites were associated with at least one phenotype (FDR 5% considering all 1344 metabolite-phenotype pairs), resulting in 132 unique metabolite-phenotype associations (Supplementary file 1d).

Mapping the transcriptome onto complex phenotypes

We applied univariable transcriptome-wide MR (TWMR [Porcu et al., 2019] Figure 1B) to identify associations between expression levels of 10,435 transcripts from the eQTLGen Consortium with ≥3 IVs (Figure 1—figure supplement 1C) measured in both exposure and outcome datasets and the same 28 UKB phenotypes described in the previous section (Supplementary file 1c). In total, 5140 transcripts associated with at least one phenotype (FDR 5% considering all 292,170 gene-phenotype pairs) resulting in 13,141 unique transcript-phenotype associations (Supplementary file 1e).

Mapping metabolome-mediated effects of the transcriptome onto complex phenotypes

The mapping of putative causal effects performed in the previous steps provides the opportunity to infer the mediating role of the metabolome in biological processes leading to transcript-phenotype associations. We combined the 133 transcript-metabolite (FDR ≤5%) and 132 metabolite-trait (FDR ≤5%) associations to pinpoint 216 transcript-metabolite-phenotype causal triplets (FDR = 1–0.95² = 9.75%) (Supplementary file 1f). Among the 37 triplets for which the transcript and metabolite had previously been linked through automated literature review, none remained after incorporating a third term for the phenotype in the search and manually removing abstracts for which the search terms were used in an erroneous context. Relaxing the search criteria by omitting the metabolite term, 13/37 (35%) triplets returned at least one match for the gene-trait association.

For each of these 216 putative mechanisms, an MVMR approach to compute the direct effect of gene expression on the phenotype was applied (see Methods; Figure 1C; Supplementary file 1f). Regressing direct effects ( $α_{d}$ ) on total effects ( $α_{T P}$ ) and accounting for regression dilution bias (see Methods; Figure 2A), it was estimated that 77% [95% CI: 70–85%] of the transcript effect on the phenotype was direct and thus not mediated by the metabolites (Figure 2B).

Figure 2

Download asset Open asset

Direct and mediated effects.

(A) Graphical representation of the multivariable Mendelian randomization (MVMR) framework for mediation analysis: DNA represents genetic instrumental variables (IVs) chosen to be directly associated with either the exposure (transcript; $β_{e Q T L}$ ) or the mediator (metabolite; $β_{m Q T L}$ ) through summary statistics. The effect of these IVs on the outcome (phenotype; $β_{G W A S}$ ) originates from genome-wide association studies (GWASs) summary statistics. Total effects $α_{T P}$ of transcripts on phenotypes are estimated by transcriptome-wide Mendelian randomization (TWMR), while direct effects $α_{d}$ are estimated by MVMR. Total effects $α_{T P}$ are assumed to equal the sum of the direct $α_{d}$ and indirect $α_{i}$ (i.e., mediated) effects, the two former being depicted in B. (B) Direct ( $α_{d}$ ; y-axis) and total ( $α_{T P}$ ; x-axis) effects for the 216 transcript-metabolite-trait causal triplets. The dashed line represents the identity, while the purple line represents the regression line with a shaded 95% confidence interval. Data related to Figure 2 panel B are available in Figure 2—source data 1.

Figure 2—source data 1 Direct and mediated effects. Total ( $α_{T P}$ ; transcript-to-phenotype effect) and direct ( $α_{d}$ ; direct_effect) effects for the 216 transcript-metabolite-trait causal triplets involving the listed transcript (Gene_ID), metabolite (Shin_ID), and complex phenotype. This file relates to Figure 2B.: https://cdn.elifesciences.org/articles/81097/elife-81097-fig2-data1-v1.xlsx
Download elife-81097-fig2-data1-v1.xlsx

Molecular mechanisms of genotype-to-phenotype associations

Dissecting causal triplets allows gaining mechanistc insights into biological pathways linking genes to phenotypes. For instance, expression of TMEM258 [MIM: 617615], FADS1 [MIM: 606148], and FADS2 [MIM: 606149], all mapping to a region on chromosome 11 (Figure 3A), were found to influence a total of 17 complex phenotypes through modulation of 1-arachidonoylglycerophosphocholine (LPC(20:4); HMDB0010395; $α_{T M E M 258 \to L P C (20 : 4)} = - 1.02$ ; P = 8.0×10^–81; $α_{F A D S 1 \to L P C (20 : 4)} = - 0.39$ ; P = 4.6×10^-15; $α_{F A D S 2 \to L P C (20 : 4)} = - 0.63$ ; P = 5.1×10^–62), 1-arachidonoylglycerophosphoethanolamine (LPE(20:4);HMDB0011517; $α_{T M E M 258 \to L P E (20 : 4)} = - 0.68$ ; P = 1.1×10^–37; $α_{F A D S 1 \to L P E (20 : 4)} = - 0.30$ ; P = 1.4×10^–07; $α_{F A D S 2 \to L P E (20 : 4)} = - 0.37$ ; P = 1.2×10^–18), and 1-arachidonoylglycerophosphoinositol (LPI(20:4); HMDB0061690; $α_{T M E M 258 \to L P I (20 : 4)} = - 0.51$ ; P = 8.2×10^–18; $α_{F A D S 2 \to L P I (20 : 4)} = - 0.28$ ; P = 6.3×10^–16) levels (Figure 3B–C). These results align with the known pleiotropy of the region (i.e., >6000 associations reported in the GWAS Catalog as of May 2022). Interestingly, involved metabolites are complex lipids synthesized from arachidonic acid, a product of the rate-limiting enzymes encoded by FADS1 and FADS2 (Figure 3B). Recently, polymorphisms affecting the expression of these genes were shown to associate with the levels of over 50 complex lipids, including the ones identified by our study (Reynolds et al., 2020). Overall, this example illustrates how our method can capture meaningful biological associations and shed light on underlying molecular pathways of pleiotropy.

Figure 3

Download asset Open asset

Molecular pleiotropy at the *FADS* locus.

(A) Genome browser (GRCh37/hg19) view of the genomic region on chromosome 11 encompassing *TMEM258*, *FADS1*, and *FADS2* (red). (B) Diagram of the mediation signals detected for *TMEM258*, *FADS1*, and *FADS2*. Two of the implicated genes encode enzymes involved in arachidonic synthesis (purple). Involved genes impact 17 phenotypes (pink) through alteration of the levels of three metabolites, 1-arachidonoylglycerophosphocholine (LPC(20:4)), 1-arachidonoylglycerophosphoethanolamine (LPE(20:4)), and 1-arachidonoylglycerophosphoinositol (LPI(20:4)) whose structure is depicted (orange). (C) Network of the 65 transcript-metabolite-trait causal triplets involving *TMEM258*, *FADS1*, and *FADS2*. Nodes represent genes (purple), metabolites (orange), or phenotypes (pink). Edges indicate the direction of the effects estimated through univariable Mendelian randomization. Width of edges is proportional to effect size and color indicates if the effect is positive (red) or negative (blue).

Power analysis

Importantly, only 42% (90/216) of the causal triplets showed a significant total transcript-to-phenotype effect (i.e., estimated by TWMR), suggesting that the method lacks power under current settings. To characterize the parameter regime where the power to detect indirect effects is larger than it is for total effects, we performed simulations using different settings for the mediated effect. In each scenario we evaluated 500 transcripts and 80 metabolites and varied two parameters characterizing the mediation:

The proportion ( $ρ$ ) of direct ( $α_{d}$ ) to total ( $α_{T P}$ ) effect (i.e., effect not mediated by the metabolite) from –2 to 2 to cover the cases where direct and mediated effect have opposite directions (51 values).
The ratio ( $σ$ ) between the transcript-to-metabolite ( $α_{T M}$ ) and the metabolite-to-phenotype ( $α_{M P}$ ) effects, exploring the range from 0.1 to 10 (51 values).

Transcripts were simulated with 6% heritability (i.e., median $h^{2}$ in the eQTLGen data) and a causal effect of 0.035 (i.e., ~65% of power in TWMR at α=0.05) on a phenotype. Each scenario was simulated 10 times and results were averaged to assess the mean difference in power (see Methods).

Simulations show that with current sample sizes (i.e., $N_{G W A S} = 300,000$ , $N_{e Q T L} = 32, 000$ , and $N_{m Q T L} = 8000$ ), when $α_{M P} > α_{T M}$ (i.e., $σ < 1$ ), TWMR has increased power to detect significant transcript-to-phenotype associations, especially when $ρ > 0$ (i.e., direct and total effect have the same direction (Figure 4A)). However, for all 216 causal triplets, we observed $σ > 1$ (Figure 4—figure supplement 1). Under this condition, and assuming that the total effect of the transcript on the phenotype is dominated by the effect mediated by the metabolite (i.e., $ρ < 0.5$ and $ρ > 1.5$ ), TWMR had less power than the approach identifying mediators (Figure 4A), confirming that significant associations were missed by TWMR due to power issues related to the proportion of mediated effect.

Figure 4 with 2 supplements see all

Download asset Open asset

Power comparison between transcriptome-wide Mendelian randomization (TWMR) and multivariable Mendelian randomization (MVMR).

Heatmap showing the difference in statistical power between TWMR and mediation analysis through MVMR at current (A; $N = 8000$ ) and realistic future (B; $N = 90, 000$ ) metabolic quantitative trait loci (mQTL) dataset sample sizes. The x-axis shows the proportion ( $ρ$ ) of direct ( $α_{d}$ ) to total ( $α_{T P}$ ) effect (i.e., effect not mediated by the metabolite) ranging from –2 to 2, arrows indicating increasing proportion of direct effect. The y-axis shows the ratio ( $σ$ ) between the transcript-to-metabolite ( $α_{T M}$ ) and the metabolite-to-phenotype ( $α_{M P}$ ) effects, ranging from 0.1 to 10. Red vs. gray indicates higher power for TWMR vs. mediation analysis, respectively, while white represents equal power between the two approaches. Data related to Figure 4 panels A and B are available in Figure 4—source data 1 and Figure 4—source data 2, respectively.

Figure 4—source data 1 Difference in statistical power between transcriptome-wide Mendelian randomization (TWMR) and mediation analysis at N = 8000 metabolic quantitative trait locus (mQTL) dataset sample size. Each cell represents the mean difference in power between TWMR and mediation analysis for a given scenario across 10 simulations. Rows reflect decreasing ratio between transcript-to-metabolite and metabolite-to-phenotype effects from 10 to 0.1 (sigma). Columns reflect increasing proportion of direct to total effect from –2 to 2 (rho). This file relates to Figure 4A.: https://cdn.elifesciences.org/articles/81097/elife-81097-fig4-data1-v1.xlsx
Download elife-81097-fig4-data1-v1.xlsx
Figure 4—source data 2 Difference in statistical power between transcriptome-wide Mendelian randomization (TWMR) and mediation analysis at N = 90,000 metabolic quantitative trait locus (mQTL) dataset sample size. Each cell represents the mean difference in power between TWMR and mediation analysis for a given scenario across 10 simulations. Rows reflect decreasing ratio between transcript-to-metabolite and metabolite-to-phenotype effects from 10 to 0.1 (sigma). Columns reflect increasing proportion of direct to total effect from –2 to 2 (rho). This file relates to Figure 4A.: https://cdn.elifesciences.org/articles/81097/elife-81097-fig4-data2-v1.xlsx
Download elife-81097-fig4-data2-v1.xlsx

Repeating the simulations with an mQTL sample size of 90,000, nearing state-of-the-art sample sizes (Lotta et al., 2021), leads to a strong shift in the above-described trends (Figure 4B). Specifically, when the effect of the transcript on the phenotype is dominated by the effect mediated by the metabolite ( $ρ < 0.3$ and $ρ > 1.7$ ), mediation analysis has more power than TWMR when $σ > 0.2$ . For larger proportions of direct effect, TWMR has increased power the more $σ$ differs from 1. In line with the increased power of mediation analysis with larger mQTL datasets, the gain in power of mediation analysis over TWMR decreases with decreasing mQTL dataset sample sizes (ranging between N = 1000 and N = 4000; Figure 4—figure supplement 2), indicating that our approach is dependent on large sample sizes to reach its full potential.

Identifying new genotype-to-phenotype associations

The 126 triplets that were not identified through TWMR due to power issues represent putative new causal relations. This is well illustrated by a proof-of concept example involving ANKH [MIM: 605145] and calcium levels, for which 48 publications were identified through automated literature review (Supplementary file 1f). While the TWMR effect of ANKH expression on calcium levels was not significant ( $α_{A N K H \to c a l c i u m} = - 0.02$ ; P=0.03), ANKH expression decreased citrate levels ( $α_{A N K H \to c i t r a t e} = - 0.30$ ; P=2.2×10^–06), which itself increased serum calcium levels ( $α_{c i t r a t e \to c a l c i u m} = 0.07$ ; P=6.5×10^–0). Mutations in ANKH have been associated with several rare mineralization disorders [MIM: 123000, 118600] (Williams, 2016) due to the gene encoding a transmembrane protein that channels inorganic pyrophosphate to the extracellular matrix, where at low concentrations it inhibits mineralization (Ho et al., 2000). Recently, a study proposed that ANKH instead exports ATP to the extracellular space (where it is then rapidly converted to inorganic pyrophosphate), along with citrate (Szeri et al., 2020). Citrate has a high binding affinity for calcium and influences its bioavailability by complexing calcium-phosphate during extracellular matrix mineralization and releasing calcium during bone resorption (Granchi et al., 2019). Together, our data support the role of ANKH in calcium homeostasis through regulation of citrate levels, connecting previously established independent links into a causal triad.

In another example, SLC6A12 [MIM: 603080], which encodes the betaine/GABA transporter-1, involved in betaine and GABA uptake (Borden et al., 1995), was identified as a negative regulator of betaine ( $α_{S L C 6 A 12 \to b e t a i n e} = - 0.37$ ; P = 8.2×10^–08). While blood betaine levels negatively impacted serum creatinine levels ( $α_{b e t a i n e \to c r e a t i n i n e} = - 0.06$ ; P = 1.7×10^–07), the effect of SLC6A12 expression on creatinine was not significant ( $α_{S L C 6 A 12 \to c r e a t i n i n e} = 0.02$ ; P = 1.5×10^–03). This observation is particularly interesting given that betaine acts as a protective renal osmolyte whose plasma and kidney tissue concentration were found to be downregulated in renal ischemia/reperfusion injury (Jouret et al., 2016; Wei et al., 2014) and whose urine levels have been proposed as a biomarker for chronic kidney disease progression (Gil et al., 2018). As both renal conditions are commonly monitored through serum creatinine levels, our data support the critical role of osmolyte homeostasis in renal health.

Discussion

In this study, we combined MR approaches integrating eQTL, mQTL, and GWAS summary statistics to explore the role of the metabolome in mediating the effect of the transcriptome on complex phenotypes. Applied to 28 medically relevant traits, our approach revealed 216 causal transcript-metabolite-phenotype triplets. Our automated literature review indicates that while some detected associations were previously reported, a large fraction, especially among the triplets, appears to be novel. It should be noted that the number of previously reported associations is likely underestimated as our approach does not account for all synonyms of a given feature and requires the terms to appear in the title or abstract of the publication. This makes it more likely for hypothesis-driven studies, inherently biased toward well-studied genes and metabolites, to be identified. Conversely, high-throughput, hypothesis-free studies that report the given association in a supplemental table are likely to be missed. Furthermore, due to its automated nature, our search is context-blind, so that some of the identified studies might report negative results, associations only under specific conditions (e.g., different organisms, experimental settings), or usage of the search term with a different meaning. To attenuate the latter, we also performed manual review of the retained abstracts for transcript-to-phenotype searches. While flawed, this rough estimate of the amount of existing evidence supporting our findings can be interpreted in combination with other lines of evidence. For instance, among the 90 signals that were also identified through TWMR, 93% showed a directionally concordant effect between the transcript-to-phenotype, transcript-to-metabolite, and metabolite-to-phenotype estimates (i.e., sign of product of the transcript-to-metabolite and metabolite-to-phenotype effects agrees with the sign of the transcript-to-phenotype effect). In these situations, dissection of causal effects provides clues as to the molecular mechanism through which involved genes modify complex phenotypes. This information is particularly valuable to identify key molecular mediators of highly pleiotropic genetic regions, such as the TMEM258/FADS1/FADS2 locus (Figure 3). While transcript levels of these genes affected eleven metabolites, three complex lipids were highlighted as strong molecular mediators of the transcript-to-phenotype effects.

Strikingly, 58% of the 216 causal transcript-metabolite-phenotype triplets were missed by TWMR – an approach that only considers gene expression and GWAS data. We highlight two novel but biologically plausible mechanisms linking ANKH to calcium levels through modulation of citrate and SLC6A12 to serum creatinine levels through regulation of the renal osmolyte betaine. Simulation analyses showed that these signals were likely missed by TWMR due to lack of power, as mediation analysis is better suited to detect associations with a low direct to total effect proportion and stronger transcript-to-metabolite than metabolite-to-phenotype effect. Promisingly, our simulations showed that mediation analysis becomes increasingly powerful over a wider range of parameter settings as the sample size of the mediator QTL study increases, highlighting the importance of generating large and publicly available molQTL datasets that can help to unravel functional gene-to-phenotype mechanisms.

As illustrated through the selected examples, a large fraction of detected mediations involves genes encoding metabolic enzymes or transporters/channels, with an enrichment for ‘secondary active transmembrane transporter activity’, for example (GO:0015291; FDR=0.021; background: 7884 genes with ≥3 IVs assessed through TWMR; STRING database). Matching the finding that the most likely effector genes of mQTLs are enriched for pathway-relevant enzymes and transporters (Smith et al., 2022), these results are not surprising given that the proteins encoded by these genes directly interact with metabolites, making it more likely that the effect of changes in their expression is mediated by metabolites. While our method is well suited to detect such effects, interpretation of discovered mediations is limited by the lack of spatial resolution of the mQTL data. Access to metabolite concentrations in different cellular compartments (e.g., extracellular space, cytosol, mitochondrial matrix, etc.) would generate more fine-tuned mechanistic hypotheses that consider the directionality of metabolite fluxes.

The observation that 77% of the transcript’s effect on the phenotype is not mediated by metabolites suggests that either true direct effects are frequent or that other unassessed metabolites or molecular layers (e.g., proteins, post-translational modifications, etc.) play a crucial role in mediation. It is to note that in the presence of unmeasured mediators or measured mediators without genetic instruments, our mediation estimates are lower bounds of the total existing mediation. In addition, unmeasured mediators sharing genetic instruments with the measured ones can modify result interpretation as some of the observed mediators may simply be correlates of the true underlying mediators. While this is a limitation of all MR methods, metabolic networks may harbor particularly large number of genetically correlated metabolite species. Similarly, owing to linkage disequilibrium and regulatory variants affecting multiple genes, transcripts from adjacent genes might appear to be involved in the same signals, as exemplified with the TMEM258/FADS1/FADS2 locus (Figure 3). While literature supports the role of the FADS genes, one cannot exclude a role for TMEM258, nor disentangle the specific function of FADS1 and FADS2. Thanks to the flexibility of the proposed framework, we expect that in the future and upon availability of ever larger and more diverse datasets, our method could be applied to estimate the relative contribution of currently unassessed mediators in translating genotypic cascades.

Another consideration is that complex phenotypes can have a stronger impact on gene expression than the opposite (Porcu et al., 2021b). Due to the lack of genome-wide trans-eQTL association summary statistics, our method does not investigate reverse causality on metabolites and gene expression, nor the role of metabolites as regulators of gene expression. Metabolites might also integrate the effect of several transcripts (i.e., multiple transcripts causally impact the levels of the same metabolite) before affecting complex phenotypes (Supplementary file 1g) or multiple metabolites may jointly mediate the impact of a single transcript. Modelling the latter phenomenon, which is beyond the scope of our current work, requires the development of structural equation models accounting for such effects and will eventually lead to a more comprehensive modelling of causal relations in complex biological networks, nuancing the interpretation of the molecular mechanisms shaping complex traits.

In conclusion, we developed a modular MR framework that has increased power over classical MR approaches to detect causal transcript-to-phenotype relationships when these are mediated by alteration of metabolite levels and is likely to become increasingly powerful upon release of larger molQTL datasets.

Methods

Univariable MR analyses

TWMR and MWMR (Porcu et al., 2019) were used to estimate the causal effects of transcript and metabolite levels (exposure) on various outcomes. For each transcript/metabolite, using inverse-variance weighted (IVW) method for summary statistics (Burgess et al., 2013), the causal effect of the molecular traits on the outcome was defined as

\hat{α} = {(β^{'} C^{- 1} β)}^{- 1} (β^{'} C^{- 1} γ)

Here, $β$ is a vector of length n containing the standardized effect size of n independent SNPs on the gene/metabolite, derived from eQTL/mQTL studies, with $β^{`}$ being the transpose of $β$ . $γ$ is a vector of length n containing the standardized effect size of each SNP on the outcome. C is the pairwise LD matrix between the n SNPs. The standardized effect sizes for molecular and outcome GWASs were obtained from Z-score of summary statistics standardized by the square root of the sample size to be on the same standard deviation scale.

IVs were selected as autosomal, non-strand ambiguous, independent (r² <0.01), and significant ( $P_{e Q T L} < 1.8 \times 10^{- 05}$ / $P_{m Q T L} < 1.0 \times 10^{- 07}$ ) eQTL/mQTLs available in the UK10K reference panel (Huang et al., 2015) using PLINK (v1.9) (Chang et al., 2015). As retained SNPs are independent, we used the identity matrix to approximate C. SNPs with larger effects on the outcome than on the exposure were removed, as these potentially indicate violation of the MR assumptions (i.e., likely reverse causality and/or confounding).

The variance of α can be calculated approximately by the Delta method

v a r (\hat{α}) = {(\frac{\partial \hat{α}}{\partial β})}^{2} * v a r (\hat{β}) + {(\frac{\partial \hat{α}}{\partial γ})}^{2} * v a r (\hat{γ}) + (\frac{\partial \hat{α}}{\partial β}) * (\frac{\partial \hat{α}}{\partial γ}) * c o v (\hat{β}, \hat{γ})

where cov( $β$ , $γ$ ) is 0 if $β$ and $γ$ are estimated from independent samples. The causal effect Z-statistic for transcript/metabolite i was defined as $\frac{{\hat{α}}_{i}}{S E ({\hat{α}}_{i})}$ , where $S E ({\hat{α}}_{i}) = \sqrt{v a r {(\hat{α})}_{i, i}}$ .

The IVW method provides an unbiased estimate under the assumption that all genetic variants are valid IVs, that is, all three MR assumption hold. However, the third assumption (no pleiotropy) is easily violated, leading to inaccurate estimates when horizontal pleiotropy occurs (Verbanck et al., 2018). To test for the presence of pleiotropy, we used Cochran’s Q test (Bowden et al., 2015; Burgess et al., 2017) to assess whether there were significant differences between the MR effects of an instrument (i.e., $α β_{i}$ ) and the estimated effect of that instrument on phenotype/metabolite levels ( $γ_{i}$ ). We defined

d_{i} = γ_{i} - α β_{i}

and its variance as

v a r (d_{i}) = v a r (γ_{i}) + {(β_{i})}^{2} * v a r (α) + v a r (γ_{i}) * {(α)}^{2} + v a r (β_{i}) * v a r (α)

Next, the deviation of each SNP was tested using the test statistic

T_{i} = \frac{d_{i}^{2}}{v a r (d_{i})} \sim χ_{1}^{2}

When p<0.05, the SNP with largest | $d_{i}$ | was removed and the test was repeated.

Mediation analysis through MVMR analyses

An MVMR approach was used to dissect the total causal effect of transcript levels on phenotypes ( $α_{T P}$ ) into a direct ( $α_{d}$ ) and indirect ( $α_{i}$ ) effect measured through a metabolite. Through inclusion of a metabolite and its associated genetic variants (r² <0.01, p_mQTL<1 × 10^–07), the direct effect of gene expression on a phenotype can be estimated using a multivariable regression model (Burgess et al., 2013) as the first element of

\hat{α} = {(B^{'} C^{- 1} B)}^{- 1} (B^{'} C^{- 1} γ)

where $B$ is a matrix with two columns containing the standardized effect sizes of n IVs on transcript levels in the first column and on the metabolite levels in the second column, $γ$ is a vector of length n containing the standardized effect size of each SNP on the phenotype, and C is the pairwise LD matrix between the n SNPs.

The proportion of direct effect $(ρ)$ is calculated by regressing direct effects ( $α_{d}$ ) on total effects ( $α_{T P}$ ) and then correcting for regression dilution bias:

ρ_{c o r r e c t e d} = \frac{ρ}{\sqrt{1 - \frac{\sum {(S E (α_{T P}))}^{2}}{\sum α_{T P}^{2}}}}

Omics and traits summary statistics

Expression QTL data originated from the eQTLGen Consortium (Võsa et al., 2021) (N = 31,684), which includes cis-eQTLs (<1 Mb from gene center, two-cohort filter) for 19,250 transcripts (16,934 with at least one significant cis-eQTL at FDR <0.05 corresponding to p<1.8 × 10^–05). mQTL data originate from Shin et al., 2014, which used ultra-high performance liquid chromatography-tandem mass spectrometry to measure 486 whole blood metabolites in 7824 European individuals. Association analyses were carried out on ~2.1 million SNPs and are available for 453 metabolites at the Metabolomics GWAS Server (http://metabolomics.helmholtz-muenchen.de/gwas/). Among these metabolites, 242 were manually annotated with the HMDB identifiers (Supplementary file 1a) and used in this study. GWAS summary statistics for 28 outcome traits measured in the UKB originate from the Neale Lab (http://www.nealelab.is/uk-biobank/). Protein interactions with metabolites were downloaded from HMDB v5.0 (https://hmdb.ca/downloads/) and were used to annotate transcript-metabolites associations.

Automated literature review

An automated literature review of all transcript-metabolite associations (Supplementary file 1b) was conducted in PubMed (September 20, 2022) following the scheme:

(G e n e [T i t l e / A b s t r a c t]) A N D ((M e t_{H M D B} [T i t l e / A b s t r a c t]) O R (M e t_{S h i n} [T i t l e / A b s t r a c t]))

With $G e n e$ being the name of the gene whose transcript is involved in the association, ${M e t}_{H M D B}$ the involved metabolite’s common name on HMDB, and ${M e t}_{S h i n}$ the involved metabolite’s name as reported in Shin et al., 2014. Returned PubMed identifiers were retrieved (Supplementary file 1b).

For transcript-metabolite associations involved in a causal triplet and for which the transcript-metabolites returned at least one publication (Supplementary file 1f), the search was extend by (i) adding an additional search term for the trait (i.e., $A N D (t r a i t [T i t l e / A b s t r a c t])$ ) and (ii) substituting the metabolite term for the trait term. Returned PubMed identifiers were retrieved and corresponding abstracts were manually curated to exclude abstracts in which the search terms were used in a meaning other than the intended one (Supplementary file 1f).

Simulation analyses

Simulation analyses were conducted to assess the gain in power upon inclusion of metabolomics data in the MR framework. In the simulated scenario, a transcript has an effect on a phenotype mediated by a metabolite. Two parameters were allowed to vary: the proportion ( $ρ$ ) of direct effect (i.e., effect not mediated by the metabolite) and the ratio ( $σ$ ) between the effect of the transcript on the metabolite ( $α_{T M}$ ) and of the metabolite on the phenotype ( $α_{M P}$ ). Other parameters were fixed, including the heritability of the transcript at $h_{T}^{2} = 0.06$ (corresponding to the median $h^{2}$ in the eQTLGen data), the number of IVs $N_{I V s}$ at 6 (corresponding to the median number of IVs used in TWMR analyses). Effect sizes $β_{e Q T L}$ are from a normal distribution $β_{e Q T L} \sim N (0, \frac{h_{T}^{2}}{N_{I V s}})$ . The causal effect of the transcript on the phenotype ( $α_{T P}$ ) was fixed to 0.035, which results in ∼65% power to detect a significant effect with TWMR. These quantities allowed to define $β_{G W A S}$ as $β_{G W A S} = α_{T P} * β_{e Q T L} + ε_{P}$ , where $ε_{P} \sim N (0, \frac{1}{N_{G W A S}})$ with $N_{G W A S} = 300,000$ to reflect the sample size of UKB GWASs. The same vector of $β_{e Q T L}$ was used to define $β_{m Q T L}$ and estimate the causal effect of the transcript on the metabolite. $β_{m Q T L}$ was defined as $β_{m Q T L} = α_{T M} * β_{e Q T L} + ε_{M}$ , where $ε_{M} \sim N (0, \frac{1}{N_{m Q T L}})$ and $N_{m Q T L} = 8000$ to reflect the sample size of the mQTL study used in this work. Simulations were also performed at $N_{m Q T L} = 90,000,$ to reflect sample size of potential future studies and $N_{m Q T L} = 1000$ , $N_{m Q T L} = 2000$ and $N_{m Q T L} = 4000$ , to compare the two approaches’ power were the developed framework to be applied on existing smaller mQTL datasets. The total effect $α_{T P}$ can be expressed as $α_{T P} = α_{T M} * α_{M P} + α_{d i r e c t}$ , where $α_{d i r e c t}$ represents the direct effect of the transcript on the phenotype and $α_{T M} * α_{M P}$ is the indirect effect mediated by the metabolite. Equivalently, $α_{T M} * α_{M P} = α_{T P} * (1 - ρ)$ where $ρ = \frac{α_{d i r e c t}}{α_{T P}} .$ To assess the ratio between the effect of the transcript on the metabolite and the effect of the metabolite on the phenotype (i.e., $σ = α_{T M} / α_{M P}$ ), $α_{T M}$ can be expressed as $α_{T M} = \sqrt{α_{T P} * (1 - ρ) * σ}$ . Similarly, to estimate the effect of the metabolite on the phenotype, a metabolite with heritability $h_{M}^{2} = 0.04$ (corresponding to the median of $h^{2}$ in the KORA +TwinsUK mQTL data) and $N_{I V s} = 5$ (corresponding to the median number of IVs used in MWMR analyses) is considered. Effect size $β_{m Q T L}$ are from a normal distribution $β_{m Q T L} \sim N (0, \frac{h_{M}^{2}}{N_{I V s}})$ . These quantities allowed to define $β_{G W A S}$ as $β_{G W A S} = \sqrt{α_{T P} * (1 - ρ) / σ} * β_{m Q T L} + ε_{P}$ , where $ε_{P} \sim N (0, \frac{1}{N_{G W A S}})$ . Ranging $ρ$ and $σ$ from –2 to 2 and from 0.1 and 10, respectively, we run each simulation for 500 transcripts measuring 80 metabolites at each run and performed TWMR and MWMR starting from above-described $β_{e Q T L}$ , $β_{m Q T L}$ , and $β_{G W A S} .$ For each MR analysis the power to detect a significant association as well as the difference in power between TWMR and the mediation analyses (i.e., ${p o w e r}_{T P} - {p o w e r}_{T M} * {p o w e r}_{M P}$ ) was calculated. Each specific scenario was repeated 10 times and the average difference in power across simulation was plotted as a heatmap.

Data and code availability

All data used in this study are publicly available. GWAS summary statistics for outcome traits measured in the UKB originate from the Neale Lab (http://www.nealelab.is/uk-biobank/). eQTL data originated from the eQTLGen Consortium (https://www.eqtlgen.org) and was published in Võsa et al., 2021. mQTL data originate from Shin et al., 2014, and are available at the Metabolomics GWAS Server (http://metabolomics.helmholtz-muenchen.de/gwas/). The HMDB was used to annotate metabolites and the v5.0 release from November 9, 2021, of the ‘All proteins’ file was downloaded to extract transcript-metabolite interactions (https://hmdb.ca/downloads). PubMed was used for the automated literature review (https://pubmed.ncbi.nlm.nih.gov). The UCSC Genome Browser (https://genome.ucsc.edu/) was used to visualize the FADS locus, while the GWAS Catalog was used to assess the number of reported GWAS signals in the region (https://www.ebi.ac.uk/gwas/). The STRING database was used for the enrichment analysis (https://string-db.org/). Produced data is available as Supplementary file 1 and Source Data. Code used to perform analyses is freely available at https://github.com/eleporcu/Gene_Metab_Pheno; (Porcu, 2022 copy archived at swh:1:rev:c6bff8d094e369ff0d399751fc85fcd5ea250134).

Data availability

All data used in this study are publicly available. GWAS summary statistics for outcome traits measured in the UK Biobank originate from the Neale Lab (http://www.nealelab.is/uk-biobank/). eQTL data originated from the eQTLGen Consortium (https://www.eqtlgen.org) and was published in Vosa et al., 2021 [3]. mQTL data originate from Shin et al. 2014 [6], and are available at the Metabolomics GWAS Server (http://metabolomics.helmholtz-muenchen.de/gwas/). The Human Metabolome Database (HMDB) was used to annotate metabolites and the v5.0 release from 2021-11-09 of the "All proteins" file was downloaded to extract transcript-metabolite interactions (https://hmdb.ca/downloads). PubMed was used for the automated literature review (https://pubmed.ncbi.nlm.nih.gov). The UCSC Genome Browser (https://genome.ucsc.edu/) was used to visualize the FADS locus, while the GWAS Catalog was used to assess the number of reported GWAS signals in the region (https://www.ebi.ac.uk/gwas/). The STRING database was used for the enrichment analysis (https://string-db.org/). Produced data is available as Supplementary File 1 and Source Data. Code used to perform analyses is freely available at https://github.com/eleporcu/Gene_Metab_Pheno; (copy archived at swh:1:rev:c6bff8d094e369ff0d399751fc85fcd5ea250134).

References

(1995) Cloning and expression of a betaine/GABA transporter from human brain
Journal of Neurochemistry 64:977–984.

https://doi.org/10.1046/j.1471-4159.1995.64030977.x
- PubMed
- Google Scholar
(2015) Mendelian randomization with invalid instruments: effect estimation and bias detection through egger regression
International Journal of Epidemiology 44:512–525.

https://doi.org/10.1093/ije/dyv080
- PubMed
- Google Scholar
(2013) Mendelian randomization analysis with multiple genetic variants using summarized data
Genetic Epidemiology 37:658–665.

https://doi.org/10.1002/gepi.21758
- PubMed
- Google Scholar
1. Burgess S
2. Bowden J
3. Fall T
4. Ingelsson E
5. Thompson SG
(2017) Sensitivity analyses for robust causal inference from mendelian randomization analyses with multiple genetic variants
Epidemiology 28:30–42.

https://doi.org/10.1097/EDE.0000000000000559
- PubMed
- Google Scholar
1. Bycroft C
2. Freeman C
3. Petkova D
4. Band G
5. Elliott LT
6. Sharp K
7. Motyer A
8. Vukcevic D
9. Delaneau O
10. O’Connell J
11. Cortes A
12. Welsh S
13. Young A
14. Effingham M
15. McVean G
16. Leslie S
17. Allen N
18. Donnelly P
19. Marchini J
(2018) The UK Biobank resource with deep phenotyping and genomic data
Nature 562:203–209.

https://doi.org/10.1038/s41586-018-0579-z
- PubMed
- Google Scholar
1. Chang CC
2. Chow CC
3. Tellier LC
4. Vattikuti S
5. Purcell SM
6. Lee JJ
(2015) Second-Generation PLINK: rising to the challenge of larger and richer datasets
GigaScience 4:7.

https://doi.org/10.1186/s13742-015-0047-8
- PubMed
- Google Scholar
(2018) A bayesian framework for multiple trait colocalization from summary association statistics
Bioinformatics 34:2538–2545.

https://doi.org/10.1093/bioinformatics/bty147
- PubMed
- Google Scholar
(2018) Increased urinary osmolyte excretion indicates chronic kidney disease severity and progression rate
Nephrology, Dialysis, Transplantation 33:2156–2164.

https://doi.org/10.1093/ndt/gfy020
- PubMed
- Google Scholar
1. Gleason KJ
2. Yang F
3. Pierce BL
4. He X
5. Chen LS
(2020) Primo: integration of multiple GWAS and omics QTL summary statistics for elucidation of molecular mechanisms of trait-associated snps and detection of pleiotropy in complex traits
Genome Biology 21:236.

https://doi.org/10.1186/s13059-020-02125-w
- PubMed
- Google Scholar
(2019) Role of citrate in pathophysiology and medical management of bone diseases
Nutrients 11:2576.

https://doi.org/10.3390/nu11112576
- PubMed
- Google Scholar
1. Gusev A
2. Ko A
3. Shi H
4. Bhatia G
5. Chung W
6. Penninx BWJH
7. Jansen R
8. de Geus EJC
9. Boomsma DI
10. Wright FA
11. Sullivan PF
12. Nikkola E
13. Alvarez M
14. Civelek M
15. Lusis AJ
16. Lehtimäki T
17. Raitoharju E
18. Kähönen M
19. Seppälä I
20. Raitakari OT
21. Kuusisto J
22. Laakso M
23. Price AL
24. Pajukanta P
25. Pasaniuc B
(2016) Integrative approaches for large-scale transcriptome-wide association studies
Nature Genetics 48:245–252.

https://doi.org/10.1038/ng.3506
- PubMed
- Google Scholar
(2000) Role of the mouse ANK gene in control of tissue calcification and arthritis
Science 289:265–270.

https://doi.org/10.1126/science.289.5477.265
- PubMed
- Google Scholar
1. Hormozdiari F
2. Kostem E
3. Kang EY
4. Pasaniuc B
5. Eskin E
(2014) Identifying causal variants at loci with multiple signals of association
Genetics 198:497–508.

https://doi.org/10.1534/genetics.114.167908
- PubMed
- Google Scholar
1. Huang J
2. Howie B
3. McCarthy S
4. Memari Y
5. Walter K
6. Min JL
7. Danecek P
8. Malerba G
9. Trabetti E
10. Zheng HF
11. Gambaro G
12. Richards JB
13. Durbin R
14. Timpson NJ
15. Marchini J
16. Soranzo N
17. UK10K Consortium
(2015) Improved imputation of low-frequency and rare variants using the UK10K haplotype reference panel
Nature Communications 6:8111.

https://doi.org/10.1038/ncomms9111
- PubMed
- Google Scholar
(2016) Nuclear magnetic resonance metabolomic profiling of mouse kidney, urine and serum following renal ischemia/reperfusion injury
PLOS ONE 11:e0163021.

https://doi.org/10.1371/journal.pone.0163021
- PubMed
- Google Scholar
(2021) Metabolic biomarker profiling for identification of susceptibility to severe pneumonia and COVID-19 in the general population
eLife 10:e63033.

https://doi.org/10.7554/eLife.63033
- PubMed
- Google Scholar
1. Lord J
2. Jermy B
3. Green R
4. Wong A
5. Xu J
6. Legido-Quigley C
7. Dobson R
8. Richards M
9. Proitsi P
(2021) Mendelian randomization identifies blood metabolites previously linked to midlife cognition as causal candidates in alzheimer’s disease
PNAS 118:e2009808118.

https://doi.org/10.1073/pnas.2009808118
- PubMed
- Google Scholar
1. Lotta LA
2. Pietzner M
3. Stewart ID
4. Wittemans LBL
5. Li C
6. Bonelli R
7. Raffler J
8. Biggs EK
9. Oliver-Williams C
10. Auyeung VPW
11. Luan J
12. Wheeler E
13. Paige E
14. Surendran P
15. Michelotti GA
16. Scott RA
17. Burgess S
18. Zuber V
19. Sanderson E
20. Koulman A
21. Imamura F
22. Forouhi NG
23. Khaw K-T
24. MacTel Consortium
25. Griffin JL
26. Wood AM
27. Kastenmüller G
28. Danesh J
29. Butterworth AS
30. Gribble FM
31. Reimann F
32. Bahlo M
33. Fauman E
34. Wareham NJ
35. Langenberg C
(2021) A cross-platform approach identifies genetic regulators of human metabolism and health
Nature Genetics 53:54–64.

https://doi.org/10.1038/s41588-020-00751-5
- PubMed
- Google Scholar
1. MacArthur J
2. Bowler E
3. Cerezo M
4. Gil L
5. Hall P
6. Hastings E
7. Junkins H
8. McMahon A
9. Milano A
10. Morales J
11. Pendlington ZM
12. Welter D
13. Burdett T
14. Hindorff L
15. Flicek P
16. Cunningham F
17. Parkinson H
(2017) The new NHGRI-EBI catalog of published genome-wide association studies (GWAS catalog)
Nucleic Acids Research 45:D896–D901.

https://doi.org/10.1093/nar/gkw1133
- PubMed
- Google Scholar
1. Min JL
2. Hemani G
3. Hannon E
4. Dekkers KF
5. Castillo-Fernandez J
6. Luijk R
7. Carnero-Montoro E
8. Lawson DJ
9. Burrows K
10. Suderman M
11. Bretherick AD
12. Richardson TG
13. Klughammer J
14. Iotchkova V
15. Sharp G
16. Al Khleifat A
17. Shatunov A
18. Iacoangeli A
19. McArdle WL
20. Ho KM
21. Kumar A
22. Söderhäll C
23. Soriano-Tárraga C
24. Giralt-Steinhauer E
25. Kazmi N
26. Mason D
27. McRae AF
28. Corcoran DL
29. Sugden K
30. Kasela S
31. Cardona A
32. Day FR
33. Cugliari G
34. Viberti C
35. Guarrera S
36. Lerro M
37. Gupta R
38. Bollepalli S
39. Mandaviya P
40. Zeng Y
41. Clarke T-K
42. Walker RM
43. Schmoll V
44. Czamara D
45. Ruiz-Arenas C
46. Rezwan FI
47. Marioni RE
48. Lin T
49. Awaloff Y
50. Germain M
51. Aïssi D
52. Zwamborn R
53. van Eijk K
54. Dekker A
55. van Dongen J
56. Hottenga J-J
57. Willemsen G
58. Xu C-J
59. Barturen G
60. Català-Moll F
61. Kerick M
62. Wang C
63. Melton P
64. Elliott HR
65. Shin J
66. Bernard M
67. Yet I
68. Smart M
69. Gorrie-Stone T
70. BIOS Consortium
71. Shaw C
72. Al Chalabi A
73. Ring SM
74. Pershagen G
75. Melén E
76. Jiménez-Conde J
77. Roquer J
78. Lawlor DA
79. Wright J
80. Martin NG
81. Montgomery GW
82. Moffitt TE
83. Poulton R
84. Esko T
85. Milani L
86. Metspalu A
87. Perry JRB
88. Ong KK
89. Wareham NJ
90. Matullo G
91. Sacerdote C
92. Panico S
93. Caspi A
94. Arseneault L
95. Gagnon F
96. Ollikainen M
97. Kaprio J
98. Felix JF
99. Rivadeneira F
100. Tiemeier H
101. van IJzendoorn MH
102. Uitterlinden AG
103. Jaddoe VWV
104. Haley C
105. McIntosh AM
106. Evans KL
107. Murray A
108. Räikkönen K
109. Lahti J
110. Nohr EA
111. Sørensen TIA
112. Hansen T
113. Morgen CS
114. Binder EB
115. Lucae S
116. Gonzalez JR
117. Bustamante M
118. Sunyer J
119. Holloway JW
120. Karmaus W
121. Zhang H
122. Deary IJ
123. Wray NR
124. Starr JM
125. Beekman M
126. van Heemst D
127. Slagboom PE
128. Morange P-E
129. Trégouët D-A
130. Veldink JH
131. Davies GE
132. de Geus EJC
133. Boomsma DI
134. Vonk JM
135. Brunekreef B
136. Koppelman GH
137. Alarcón-Riquelme ME
138. Huang R-C
139. Pennell CE
140. van Meurs J
141. Ikram MA
142. Hughes AD
143. Tillin T
144. Chaturvedi N
145. Pausova Z
146. Paus T
147. Spector TD
148. Kumari M
149. Schalkwyk LC
150. Visscher PM
151. Davey Smith G
152. Bock C
153. Gaunt TR
154. Bell JT
155. Heijmans BT
156. Mill J
157. Relton CL
(2021) Genomic and phenotypic insights from an atlas of genetic effects on DNA methylation
Nature Genetics 53:1311–1321.

https://doi.org/10.1038/s41588-021-00923-x
- PubMed
- Google Scholar
(2017) Estimating the causal tissues for complex traits and diseases
Nature Genetics 49:1676–1683.

https://doi.org/10.1038/ng.3981
- PubMed
- Google Scholar
(2019) Mendelian randomization integrating GWAS and eqtl data reveals genetic determinants of complex and clinical traits
Nature Communications 10:3300.

https://doi.org/10.1038/s41467-019-10936-0
- PubMed
- Google Scholar
1. Porcu E
2. Gilardi F
3. Darrous L
4. Yengo L
5. Bararpour N
6. Gasser M
7. Marques-Vidal P
8. Froguel P
9. Waeber G
10. Thomas A
11. Kutalik Z
(2021a) Triangulating evidence from longitudinal and mendelian randomization studies of metabolomic biomarkers for type 2 diabetes
Scientific Reports 11:6197.

https://doi.org/10.1038/s41598-021-85684-7
- PubMed
- Google Scholar
1. Porcu E
2. Sadler MC
3. Lepik K
4. Auwerx C
5. Wood AR
6. Weihs A
7. Sleiman MSB
8. Ribeiro DM
9. Bandinelli S
10. Tanaka T
11. Nauck M
12. Völker U
13. Delaneau O
14. Metspalu A
15. Teumer A
16. Frayling T
17. Santoni FA
18. Reymond A
19. Kutalik Z
(2021b) Differentially expressed genes reflect disease-induced rather than disease-causing changes in the transcriptome
Nature Communications 12:5647.

https://doi.org/10.1038/s41467-021-25805-y
- PubMed
- Google Scholar
Software
1. Porcu E
(2022) Gene_Metab_Pheno, version swh:1:rev:c6bff8d094e369ff0d399751fc85fcd5ea250134
Software Heritage.

https://archive.softwareheritage.org/swh:1:dir:d314533b7c8b4abe22b4fdf9be7499279fdc1e6e;origin=https://github.com/eleporcu/Gene_Metab_Pheno;visit=swh:1:snp:ab86f3b131ad58088c7db0744b97590acbde2121;anchor=swh:1:rev:c6bff8d094e369ff0d399751fc85fcd5ea250134
1. Qian L
2. Fan Y
3. Gao F
4. Zhao B
5. Yan B
6. Wang W
7. Yang J
8. Ma X
(2021) Genetically determined levels of serum metabolites and risk of neuroticism: a Mendelian randomization study
The International Journal of Neuropsychopharmacology 24:32–39.

https://doi.org/10.1093/ijnp/pyaa062
- PubMed
- Google Scholar
1. Reynolds LM
2. Dutta R
3. Seeds MC
4. Lake KN
5. Hallmark B
6. Mathias RA
7. Howard TD
8. Chilton FH
(2020) Fads genetic and metabolomic analyses identify the ∆5 desaturase (FADS1) step as a critical control point in the formation of biologically important lipids
Scientific Reports 10:15873.

https://doi.org/10.1038/s41598-020-71948-1
- PubMed
- Google Scholar
1. Sabatine MS
2. Liu E
3. Morrow DA
4. Heller E
5. McCarroll R
6. Wiegand R
7. Berriz GF
8. Roth FP
9. Gerszten RE
(2005) Metabolomic identification of novel biomarkers of myocardial ischemia
Circulation 112:3868–3875.

https://doi.org/10.1161/CIRCULATIONAHA.105.569137
- PubMed
- Google Scholar
1. Sadler MC
2. Auwerx C
3. Porcu E
4. Kutalik Z
(2022) Quantifying the role of transcript levels in mediating DNA methylation effects on complex traits and diseases
Nat Commun 13:7559.

https://doi.org/10.1038/s41467-022-35196-3
- Google Scholar
1. Shah SH
2. Bain JR
3. Muehlbauer MJ
4. Stevens RD
5. Crosslin DR
6. Haynes C
7. Dungan J
8. Newby LK
9. Hauser ER
10. Ginsburg GS
11. Newgard CB
12. Kraus WE
(2010) Association of a peripheral blood metabolic profile with coronary artery disease and risk of subsequent cardiovascular events
Circulation. Cardiovascular Genetics 3:207–214.

https://doi.org/10.1161/CIRCGENETICS.109.852814
- PubMed
- Google Scholar
1. Shin SY
2. Fauman EB
3. Petersen AK
4. Krumsiek J
5. Santos R
6. Huang J
7. Arnold M
8. Erte I
9. Forgetta V
10. Yang TP
11. Walter K
12. Menni C
13. Chen L
14. Vasquez L
15. Valdes AM
16. Hyde CL
17. Wang V
18. Ziemek D
19. Roberts P
20. Xi L
21. Grundberg E
22. Waldenberger M
23. Richards JB
24. Mohney RP
25. Milburn MV
26. John SL
27. Trimmer J
28. Theis FJ
29. Overington JP
30. Suhre K
31. Brosnan MJ
32. Gieger C
33. Kastenmüller G
34. Spector TD
35. Soranzo N
36. Multiple Tissue Human Expression Resource Consortium
(2014) An atlas of genetic influences on human blood metabolites
Nature Genetics 46:543–550.

https://doi.org/10.1038/ng.2982
- PubMed
- Google Scholar
(2022) Integrative analysis of metabolite GWAS illuminates the molecular basis of pleiotropy and genetic correlation
eLife 11:e79348.

https://doi.org/10.7554/eLife.79348
- PubMed
- Google Scholar
1. Sun BB
2. Maranville JC
3. Peters JE
4. Stacey D
5. Staley JR
6. Blackshaw J
7. Burgess S
8. Jiang T
9. Paige E
10. Surendran P
11. Oliver-Williams C
12. Kamat MA
13. Prins BP
14. Wilcox SK
15. Zimmerman ES
16. Chi A
17. Bansal N
18. Spain SL
19. Wood AM
20. Morrell NW
21. Bradley JR
22. Janjic N
23. Roberts DJ
24. Ouwehand WH
25. Todd JA
26. Soranzo N
27. Suhre K
28. Paul DS
29. Fox CS
30. Plenge RM
31. Danesh J
32. Runz H
33. Butterworth AS
(2018) Genomic atlas of the human plasma proteome
Nature 558:73–79.

https://doi.org/10.1038/s41586-018-0175-2
- PubMed
- Google Scholar
(2020) The membrane protein ANKH is crucial for bone mechanical performance by mediating cellular export of citrate and ATP
PLOS Genetics 16:e1008884.

https://doi.org/10.1371/journal.pgen.1008884
- PubMed
- Google Scholar
1. Taylor DL
2. Jackson AU
3. Narisu N
4. Hemani G
5. Erdos MR
6. Chines PS
7. Swift A
8. Idol J
9. Didion JP
10. Welch RP
11. Kinnunen L
12. Saramies J
13. Lakka TA
14. Laakso M
15. Tuomilehto J
16. Parker SCJ
17. Koistinen HA
18. Davey Smith G
19. Boehnke M
20. Scott LJ
21. Birney E
22. Collins FS
(2019) Integrative analysis of gene expression, DNA methylation, physiological traits, and genetic variation in human skeletal muscle
PNAS 116:10883–10888.

https://doi.org/10.1073/pnas.1814263116
- PubMed
- Google Scholar
1. Vandiedonck C
(2018) Genetic association of molecular traits: a help to identify causative variants in complex diseases
Clinical Genetics 93:520–532.

https://doi.org/10.1111/cge.13187
- PubMed
- Google Scholar
1. Verbanck M
2. Chen C-Y
3. Neale B
4. Do R
(2018) Detection of widespread horizontal pleiotropy in causal relationships inferred from Mendelian randomization between complex traits and diseases
Nature Genetics 50:693–698.

https://doi.org/10.1038/s41588-018-0099-7
- PubMed
- Google Scholar
1. Visscher PM
2. Wray NR
3. Zhang Q
4. Sklar P
5. McCarthy MI
6. Brown MA
7. Yang J
(2017) 10 years of GWAS discovery: biology, function, and translation
American Journal of Human Genetics 101:5–22.

https://doi.org/10.1016/j.ajhg.2017.06.005
- PubMed
- Google Scholar
1. Võsa U
2. Claringbould A
3. Westra H-J
4. Bonder MJ
5. Deelen P
6. Zeng B
7. Kirsten H
8. Saha A
9. Kreuzhuber R
10. Yazar S
11. Brugge H
12. Oelen R
13. de Vries DH
14. van der Wijst MGP
15. Kasela S
16. Pervjakova N
17. Alves I
18. Favé M-J
19. Agbessi M
20. Christiansen MW
21. Jansen R
22. Seppälä I
23. Tong L
24. Teumer A
25. Schramm K
26. Hemani G
27. Verlouw J
28. Yaghootkar H
29. Sönmez Flitman R
30. Brown A
31. Kukushkina V
32. Kalnapenkis A
33. Rüeger S
34. Porcu E
35. Kronberg J
36. Kettunen J
37. Lee B
38. Zhang F
39. Qi T
40. Hernandez JA
41. Arindrarto W
42. Beutner F
43. BIOS Consortium
44. i2QTL Consortium
45. Dmitrieva J
46. Elansary M
47. Fairfax BP
48. Georges M
49. Heijmans BT
50. Hewitt AW
51. Kähönen M
52. Kim Y
53. Knight JC
54. Kovacs P
55. Krohn K
56. Li S
57. Loeffler M
58. Marigorta UM
59. Mei H
60. Momozawa Y
61. Müller-Nurasyid M
62. Nauck M
63. Nivard MG
64. Penninx BWJH
65. Pritchard JK
66. Raitakari OT
67. Rotzschke O
68. Slagboom EP
69. Stehouwer CDA
70. Stumvoll M
71. Sullivan P
72. ’t Hoen PAC
73. Thiery J
74. Tönjes A
75. van Dongen J
76. van Iterson M
77. Veldink JH
78. Völker U
79. Warmerdam R
80. Wijmenga C
81. Swertz M
82. Andiappan A
83. Montgomery GW
84. Ripatti S
85. Perola M
86. Kutalik Z
87. Dermitzakis E
88. Bergmann S
89. Frayling T
90. van Meurs J
91. Prokisch H
92. Ahsan H
93. Pierce BL
94. Lehtimäki T
95. Boomsma DI
96. Psaty BM
97. Gharib SA
98. Awadalla P
99. Milani L
100. Ouwehand WH
101. Downes K
102. Stegle O
103. Battle A
104. Visscher PM
105. Yang J
106. Scholz M
107. Powell J
108. Gibson G
109. Esko T
110. Franke L
(2021) Large-Scale cis- and trans-eqtl analyses identify thousands of genetic loci and polygenic scores that regulate blood gene expression
Nature Genetics 53:1300–1310.

https://doi.org/10.1038/s41588-021-00913-z
- PubMed
- Google Scholar
1. Wang TJ
2. Larson MG
3. Vasan RS
4. Cheng S
5. Rhee EP
6. McCabe E
7. Lewis GD
8. Fox CS
9. Jacques PF
10. Fernandez C
11. O’Donnell CJ
12. Carr SA
13. Mootha VK
14. Florez JC
15. Souza A
16. Melander O
17. Clish CB
18. Gerszten RE
(2011) Metabolite profiles and the risk of developing diabetes
Nature Medicine 17:448–453.

https://doi.org/10.1038/nm.2307
- PubMed
- Google Scholar
1. Wang-Sattler R
2. Yu Z
3. Herder C
4. Messias AC
5. Floegel A
6. He Y
7. Heim K
8. Campillos M
9. Holzapfel C
10. Thorand B
11. Grallert H
12. Xu T
13. Bader E
14. Huth C
15. Mittelstrass K
16. Döring A
17. Meisinger C
18. Gieger C
19. Prehn C
20. Roemisch-Margl W
21. Carstensen M
22. Xie L
23. Yamanaka-Okumura H
24. Xing G
25. Ceglarek U
26. Thiery J
27. Giani G
28. Lickert H
29. Lin X
30. Li Y
31. Boeing H
32. Joost H-G
33. de Angelis MH
34. Rathmann W
35. Suhre K
36. Prokisch H
37. Peters A
38. Meitinger T
39. Roden M
40. Wichmann H-E
41. Pischon T
42. Adamski J
43. Illig T
(2012) Novel biomarkers for pre-diabetes identified by metabolomics
Molecular Systems Biology 8:615.

https://doi.org/10.1038/msb.2012.43
- PubMed
- Google Scholar
1. Wei Q
2. Xiao X
3. Fogle P
4. Dong Z
(2014) Changes in metabolic profiles during acute kidney injury and recovery following ischemia/reperfusion
PLOS ONE 9:e106647.

https://doi.org/10.1371/journal.pone.0106647
- PubMed
- Google Scholar
1. Williams CJ
(2016) The role of ANKH in pathologic mineralization of cartilage
Current Opinion in Rheumatology 28:145–151.

https://doi.org/10.1097/BOR.0000000000000247
- PubMed
- Google Scholar
1. Wishart DS
2. Guo A
3. Oler E
4. Wang F
5. Anjum A
6. Peters H
7. Dizon R
8. Sayeeda Z
9. Tian S
10. Lee BL
11. Berjanskii M
12. Mah R
13. Yamamoto M
14. Jovel J
15. Torres-Calzada C
16. Hiebert-Giesbrecht M
17. Lui VW
18. Varshavi D
19. Varshavi D
20. Allen D
21. Arndt D
22. Khetarpal N
23. Sivakumaran A
24. Harford K
25. Sanford S
26. Yee K
27. Cao X
28. Budinski Z
29. Liigand J
30. Zhang L
31. Zheng J
32. Mandal R
33. Karu N
34. Dambrova M
35. Schiöth HB
36. Greiner R
37. Gautam V
(2022) HMDB 5.0: the human metabolome database for 2022
Nucleic Acids Research 50:D622–D631.

https://doi.org/10.1093/nar/gkab1062
- PubMed
- Google Scholar
1. Wu Y
2. Zeng J
3. Zhang F
4. Zhu Z
5. Qi T
6. Zheng Z
7. Lloyd-Jones LR
8. Marioni RE
9. Martin NG
10. Montgomery GW
11. Deary IJ
12. Wray NR
13. Visscher PM
14. McRae AF
15. Yang J
(2018) Integrative analysis of omics summary data reveals putative mechanisms underlying complex traits
Nature Communications 9:918.

https://doi.org/10.1038/s41467-018-03371-0
- PubMed
- Google Scholar
1. Yin X
2. Bose D
3. Kwon A
4. Hanks SC
5. Jackson AU
6. Stringham HM
7. Welch R
8. Oravilahti A
9. Fernandes Silva L
10. FinnGen
11. Locke AE
12. Fuchsberger C
13. Service SK
14. Erdos MR
15. Bonnycastle LL
16. Kuusisto J
17. Stitziel NO
18. Hall IM
19. Morrison J
20. Ripatti S
21. Palotie A
22. Freimer NB
23. Collins FS
24. Mohlke KL
25. Scott LJ
26. Fauman EB
27. Burant C
28. Boehnke M
29. Laakso M
30. Wen X
(2022) Integrating transcriptomics, metabolomics, and GWAS helps reveal molecular mechanisms for metabolite levels and disease risk
American Journal of Human Genetics 109:1727–1741.

https://doi.org/10.1016/j.ajhg.2022.08.007
- PubMed
- Google Scholar
1. Zhu Z
2. Zhang F
3. Hu H
4. Bakshi A
5. Robinson MR
6. Powell JE
7. Montgomery GW
8. Goddard ME
9. Wray NR
10. Visscher PM
11. Yang J
(2016) Integration of summary data from GWAS and eQTL studies predicts complex trait gene targets
Nature Genetics 48:481–487.

https://doi.org/10.1038/ng.3538
- PubMed
- Google Scholar

Article and author information

Author details

Chiara Auwerx
1. Center for Integrative Genomics, University of Lausanne, Lausanne, Switzerland
2. Swiss Institute of Bioinformatics, Lausanne, Switzerland
3. University Center for Primary Care and Public Health, Lausanne, Switzerland
4. Department of Computational Biology, University of Lausanne, Lausanne, Switzerland
Contribution
Writing – original draft, Investigation, Visualization

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-3613-8450
Marie C Sadler
1. Swiss Institute of Bioinformatics, Lausanne, Switzerland
2. University Center for Primary Care and Public Health, Lausanne, Switzerland
3. Department of Computational Biology, University of Lausanne, Lausanne, Switzerland
Contribution
Formal analysis, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0002-2599-9207
Tristan Woh

Department of Computational Biology, University of Lausanne, Lausanne, Switzerland

Contribution
Formal analysis, Writing – review and editing

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-6916-0174
Alexandre Reymond

Center for Integrative Genomics, University of Lausanne, Lausanne, Switzerland

Contribution
Writing – review and editing

Competing interests
No competing interests declared
Zoltán Kutalik
1. Swiss Institute of Bioinformatics, Lausanne, Switzerland
2. University Center for Primary Care and Public Health, Lausanne, Switzerland
3. Department of Computational Biology, University of Lausanne, Lausanne, Switzerland
Contribution
Conceptualization, Writing – original draft

For correspondence
zoltan.kutalik@unil.ch

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0001-8285-7523
Eleonora Porcu
1. Center for Integrative Genomics, University of Lausanne, Lausanne, Switzerland
2. Swiss Institute of Bioinformatics, Lausanne, Switzerland
3. University Center for Primary Care and Public Health, Lausanne, Switzerland
Contribution
Conceptualization, Formal analysis, Writing – original draft

For correspondence
eleonora.porcu@unil.ch

Competing interests
No competing interests declared

"This ORCID iD identifies the author of this article:" 0000-0003-2878-7485

Funding

Swiss National Science Foundation (310030-189147)

Zoltán Kutalik

Swiss National Science Foundation (31003A_182632)

Alexandre Reymond

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

Computations were carried out on the high-performance cluster of the Lausanne University Hospital (CHUV). This work was supported by funding from the Department of Computational Biology (ZK) and the Center for Integrative Genomics (AR) from the University of Lausanne, as well as funding from the Swiss National Science Foundation (310030-189147 to ZK and 31003A_182632 to AR).

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.