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

Curation statements for this article:
  • Curated by eLife

    eLife logo

    eLife assessment

    Auwerx and colleagues take a new approach to mine large datasets of the intermediary molecular data between GWAS and phenotype, touncover molecular mechanisms that lead from a GWAS hit to a phenotypic effect. The approach should be of great use to all (human) geneticists. Revisions are necessary to ensure that the significant findings from this approach are understood by the bioinformatic community and that these methods can be applied generally, given that the paper's main novelty is in its approach to mine large datasets, rather than a specific, key molecular finding.

This article has been Reviewed by the following groups

Read the full article See related articles

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.

Article activity feed

  1. Author Response

    Reviewer #1 (Public Review):

    Auwerx et al. have taken a new approach to mine large existing datasets of intermediary molecular data between GWAS and phenotype, with the aim of uncovering novel insight into the molecular mechanisms which lead a GWAS hit to have a phenotypic effect. The authors show that you can get additional insight by integrating multiple omics layers rather than analyzing only a single molecular type, including a handful of specific examples, e.g. that the effect of SNPs in ANKH on calcium are mediated by citrate. Such additional data is necessary because, as the authors' point out, while we have thousands of SNPs with significant impact on phenotypes of interest, we often don't know at all the mechanism, given that the majority of significant SNPs found through GWAS are in non-coding (and often intergenic) regions.

    This paper shows how one can mine large existing datasets to better estimate the cellular mechanism of significant, causal SNPs, and the authors have proven that by providing insight into the links between a couple of genes (e.g. FADS2, TMEM258) and metabolite QTLs and consequent phenotypes. There is definitely a need and utility for this, given how few significant SNPs (and even fewer recently-discovered ones) hit parts of the DNA where the causal mechanism is immediately obvious and easily testable through traditional molecular approaches.

    I find the paper interesting and it provides useful insight into a still relatively new approach. However, I would be interested in knowing how well this approach scales to the general genetics community: would this method work with a much smaller N (e.g. n = 500)? Being able to make new insights using cohorts of nearly 10,000 patients is great, but the vast majority of molecular studies are at least an order of magnitude smaller. While sequencing and mass spectrometry are becoming exponentially cheaper, the issue of sample size is likely to remain for the foreseeable future due to the challenges and expenses of the initial sample collection.

    We thank the reviewer for his assessment and have now addressed – in the revised version of the manuscript, as well as in the below point-by-point reply – his specific comments/questions.

    Reviewer #2 (Public Review):

    Auwerx et al. present a framework for the integration of results from expression quantitative trait loci (eQTL), metabolite QTL (mQTL) and genome-wide association (GWA) studies based on the use of summary statistics and Mendelian Randomization (MR). The aim of their study is to provide the field with a method that allows for the detection of causal relationships between transcript levels and phenotypes by integrating information about the effect of transcripts on metabolites and the downstream effect of these metabolites on phenotypes reported by GWA studies. The method requires the mapping of identical SNPs in disconnected mQTL and eQTL studies, which allows MRbased inference of a causal effect from a transcript to a metabolite. The effect of both transcripts and metabolites on phenotypes is evaluated in the same MR-based manner by overlaying eQTL and mQTL SNPs with SNPs present in phenotypic GWA studies.

    The aim of the presented approach is two-fold: (1) to allow identification of additional causal relationships between transcript levels and phenotypes as compared to an approach limited to the evaluation of transcript-to-phenotype associations (transcriptome-wide MR, TWMR) and (2) to provide information about the mechanism of effects originating from causally linked transcripts via the metabolite layer to a phenotype.

    The study is presented in a very clear and concise way. In the part based on empirical study results, the approach leads to the identification of a set of potential causal triplets between transcripts, metabolites and phenotypes. Several examples of such causal links are presented, which are in agreement with literature but also contain testable hypotheses about novel functional relationships. The simulation study is well documented and addresses an important question pertaining to the approach taken: Does the integration of mQTL data at the level of a mediator allow for higher power to detect causal transcript to phenotype associations?

    We thank the reviewer for his/her assessment and have now addressed – in the revised version of the manuscript, as well as in the below point-by-point reply – his/her specific comments/questions.

    Major Concerns

    1. Our most salient concern regarding the presented approach is the presence of multiple testing problems. In the analysis of empirical datasets (p. 4), the rational for setting FDR thresholds is not clearly stated. While this appears to be a Bonferroni-type correction (p-value threshold divided by number of transcripts or metabolites tested), the thresholds do not reflect the actual number of tests performed (7883 transcripts times 453 metabolites for transcript-metabolite associations, 87 metabolites or 10435 transcripts times 28 complex phenotypes). The correct and more stringent thresholds certainly decrease the overlap between causal relationships and thus reduce the identifiable number of causal triplets. Furthermore, we believe that multiple testing has to be considered for correct interpretation of the power analysis. The study compares the power of a TWMR-only approach to the power of mediation-based MR by comparing "power(TP)" against "power(TM) * power(MP)" (p. 12). This comparison is useful in a hypothetical situation given data on a single transcript affecting a single phenotype, and with potential mediation via a single metabolite. However, in an actual empirical situation, the number of non-causal transcript-metabolite-phenotype triplets will exceed the number of non-causal transcript-phenotype associations due to the multiplication with the number of metabolites that have to be evaluated. This creates a tremendous burden of multiple testing, which will very likely outweigh the increase in power afforded by the mediation-based approach in the hypothetical "single transcript-metabolite-phenotype" situation described here. Thus, for explorative detection of causal transcript-phenotype relationships, the TWMR-only method might even outperform the mediation-based method described by the authors, simply because the former requires a smaller number of hypotheses to be tested compared to the latter. The presented simulation would only hold in cases where a single path of causality with a known potential mediator is to be tested.

    We thank the reviewer for pointing out the multiple testing issue. Based on this comment, we have revised our approach by mainly implementing two major modifications to our approach.

    First, we reduce the number of assessed metabolites to 242 compounds for which we were able to identify a Human Metabolome Database (HMDB) identifier through manual curation. This was triggered by the suggestion of reviewer #1 to facilitate the database/literature-based follow-up of our discoveries. The motivation is to only test metabolites that if found to be significantly associated would yield interpretable results, thereby reducing the number of tests to be performed. This modification is described in the revised manuscript:

    Results: “Summary statistics for cis-eQTLs stem from the eQTLGen Consortium metaanalysis of 19,942 transcripts in 31,684 individuals [3], while summary statistics for mQTLs originate from a meta-analysis of 453 metabolites in 7,824 individuals from two independent European cohorts: TwinsUK (N = 6,056) and KORA (N = 1,768) [6]. After selecting SNPs included in both the eQTL and mQTL studies, our analysis was restricted to 7,884 transcripts with ≥ 3 instrumental variables (IVs) (see Methods, Supplemental Figure 1) and 242 metabolites with an identifier in The Human Metabolome Database (HMDB) [28] (see Methods, Supplemental Table 1).”

    Methods: “mQTL data originate from Shin et al. [6], which used ultra-high performance liquid chromatography-tandem mass spectrometry (UPLC-MS/MS) to measure 486 whole blood metabolites in 7,824 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 Human Metabolome Database (HMDB) identifiers (Supplemental Table 1) and used in this study.”

    Second, to account for all remaining tests, we now select significant causal effects based on FDR < 5% in all performed univariable MR analyses. With 5% FDR on both the transcript-to-metabolite and metabolite-to-phenotype effects, the FDR for triplets is slightly inflated to 9.75% (= 1-0.952), a consideration that we now explicitly describe. Note that selecting triplets based on transcript-tometabolite and metabolite-to-phenotype effects FDR < 2.5%, result in a FDR < 5% (1-0.9752) for the triplets. This more stringent threshold identifies 135 causal triplets, 39 of which would be missed by TWMR. Overall, Results and Supplemental Tables have been updated and now read as follow:

    “Mapping the transcriptome onto the metabolome […] 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 transcriptmetabolite associations (FDR 5% considering all 1,907,690 instrumentable gene-metabolite pairs Supplemental Table 2) […].

    Mapping the metabolome onto complex phenotypes […] Overall, 34 metabolites were associated with at least one phenotype (FDR 5% considering all 1,344 metabolite-phenotype pairs), resulting in 132 unique metabolitephenotype associations (Supplemental Table 4).

    Mapping the transcriptome onto complex phenotypes […] In total, 5,140 transcripts associated with at least one phenotype (FDR 5% considering all 292,170 gene-phenotype pairs) resulting in 13,141 unique transcript-phenotype associations (Supplemental Table 5).

    Mapping metabolome-mediated effects of the transcriptome onto complex phenotypes […] 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.952 = 9.75%) (Supplemental Table 6).”

    In the simulations performed for the power analysis, we used a Bonferroni correction. We ran each simulation for 500 transcripts, measuring 80 metabolites at each run and performed TWMR and MWMR. The power of TWMR was calculated by counting how many times we obtain p-values ≤ 0.05/500. The power of the mediation analysis was calculated as 𝑝𝑜𝑤𝑒𝑟"$ ∗ 𝑝𝑜𝑤𝑒𝑟$#, where 𝑝𝑜𝑤𝑒𝑟"$ was calculated by counting how many times we obtain p-values ≤ 0.05/(500*80), and 𝑝𝑜𝑤𝑒𝑟$# was calculated by counting how many times we obtain p-values ≤ 0.05/80. In the revised manuscript, we additionally repeated each simulated scenario 10 times to increase robustness of results. This has been clarified in both the Methods and Results sections of the revised manuscript:

    Methods: “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 𝛽7<"=, 𝛽4<"= and 𝛽>?,(. For each MR analysis we calculated the power to detect a significant association as well as the difference in power between TWMR and the mediation analyses (i.e., 𝑝𝑜𝑤𝑒𝑟"# − 𝑝𝑜𝑤𝑒𝑟"$ ∗ 𝑝𝑜𝑤𝑒𝑟$#). Each specific scenario was repeated 10 times and the average difference in power across simulation was plotted as a heatmap.”

    Results: “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: a. the proportion (𝜌) of direct (𝛼!) to total (𝛼"#) 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); b. the ratio (𝜎) between the transcript-to-metabolite (𝛼"$) and the metabolite-to-phenotype (𝛼$#) effects, exploring the range from 0.1 to 10 (51 values).
    Transcripts were simulated with 6% heritability (i.e., median ℎ@ in the eQTLGen data) and a causal effect of 0.035 (i.e., ~65% of power in TWMR at a = 0.05) on a phenotype. Each scenario was simulated 10 times and results were averaged to assess the mean difference in power (see Methods).”

    1. A second concern regards the interpretation of the results based on the empirical datasets. For the identified 206 transcript-metabolite-phenotype causal triplets, the authors show a comparison between TWMR-based total effect of transcripts on phenotypes and the calculated direct effect based on a multivariable MR (MVMR) test (Figure 2B), which corrects for the indirect effect mediated by the metabolite in the causal triplet. The comparison shows a strong correlation between direct and total effect. A thorough discussion of the potential reasons for deviation (in both negative and positive directions) from the identity line is missing.

    Deviation from the identity line, as observed in Figure 2B, indicates that while there is a strong correlation between direct and total effect, it is not perfect, and part of the total effect is due to an indirect effect mediated by metabolites. This is explained and discussed in the Results and Discussion section:

    Results: “Regressing direct effects (𝛼!) on total effects (𝛼"#) on (Figure 2A), we estimated that for our 216 mediated associations, 77% [95% CI: 70%-85%] of the transcript effect on the phenotype was direct and thus not mediated by the metabolites (Figure 2B).”

    Discussion: “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 such 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. […] 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.”

    Furthermore, no test of significance for potential cases of mediation is presented. Due to the issues of multiple testing discussed above, the significance of the inferred cases of mediation is drawn into question. The examples presented for causal triplets (involving the ANKH and SLC6A12 transcripts) feature transcripts with low total effects and a small ratio between direct and total effect, in line with the power analysis. However, in these examples, the total effects are also quite low. Its significance has to be tested with an appropriate statistical test, incorporating multiple testing correction.

    Following the reviewer’s suggestion, we have modified our criteria to call significant associations to account for multiple testing (see extensive reply to major concern #1). With 5% FDR on both the transcript-to-metabolite and metabolite-to-phenotype effects, the FDR for triplets is slightly inflated to 9.75% (= 1-0.952). We mention this limitation in the revised manuscript:

    “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.952 = 9.75%) (Supplemental Table 6).”

    All examples presented in the original manuscript remained significant. The fact that the total effect in these examples is low makes them particularly interesting as it highlights how our approach can detect biologically plausible associations between a transcript and a phenotype that only show mild evidence through TWMR but are strongly supported when accounting for metabolites that mediate the transcript-phenotype relation, showcasing situations in which our method can provide a true advantage over classical approaches such as TWMR. Such examples may emerge due to opposite signed direct and indirect effects, which cancel each other out when it comes to testing total effects. What is key that we do not claim the total and the mediated effects to be different (as we would have very limited power to do so), but simply point out that under certain settings we are better powered to detect mediated effects than total ones. In the ANKH example (more details below), the total ANKH-calcium effect is almost exactly the same as the product of the 𝛼,-.%→056157 and 𝛼056157→0120*34 effects, simply the latter ones are detectable, while the total effect is not.

    In the revised manuscript the case for our selected examples is made even stronger thanks to an analysis proposed by Reviewer #1 that aimed at estimating the proportion of previously reported associations through automated literature review. For instance, while our literature review found previously reported evidence of the ANKH-calcium link and of the ANKH-citrate link, we did not identify any publication mentioning all 3 terms in combination in the abstract and/or title, illustrating how our approach can establish bridges between knowledge gaps. We revised the Results section describing the ANKH example accordingly:

    “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 (Supplemental Table 6). While the TWMR effect of ANKH expression on calcium levels was not significant (𝛼,-.%→012034 = −0.02; 𝑃 = 0.03), we observed that ANKH expression decreased citrate levels (𝛼,-.%→056157 = −0.30; 𝑃 = 2.2 × 1089:), which itself increased serum calcium levels (𝛼056157→012034 = 0.07; 𝑃 = 6.5 × 108;9). Mutations in ANKH have been associated with several rare mineralization disorders [MIM: 123000, 118600] [32] due to the gene encoding a transmembrane protein that channels inorganic pyrophosphate to the extracellular matrix, where at low concentrations it inhibits mineralization [33]. Recently, a study proposed that ANKH instead exports ATP to the extracellular space (which is then rapidly converted to inorganic pyrophosphate), along with citrate [34]. 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 [35]. 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.”

    Furthermore, the analysis of the empirical data indicates that the ratio between direct and indirect effect of a transcript on a phenotype is in most cases close to identity, except for triplets with low total effects. This fact should be considered in the power analysis, which assigned the highest gain in power by the mediation analysis to cases of low direct to total effect ratio. The empirical data indicate that these cases might be rare or of minor relevance for the tested phenotypes.

    As our previous power analyses did not fully reflect scenarios observed from empirical data, we extended the range of covered 𝜌 (i.e., the ratio between direct and total effect), so that it mimics more closely the observed range of 𝜌. In the revised manuscript, 𝜌 varies from -2 to 2, so that we also consider configurations where direct and total effects have opposite direction. To provide the readers with a rough idea how frequent the different parameter combinations occur in real data, we now provide another heatmap indicating the density of detected associations in those parameter regimes as Supplemental Figure 4.

    This map can be brought in perspective of Figure 4A that illustrates the power of TWMR vs. mediation analysis over the same range of parameter settings.

    It becomes apparent from Supplemental Figure 4 that in real data, 𝜎 is always larger than 1 and often exceeds 10. Note, however, that this heatmap must be interpreted with care, since the “detected” density will be low in regions where both methods have low power.

    1. Related to the interpretation of causal links: horizontal pleiotropy needs to be considered. The authors report the identification of causal links between TMEM258, FADS1 and FADS2, arachidonic acid-derived lipids and complex phenotypes. However, they also mention the high degree of pleiotropy due to linkage disequilibrium at the underlying eQTL and mQTL region as well as the network of over 50 complex lipids known to be associated with the expression of the above transcripts. Thus, it seems possible that the levels of undetected lipid species may be more important for the phenotypic effect of variation in these transcripts and that the reported "mediators" are rather covariates. Such horizontal pleiotropy would violate a basic assumption of the MR approach. While we think that this does not invalidate the approach altogether, it does affect the interpretation of specific metabolites as mediators. This is aggravated by the fact that metabolic networks are more tightly interconnected than macromolecular interaction networks (assortative nature of metabolic networks) and that single point-measurements of metabolites may not be generally informative about the flux through a specific metabolic pathway.

    This is a valid point and we discuss this limitation in the revised Discussion:

    “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.”

  2. eLife assessment

    Auwerx and colleagues take a new approach to mine large datasets of the intermediary molecular data between GWAS and phenotype, touncover molecular mechanisms that lead from a GWAS hit to a phenotypic effect. The approach should be of great use to all (human) geneticists. Revisions are necessary to ensure that the significant findings from this approach are understood by the bioinformatic community and that these methods can be applied generally, given that the paper's main novelty is in its approach to mine large datasets, rather than a specific, key molecular finding.

  3. Reviewer #1 (Public Review):

    Auwerx et al. have taken a new approach to mine large existing datasets of intermediary molecular data between GWAS and phenotype, with the aim of uncovering novel insight into the molecular mechanisms which lead a GWAS hit to have a phenotypic effect. The authors show that you can get additional insight by integrating multiple omics layers rather than analyzing only a single molecular type, including a handful of specific examples, e.g. that the effect of SNPs in ANKH on calcium are mediated by citrate. Such additional data is necessary because, as the authors' point out, while we have thousands of SNPs with significant impact on phenotypes of interest, we often don't know at all the mechanism, given that the majority of significant SNPs found through GWAS are in non-coding (and often intergenic) regions.

    This paper shows how one can mine large existing datasets to better estimate the cellular mechanism of significant, causal SNPs, and the authors have proven that by providing insight into the links between a couple of genes (e.g. FADS2, TMEM258) and metabolite QTLs and consequent phenotypes. There is definitely a need and utility for this, given how few significant SNPs (and even fewer recently-discovered ones) hit parts of the DNA where the causal mechanism is immediately obvious and easily testable through traditional molecular approaches.

    I find the paper interesting and it provides useful insight into a still relatively new approach. However, I would be interested in knowing how well this approach scales to the general genetics community: would this method work with a much smaller N (e.g. n = 500)? Being able to make new insights using cohorts of nearly 10,000 patients is great, but the vast majority of molecular studies are at least an order of magnitude smaller. While sequencing and mass spectrometry are becoming exponentially cheaper, the issue of sample size is likely to remain for the foreseeable future due to the challenges and expenses of the initial sample collection.

  4. Reviewer #2 (Public Review):

    Auwerx et al. present a framework for the integration of results from expression quantitative trait loci (eQTL), metabolite QTL (mQTL) and genome-wide association (GWA) studies based on the use of summary statistics and Mendelian Randomization (MR). The aim of their study is to provide the field with a method that allows for the detection of causal relationships between transcript levels and phenotypes by integrating information about the effect of transcripts on metabolites and the downstream effect of these metabolites on phenotypes reported by GWA studies. The method requires the mapping of identical SNPs in disconnected mQTL and eQTL studies, which allows MR-based inference of a causal effect from a transcript to a metabolite. The effect of both transcripts and metabolites on phenotypes is evaluated in the same MR-based manner by overlaying eQTL and mQTL SNPs with SNPs present in phenotypic GWA studies.

    The aim of the presented approach is two-fold: (1) to allow identification of additional causal relationships between transcript levels and phenotypes as compared to an approach limited to the evaluation of transcript-to-phenotype associations (transcriptome-wide MR, TWMR) and (2) to provide information about the mechanism of effects originating from causally linked transcripts via the metabolite layer to a phenotype.

    The study is presented in a very clear and concise way. In the part based on empirical study results, the approach leads to the identification of a set of potential causal triplets between transcripts, metabolites and phenotypes. Several examples of such causal links are presented, which are in agreement with literature but also contain testable hypotheses about novel functional relationships. The simulation study is well documented and addresses an important question pertaining to the approach taken: Does the integration of mQTL data at the level of a mediator allow for higher power to detect causal transcript to phenotype associations?

    Major Concerns
    1. Our most salient concern regarding the presented approach is the presence of multiple testing problems. In the analysis of empirical datasets (p. 4), the rational for setting FDR thresholds is not clearly stated. While this appears to be a Bonferroni-type correction (p-value threshold divided by number of transcripts or metabolites tested), the thresholds do not reflect the actual number of tests performed (7883 transcripts times 453 metabolites for transcript-metabolite associations, 87 metabolites or 10435 transcripts times 28 complex phenotypes). The correct and more stringent thresholds certainly decrease the overlap between causal relationships and thus reduce the identifiable number of causal triplets. Furthermore, we believe that multiple testing has to be considered for correct interpretation of the power analysis. The study compares the power of a TWMR-only approach to the power of mediation-based MR by comparing "power(TP)" against "power(TM) * power(MP)" (p. 12). This comparison is useful in a hypothetical situation given data on a single transcript affecting a single phenotype, and with potential mediation via a single metabolite. However, in an actual empirical situation, the number of non-causal transcript-metabolite-phenotype triplets will exceed the number of non-causal transcript-phenotype associations due to the multiplication with the number of metabolites that have to be evaluated. This creates a tremendous burden of multiple testing, which will very likely outweigh the increase in power afforded by the mediation-based approach in the hypothetical "single transcript-metabolite-phenotype" situation described here. Thus, for explorative detection of causal transcript-phenotype relationships, the TWMR-only method might even outperform the mediation-based method described by the authors, simply because the former requires a smaller number of hypotheses to be tested compared to the latter. The presented simulation would only hold in cases where a single path of causality with a known potential mediator is to be tested.

    2. A second concern regards the interpretation of the results based on the empirical datasets. For the identified 206 transcript-metabolite-phenotype causal triplets, the authors show a comparison between TWMR-based total effect of transcripts on phenotypes and the calculated direct effect based on a multivariable MR (MVMR) test (Figure 2B), which corrects for the indirect effect mediated by the metabolite in the causal triplet. The comparison shows a strong correlation between direct and total effect. A thorough discussion of the potential reasons for deviation (in both negative and positive directions) from the identity line is missing. Furthermore, no test of significance for potential cases of mediation is presented. Due to the issues of multiple testing discussed above, the significance of the inferred cases of mediation is drawn into question. The examples presented for causal triplets (involving the ANKH and SLC6A12 transcripts) feature transcripts with low total effects and a small ratio between direct and total effect, in line with the power analysis. However, in these examples, the total effects are also quite low. Its significance has to be tested with an appropriate statistical test, incorporating multiple testing correction. Furthermore, the analysis of the empirical data indicates that the ratio between direct and indirect effect of a transcript on a phenotype is in most cases close to identity, except for triplets with low total effects. This fact should be considered in the power analysis, which assigned the highest gain in power by the mediation analysis to cases of low direct to total effect ratio. The empirical data indicate that these cases might be rare or of minor relevance for the tested phenotypes.

    3. Related to the interpretation of causal links: horizontal pleiotropy needs to be considered. The authors report the identification of causal links between TMEM258, FADS1 and FADS2, arachidonic acid-derived lipids and complex phenotypes. However, they also mention the high degree of pleiotropy due to linkage disequilibrium at the underlying eQTL and mQTL region as well as the network of over 50 complex lipids known to be associated with the expression of the above transcripts. Thus, it seems possible that the levels of undetected lipid species may be more important for the phenotypic effect of variation in these transcripts and that the reported "mediators" are rather covariates. Such horizontal pleiotropy would violate a basic assumption of the MR approach. While we think that this does not invalidate the approach altogether, it does affect the interpretation of specific metabolites as mediators. This is aggravated by the fact that metabolic networks are more tightly interconnected than macromolecular interaction networks (assortative nature of metabolic networks) and that single point-measurements of metabolites may not be generally informative about the flux through a specific metabolic pathway.