The association between body mass index and metabolite response to a liquid mixed meal challenge: a Mendelian randomization study

Background Metabolite abundance is a dynamic trait that varies in response to environmental stimuli and phenotypic traits, such as food consumption and body mass index (BMI, kg/m2). Objectives In this study, we used the Netherlands Epidemiology of Obesity (NEO) study data to identify observational and causal associations between BMI and metabolite response to a liquid meal. Methods A liquid meal challenge was performed, and Nightingale Health metabolite profiles were collected in 5744 NEO participants. Observational and one-sample Mendelian randomization (MR) analysis were conducted to estimate the effect of BMI on metabolites (n = 229) in the fasting, postprandial, and response (or change in abundance) states. Results We observed 473 associations with BMI (175 fasting, 188 postprandial, and 110 response) in observational analyses. In MR analyses, we observed 20 metabolite traits (5 fasting, 12 postprandial, and 3 response) to be associated with BMI. MR associations included the glucogenic amino acid alanine, which was inversely associated with BMI in the response state (β: −0.081; SE: 0.023; P = 5.91 × 10−4), suggesting that as alanine increased in postprandial abundance, that increase was attenuated with increasing BMI. Conclusions Overall, this study showed that MR estimates were strongly correlated with observational effect estimates, suggesting that the broad associations seen between BMI and metabolite variation has a causal underpinning. Specific effects in previously unassessed postprandial and response states are detected, and these may likely mark novel life course risk exposures driven by regular nutrition.


Annotation of Supplementary Tables
1. Supplementary Table 1 a. metabolite annotation 2. Supplementary Table 2 a. annotation of other study (co)variables 3. Supplementary Table 3 a.observational and MR summary statistics for alternative response traits 4. Supplementary Table 4 a. primary study results -observational and MR summary statistics for the wNEO framework 5. Supplementary Table 5 a. observational and MR summary statistics for the wNEO sensitivity (additional covariables or confounders) analyses 6. Supplementary Table 6 a. observational and MR summary statistics for the wNEO framework but with no transformation of the outcome (metabolite) traits 7. Supplementary Table 7 a.observational and MR summary statistics for NEO cohort but with no weights 8. Supplementary Table 8 a. observational and MR summary statistics for the Leiderdorp NEO sub-sample 9. Supplementary Table 9 a. observational and MR summary statistics for the Leiden NEO sub-sample 10.Supplementary Table 10 a. observational and MR summary statistics for females using the wNEO framework 11.Supplementary Table 11 a. observational and MR summary statistics for males using the wNEO framework 12. Supplementary Table 12 a.observational and MR summary statistics for females in the Leiderdorp sub-sample 13.Supplementary Table 13 a.observational and MR summary statistics for males in the Leiderdorp sub-sample 14.Supplementary Table 14 a. fasting -v-postprandial metabolite paired t-tests 15.Supplementary Table 15 a. associations and variance explained between study variables (Table 2) and BMI and the PGS

Metabolite data quality control
The initial NEO metabolite data set contained data for 5,744 individuals, 229 fasting metabolite, 229 postprandial metabolites, and 148 previously derived (38) ornls response metabolite traits.We used the R package metaboprep to perform quality control on the data prior to analysis.The following parameter values were used when running metaboprep: (1) feature missingness 0.2, (2) sample missingness 0.2, (3) total peak area standard deviation 5, (4) principal component (PC) outlier standard deviations 5, (5) tree cut height 0.5, (6) derived variable exclusion FALSE.Metaboprep performs an initial exclusion of samples with missingness greater than or equal to 80% (n = 9), followed by exclusion of features with missingness greater than equal to 80% (n = 0).After this initial filtering, the data set is filtered using the parameters defined above.This will filter (1) samples (n = 217) and then (2) features (n = 3) with greater than or equal to 20% missingness, (3) then filter samples with a total sum abundance (TSA) at complete features (no missingness) that fall beyond 5 standard deviations from the mean of the observed TSA distribution (n = 0), and then (4) filter samples (n = 0) that are 5 standard deviations from the mean of the first 'i' informative PCs, here defined by the Cattel's Scree test acceleration factor (i = 2).Parameter 5 defines the clustering dendrogram tree cut height which was used to estimate the number of clusters, and then used to identify the representative or principal variables in the data set.These principal variables were used in the data filter PC analysis used by parameter 4 above.Derived variables or features such as ratios that are derived from multiple features in the data set were (6) not excluded from any analyses.Given the defined metaboprep parameters 226 samples and 3 metabolites, all ornls response traits, were filtered from the data set.The metaboprep log file and metaboprep report, as generated by metaboprep, are available at the end of this Supplementary Materials document.
Following metaboprep quality control, additional quality control steps were taken.First, all zero values were turned into NAs.Second, for each metabolite (in the fasting and postprandial state, individually) any sample with a value 10 interquartile distances from the median was turned into NAs.This QC-step removes all extreme observations that are unrealistic given the empirical distribution of a single metabolite trait.Third, for each metabolite we used the expectation that fasting and postprandial data are correlated to estimate a delta value (postprandial -fasting), from which all samples that fell five interquartile distances from the median of that metabolite's delta distribution were turned into NAs -in both the fasting and postprandial data.This QC-step removes all extreme observations that are unrealistic given the observed, empirical, correlated nature of this bivariate data.See Supplementary Figure 3 for an example illustration of these QC steps.

Effective number of tested metabolites
There are 687 metabolite traits in the data set (fasting = 229, postprandial = 229, response = 229) but these traits are not strictly independent.First, given the presence of both fasting and postprandial data, and a response trait derived from the two, we would expect there to be a degree of dependency.Second, given the focus on lipids and lipoproteins on this platform we also expect an abundance of inter-correlated structure (Supplementary Figure 4).As such we used the R package iPVs (https://github.com/hughesevoanth/iPVs) to estimate the effective number of metabolites in the data set.The effective number represents an estimate of the number of representative or independent traits present in the data (47-49).The method implemented by iPVs provides an estimate by the iterative construction of a hierarchical clustering dendrogram followed by a tree cut to identify clusters of correlated variables and identification of a principal variable (PVs) for each cluster.The iPVs() function parameters were set to: "spearman" for the correlation matrix, "R" for the distance matrix, "complete" for the hierarchical cluster method, and 0.5 for the tree cut height.The distance matrix produced is equal to one minus the absolute Spearman's rho.At a tree cut height of 0.5, metabolites with a Spearman's rho greater than 0.5 cluster together.If the tree cut height was, for example, set to a value of 0.8, then metabolites with a Spearman's rho greater than 0.2 would cluster together.In total 43 clusters or representative variables in the NEO metabolite data set were identified.Resultantly our data reduced, study-wide Bonferroni (BF) corrected p-value was set to 0.05/43 or 1.163x10 -3 .

Identification of possible confounders
To identify possible confounders, or variables that may violate MR IV assumption number two (independence) we tested, in a univariate fashion, for an association between (a) sampling date and BMI and BMI-PGS, (b) study covariables and BMI, (c) study covariables and BMI-PGS and (d) study covariables and metabolite traits (Supplementary Table 15).An ANOVA analysis of the univariate linear models was then used to partition the sums of squares and estimate the variance explained, in the form of an eta-squared statistic (h 2 ).Eta-squared (h 2 ) is a measure of the variance explained and can be derived from sums of squares or deviances extracted from ANOVAs.In addition, an ANOVA F-test was used to estimate a p-value for each association.These analyses were carried out for the NEO cohort as a whole and for each of the two NEO sub-samples and provides a means to identify possible confounders that may invalidate the MR framework.
The association between BMI-PGS and principal component three is partially driven by its association with sub-population structure (h 2 = 0.032, P = 1.71x10 -41 , Supplementary Figure 13).All 17 covariables associated the BMI-PGS were also tested for an association with each metabolite trait in a univariable linear model.Results illustrated broad association across all traits and covariables, defining each covariable as a confounder in MR analysis (Supplementary Figure 10).

Sub-population analyses
All observational and MR analyses were repeated in each of the two sub-populations and in the NEO cohort without the inclusion of weights.The (i) randomly sampled Leiderdorp subpopulation (n=1406) provides a means to verify the point estimates derived from the primary weighted NEO (wNEO) results presented above.The (ii) unweighted NEO cohort provides an evaluation of the weights, and the (iii) Leiden sub-population sample (n=4111) provides an evaluation of effect estimates derived from a sampling biased for elevated BMI.
First, observational effect estimates from the wNEO analysis strongly correlate with those from the Leiderdorp sub-sample, with a Pearson's r of 0.983.The correlation does reduce when compared to the un-weighted NEO analysis (r = 0.96) and the biased Leiden sample (r = 0.91, Supplementary Figure 12).Second, MR effect estimates remain strongly correlated between the wNEO and Liederdorp sample (r = 0.855), but as the shift in sample population BMI increases the correlation reduces.When compared to the un-weighted NEO sample the Pearson's r is 0.445 (P = 1.10x10 -34 ) and even becomes negative when compared to the Leiden sample (Pearson's r = -0.107,P = 4.84x10 -03 ; Supplementary Figure 12).The comparison to the randomly sample Leiderdorp sample would suggest that the weights used in the wNEO analysis did not introduce any strong error in effect estimates.Yet, the decrease in congruency between analyses as the sample population mean BMI shifts would indicate that there are either un-accounted for confounders influencing the results or that the relationship between BMI and metabolite trait variation are not linear.

Sensitivity analyses
To evaluate the influence of confounders on MR effect estimates we reran the association analysis in the wNEO data set.We assumed that the association between BMI-PGS and hip circumference, weight, waist circumference, basal metabolic rate, waist-to-hip ratio, and mean CO 2 production were the product of biological similarities with our exposure of interest, BMI, and not confounders.We did include a smoking variable (packyears), a diet variable (on a weight loss diet), and PC3 as new covariates in the MR analysis.We also excluded the 101 samples with a peroxidation flag on their samples.Overall primary (wNEO) observational (Pearson's r = 0.998) and MR (Pearson's r = 0.963) effect estimates correlate strongly with those from the sensitivity analysis (Supplementary Table 5).Further no effect estimates differ between the primary and sensitivity models (z-test, P< 0.05) indicating that the inclusion of the additional confounders had minimal effect on the association statistics (Supplementary Figure 14 and 15).(Right) is a scatter plot for fasting (x-axis) and postprandaial (y-axis) total lipids in small HDL (S-HDL-L) after quality control.An an equivalency line (x = y) is illustrated as a dashed grey line.In addition, six linear models were fit to the data and are plotted.They are a (black) a Deming regression, (orange) a linear model, (blue) a linear model were the independent term or fasting data was fit as a quadratic term, (purple) a linear model where the fasting data was fit as a cubic term, (red) a median regression, and (green) a generalized additive model or GAM where the fasting data was fit as a smooth.Each of these models were fit to evaluate how much variability was present among each model type.

Supplementary
and response).Along the x-axis are the point or effect estimate and along the y-axis are the - Tiles with an effect estimate provided in text are those with a p-value smaller than 0.05.The lipoproteins are labeled and color corradiated the x-axis, and the component or ratio being measured is along the y-axis.Lower: A dot plot or profile of postprandial observational effect estimates for lipoproteins (x-axis) ordered by lipoprotein size or density is provided to illustrate the correlation between effect estimates (y-axis) within a lipoprotein and the structure of estimates between lipoproteins by size.The component or measurement of each lipoprotein are defined by the color as described in the key.All effect estimates can be found in Supplementary Table 4.

Supplementary Figure 18: Fasting and Postprandial MR profiles. (A) Fasting MR and (B)
Postprandial MR.Upper: A tile plot of MR effect estimates for lipoproteins and lipoprotein ratios in the fasting state.Tiles with an effect estimate provided in text are those with a pvalue smaller than 0.05.The lipoproteins are labeled and color corradiated the x-axis, and the component or ratio being measured is along the y-axis.Lower: A dot plot or profile of fasting MR effect estimates for lipoproteins (x-axis) ordered by lipoprotein size or density is provided to illustrate the correlation between effect estimates (y-axis) within a lipoprotein and the structure of estimates between lipoproteins by size.The component or measurement of each lipoprotein are defined by the color as described in the key.All effect estimates can be found in Supplementary Table 4.The association between body mass index and metabolite response to a liquid mixed meal challenge: a Mendelian randomization study; Hughes et al

BMI-PGS by sampling month
The association between body mass index and metabolite response to a liquid mixed meal challenge: a Mendelian randomization study; Hughes et al Supplementary Figure 3 The association between body mass index and metabolite response to a liquid mixed meal challenge: a Mendelian randomization study; Hughes et al      The association between body mass index and metabolite response to a liquid mixed meal challenge: a Mendelian randomization study; Hughes et al

A B
The association between body mass index and metabolite response to a liquid mixed meal challenge: a Mendelian randomization study; Hughes et al  The association between body mass index and metabolite response to a liquid mixed meal challenge: a Mendelian randomization study; Hughes et al The association between body mass index and metabolite response to a liquid mixed meal challenge: a Mendelian randomization study; Hughes et al The association between body mass index and metabolite response to a liquid mixed meal challenge: a Mendelian randomization study; Hughes et al The association between body mass index and metabolite response to a liquid mixed meal challenge: a Mendelian randomization study; Hughes et al The association between body mass index and metabolite response to a liquid mixed meal challenge: a Mendelian randomization study; Hughes et al -Estimating the number of indpendent features.
-Performing tree cut.Cut height defined at 0.5 -Identifying independent features.-Writing feature summary statistics to file.c.Performing principle component analysis and identifying outliers.
-Re-Writing sample summary statistics to file to include PCs.
-Writing PC statistics to file.
V. Performing data filtering.a. Performing data filtering.(1) Samples with missingness >=80% were excluded.(2) features with missingness >=80% were excluded (xenobiotics are not included in this step).( 3) sample exclusions based on the user defined threshold were excluded.(4) feature exclusions based on user defined threshold were excluded (xenobiotics are not included in this step).( 5) samples with a total-peak-area or total-sum-abundance that is >= N standard deviations from the mean, where N was defined by the user, were excluded.(6) samples that are >= N standard deviations from the mean on principal component axis 1 and 2, where N was defined by the user, were excluded.

Metabolite or feature reduction and principal components
A data reduction was carried out to identify a list of representative features for generating a sample principal component analysis.This step reduces the level of inter-correlation in the data to ensure that the principal components are not driven by groups of correlated features.The data reduction table on the left presents the number of metabolites at each phase of the data reduction (Spearman's correlation distance tree cutting) analysis.On the right principal components 1 and 2 are plotted for all individuals, using the representitve features identified in the data reduction analysis.The red vertical and horizontal lines indicate the standard deviation (SD) cuto s for identifying individual outliers, which are plotted in red.The standard deviations cutto were defined by the user.
2 Filtered data

N
• The number of samples in data = 5518 • The number of features in data = 603 metaboprep report derived from the representative metabolites.The Scree plot also identifies the number of PCs estimated to be informative (vertical lines) by the Cattel's Scree Test acceleration factor (red, n = 2) and Parallel Analysis (green, n = 18).Individuals in the PC plot were clustered into 4 kmeans (k) clusters, using data from PC1 and PC2.The kmeans clustering and color coding is strictly there to help provide some visualization of the major axes of variation in the sample population(s).Analysis details: Of the 603 features in the data 0 features were excluded from this analysis because of no variation or too few observations (n < 40).Of the remaining 603 metabolite features, a total of 271 may be considered normally distributed given a Shapiro W-statistic >= 0.95.

Distributions
A pdf report is being written to NEO_outlier_detection_pre_filtering.pdf that contains dotplot, histogram and distribution summary statistics for each metabolite in your data set, providing an opportunity to visually inspect all your metabolites feature data distributions.

Notes on outlying samples at each metabolite|feature
There may be extreme outlying observations at individual metabolites|features that have not been accounted for.You may want to: 1. Turn these observations into NAs.

Filtered data sample missingness: influenced by possible explanatory variables
The figure provides an illustrative evaluation of the proportion of sample missigness as a product of sample batch variables provided by your supplier.This is the univariate influence of batch e ects on sample missingness.

Multivariate evaluation: batch variables
## [1] " --No sample level batch variables were provided or all were invariable --" Table Legend: TypeII ANOVA: the eta-squared (eta-sq) estimates are an estimation of the percent of variation explained by each independent variable, after accounting for all other variables, as derived from the sum of squares.This is a multivariate evaluation of batch variables on *sample missingness*.9 metaboprep report 4 Sample Total Peak|Abundance Area (TPA): Total peak|abundance area (TPA) is simply the sum of the abundances measured across all features.TPA is one measure that can be used to identify unusual samples given their entire profile.However, the level of missingness in a sample may influence TPA.To account for this we: 1. Evaluate the correlation between TPA estimates across all features with TPA measured using only those features with complete data (no missingness).
2. Determine if the batch e ects have a measurable impact on TPA.

Relationship with missingness
Correlation between total peak area (at complete features) and missingness

Univariate evaluation: batch e ects
The figure below provides an illustrative evaluation of the total peak area as a product of sample batch variables provided by your supplier.Table Legend: TypeII ANOVA: the eta-squared (eta-sq) estimates are an estimation on the percent of variation explained by each independent variable, after accounting for all other variables, as derived from the sum of squares.This is a multivariate evaluation of batch variables on *total peak|abundance area* at complete features.

Power analysis
Exploration for case/control and continuous outcome data using the filtered data set Analytical power analysis for both continuous and imbalanced presence/absence correlation analysis.estimates of power for continuous traits with the total sample size on the x-axis and the estimated power on the y-axis.Figure (B) provides estimates of power for presence/absence (or binary) traits in an imbalanced design.The estimated power is on the y-axix.The total sample size is set to 5518 and the x-axis depicts the number of individuals present (or absent) for the trait.The e ects sizes illustrated here were chosen by running an initial set of simulations which identified e ects sizes that would span a broad range of power estimates given the sample population's sample size.

Figure 2 :Supplementary Figure 3 :
BMI and PGS by visit date.Box-plots of NEO sampling dates (month-year) along the x-axis with the mean (box vertical line), 1 st and 3 rd quartile (box) and 1.5*IQR of box (whiskers) illustrated.The red vertical line illustrates the mean value of the total sample.Upper plot illustrates variation in BMI.The lower plot illustrates the variation in the instrumental variable -BMI-PGS.Defined in each sub-header is the proportion of variation in BMI (upper) or BMI-PGS (lower) explained by sample date.Raw and QC'd fasting and postprandial scatter plots for S-HDL-C.(Left) A scatter plot for fasting (x-axis) and postprandaial (y-axis) cholesterol in small HDL (S-HDL-C).An equivalency line (x = y) is illustrated as a dashed grey line.Each dot is a sample, those colored blue had a zero value and were turned into NA, those colored green were identified as outliers (10 interquartile distances from the median) in one of the two dietary states and turned into NA in that dietary state, those colored red were identified as delta outliers (5 interquartile distances from the median).Delta was estimates as the postprandial value minus the fasting value.Those samples colored as grey were not altered.
sub-sample populationsVariation in BMI by sub−populationSupplementary Figure1 body mass index and metabolite response to a liquid mixed meal challenge: a Mendelian randomization study; Hughes et al body mass index and metabolite response to a liquid mixed meal challenge: a Mendelian randomization study; Hughes et al figer_g height cancer _ r a l a _ f a l a _ p a l a _ r a l b _ f a l b Missingness is evaluated across samples and features using the original/raw data set.1.1.1Visual structure of missingness in your raw data set.

Figure Legend :
Figure Legend: Missingness structure across the raw data table.White cells depict missing data.Individuals are in rows, metabolites are in columns.

−
Figure Legend:The data reduction table on the left presents the number of metabolites at each phase of

Figure
Figure Legend:A matrix plot of the top five principal components including demarcations of the 3rd (grey), 4th (orange), and 5th (red) standard deviations from the mean.Samples are color coded as in the summary PC plot above using a kmeans analysis of PC1 and PC2 with a k (number of clusters) set at 4. The choice of k = 4 was not robustly chosen it was a choice of simplicity to help aid visualize variation and sample mobility across the PCs.

Figure Legend :
Figure Legend: Histogram plots of Shapiro W-statistics for raw (left) and log transformed (right) data distributions.A W-statistic value of 1 indicates the sample distribution is perfectly normal and value of 0 indicates it is perfectly uniform.Please note that log transformation of the data *may not* improve the normality of your data.

2 . 3 Influence of batch variables on filtered data 3 . 1
Figure Legend: Box plot illustration(s) of the relationship that available batch and biological variables have with feature missingness.

Figure Legend :
Figure Legend: Box plot illustration(s) of the relationship that available batch variables have with sample missingness.

Figure Legend :
Figure Legend: Relationship between total peak area at complete features (x-axis) and sample missingness (y-axis).

Figure Legend :4. 2 . 1
Figure Legend: Violin plot illustration(s) of the relationship between total peak area (TPA) and sample batch variables that are available in your data.

Figure Legend :
Figure Legend: Simulated e ect sizes are illustrated by their color in each figure.Figure (A) provides Figure Legend: Simulated e ect sizes are illustrated by their color in each figure.Figure (A) provides

Supplementary Figure 16: NEO and Wurtz et al MR estimate scatter plot. Each
metabolites are labelled with their metabolite ID in the plot.See Supplementary Table1for further annotation.Supplementary Figure17: Postprandial observational profile.Upper: A tile plot of observational effect estimates for lipoproteins and lipoprotein ratios in the postprandial state.
The association between body mass index and metabolite response to a liquid mixed meal challenge: a Mendelian randomization study; Hughes et al −− NEO: fasting, postprandial, and response metabolite trait dendrogram −−

Table Legend :
Six primary data filtering exclusion steps were made during the preparation of the data.

Table Legend :
The table reports the number of point estimates for the minimum (0%) median (50%) and maximum (100%) number of outlying features across samples and the number of outlying samples across features.