Integrated analyses of growth differentiation factor-15 concentration and cardiometabolic diseases in humans

Growth differentiation factor-15 (GDF15) is a stress response cytokine that is elevated in several cardiometabolic diseases and has attracted interest as a potential therapeutic target. To further explore the association of GDF15 with human disease, we conducted a broad study into the phenotypic and genetic correlates of GDF15 concentration in up to 14,099 individuals. Assessment of 772 traits across 6610 participants in FINRISK identified associations of GDF15 concentration with a range of phenotypes including all-cause mortality, cardiometabolic disease, respiratory diseases and psychiatric disorders, as well as inflammatory markers. A meta-analysis of genome-wide association studies (GWAS) of GDF15 concentration across three different assay platforms (n=14,099) confirmed significant heterogeneity due to a common missense variant (rs1058587; p.H202D) in GDF15, potentially due to epitope-binding artefacts. After conditioning on rs1058587, statistical fine mapping identified four independent putative causal signals at the locus. Mendelian randomisation (MR) analysis found evidence of a causal relationship between GDF15 concentration and high-density lipoprotein (HDL) but not body mass index (BMI). Using reverse MR, we identified a potential causal association of BMI on GDF15 (IVW pFDR = 0.0040). Taken together, our data derived from human population cohorts do not support a role for moderately elevated GDF15 concentrations as a causal factor in human cardiometabolic disease but support its role as a biomarker of metabolic stress.


Introduction
Obesity accounts for an estimated 2.8 million deaths worldwide and 2.3% of global disability-adjusted life years (World Health Organization, 2008). Obesity has been causally linked to a variety of cardiometabolic risk factors and diseases, including fasting insulin, systolic blood pressure, and type 2 diabetes (Holmes et al., 2014). Therefore, interventions to reduce obesity are recommended in the management of these diseases (Yumuk and Tsigos, 2016). Currently, the most effective obesity intervention is bariatric surgery (Jammah, 2015); however, its clinical utility is limited as it is highly invasive and accompanied by a high risk of long-term complications such as gastroesophageal reflux disease and nutritional deficiencies (Mesureur and Arvanitakis, 2017). New therapies such as GLP-1 receptor agonists aim to target obesity via appetite suppression but are associated with common side effects including nausea and diarrhoea (Wilding et al., 2021). The development of additional therapeutic strategies targeting obesity as an underlying cause and pathology of cardiometabolic disease represents an area of unmet clinical need.
GDF15, a distant member of the transforming growth factor-β family (Bootcov et al., 1997;Emmerson et al., 2018), is upregulated in response to cellular stress and several diseases, including cancer. It has been suggested that GDF15 functions as a stress response agent implicated in organ injury (Kempf et al., 2006;Tsai et al., 2018;Zimmers et al., 2005). GDF15 was originally identified to play a role in tumour-induced anorexia/cachexia in mice through appetite regulation . The body weight reduction is considered to be mainly driven by food intake inhibition mediated by its receptor GDNF-family receptor alpha like (GFRAL) in distinct areas of the brainstem (Emmerson et al., 2017;Hsu et al., 2017;Mullican et al., 2017;Yang et al., 2017). GDF15 plasma levels are elevated in animal models of obesity (Lockhart et al., 2020;Xiong et al., 2017); however, administration of recombinant protein robustly lowers body weight in obese and diabetic animals, including non-human primates (Mullican et al., 2017;Xiong et al., 2017). Similarly, genetic overexpression of Gdf15 results in decreased body weight and increased resistance to obesity associated with high-fat diets, whereas Gdf15 and Gfral deficient mice are more susceptible to diet-induced obesity (Mullican et al., 2017;O'Rahilly, 2017;Tran et al., 2018). Recent preclinical work has identified GDF15's role as a sentinel protein, upregulated in response to various ingested toxins triggering nausea-related behaviour in rodents assessed by pica and CTA (Conditioned Taste Aversion: Borner et al., 2020b;Patel et al., 2019). It is hypothesised that this contradictory finding of upregulated levels in obesity but weight loss upon overexpression may be explained by GDF15 acting as a signal to the brain to induce weight loss (Villanueva, 2017).
Currently, there is a lack of data supporting GDF15 as a potential therapeutic target in humans. However, observational studies have revealed strong positive correlations of GDF15 plasma levels with body mass index (BMI), insulin resistance, age, and mean arterial blood pressure (in obese individuals: Tsai et al., 2015;Vila et al., 2011). Higher GDF15 levels have also been associated with all-cause mortality as well as mortality associated with heart failure and acute myocardial infarction, cancer, advanced heart failure, and end-stage chronic kidney disease (Adela and Banerjee, 2015;Khan et al., Lemmelä 2009; Nair et al., 2017;Wiklund et al., 2010). Metformin therapy, a type 2 diabetes treatment, has been shown to depend on GDF15 to lower body weight in mice and plasma GDF15 levels increased up to 40% upon metformin therapy in humans (Coll et al., 2020;Day et al., 2019). In fact, an analysis of changes induced by Metformin treatment across 237 biomarkers showed GDF15 levels to be most significantly altered, independent of glucose and glycosylated haemoglobin (Gerstein et al., 2017). Recent data also suggest that GDF15 is involved in hyperemesis gravidarum (Fejzo et al., 2018a), supporting the hypothesis that GDF15 triggers anorexia and subsequent weight loss, at least partly, through the induction of malaise (Borner et al., 2020a). Mendelian randomisation (MR) has previously been applied to explore the causal relationship of GDF15 levels and cardiometabolic diseases, with causal associations reported for high-density lipoprotein (HDL) cholesterol and bone mineral density (BMD: Cheung et al., 2019;Folkersen et al., 2020). However, these analyses were based on small genetic studies for GDF15 and have not been replicated. Opposing results have also been reported regarding the association of GDF15 on BMI, with one study finding a causal relationship (Karhunen et al., 2021) and another reporting no significant effect (Au Yeung et al., 2019).
In this study, we used data from several large biobanks to conduct a systematic and extensive phenotypic and genotypic analysis of GDF15 with cardiometabolic traits and diseases, and to ascertain the causal relationship between GDF15 levels and cardiometabolic traits using MR and proteintruncating variant (PTV) analysis.

Association of GDF15 plasma levels with 676 disease outcomes
To assess systematically the relationship between GDF15 plasma levels and clinical phenotypes, we utilised a large Finnish biobank, FINRISK. The FINRISK cohort comprises a cross-sectional population survey carried out over a 40-year period in Finland. GDF15 plasma levels were available for 6610 participants from the 1997 recruitment cohort linked to 676 disease outcomes and were included in this analysis. The median GDF15 concentration (measured using an immunoluminometric assay [ILMA]) was 796 ng/L (interquartile range [IQR] = 664-986). First, we examined the association of baseline characteristics with GDF15 plasma levels (Supplementary file 1a). Age explained 28% of the observed GDF15 variance, with current smoking and BMI accounting for 1.9% and 0.08% of the variance, respectively (Supplementary file 1b). Gender did not show a significant association. All subsequent analyses have been corrected for these covariates (age, gender, smoking, and BMI). Analyses corrected for age and gender alone are also presented for comparison , although results were similar.
We then investigated potential associations of GDF15 plasma levels with a range of disease phenotypes (both prevalent and incident), focusing on the phenotypes defined by the FinnGen consortium (Tuomo et al., 2020), as these endpoints have undergone additional clinical validation (see Materials and methods). A total of 676 disease endpoints were examined as dependent variables in association analyses with GDF15 plasma levels (the independent variable). After multiple-testing correction using false discovery rate (FDR, p FDR < 0.05), GDF15 was significantly associated with 80 disease endpoints ( Table 1 and Supplementary file 2c).

GDF15 plasma concentration associates with prevalent and incident disease and independently predicts all-cause mortality and cardiometabolic diseases
To determine the prognostic potential of GDF15 in predicting incident disease, we investigated its plasma levels in relation to prevalent and incident diseases in the FINRISK cohort. We selected disease endpoints showing association to the GDF15 plasma levels, namely type 2 diabetes, atherosclerosis (excluding cerebral, coronary, and PAD), COPD, psychiatric disorders, and malignant neoplasm of respiratory system and intrathoracic organs (Supplementary file 2d).
Other tested endpoints (psychiatric disorders; malignant neoplasm of respiratory system and intrathoracic organs) also showed similar results for which plasma GDF15 levels seemed to predict the disease (Supplementary file 2d). Together, these results of several disease endpoints indicate that increased GDF15 plasma levels could represent an early biomarker in individuals before disease diagnosis.

Association of GDF15 plasma levels with 96 quantitative biomarkers
We then performed an association analysis of GDF15 plasma levels with 96 quantitative biomarkers in FINRISK. After multiple-testing correction (p FDR < 0.05), significant associations were observed with 45 biomarkers, most notably, with blood markers known to be associated with inflammation (Supplementary file 2f).
Genome-wide association studies of plasma levels of GDF15 quantified using three different assays in two independent cohorts We performed GWAS to identify the genetic determinants of GDF15 plasma levels using two independent cohorts (FINRISK and INTERVAL) and three GDF15 quantification methodologies (ILMA in FINRISK and Olink and SomaScan proteomic assays in INTERVAL). Results from individual GWAS analyses are reported in Supplementary file 3a-c. In FINRISK GWAS was conducted on 5817 individuals with available GDF15 plasma levels, measured using an ILMA (Sinning et al., 2017;Supplementary file 3a). The most significant associations were found within or in proximity to the GDF15 gene on chromosome 19, comprising 159 significant variants (p-value < 5 × 10 -8 ). In INTERVAL we assayed plasma GDF15 in 3301 participants using SomaScan (an aptamer-based multiplex protein assays) and in 4998 participants using a proximity extension-based antibody assay (Olink). The most significant associations were found around the GDF15 gene on chromosome 19, comprising 134 significant variants for INTERVAL-SomaScan and 72 significant variants for INTERVAL-Olink (p-value < 5 × 10 -8 :Supplementary file 3b-c).
We then performed a meta-analysis of the genome-wide significant variants identified in the GWAS performed in the FINRISK, INTERVAL-SomaScan, and INTERVAL-Olink cohorts to generate a combined statistic and heterogeneity value. The analysis identified four genetic variants (rs1059369, rs1054221, rs1227734, rs189593084) at the GDF15 locus to be associated with GDF15 plasma levels (p-value < 5 × 10 -8 : Supplementary file 3d).
Two variants had a strengthened signal in the combined meta-analysis, rs1054221 and rs1227734 (Supplementary file 3d), However, five genetic associations, which were genome-wide significant for two of the three GWAS, displayed substantial heterogeneity, for example, rs16982345 (heterogeneity I 2 =99.8%, heterogeneity p-value = 7.6 × 10 -180 ; Supplementary file 3d). By exploring the LD structure between the heterogeneous variants (Supplementary file 3e), we found that these variants were all located within the same LD block that included the missense variant rs1058587 (p.H202D). Heterogeneity in GWAS results of GDF15 levels has been observed in other studies Pietzner et al., 2020) and significant variations in GDF15 antibody binding affinity dependent on rs1058587 have been previously reported (Fairlie et al., 2001).

Mitigation of epitope binding effects due to the rs1058587 missense variant
To remove any potential epitope effects, we applied conditional analysis on rs1058587 when linking plasma GDF15 concentrations with disease phenotypes and quantitative traits. A total of 59 significant disease associations and 43 biomarker associations were identified. These associations largely did not differ from the unconditioned results, although 21 fewer disease associations reached significance .
To minimise confounding of the GWAS results from rs1058587 epitope effects, we conditioned the association signal on chromosome 19 on rs1058587 in both FINRISK (Supplementary file 3f) and INTERVAL-SomaScan (Supplementary file 3g). INTERVAL-Olink did not have a significant association at this LD block. We then performed a meta-analysis across the genome for all three cohorts. Only signals within 1 Mb of the GDF15 gene reached significance in the meta-analysis. Manhattan and Quantile-Quantile (QQ) plots of this GWAS meta-analysis are shown in Table 2, Figure 2, and Supplementary file 3h. We found 146 significant associations in this region, with rs1227734 identified as the strongest association (beta = 0.50, p-value = 2.5 × 10 -187 ).

MR analysis suggests a causal relationship between GDF15 plasma levels and HDL
Observed associations of GDF15 plasma levels with obesity-related diseases have been presented here and found in previous reports (Tsai et al., 2015;Vila et al., 2011) and variants in the GDF15 gene region have been associated with cardiovascular traits, cholesterol, WHR, and BMI (Wang et al., 2021a). We therefore applied MR to assess the relationship between genetically determined GDF15 plasma levels and BMI, WHR, glucose, and type 2 diabetes. We also included HDL cholesterol and estimated BMD (eBMD) as additional outcomes in our MR analysis due to the positive MR results in previous studies (Cheung et al., 2019;Folkersen et al., 2020).
We applied two-sample MR (Bowden et al., 2015) using LD clumping (R 2 < 0.01, within 1.5 Mb of either side of the lead SNP, see Materials and methods) and identified five independent genetic instruments from the conditioned meta-analysis on chromosome 19: rs111527728, rs113700483, rs1227734, rs150286074, and rs73923175 (regional heritability = 0.061, F-statistics reported in Supplementary file 4a). Association statistics between the conditioned instrumental variables (IVs) and exposure (GDF15) were taken from the above meta-analysis, and the association statistics between the IVs and outcome were extracted from publicly available GWAS summary statistics (see Materials and methods). These values were applied to a random-effect inverse variance weighted (IVW) MR, as well as MR-Egger, weighted median MR, and MR-PRESSO methods (Verbanck et al., 2018), in order to detect horizontal pleiotropy (see Materials and methods for further details on the underlying assumptions of these methods). MR using the conditioned variants found two significant associations of genetically determined GDF15 plasma levels with HDL (IVW estimate = −0.0085, p FDR mass index; GDF15, growth differentiation factor-15; CHD, coronary heart disease; STR, stroke; HDL, high-density lipoprotein. The online version of this article includes the following figure supplement(s) for figure 1:    = 0.0024) and WHR (IVW estimate = 0.017, p FDR = 0.0013). The effects sizes for both of these were very small and the WHR result was not robust to multiple-testing correction in MR-PRESSO (Table 3 and Supplementary file 4b).
Horizontal pleiotropy was not detected in any of the MR analyses. We did not replicate findings of causality of genetically determined GDF15 plasma levels with eBMD but we were able to provide additional validation of the causal role of genetically determined GDF15 in HDL cholesterol (Cheung et al., 2019;Folkersen et al., 2020).
Reverse MR analysis identifies BMI as a causal factor for GDF15 plasma levels We applied reverse two-sample MR (as described above) using GDF15 plasma levels as the outcome variable to assess associations with BMI, WHR, glucose, diabetes, HDL cholesterol, and eBMD as the MR exposures. GDF15 GWAS summary statistics were taken from the conditioned meta-analysis     of GDF15 levels in FINRISK and INTERVAL in all chromosomes. LD clumping was used to identify the genetic instruments (see Materials and methods, F-statistics are reported in Supplementary file 4c-h). We found a significant association between higher genetically predicted BMI and higher GDF15 plasma levels (IVW estimate = 0.097, p FDR = 0.0040) but not any other tested trait ( Carriers of GDF15 PTVs present an opportunity to explore the phenotypic consequences of predicted loss-of-function (LOF) of GDF15 on human disease. We analysed whole-exome sequencing data from 302,388 participants from the UK Biobank, of whom 109 carried GDF15 PTVs in the heterozygous state. We assessed differences in BMI, WHR, glucose, eBMD, and HDL cholesterol between carriers and non-carriers using the Mann-Whitney U test. We were not able to carry out this analysis on type 2 diabetes given too few patients in this subset (n=2). The analysis was restricted to unrelated individuals of European ancestry, leaving 91 carriers of GDF15 PTVs (Supplementary file 4j, male N=42, female N=49, mean age = 56.3) and 40,000 randomly selected European-ancestry non-carriers (male N=20,000, female N=20,000, mean age = 56.9). In line with the MR analyses, we observed no differences in BMI (mean difference = 0.1, p-value = 0.90), WHR (male: mean difference = 0.01, p-value = 0.64, female: mean difference = 0.01, p-value = 0.55), glucose (mean difference = 0.14, p-value = 0.32), eBMD (mean difference = 0.02, p-value = 0.82), or HDL cholesterol (mean difference = 0.08, p-value = 0.08), suggesting that mono-allelic GDF15 LOF does not have a strong effect on these traits. It is possible that the other allele is compensating for the LOF but given that no homozygous GDF15 LOF was found in over 300,000 patients we were not able to test this hypothesis. It may be worthwhile for future analysis to experimentally explore the effect of heterozygous LOF on GDF15 protein levels.

Discussion
In this study, we present a systematic phenotypic and genotypic analysis of GDF15 with a wide range of health outcomes and biomarkers. In line with previous findings, our analysis confirmed strong correlations of GDF15 plasma levels with a range of clinical parameters (e.g. age, smoking, and BMI) and human diseases (e.g. diabetes, cardiovascular and respiratory disease), as well as several inflammatory biomarkers. As preclinical data in mice strongly implicated Gdf15 in the aetiology of obesity and glucose tolerance, we specifically investigated whether human genetic evidence supports these findings. Using data from large biobanks, we found that neither MR nor GDF15 PTV analyses supported a causal role for GDF15 plasma levels in influencing BMI, glucose, diabetes, or eBMD in humans. Instead, we found that higher BMI may cause increases in GDF15 plasma levels highlighting the role of GDF15 as a likely marker of metabolic stress in humans. Additionally, we replicated a causal association of GDF15 with HDL cholesterol and identified a nominal association with WHR. Access to the FINRISK cohort provided a valuable opportunity to explore GDF15 plasma level associations with multiple phenotypes within a single large population. We found the strongest associations with all-cause mortality and cardiometabolic diseases (cardiovascular disease and diabetes), as well as cardiometabolic risk factors (e.g. hypertension, serum triglycerides, BMI), as previously reported (Bonaca et al., 2011;Brown et al., 2002;Daniels et al., 2011;Ho et al., 2012;Kempf et al., 2007;Khan et al., 2009;Lind et al., 2009;Rohatgi et al., 2012;Wiklund et al., 2010;Wollert et al., 2007). Our analyses also found that GDF15 plasma levels are a strong predictor of incident diabetes and an independent predictor of all-cause mortality, cardiovascular disease, and diabetes morbidity, suggesting its potential use as a pre-diabetic prognostic biomarker. We found associations with neoplasms of the lung and digestive system but not with other types of cancer. Our data further support the association of GDF15 plasma levels with inflammatory phenotypes and biomarkers (CRP, mid-regional pro-adrenomedullin: Brown et al., 2007) and uncovered less well-described associations with respiratory disease and psychiatric disorders. Despite previous observations suggesting a role for GDF15 in anorexia (Borner et al., 2020a), we did not identify any association of GDF15 plasma levels with this trait. However, we note that statistical power may have been limited due to the low number of cases (Supplementary file 2c). To identify phenotype associations in an unbiased way, we performed analyses encompassing all phenotypes available applying covariates uniformly across all analyses. We note that covariates specific to a particular disease might have been omitted introducing bias or leading to inflation or deflation of statistics. We detected a strong association between GDF15 levels and smoking, which was recorded and adjusted in our analyses as a binary phenotype. Smoking is an important contributing factor in multiple diseases and it is therefore likely that adjustment as a continuous variable would have been able to better explain the contribution of smoking to the phenotypic association signals, especially in respiratory disease, lung cancer, and some psychiatric disorders. These findings demonstrate that GDF15 plasma levels are a general marker for risk in multiple diseases and are associated with non-specific biomarkers such as CRP.
To define the genetic architecture of plasma GDF15 levels, we performed a GWAS meta-analysis in two large cohorts, FINRISK and INTERVAL, across three different assay platforms. A striking finding of this analysis was the substantial heterogeneity between variants across these studies, consistent with previous reports . Exploring the LD between the variants displaying heterogeneity identified that a large proportion resided within an LD block encompassing a common missense variant, rs1058587 (p.H202D). This variant has been previously associated with hyperemesis gravidarum (Fejzo et al., 2018b) but the presence of this missense variant within this locus raises the possibility of epitope artefacts, indeed a previous study has identified epitope effects in this region (Fairlie et al., 2001). Therefore, it is highly likely that the inconsistency in GWAS results across this locus is driven by differences in the properties of the GDF15 assay used in each study (ILMA in FINRISK, SomaScan, and Olink assays in INTERVAL). Therefore, it may be beneficial to adopt a more consistent method for measuring GDF15 levels within the research environment. Approaches to assay protein levels that do not involve binding to epitopes, such as mass spectrometry, will be informative for avoiding potential artefacts due to protein-altering variants in the future.
Our analyses did not support a causal role for normal human GDF15 plasma levels with BMI, diabetes, glucose, or eBMD. These results were further supported by the lack of association between these phenotypes and GDF15 PTV carriers in the UK Biobank dataset, albeit only heterozygous LOF carriers could be assessed. A meta-analysis completed by Cheung et al. (Ho et al., 2012) identified significant associations of genetically predicted GDF15 levels with HDL cholesterol (beta = −0.048, p-value = 0.001) and eBMD (beta = 0.026, p-value <0.001). Whilst we did not replicate the association with BMD, we do replicate the association of GDF15 with HDL. We additionally found a causal association of genetically determined GDF15 with WHR, although this did not survive multiple-testing correction in MR-PRESSO. This finding may be worth further study although the estimated effect is weak (beta = 0.017). We do not report a significant causal effect of genetically determined GDF15 with BMI consistent with a recent large-scale, pan-ancestry exome-sequencing study of over 640,000 individuals exploring the associations of rare coding variants with BMI that did not report an association with GDF15 (Akbari et al., 2021). Similarly a further study using genetic instruments weakly correlated to rs1058587 (strongest R 2 ≤ 0.02) found no significant causal effect of GDF15 levels on BMI or diabetes traits (Au Yeung et al., 2019). In contrast to these results, Karhunen et al., 2021 recently reported a significant causal association of GDF15 levels with BMI; the study used instruments that included SNPs in strong LD with rs1058587 the variant we identified to lead to a binding artefact in the assays, highlighting the importance of exploring assay bias in genetic studies. Additionally, the larger sample size of the study presented here (n=14,099) offers considerably improved statistical power over previous studies.
The absence of genetic data supporting a causal role for GDF15 plasma levels raises the possibility that GDF15 may behave differently in humans compared to rodents and non-human primates, where there are robust data demonstrating that GDF15 potently reduces body weight. GDF15 levels have also been demonstrated to provoke nausea-related pica and CTA responses in rat and mice, respectively (Borner et al., 2020a;Patel et al., 2019), as well as vomiting in shrews, Suncus murinus, but studies in non-human primates report no signs of nausea (Xiong et al., 2017). Several clinical trials testing the safety and efficacy of GDF15 therapy have been planned or initiated but no data have yet been disclosed. Nevertheless, body weight loss in pregnant women suffering from hyperemesis gravidarium has been linked with higher plasma GDF15 levels (Fejzo et al., 2018a). In FINRISK, information on nausea was not available but diseases with nausea as the main symptom were available (such as gastroesophageal reflux disease) and no significant association with GDF15 plasma levels was found (OR = 1.1, p=0.63). In pregnancy GDF15 levels are found at much higher serum levels than normal (Moore et al., 2000) and it is therefore possible that at higher levels GDF15 may induce nausea and weight loss. Further studies into the effect of higher-concentration GDF15 in humans and the outcomes of clinical trials will elucidate its role further.
As our analyses did not support a causal role for normal human GDF15 plasma levels in the phenotypes examined, we sought to investigate if genetics supported the role for GDF15 as a consequence of these conditions. We conducted reverse MR using GDF15 plasma levels as the outcome and the previously assessed cardiometabolic traits and diseases as the exposures. We identified a significant causal association of BMI on GDF15 plasma levels, suggesting that higher GDF15 levels could be a consequence of higher BMI, supporting the recently reported role of GDF15 in response to stress (Patel et al., 2019). However, the contradictory findings of horizontal pleiotropy between MR-Egger and MR-PRESSO lead to residual uncertainty in the interpretation of the MR analysis. It is possible that the heterogeneity found in our study caused by the use of different assays is driving this pleiotropic finding. A recent report in a MR study of GDF15 levels performed by the SCALLOP consortium also found a significant causal effect of BMI (IVW effect = 0.20, p-value = 1 × 10 -7 : Folkersen et al., 2020). Our replication of this finding strengthens the evidence that GDF15 levels are likely a consequence of differences in BMI. Further research will be required to determine conclusively the role of horizontal pleiotropy in this relationship.
Using valid genetic instruments in MR is of key importance as an invalid instrument would violate the assumptions of the model and lead to unreliable results. To minimise bias from differential assay binding, we conditioned on rs1058587 in our analyses, which could be overcautious. The association of rs1058587 with adiposity traits and hyperemesis gravidarum raises the possibility that GDF15 exerts a causal effect that is not mediated by altered plasma levels, such as via altered proteoforms and/or altered GFRAL binding. With the use of a wide range of GDF15 assays and the possibility that this missense variant is impacting the assay functionality in different assays to varying degrees, the generation of functional data quantifying assay performance is essential for understanding the contribution of this association to GDF15 plasma levels and improving the power of MR analyses. Similarly, as our GDF15 PTV data is based solely on heterozygous carriers, it is possible that the other allele is compensating for the loss and further assessment of the impact of heterozygote GDF15 PTVs on the presence of the protein will be required to validate our findings.
Taken together, our study provides a systematic and broad investigation of GDF15 phenotypic and genotypic associations, identifying possible epitope artefacts that could affect the data validity of GDF15 assays and introduce systematic bias in genetic analyses. Taking into account these biases, our genetic analyses did not support a causal association between normal human GDF15 plasma levels and obesity and diabetes, albeit a nominal finding of a causal association of WHR may warrant further investigation. Conversely, we found that increased GDF15 plasma levels may be a consequence of higher BMI, suggesting that GDF15 may act as a stress-induced biomarker. However, it is conceivable that at normal human plasma levels we are underpowered to detect GDF15's impact on BMI and elevated levels may induce a different or more pronounced effect. Additionally, whilst we have controlled for the assay bias caused by rs1058587, it is possible that this adjustment is too conservative. We caution that future studies with unbiased GDF15 measurements should reassess the impact of the variant using further Mendelian randomisation analyses. Our results do not support GDF15 plasma levels as a causal factor at normal human plasma levels for obesity or its related cardiometabolic diseases.

Study population and phenotypes FINRISK
This study was carried out in accordance with the recommendations of the Declaration of Helsinki. All participants of studies have given written informed consent. FINRISK study was approved by the Ethics Committee of Helsinki and Uusimaa Hospital District.
The FINRISK cohort comprises a cross-sectional population survey carried out over a 40-year period from the year 1972 across multiple regions in Finland. The study aimed to assess the risk factors of chronic diseases and health behaviours in a working age population (Borodulin et al., 2018;Vartiainen et al., 2000) and consisted of 6000-8800 individuals per survey. Measurement of GDF15 plasma levels using an ILMA was included for participants from the 1997 recruitment cohort. Participants were additionally matched to their electronic health records giving access to longitudinal prescription records, death records, and diagnosis history. The cohort characteristics are summarised in Supplementary file 1a. In total, 6610 Finnish individuals from the FINRISK 1997 cohort with available GDF15 plasma concentrations and up to 676 disease outcomes were included in this analysis. Nurse assessment/interview, self-report data, and blood samples were all collected at various time points during the study. A wealth of quantitative biomarkers (e.g. blood lipids, blood sugar markers, inflammatory biomarkers, cytokines, blood count, fatty acids, metabolome) and physiological measures (e.g. anthropometrics, cardiovascular physiology, body composition) were collected from the study participants. In addition, participants who consented were matched to their electronic health records giving access to hospital discharge registry (years 1969-2015), hospital discharge registry of specialist health care operations (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), death registry , drug reimbursement registry , drug purchase registry , and cancer registry . Gender was obtained from social security numbers. A group of clinicians working in the FinnGen consortium constructed 676 disease endpoints combining information from multiple health registries. ICD8, ICD9, and ICD10 codes were utilised in the disease endpoint definitions, for more information see https://www.finngen.fi/en/researchers/clinical-endpoints. Genotyping was carried out for in 6538 individuals who were recruited in 1997, for more information see Supplementary materials.

INTERVAL
The INTERVAL study is a prospective cohort study comprising approximately 50,000 participants nested within a pragmatic randomised trial of blood donors (Jammah, 2015). Between 2012 and 2014, blood donors aged 18 years and older were recruited at 25 NHSBT (National Health Service Blood and Transplant) donor centres across England. Participants were generally healthy because people with a history of major diseases (such as myocardial infarction, stroke, cancer, HIV, and hepatitis B or C) and those who have had recent illness or infection were ineligible to donate blood. Participants completed an online questionnaire which included questions about demographic characteristics (e.g. age, sex, and ethnicity), anthropometry measures (height, weight), and lifestyle information (alcohol intake, smoking, physical activity, and diet).
Informed consent was obtained from all participants and the INTERVAL study was approved by the National Research Ethics Service (11/EE/0538).
For the SomaScan assays, two non-overlapping subcohorts of 2731 and 831 participants were randomly selected, of which 3301 participants (2481 and 820 in the two subcohorts) remained after genetic quality control.
For the Olink assay, protein measurements were conducted in an additional subcohort of 4998 INTERVAL participants aged over 50 at baseline; 4987 samples passed quality control for this panel and were included in the analyses.
For information on genotyping in the INTERVAl cohort, see Suplementary materials.

UK Biobank
UK Biobank is a population-based cohort consisting of ~500,000 individuals with participants recruited between 2006 and 2010 and aged between 40 and 69. Recruitment and cohort information has been previously described (Sudlow et al., 2015). Electronic health records, self-report questionnaires, diet and lifestyle information, and biomarker data are available on this cohort. Here, we utilised BMI and WHR, which were measured during a nurse assessment, as well as diabetes information. Touchscreen questionnaire diabetic patients were applied in this study rather than ICD10 diagnosed in order to increase patient numbers, given the small subgroup of individuals being examined.

Laboratory methods for GDF15 measurement FINRISK
Blood samples were collected after an advisory 4 hr fast, immediately centrifuged and then stored at -70°C until GDF15 measurement which was undertaken using an ILMA with a limit of detection of 24 ng/L and a linear range from 200 to 50,000 ng/L . The ILMA is technically identical to immunoradiometric assay (IRMA: Horn et al., 1996) except that the GDF15 detection antibody was labelled with acridinium ester and assay results were quantified in a luminometer (Berthold). ILMA and IRMA use an antibody that binds over a sequence of Ala197-Ile308. GDF-15 concentrations measured with the ILMA and IRMA are very similar [(ILMA concentration/ILMA concentration)×100% = 97.8% ± 1.3%] and closely correlated (r 2 =0.99, p<0.001, n=31 samples: H). Head-to-head comparison of IRMA with the clinical 'gold standard' Roche assay have been previously conducted and reported GDF15 levels were comparable .

INTERVAL
Blood sample collection procedures have been described previously (Moore et al., 2014). Blood samples were collected in 6 mL EDTA tubes and transferred at ambient temperature to UK Biocentre (Stockport, UK) for processing. Plasma was extracted by centrifugation into two 0.8 mL plasma aliquots and stored at -80°C before use. The procedures for obtaining protein measurements using the SomaScan assay have also been described previously (Sun et al., 2018). Briefly, SomaScan is a multiplexed, aptamer-based approach that was used to measure the relative concentrations of 3622 plasma proteins or protein complexes using modified aptamers ('SOMAmer reagents', hereafter referred to as SOMAmers). Aliquots of 150 µL of plasma from INTERVAL baseline samples were sent on dry ice to SomaLogic Inc (Boulder, CO) for protein measurement. Modified single-stranded DNA SOMAmers were used to bind to specific protein targets and these are then quantified using a DNA microarray. Protein concentrations are quantified as relative fluorescent units. Quality control (QC) was completed at the sample and SOMAmer levels.
The Olink assay uses pairs of monoclonal antibodies (or a single polyclonal antibody) to bind each protein target, and then uses proximity extension followed by qPCR to quantify protein abundance. Aliquots of plasma from samples taken from INTERVAL participants at the 2-year follow-up survey were shipped on dry ice to Olink Proteomics (Uppsala, Sweden) for assay. Protein concentrations were recorded as normalised relative protein abundances ('NPX').

Genotyping and imputation FINRISK
A total of 26,404 FINRISK samples were genotyped using several arrays: the HumanCoreExome BeadChip, the Human610-Quad BeadChip, the Affymetrix6.0, and the Infinium HumanOmniExpress (Illumina Inc, San Diego, CA). The present study, using samples taken in FINRISK in 1997, consisted of 6538 individuals which were genotyped using three genotyping arrays: the HumanCoreExome BeadChip, the Human610-Quad BeadChip, and the Infinium HumanOmniExpress (Illumina Inc, San Diego, CA). Genotype calls were generated together with other available datasets using zCall at the Institute for Molecular Medicine Finland (FIMM). After sample-wise QC (exclude samples with ambiguous gender, missingness [>5%], excess heterozygosity [±4 SD], non-European ancestry) and variantwise QC (exclude SNPs with high missingness [>2%]), low HWE p-value (<1 x 10 -6 ), minor allele count (MAC) < 3 (in case Zcall'ed chip data) or MAC < 10 (chip data called using Illumina GenCall) steps, the samples were pre-phased using Eagle2 (version 2.3). Genotype imputation was carried out by using a Finnish population-specific reference panel consisting of 2690 high-coverage WGS and 5092 WES samples with IMPUTE2 (version 2.3.2: Howie et al., 2009) that allows the usage of two panels at the same time (the 'merge_ref_panels' option). Post-imputation QC involved excluding variants imputed with imputation INFO < 0.7.

INTERVAL
The genotyping protocol and QC for the INTERVAL samples have previously been described in detail (Astle et al., 2016). DNA was extracted from buffy coat at LGC Genomics (UK) and was used to assay approximately 830,000 variants on the UK Biobank Affymetrix Axiom genotyping array at Affymetrix (Santa Clara, CA). Genotyping was performed in batches of approximately 4800 samples. Variants were excluded from a batch if they strongly deviated from HWE (p-value < 5 x 10 -6 ) or had a within-batch call rate < 0.97. Sample QC included removing duplicate samples and samples with non-European ancestry, missing phenotypic sex and sex mismatches, and extreme heterozygosity (±3 SD). Relatedness was removed by excluding one participant from each pair of close (first-or second-degree) relatives, defined as π > 0.187.
Additional variant QC steps were performed prior to imputation to establish a high-quality imputation scaffold. This included imposing a global HWE filter of p-value < 5 x 10 -6 , a call rate filter of 99% over the INTERVAL genotyping batches that a variant was not failed in, and a global call rate filter of 75% across all INTERVAL genotyping batches. All monomorphic variants, non-autosomal and multi-allelic variants were removed and 654,966 high-quality variants remained to be used for imputation. Phasing was conducted using SHAPEIT3 and variants were imputed using a combined 1000 Genomes Phase 3-UK10K imputation panel. Imputation was performed on the Sanger Imputation Server (https://imputation.sanger.ac.uk), resulting in 87,696,888 imputed variants.

Genome-wide association analysis FINRISK
For the FINRISK cohort genome-wide association analyses were performed for 5817 individuals with plasma GDF15 concentrations available. A normal distribution of GDF15 plasma concentrations was achieved through applying an inverse normal transformation. Multi-dimensional scaling was done for genetic data of both studies using PLINK version 1.07 (Purcell et al., 2007). Only good quality autosomal markers passing the following criteria: imputation informativeness > 0.7, and MAF > 0.001, were included in further evaluations. Frequentist test with method 'expected', assuming an additive genetic model, was performed using SNPTEST version 2 (Marchini et al., 2007). Results were adjusted for the first five principal components of the genetic data to account for the population stratification and for gender, age, and genotyping set. The QQ and Manhattan plots and regional plots were created using R-2.11 to visualise genome-wide association results. The genomic positions indicated throughout this study are based on NCBI human genome build 37.

INTERVAL
The GWAS for INTERVAL-SomaScan has previously been described in detail (Sun et al., 2018). Relative protein abundances were first natural log-transformed within each subcohort. Log-transformed protein levels were then adjusted in a linear regression for age, sex, duration between blood draw and processing (binary, ≤1 day/>1 day) and the first three principal components of ancestry from multi-dimensional scaling. A normal distribution of the protein residuals from this linear regression was achieved through rank-inverse normal transformation and these were used as phenotypes for association testing. Simple linear regression assuming an additive genetic model was used to test for genetic associations. Genotype uncertainty was accounted for by carrying out association tests on allelic dosages ('method expected' option) using SNPTEST version 2.5.2 (Mullican et al., 2017).
For the INTERVAL-Olink GWAS, normalised protein levels ('NPX') were first regressed on age, sex, plate, time from blood draw to processing (in days), and season (categorical: 'Spring', 'Summer', 'Autumn', 'Winter'). The residuals from this linear regression were rank-inverse normalised. The rankinverse normalised residuals were then fitted in a linear regression model adjusted for ancestry by including the first three components of multi-dimensional scaling as covariates. GWAS was conducted using SNPTEST version 2.5.2.

Fine mapping
Fine mapping was performed using the FINEMAP program (Benner et al., 2016). FINEMAP can potentially identify sets of variants with more evidence of being causal than those highlighted by a stepwise conditional analysis. The FINEMAP provides (1) a list of potential causal configurations together with their posterior probabilities and Bayes factors and, (2) for each variant, the posterior probability and Bayes factor of being causal. FINEMAP was applied with its default settings. LD between SNPs for the conditional analysis was combined using the formula: (N1 * LD1 + N2 * LD2 + N3 * LD3)/(N1 + N2 + N3), where N represents the number of individuals in the GWAS in each cohort, LD represents the R 2 value, and 1, 2, and 3 represent the three different cohorts: FINRISK, INTERVAL-SomaScan, and INTERVAL-Olink.

Variant annotation
The functional significance of the fine mapped variants were explored using public databases and browsers, including GTEx, Ensmbl, SNiPA, RegulomeDB (Võsa et al., 2021).

Meta-analysis
Meta-analysis was completed on the chromosome 19 locus around the GDF15 gene in METAL (Willer et al., 2010) using the IVW method. A total of 3,799 SNPs were found in common across FINRISK, INTERVAL-SomaScan, and INTERVAL-Olink.

Statistical analysis
Association study between GDF15 levels and disease endpoints and quantitative biomarkers To examine the association of GDF15 with disease endpoints and quantitative biomarkers, logistic and linear regression models were used, respectively. Inverse variance transformed GDF15 levels were used in the analyses due to right-skewed distribution. Association of GDF15 with disease endpoints and quantitative biomarkers were examined in (1) age and sex, (2) age, sex, and BMI, (3) age, sex, and smoking-adjusted models. Analyses using linear and logistic regression models were performed using R2.11 (see URLs). FDR was applied to test for multiple test correction.

Survival analysis
All analyses were performed in R version 3.4.0. Cox proportional hazard regression model was conducted to identify predictors of outcomes during 10 years. The survival curves for a Cox proportional hazards model were used to illustrate the timing of the death, type 2 diabetes, and CVD events during 10 years in relation to GDF15 quartiles, and statistical assessment between upper quartile and other quartiles was preformed using Cox proportional hazard regression model. To inspect the validity of the Cox model, the test of the proportional hazards assumption for a Cox regression model was used. Hosmer-Lemeshow goodness of fit test was used to determine whether models were calibrated with similar expected and observed event rates in both low-and high-risk individuals.

Wilcoxon test for prevalent and incident disease
Wilcoxon rank sum test (Mann-Whitney U test) was performed, in R using function ' wilcox. test' from stats package, to test whether the distributions of the incident disease case group, prevalent disease case group, and control group were systematically different from one another. Wilcoxon rank sum test was used for the prevalent and incident analysis given non-transformed GDF15 plasma levels (that were not normally distributed) were used.

Type 2 diabetes endpoint used for prevalent and incident analysis
The ICD codes used for type 2 diabetes endpoint definition in the incident prevalent analysis were: ICD9 codes: 2502A, 2501A, 2503A, 2504A, 2505A, 2506A, 2507A, 2508A, 2500A (diabetes mellitus, type 2); ICD10 codes: E110, E111, E112, E113, E114, E115, E11[6-8], E119 (adult type diabetes); ATCcode A10B (blood glucose lowering drugs, excluding insulins) from kela drug purchase register as well as ICD10 code E11 (adult type diabetes) from kela reimbursement register. During the ICD8 era, type 1 and type 2 diabetes were not separated in Finnish ICD8 diagnosis codes, and thus ICD8 diabetes diagnosis codes were not included to the type 2 diabetes endpoint. To ensure that type 1 diabetes cases are not mixed with type 2 diabetes cases, type 1 diabetes endpoint cases were removed from the type 2 diabetes endpoint cases.

Mendelian randomisation
In brief, MR is a method that explores whether an intermediate trait has a causal relationship with an endpoint by using genetic instruments. The use of genetic data means there is less environmental bias. Variants identified as significantly associated with the intermediate trait are applied and summary results from intermediate trait -variants and endpoint -variants are applied in the model. MR-Egger regression does not constrain the intercept, therefore it is not biased by invalid IV but has reduced power. MR relies on three main assumptions that (1) the instruments are associated with the exposure, (2) the instruments are not associated with any confounders, and (3) the instruments are associated with the outcome only through the exposure. Here, we applied random effect IVW analysis (which has the most statistical power but can be biased in the presence of horizontal pleiotropy), MR-Egger (which relies on the InSIDE assumption and has less statistical power but can be used to identify pleiotropy), weighted median MR (which is robust to outliers but sensitive to addition/removal of IVs), and MR-PRESSO (which is able to remove outliers but is prone to false positives in the presence of multiple invalid instruments). The strengths and weaknesses listed here are summarised in this paper (Burgess et al., 2019).
We applied MR (Bowden et al., 2015) to explore causality of GDF15 with BMI, WHR, glucose, type 2 diabetes, HDL cholesterol, and eBMD. We additionally explored these relationship in reverse: with GDF15 as an outcome and other traits as the exposure. GWAS summary statistics for GDF15 were taken from the conditional meta-analysis of FINRISK and INTERVAL for the forward MR and from FINRISK only for the reverse MR. GWAS in INTERVAL being completed on chromosome 19 only whereas GWAS in FINRISK were completed on all chromosomes. Summary statistics for the other traits were obtained from large public data resources: BMI (n=590,827: Yengo et al., 2018), WHR (n=485,486: Shungin et al., 2015), glucose (n=314,916: http://www.nealelab.is/uk-biobank), type 2 diabetes (n cases = 74,124, n controls = 824,006: Mahajan et al., 2018), HDL cholesterol (n=315,133: http://www.nealelab.is/uk-biobank), and eBMD (n=426,824: Morris et al., 2019). For BMI, MR with GDF15 as the exposure was applied using summary statistics from UK Biobank only (n=361,194) as the LD clumped variants were not available in the larger meta-analysis, that was utilised for reverse MR. We utilised LD clumping to identify independent genetic variants in summary statistics and applied these as genetic instruments in MR analyses. LD clumping was carried out in PLINK (version 1.9) with the following thresholds applied: R 2 < 0.01, p-value < 5 × 10 -8 and clump radius < 3000 kb. MR analysis was run in R (version 4.0.2) using the package 'MendelianRandomization'. The IVW method from this package was applied to test for causality and a sensitivity analysis using MR-Egger to testing for horizontal pleiotropy. A finding of significant causality is indicated with the IVW p-value and a significant MR-Egger intercept p-value indicates a finding of horizontal pleiotropy. Weighted median MR and MR-PRESSO (Mendelian randomisation pleiotropy residual sum and outlier: Verbanck et al., 2018) were additionally applied as a further sensitivity analysis and to test for horizontal pleiotropy and was run in R using the 'MRPRESSO' package. The Global p-value indicates a significant finding of horizontal pleiotropy and the software also tests for outliers, the finding of which instigates removal of outliers from the analysis. MR-PRESSO then re-runs IVW testing and gives a significant distortion p-value to indicate if the removal of outliers reduced the horizontal pleiotropy found in the analysis.
GDF15 PTV analysis UK Biobank exome sequencing UK Biobank was whole exome sequenced (paired-end 75 bp) at Regeneron Pharmaceuticals using the IDT xGen version 1 capture kit and NovaSeq6000 for more information, see previous publication (Wang et al., 2021b). Data was available on 302,362 individuals of which >95% of CCDS had at least ×10 coverage and average coverage of CCDS was ×59. Alignment to GRCh38 and SNV and indel calling were completed utilising Illumina DRAGEN Bio-IT Platform Germline Pipeline version 3.07 on a custom-built Amazon Web Services (AWS) cloud platform. Annotation of SNVs and indels was performed using SnpEFF v4.3 against Ensembl Build 38.92 and MAPQ < 30 were excluded. Related individuals were identified using KING and first-and second-degree relatives were excluded from downstream analysis. Sex and ancestry checks were completed using PEDDY (Pedersen and Quinlan, 2017) and individuals with gender mismatch were excluded.

Quality control
Individuals with GDF15 PTVs were identified by extracting variants annotated as PTVs (including exon loss variants, frameshift variants, start lost, stop gained, stop lost, splice acceptor variants, splice donor variant, rare amino acid variant, transcript ablation, gene fusion, and bidirectional gene fusion). Supplementary file 4j, demonstrates the frequency of LOF variants in each ancestry. From the individuals that did not carry GDF15 PTV variants 20,000 males and 20,000 females were randomly extracted as controls.

Mann-Whitney U test
Analysis was completed in R (version 3.2.4) using ' wilcox. test' and models compared BMI, WHR (in males and females separately), diabetic status (self-reported), eBMD, and HDL in GDF15 PTV carriers and non-carriers.
has ongoing research collaboration with Bayer Ltd (all outside this work). Adam S Butterworth: reports grants outside of this work from Biogen, BioMarin, Bioverativ, Merck, Novartis, Pfizer and Sanofi and personal fees from Novartis. Athena Matakidou: is an employee of AstraZeneca and currently an employee of GlaxoSmithKline (although not an employee of GlaxoSmithKline when this work was carried out). The other authors declare that no competing interests exist. The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.