There is increasing evidence that the human gut microbiome plays a role in immune function and metabolic disease^1,4,5. Manipulation of the gut microbiome offers an alternative to pharmacological interventions provided it can be demonstrated that altering microbiota composition and/or function (e.g. through personalized nutrition) has clinical benefit. To demonstrate this, it is essential to discriminate between microbiome features that are causal for disease, from those that are a consequence of disease or its treatment, and those that show statistical correlation due to confounding or pleiotropy.

Animal studies support a causal role for the gut microbiome in the development of type 2 diabetes (T2D), insulin resistance and obesity^6,7, but translating these findings to humans and identifying the specific bacterial species responsible has proven challenging⁸. Cross-sectional studies have confirmed that gut microbiota composition is altered in subjects with pre-diabetes or T2D compared to controls, while fecal transplantation studies have shown that insulin sensitivity increases in obese subjects with metabolic syndrome after the transfer of gut microbiota from lean donors^4,5,9,10. Whilst the specific microbiome features identified as responsible for these effects have differed between studies, one consistent finding in T2D subjects is a shift in microbiome composition away from species able to produce butyrate. Butyrate and other short-chain fatty acids (SCFAs), such as acetate and propionate, are produced by gut bacterial fermentation of undigested food components. Following absorption by the colonocytes, these SCFAs are either used locally as fuel for colonic mucosal epithelial cells or they enter the portal bloodstream¹¹. While the bulk of evidence suggests that increased SCFA production benefits the host by exerting anti-obesity and anti-diabetic effects^4,10,12–14, some in vitro and in vivo studies have indicated that over-production or accumulation of SCFAs in the bowel may also lead to obesity due to increased energy accumulation^15,16. Resolution of these conflicting data requires a detailed understanding of the causal relationships between gut microbiome composition, SCFA abundance and host energy metabolism.

Using a Mendelian Randomization (MR) approach³, we set out to identify if any bacterial species or pathways, i.e. sets of species grouped according to the specific functions they play in the gut, have a causal effect on metabolic traits. We and others have recently shown that it is possible to detect variants in the host genome that influence the composition of the gut microbiota^2,17,18. This allowed us to deploy a MR approach to infer causal relationships by asking whether genetic predictors of microbiome content influence metabolic traits—or the reverse. This formulation holds even though the quantitative contribution of host genetics to variation in microbiome composition may be limited¹⁹.

We assembled genome-wide genetic data, gut metagenomic-sequencing, measurements of fecal SCFAs, and clinical phenotypes for 952 normo-glycemic individuals from the LifeLines-DEEP (LL-DEEP) cohort. From consortium websites (GIANT, MAGIC and DIAGRAM, see URLs), we also gathered publically-available genome-wide association (GWA) summary statistics for 17 anthropometric and glycemic traits^20–27 (Supplementary Table 1). We focused our analyses on 245 microbiome features (2 fecal SCFA levels, 57 unique taxa, 186 pathways) that were, in LL-DEEP, correlated (false discovery rate (FDR) < 0.1) with at least one of the measured anthropometric and metabolic traits (Methods, Supplementary Table 2 and 3).

For each of these features, we sought genetic predictors -- independent genetic variants (r² ≤ 0.1), associated (P < 1 × 10⁻⁵) with the respective features – using GWAS data from LL-DEEP, reprocessed from our previous study² (Methods). The threshold P < 1 × 10⁻⁵ for variant inclusion was identified by maximizing the amount of genetic variance explained by the genetic predictors in 445 independent normo-glycemic individuals (the 500FG cohort)²⁸(Methods, Supplementary Figure 1), and designed to capture sets of variants likely to be enriched for association. On average, in LL-DEEP the identified genetic predictors explained 13% (range 2%−30%) of variance in their respective microbiome features. The average F-statistic, another measure of the strength of these genetic predictors, was 21.7 (range 15.3 – 25.5); an F-statistic >10 is considered sufficiently informative for MR analyses²⁹.

We used the inverse-variance weighted (IVW) test to identify causal relationships between the 245 microbiome features and the 17 traits of interest in a two-sample bidirectional MR analysis using pairs of GWAS summary statistics (one from a microbiome feature and one from a metabolic/anthropometric trait)²⁹. Based on principal component analysis (PCA) and cluster analyses conducted on the microbiome and metabolic and anthropometric traits (Methods, Supplementary Figure 2), we adopted a conservative multiple testing adjusted threshold of P < 1.3 × 10⁻⁴ to declare a causal relationship significant. Because the presence of horizontal pleiotropy (where a genetic predictor has independent effects on the diseases through multiple traits) could bias the MR estimates, we investigated the robustness of significant findings to pleiotropy by using three additional MR tests: the MR-PRESSO³⁰, the weighted median test³¹, and the MR-Egger³². We formally examined the presence of horizontal pleiotropy using the MR-PRESSO Global test³⁰ and the modified Rücker’s Q’ test^33,34. Finally, we sought to validate these causal relationships in an independent cohort (UK Biobank)³⁵ (Figure 1).

Open in a separate window

Figure 1.

Schematic representation of the study

Figure 1 is a schematic representation of our study, highlighting for each step the research question we want to answer, the analysis workflow, and the data used. We first aimed to identify which microbiome feature (taxa, microbiome pathway or short-chain fatty acid (SCFA)) correlated with metabolic traits in the LifeLines-DEEP (LL-DEEP) cohort (Step 1). We then performed genome-wide association (GWA) analysis in LL-DEEP to identify genetic predictors of those microbiome features (Step 2), and used the genetic predictors to estimate causal relationships through bidirectional Mendelian Randomization analysis and effect sizes for metabolic traits extracted from summary statistics of large GWA studies (Step 3). Finally, we validated our causality results using the UK Biobank (Step 4).

We observed a significant causal influence for one specific microbiome feature, a microbial pathway involved in 4-aminobutanoate (GABA) degradation (MetaCyc designation PWY-5022: 4-aminobutanoate degradation V) on increased insulin secretion, and specifically the ratio of the areas under the curve for insulin and glucose, AUC_insulin/AUC_glucose, measured during an oral glucose tolerance test (oGTT) (Figure 2a). Using nine genetic predictors (variance explained = 16%; F-statistic = 21, Supplementary Table 4), we estimated that each standard deviation (SD) increase in the abundance of PWY-5022 generated a 0.16 mU/mmol increase in AUC_insulin/AUC_glucose (P = 9.8 × 10⁻⁵) (Supplementary Table 5, Supplementary Figure 3). This causal relationship was robust when additional MR tests were performed (P_MR-PRESSO = 0.02, P_{Weighted-Median} = 0.02 and P_MR-EGGER = 0.02), and there was no evidence for horizontal pleiotropy (P_{MR-PRESSOglobal} = 0.18 and P_{RückerQ’(modified)} = 0.77) (Supplementary Figure 4). The reverse MR analysis (testing the relationship between genetic predictors of AUC_Insulin/AUC_glucose and PWY-5022 abundance) was not significant (P > 0.1, Supplementary Table 6). There was no evidence of causality with seven metabolic and anthropometric traits (body-mass index (BMI), body fat %, waist-hip ratio (WHR), visceral adipose tissue, abdominal subcutaneous adipose tissue, obesity and T2D) in a MR analyses that used UK Biobank summary statistics (Supplementary Table 7); insulin secretion phenotypes after oGTT were not available. We also found supportive evidence (P < 0.05) for a causal impact of this pathway on other insulin response parameters (Figure 2b). Though other types of causal relationship are possible, these data are consistent with a model whereby host genetic variation that influences gut microbiome composition so as to modulate GABA degradation activity results in improvements in the capacity of the pancreatic islets to secrete insulin in response to a physiological glucose challenge.

Open in a separate window

Figure 2.

Causal effect of butyrate-producing activity of the gut on glucose-stimulated insulin response

a) Schematic representation of the Mendelian Randomization analysis results: genetic predisposition to higher abundance of butyrate-producing microbiome pathway PWY-5022 (4-aminobutanoate degradation V pathway) is associated with insulin response after glucose challenge. The causal effect of PWY-5022 was also seen on other insulin response parameters, and the forest plot in panel (b) represents the magnitude of the effect on each parameter per one-standard-deviation increase in pathway abundance, as estimated in the inverse-variance weighted Mendelian Randomization (MR) analysis. MR analysis was carried out using up to nine genetic predictors and their effect size from LL-DEEP (952 samples) and MAGIC summary statistics (trait specific sample sizes are: AUC_insulin/AUC_glucose = 4213; insulin at 30 min = 4,409; AUC_insulin = 4,324; correct insulin response = 4,789; insulin increase at 30 min = 4,447; Disposition index = 5,130) (Methods, Supplementary Table 4 and 5). Corresponding two-sided P-values are given in the annexed table. (c) Correlation plots with PWY-5022 abundance and the bacteria correlating the most with it in 950 LL-DEEP samples (subset of the 952 normo-glycemic samples for which presence of those bacteria was detected). The Spearman correlation coefficient ρ is given in blue in each panel.

Butyrate and acetate are products of GABA degradation. In our taxonomic analyses, the bacterial species most correlated with the abundance of PWY-5022 were Eubacterium rectale and Roseburia intestinalis (Spearman ρ = 0.52 and 0.30, respectively, Figure 2c), both well-known butyrate-producing bacteria^36,37. Plasma butyrate levels were not measured in our study; current assays are challenging to perform and provide unreliable estimates³⁸. Whilst we consider the abundance of the PWY-5022 pathway to act as a proxy for butyrate production in the gut, we were unable to directly link PWY-5022 abundance to the amount of butyrate absorbed by the host. The abundance of PWY-5022 was poorly correlated with fecal butyrate levels (Spearman ρ = 0.1), and we did not detect any causal relationships between this SCFA and the 17 traits (P > 0.05), indicating that fecal butyrate is a poor proxy for butyrate production and absorption.

These results point towards a causal role for gut-produced butyrate that is focused on the dynamic insulin response to food ingestion, rather than the homeostatic mechanisms involved in the maintenance of glucose metabolism in the fasted state. Independent clinical studies support this hypothesis. For example, an intervention study evaluating the role of Bifidobacteria-increasing prebiotics (fructo-oligosaccharides) in 35 healthy individuals showed that prebiotics decrease levels of butyrate-producing bacteria and result in adverse effect on glucose metabolism following an oGTT³⁹.

The PWY-5022 finding led us to consider the role of other SCFAs in metabolic and anthropometric traits. In our cross-sectional analysis within LL-DEEP, we had detected associations between fecal propionate levels and BMI (FDR < 0.1). Propionate is produced by different bacteria than those producing butyrate⁴⁰, and its three genetic predictors (variance explained = 6.3%, F-statistic = 21) were independent of those implicated in PWY-5022 abundance (Supplementary Table 4, 8). In MR analyses for the 17 traits of interest, we found that each SD increase in fecal propionate levels was causally associated with a 0.03 SD increase in BMI (P = 0.0068) and an odd’s ratio (OR) = 1.15 for T2D (P = 0.004) (Supplementary Table 9), although these did not pass our significance threshold. No associations were evident in the reverse MR analysis testing the effect of T2D and BMI on fecal propionate levels (P > 0.1, Supplementary Table 10).

Of the two observed effects of fecal propionate, on BMI and T2D, the latter was more robust. The causal relationship on increased T2D risk was robust when other MR tests were performed (P_MR-PRESSO = 0.03, P_{Weighted-Median} = 0.03) and there was no evidence for pleiotropy (P_{MR-PRESSOglobal} 0.75, P_{RückerQ’(modified)} = 0.50) (Supplementary Figure 5). By contrast, the effect of propionate on increased BMI was not significant when using other MR tests and there was also evidence for pleiotropy (P_{MR-PRESSOglobal} = 2.0 × 10⁻³, P_{RückerQ’(modified)} = 9.2 × 10⁻⁴)(Supplementary Table 9, Supplementary Figure 6). The pleiotropy in the BMI effect could be accounted for by SNP rs7142308 (NC_000014.8:g.79482379A>G) (P_{MR-PRESSOOutlierTest} = 0.01), located within a BMI-associated locus²⁰ but independent of the lead variant (rs7141420 (NC_000014.8:g.79899454C>T), r² = 0.01 with rs7142308 in 1000 Genomes Europeans).

Applying MR analyses to UK Biobank summary statistics, we replicated the relationship between fecal propionate levels and increased risk of T2D (P_IVW = 0.01, P_MR-PRESSO = 0.007, P_{Weighted-Median} = 0.04; P_IVWcombined = 4 × 10⁻⁵, Figure 3), and there was no evidence of pleiotropy (P_{MR-PRESSOglobal} = 0.97, P_{RückerQ’(modified)} = 0.99). The relationship between fecal propionate and BMI was again not robust to pleiotropy, highlighting the need for caution in interpreting this effect as causal (Supplementary Table 11).

Open in a separate window

Figure 3.

Causal effect of fecal propionate on type 2 diabetes (T2D)

a) Schematic representation of the Mendelian Randomization analysis results: genetic predisposition to higher fecal propionate levels is associated with increased risk of T2D. b) A forest plot depicts the magnitude of the causal effect on T2D per each one-standard-deviation increase in fecal butyrate levels, as estimated by the inverse-variance weighted Mendelian Randomization (MR) analysis. MR analysis was carried out using the three genetic predictors derived in LifeLines-DEEP (LL-DEEP) and their effects in the discovery data set (DIAGRAM; 26,676 T2D cases and 132,532 controls) and in the replication cohort (UK Biobank; 19,119 T2D cases and 423,698 controls). The effect derived combining the two causal effects (from discovery and replication) with an inverse-variance weighted meta-analysis approach is also given. Corresponding two-sided P-values are listed in the annexed table. OR, odd’s ratio.

Over 95% of gut-produced SCFAs are absorbed by the host⁴¹, such that increases in fecal propionate levels could be the consequence of either increased production or reduced absorption. The latter (which would link increased fecal propionate to reduced circulating levels) would be more consistent with the preponderance of evidence that indicates that SCFAs have a largely beneficial effect on energy balance and metabolic homeostasis^4,10,12–14. As with plasma butyrate, plasma propionate levels were not measured in our cohorts. Further studies are warranted to explore the mechanisms underlying this relationship between fecal propionate levels and T2D.

In summary, these data are consistent with a causal role for gut-produced SCFAs, specifically butyrate and propionate, with respect to energy balance and glucose homeostasis in man. We have shown that a genetically-influenced shift in the gut microbiome towards increased production of butyrate has beneficial effects on beta-cell function, though no impact on T2D risk could be detected. We have also demonstrated that host genetic variation that results in increased fecal propionate levels (reflecting some combination of increased production or impaired absorption) has impact on T2D-risk.

Although the LL-DEEP cohort represents the largest population study on the genetics of microbiome^2,17,18, it is still underpowered to capture the limited genetic component that has been estimated for microbiome features¹⁹. The results from this and other microbiome GWAS^2,17,18 showed only limited direct overlap, highlighting the need for standardized protocols for data analyses and larger sample size⁴². This will be crucial also in the context of MR analyses, as expanded GWAS will deliver more robust genetic predictors⁴³. A better understanding of the complex interplay between gut microbiome and host metabolism will require expansion of current analyses and the ability to fold in measures of circulating SCFAs. Nevertheless, this study demonstrates that microbiome GWAS provide a route to causal inference that can guide and complement more direct experimental approaches, such as those based around fecal transplantation and animal models. We envisage that with expanded microbiome-genetic studies, for example the MiBioGen consortium⁴⁴, MR will become a standard tool for systematically screening a large number of hypotheses generated in current and future microbiome-wide association studies.