A multi-variant recall-by-genotype study of the metabolomic signature of body mass index

Objective We estimated the effect of body mass index (BMI) on circulating metabolites in young adults using a recall-by-genotype (RbG) study design. Methods An RbG study was implemented in the Avon Longitudinal Study of Parents and Children. Samples from 756 participants were selected for untargeted metabolomics analysis based on low/high genetic liability for higher BMI defined by a genetic score (GS). Regression analyses were performed to investigate the association between BMI GS groups and relative abundance of 973 metabolites. Results After correction for multiple testing, 29 metabolites were associated with BMI GS group. Bilirubin was amongst the most strongly associated metabolites with reduced levels measured in individuals with the highest BMI GS (beta=-0.32, 95% confidence interval (CI): -0.46, -0.18, Benjamini-Hochberg (BH) adjusted p=0.005). We observed associations between BMI GS group and levels of several potentially diet-related metabolites including hippurate which had lower mean abundance in individuals in the high BMI GS group (beta=-0.29, 95% CI: -0.44, -0.15, BH adjusted p=0.008). Conclusions Together with existing literature our results suggest a genetic predisposition to higher BMI captures differences in metabolism leading to adiposity gain. In the absence of prospective data, separating these effects from the downstream consequences of weight gain is challenging. Study Importance questions What is already known about this subject? Metabolomics, defined as the measurement and study of circulating small molecules that are the substrates and products of cellular metabolism, is increasingly used by epidemiologists to provide a functional read-out of bulk cellular activity and a proxy to individual current health. This approach also provides insight into biological pathways linking exposures and disease. In observational studies, elevated body mass index (BMI) has been associated with a wide range of circulating metabolites. Researchers are now looking to genetic epidemiological methods, such as Mendelian randomization, to offer insight into potential causal relationships. What are the new findings in your manuscript? We identified 29 metabolites whose relative abundance varies with a genetic predisposition to higher BMI. Bilirubin, a key component of the heme catabolic pathway and a potent antioxidant, showed the strongest association with BMI score group. How might your results change the direction of research or the focus of clinical practice? Results of both Mendelian randomization and recall-by-genotype studies need to be combined with alternative study designs to distinguish between biomarkers that are intermediates on the pathway to BMI from those reflective of metabolic changes that result from increased adiposity. Separating causal biomarkers from non-causative biomarkers of adiposity is important since only the former are relevant to treatment and prevention, whilst both could be informative with respect to prediction and the downstream consequences of high BMI.


Introduction
Despite the extensive focus in the literature, the full downstream impact of high body mass index (BMI) and the potential causal mechanisms by which BMI impacts a large number of non-communicable diseases remains unclear (1). Through a combination of large-scale observational studies and intervention designs, BMI has been established as a major risk factor for many common complex diseases, including type 2 diabetes mellitus, hypertension, myocardial infarction, stroke and several types of cancer (2)(3)(4). Whilst the precision of effect estimates describing the association between BMI and disease (and our confidence in them) has increased in line with greater sample sizes and independent replication, observational studies are limited by confounding, bias and reverse causation. Meanwhile, intervention studies (e.g., weight change protocols) designed to circumvent these conventional limitations have their own challenges -notably, a limited ability to alter BMI to the extent required to quantify an effect, and the necessarily short-term and small-scale nature of such interventions.
In response to these challenges and following developments in understanding genetic contributions to adiposity/BMI, methods from within the field of applied genetic epidemiology are now in use by researchers interested in dissecting the relationship between BMI and health. One such approach in which genetic variants act as an approximation to instrumental variables to evaluate the causal effect of an exposure (adiposity) on an outcome (disease), is Mendelian randomization (MR). In MR, genetic variation fulfils the role of an instrumental variable (5) where the presence of variance in BMI explained by genotype is (in principle) orthogonal to confounding factors and where genotype is assumed to exert an effect on health outcome only through BMI. Whilst validating the likely causal nature of the relationship observed between BMI and many common diseases, these studies do little to explain how the risk is delivered.
The study of circulating metabolites uses various techniques to detect and measure low molecular weight metabolites across a range of body fluids and tissues and can be used to provide a functional read-out of an individual's current health. The use of metabolomics in epidemiological studies is increasing and has the potential to help elucidate the mechanisms linking obesity and associated co-morbidities, as well as identifying biomarkers to facilitate intervention and treatment. To date, studies have shown BMI-associated changes across a range of metabolite classes, including sex steroids, branched-chain and aromatic amino acids, acylcarnitines and lipids (6). But with much of the existing literature on the metabolomic impact of adiposity being based on observational epidemiological analyses, gaps remain in our understanding of the biology underpinning the development and direct pathophysiological consequences of obesity.
We aimed to integrate the use of genetic predictors for BMI with in-depth intermediate phenotyping to explore the relationship between BMI and metabolic health. Recall by genotype (RbG) is a study design in which participants (or samples) are selected from a pre-existing cohort based on genetic variation either at single variants or in the form of a genetic score (GS) (7). In this way, RbG exploits the concept of MR (i.e., the random assortment of genetic variants in offspring), enables greater power for a given number of samples analysed as compared to random selection, facilitates deep phenotyping and is typically less prone to confounding and reverse causality (7,8). The aim of this RbG study was to examine the effect of BMI (the exposure) on circulating metabolites (the outcome) using a GS describing a high versus low predisposition for higher BMI.

Methods
In this study, metabolomics data were derived from plasma samples collected during a routine research clinic conducted as part of the Avon Longitudinal Study of Parents and Children (ALSPAC). Individuals were selected for inclusion in the study based on a GS for BMI. The primary statistical analysis sought to identify metabolites whose levels were associated with GS group (low/high). Sensitivity analyses were performed to further characterise the associations observed. An overview of the study design is shown in Figure 1.

Study participants
ALSPAC is a prospective birth cohort of 14

Phenotype data collection
In ALSPAC, regular clinic visits of subsets of G1 were carried out from 4 months to 24 years old, including assessment of their basic anthropometric measures. For this study, data were extracted for several variables to characterise the GS-derived groups. To validate the performance of the GS, data were extracted from all available timepoints for BMI (kg/m 2 ) (calculated as (weight (kg) / height squared (m 2 )). Data were also extracted for other measures of adiposity; this included weight (kg) and variables derived from a dual energy x-ray absorptiometry (DXA) scan of body composition, specifically, total body fat mass (kg) and total body lean mass (kg).
According to the theory of MR and analogous to the situation in a randomised controlled trial, creating groups based on a GS for BMI (as opposed to using BMI itself), should ensure those groups do not differ in any other respect, meaning any downstream analyses should be free from confounding. However, in the presence of (unmeasured) population structure, associations can be induced between GS (and therefore score group) and some of the traditional confounders of the relationship between BMI (as an exposure) and selected outcomes (15). Therefore, data were extracted for several phenotypic correlates of observed BMI to check for associations with score group and to evaluate the potential for them to act as confounders in the primary analysis. Full details of all phenotypic variables can be found in Supplementary Methods.

Statistical analysis
An overview of the statistical analysis is provided in Figure 2. In all analyses, the low BMI score group was treated as the reference such that estimated effects represent the difference in the high BMI group relative to the low BMI group. All analyses were conducted in RStudio (16) using R v4.0.2 (17).

Metabolite processing
We processed the raw (original scale) data received from Metabolon (N=760 samples) in preparation for statistical analysis using an in-house pipeline developed in R (17). Data were filtered based on a series of quality metrics (e.g., missing data).
Except for those in the xenobiotic class (as annotated by Metabolon) metabolites with more than 20% missing values were excluded from the analysis. Following these exclusions, missing data were imputed using a random-forest based method and imputed data transformed using a rank-based normal transformation (RNT). In the case of xenobiotics (metabolites not produced by the human body), those with >20% missing were transformed to presence/absence (P/A) binary phenotypes such that missing values were replaced with 0 and recorded (non-zero) abundance measures replaced with 1. Xenobiotics were treated in this way because of their typically high level of missingness (or absence), which is both expected and biologically relevant given their (predominantly) exogenous origins. Xenobiotics present in <11 samples were excluded from downstream analyses on the basis that any statistical analyses would not be robust. A detailed description of this preanalysis processing can be found in Figure 2 and in Supplementary Methods. This procedure left 973 metabolites for analysis, including 905 continuous metabolites and 68 xenobiotics transformed to P/A.

Characterisation of recall groups
Between-group differences in our phenotype of interest, BMI (kg/m 2 ), the previously described adiposity traits and potential confounders, including technical covariates, were assessed. The normal distribution of continuous variables was checked by a Shapiro-Wilk test and between-group differences assessed using a Student's (twosample, two-sided) t-test assuming unequal variance. In the case of BMI and weight, between-group differences were evaluated at different ages, ranging from four months to 24 years (when the samples used in this study were collected). Betweengroup differences in body composition measures were assessed based on DXA scans conducted during the age 24 years clinic. Between-group differences of continuous variables that could confound the primary analysis were also assessed by t-test. A Fisher's exact test for count data was applied to test for possible association of group with binary and categorical variables.

Primary analysis -association of metabolites with recall group
To identify metabolite levels that differed between the low BMI and high BMI score groups, mean abundance was compared between groups using regression models.
In Model 1, the post-imputation RNT metabolites were analysed within a linear regression framework [metabolite ~ BMI.score.group]. The R 2 from the model was used as an indication of the variance explained by score group. Log 2 median fold change calculated as the ratio of median abundance (untransformed and unimputed) in the high BMI score group divided by median abundance in the low BMI score group was used to indicate relative effect sizes.
In Model 2, metabolites in the xenobiotic class with high levels of missingness and previously transformed to P/A traits, were analysed within a logistic regression framework [metabolite ~ BMI.score.group]. In this case, the variance explained by the model was estimated using the 'rsq' function in the R package of the same name (18). A Benjamini-Hochberg (BH) correction was applied to adjust the p-values obtained from each of these analyses (Model 1 and Model 2) for multiple testing.

Sensitivity analyses
Several sensitivity analyses were carried out to aid interpretation and to further characterise the associations observed in the primary analysis. These analyses were restricted to the subset of BMI score group associated metabolites (BH p<0.05) output from Model 1. A full description of the methods used can be found under corresponding headings in Supplementary Methods.
Extension of primary association analyses: Model 1 was extended to a multivariate model in which any potential confounder that had previously been shown to be associated with score group was fitted as an independent fixed effect alongside score group.
Metabolite correlation analysis: A hierarchical clustering approach was applied to the subset of associated metabolites to identify redundancy in the data (i.e., where associated metabolites were highly correlated and likely representing the same biological signal). A reduced set of 'representative' metabolites was derived forming the focus for the next steps.

Results
After filtering based on a series of pre-defined quality metrics, the study sample  Table 1. Both recall groups were consistent with the overall cohort in terms of age and sex distribution, whilst adiposity traits showed expected differences (see below for details).

Characterisation of recall groups
At the age 24 years clinic when the samples were collected, the mean (standard deviation, SD) BMI of individuals in the low BMI score group was 23.4 (3.7) kg/m 2 , falling within the 'normal weight' range as defined by the World Health Organisation (18.5kg/m 2 to <25kg/m 2 ). In contrast, the mean BMI of individuals in the high BMI score group, 26.1 (5.2) kg/m 2 , fell within the 'pre-obesity' range (25kg/m 2 to <30kg/m 2 ). Differences were also observed in weight and in both total fat mass and total lean mass at this timepoint ( Table 1). Temporal analyses showed that the between-group differences in BMI emerged at about four years of age then increased rapidly until participants reached around 13 years of age and somewhat more slowly thereafter (Figure 3 and Supplementary Table S1); a similar pattern was observed in weight ( Figure S1 and Supplementary Table S1). There was little evidence for an association between BMI score group and most potential confounders tested (Supplementary Table S2). BMI score group showed modest association with parental (mothers and mother's partners) social class (Supplementary Table S2) but with no clear direction of effect across categories (Supplementary Figure S2).

Association of metabolites with recall group
Overall, we observed relatively small differences across a wide range of molecules with median log 2 fold changes typically in the range -0.5 to 0.5 and a slight bias towards decreased abundances in the high BMI score group (Figure 4). Of

Sensitivity analyses
Characterisation of the score groups indicated some association with mother's and mother's partner's social class. Therefore, in sensitivity analyses, these variables were fitted alongside score group in a multivariate model for the 29 Figure S4).
Lower plasma levels of bilirubin degradation product (C16H18N2O5 (1) Plots for the representative metabolites not shown in Figure 5 can be found in Supplementary Figure S5.

Discussion
In this study, we characterised the metabolic profile associated with low/high genetic liability for higher BMI using an RbG framework. The mean difference in BMI between the low and high BMI score groups increased from early childhood, reaching a maximum of 2.8 kg/m 2 (95% CI: 2.1, 3.4) at time of sampling when individuals were on average 24.5 years of age. This observation reflects differences in the ability of the GS to capture variation in BMI at different ages as shown previously in the same cohort (19). We identified 29 metabolites associated with BMI GS group allocation. In most cases, associated metabolites were seen at lower levels in the high BMI score group, with the largest effects seen for bilirubin, hippurate and tridecenedioate. Two sphingomyelin metabolites were seen at increased abundance in the high BMI score group. The potential relevance of a selection of these metabolites to health and disease is explored in Supplementary   Table S7.
Conventionally, MR and RbG approaches are an attempt at isolating the causal contribution of modifiable exposures, such as BMI, to chosen outcomes. However, when the outcome is also a biological intermediate that may itself be directly proxied by (elements of) the genetic predictor used to capture variance in the exposure, as is the case with metabolites, the assumptions underpinning the causal inference no longer hold. Holmes and Davey Smith have previously described seven scenarios which could underpin observed associations between a genetic risk score for disease and potential biomarkers (20). The scenarios, all of which we propose could equally apply to the present study, include both 'real' GRS-to-trait associations (which are therefore informative with respect to underlying biology) and associations that represent potential artifacts. Whilst others have concluded that the majority of metabolic perturbations seen in obesity are a response to increased adiposity itself rather than shared genetic mechanisms (21), our results and those of Hsu et al. (22), suggest forward causal connections (i.e., metabolites involved in pathways leading to changes in BMI) may exist. Whilst both sets of metabolic pathways, cause and effect, may be informative with respect to predicting the risk of developing obesityassociated comorbidities, only the former is likely to be of therapeutic relevance for the prevention of obesity. It is within this context that we go on to discuss our findings in more detail.
Bilirubin showed the strongest association with BMI score group allocation and, together with its associated degradation products, formed the largest cluster of associated metabolites. Bilirubin is a key component of the heme catabolic pathway and a potent antioxidant, and in this study was found in lower abundance in the high BMI score group. Although bilirubin is not amongst the metabolites most commonly associated with BMI and adiposity (23)  However, the apparent instrument-dependent nature of the findings in the latter (22) point to heterogeneity in the underlying biology. Levels of circulating BCAA are also known to be influenced by dietary intake (39) and a link has been proposed between the obesity-related rise in circulating BCAA and a decline in their catabolism in adipose tissue (40,41), with further evidence suggesting that this could be tissue specific (42). Moreover, 3-hydroxy-2-ethylproprionate (a metabolite annotated to the 'Leucine, isoleucine and valine metabolism' sub-pathway and which is a product of isoleucine catabolism (43)) was observed to associate with BMI GS group but not with measured BMI. One potential explanation for this given previously reported associations of 3-hydroxy-2-ethylproprionate with muscle cross-sectional area (i.e., with body composition) (44), is that the association seen with score group is underpinned by differences in lean mass between the groups. However, we are not well-powered to investigate this hypothesis within the current study.
We observed associations between BMI score group and the levels of potentially diet-related metabolites, including hippurate and PFOS. Hippurate, or hippuric acid, a glycine conjugate of benzoic acid, is synthesised in the liver and kidney (45). The benzoic acid component is derived mainly via microbial and mammalian cometabolism of large polyphenolic molecules contained in, for example, fruits and vegetables, to a range of simpler aromatics that are then further metabolised to benzoic acids (45). We observed lower levels of hippurate in individuals in the high BMI GS group, concordant with a previous identified associations between hippurate levels and visceral body fat mass (27,46). Previous literature combining data on diet intake, visceral fat mass and gut microbial profiling, suggests the association of circulating (and urine excreted) levels of hippurate with adiposity and related health outcomes (45), is likely to be the result of a complex interaction between diet intake, microbiome diversity and composition, and adipose tissue function (27,47).
Individuals in the high BMI score group also had lower plasma PFOS levels. PFOS, as an anthropogenic organic pollutant with chemical and thermal stability, has been detected in drinking water and the diet (especially in fish and crustaceans) in multiple countries and areas across the globe and has a global toxic effect on human health (see Supplementary Table S7). Furthermore, several of the xenobiotics which appear to be present at different rates in the two groups (albeit not meeting our stringent threshold for association) may also be biomarkers of food consumption (e.g., acesulfame, betonicine, theanine).
The associations observed between BMI score groups and these metabolites suggest that at least some of the genetic predisposition to increased BMI may be conveyed either via dietary choices or through differences in nutrient metabolism.
Others have taken the fact that many genes associated with high BMI appear to be highly expressed in the central nervous system (48) (with some also having been linked to appetite regulation), as evidence that genetic susceptibility to obesity is partly attributable to appetitive phenotypes (49). Studies of the association between BMI-associated GS and appetite and satiety traits (50,51) along with results of the current study provide some support for this hypothesis. However, behavioural traits such as these are known to be particularly at risk from bias even in an MR (and likely RbG) setting (52), where population stratification (53) or complex genetic effects not accounted for in these study designs, for example, dynastic effects (54) can be problematic given the underlying assumptions of these methods.
The potential for reintroduction of such confounding effects (and pleiotropy) appears to increase with the number of variants included in the GS, either through the lowering of the p-value threshold used for selecting SNPs or through increasing the power of association analyses such that ever smaller effects are detected with sufficient statistical certainty (15). In this study, the weak correlation between score group and social class suggests some residual confounding (due to population stratification) may be present. However, given the consistency in the effect estimates after adjusting for maternal and paternal social class, we believe the effects of any such confounding on our results to be small. In this study, where metabolites showed an association with BMI score group, we also estimated their association with measured BMI. Whilst in general we saw good concordance of effects estimated in the two analyses, given our sampling frame, the observational estimates cannot necessarily be extrapolated to the wider (unselected) population. Where we see inconsistency, this could be the result of different key sources of bias effecting the different analytical strategies.
In conclusion, we used an innovative RbG study design to identify metabolites whose relative abundance varies with a genetic predisposition to increased BMI. We hypothesise that these differences may reflect gene-derived perturbations to biological pathways relevant to weight gain or be consequences of higher BMI itself.
To provide a definitive answer to the question of the role of metabolites (and related biological pathways) in obesity, results from several different approaches with unrelated sources of bias, including challenge/intervention studies, need to be integrated. In doing so, we can begin to understand the role of different metabolic pathways in weight gain and related morbidity and partition metabolites according to whether they are likely to be causative in nature and/or reflective of metabolic changes that occur after increased adiposity. This study involves the first-generation offspring in the Avon Longitudinal Study of Parents and Children (ALSPAC) multi-generational cohort, in which 14,541 pregnant women resident in the South West of England were recruited in the 90s. Firstly, we constructed genetic scores (GS) for body mass index (BMI) for all first-generation offspring. Under the recall-by-genotype study design, we recalled the plasma samples (collected at the age 24 years clinic) of individuals with a low (yellow) or high (blue) BMI GS for further analysis. Then, metabolites in those plasma samples were quantified by Metabolon. Finally, we performed statistical analysis to compare the metabolites levels between the two BMI GS groups. Our results are relevant to understanding the role of metabolites both as intermediates on the pathway to BMI and from BMI to disease.

Figure 2. Overview of statistical analysis.
'Raw data' is the original scale data normalized in terms of raw area counts (as supplied by Metabolon). Data were prepared for statistical analysis by first filtering samples and metabolites based on a series of quality metrics and then applying imputation and re-scaling procedures as appropriate. SD = standard deviations; QC = quality control; BMI = body mass index; GS = genetic score.

Figure 3. Mean differences in BMI between the high and low BMI score groups.
Error bars represent the standard errors of mean differences in weight. Sample size ranges from 108 (at age 31 months) to 743 (at age 24 years). Test results are given for a Students (two-sample, two-sided) t-test. ***: p-value < 0.001; **: p-value < 0.01; *: p-value < 0. For full results see Supplementary Table S1. Points are coloured by super-pathway. Log 2 median fold change calculated as the ratio of median abundance (untransformed and unimputed) in the high BMI score group divided by median abundance in the low BMI score group. P-values used to derive -log 10 (p) are those from the linear regression analysis. All points above the dashed line have a Benjamini-Hochberg adjusted p-value <0.05. Solid grey lines indicate the density of points. A representative selection of metabolites of known identity are labelled.
*: Indicates a compound that has not been confirmed based on a standard.  ). In the plots, solid lines denote the predicted univariate within score group relationship between BMI and metabolite with a 95% confidence interval denoted by shading.

Post hoc analysis
With z-scored metabolite levels, further examine the association between metabolites and measured BMI Age at clinic visit (years) Mean differences in BMI (kg/m2)