Multi-ancestry genome-wide association study of gestational diabetes mellitus highlights genetic links with type 2 diabetes

1Gestational diabetes mellitus (GDM), defined as hyperglycaemia with onset or first recognition during pregnancy, is associated with increased risk of pregnancy complications and adverse perinatal outcomes, including pre-eclampsia, stillbirth, large for gestational age, neonatal hypoglycaemia, preterm birth, low Apgar scores and admission to neonatal intensive care. Whilst hyperglycaemia commonly resolves postpartum, GDM often reoccurs and is associated with subsequent diagnosis of type 2 diabetes (T2D) and coronary heart disease. Although the global prevalence of GDM is increasing, it varies according to population characteristics (such as maternal age, ancestry and obesity rates) and the criteria used for screening and diagnosis.

2GDM and T2D share both genetic and non-genetic risk factors, including obesity, poor diet and sedentary lifestyle. Family studies have demonstrated that women with GDM have 30.1% probability of having at least one parent with T2D, compared to just 13.2% for pregnant women with normal glucose tolerance. Furthermore, women with a history of GDM appear to have a nearly 10-fold higher risk of developing T2D than those with a normoglycaemic pregnancy. Taken together, these observations support the hypothesis that the two diseases are part of the same underlying pathology, with pregnancy potentially acting as a stress test that reveals women at increased risk of GDM and/or T2D.

3There have been considerable advances in our understanding of the genetic contribution to T2D through large-scale genome-wide association studies (GWAS) across diverse populations. In contrast, despite the observed familial clustering of GDM, most genetic association studies of the disease have focussed on evaluating the impact of previously reported loci for T2D and glycaemic traits in modest sample sizes. The most comprehensive systematic review of genetic susceptibility to GDM (from 23 studies) revealed association with T2D risk variants from seven loci, of which six are related to insulin secretion and one to insulin resistance. A genetic risk score (GRS) of risk variants across 34 loci associated with T2D and/or fasting glucose was significantly associated with GDM and improved predictive power over a model including only clinical variables. Variants associated with both insulin secretion and insulin resistance have also been used to construct an aggregated GRS that was shown to predict GDM risk, with and without adjustment for body mass index (BMI), maternal age and gestational age, although this score was not compared with established clinical predictors. To date, the largest GWAS of GDM has been undertaken in women from a Korean population, including 468 cases and 1242 non-diabetic controls in the discovery stage, with an additional 931 cases and 783 non-diabetic controls in the follow-up stage. Two loci were associated with GDM at genome-wide significance (P < 5 × 10−8), mapping near MTNR1B and CDKAL1, both of which have also been previously implicated in T2D risk.

4To gain novel insight into the genetic architecture of GDM, the GENetics of Diabetes In Pregnancy (GenDIP) Consortium assembled GWAS of diverse ancestry in a total of 5485 women with GDM and 347 856 women without GDM. With these resources, we aimed to improve our understanding of the aetiological factors and molecular processes driving the occurrence of GDM, including the extent to which these overlap with T2D pathophysiology, and investigate the effects of potential causal metabolic risk factors on the disease through Mendelian randomization (MR).

5We began by aggregating GDM association summary statistics across GWAS through multi-ancestry meta-analysis: the effective sample size was 72.2% European, 13.4% East Asian, 9.9% South Asian, 2.8% Hispanic/Latino and 1.7% African (Supplementary Material, Tables S1 and S2). To maximize sample size, we used a phenotype definition that makes best use of the information available in each study, including data from health records, oral glucose tolerance tests and self-report (Supplementary Material, Table S1). Each GWAS was imputed to reference panels from the 1000 Genomes Project, Haplotype Reference Consortium or population-specific whole-genome sequence data (Supplementary Material, Table S3). Within each GWAS, GDM association summary statistics were derived for all single nucleotide variants (SNVs) passing quality control after appropriate adjustment to account for population structure (Supplementary Material, Table S3).

6The most powerful methods for multi-ancestry meta-analysis allow for potential allelic effect heterogeneity on disease between population groups that cannot be accommodated in a fixed-effects model. Our primary analysis used MR-MEGA, which models heterogeneity between GWAS by including axes of genetic variation that represent ancestry as covariates in a meta-regression model. We considered three axes of genetic variation that separated the five ancestry groups, but which also revealed finer-scale genetic differences between GWAS of the same ancestry (Supplementary Material, Fig. S1). We also conducted multi-ancestry and ancestry-specific fixed-effects meta-analyses. We identified five loci at genome-wide significance in the multi-ancestry meta-regression (Table 1, Supplementary Material, Figs S2 and S3), including the previously reported associations from GDM GWAS at MTNR1B (rs10830963, P = 4.3 × 10−54) and CDKAL1 (rs9348441, P = 1.6 × 10−14). The remaining three loci for GDM mapped to/near TCF7L2 (rs7903146, P = 4.0 × 10−16), CDKN2A-CDKN2B (rs10811662, P = 4.1 × 10−9) and HKDC1 (rs9663238, P = 2.9 × 10−8). Through approximate conditional analyses, conducted using ancestry-matched linkage disequilibrium (LD) reference panels for each GWAS, we observed no evidence for multiple distinct association signals at genome-wide significance at any of the five GDM loci (Supplementary Material, Fig. S4).

7We next sought to investigate the impact of differences in ancestry and phenotype definition between GWAS on heterogeneity in allelic effects at GDM loci. To do this, we extended the MR-MEGA meta-regression model to include an additional covariate to represent whether GDM status in the study was confirmed via ‘a universal blood-based test’ (Supplementary Material, Table S1). Here, we use this term to refer to a blood-based test that was applied to all participants, including a diagnostic oral glucose tolerance test (OGTT) or a screening glucose challenge or fasting glucose test, in contrast to clinician decision, risk factor screening or a lack of clarity on what basis women did or did not have a diagnostic OGTT. This model enables partitioning of heterogeneity into three components (Table 2). The first component captures heterogeneity that is correlated with genetic ancestry (that can be explained by the three axes of genetic variation), which can occur because of differences in the structure of LD between ancestry groups or interactions with lifestyle factors that vary across populations. The second component measures heterogeneity that can be explained by the use of a universal blood-based test to screen for or diagnose GDM. The final component reflects residual heterogeneity due to study design that cannot be explained by the first two components. The greatest evidence of ancestry-correlated heterogeneity (after accounting for the use of a universal blood-based test) was observed at the CDKAL1 locus (PHET = 3.4 × 10−5), where the lead SNV demonstrated stronger effects on GDM in GWAS of East Asian ancestry than in other populations, despite the risk allele being common in all ancestry groups (Supplementary Material, Fig. S5, Table S4). A similar pattern of ancestry-correlated heterogeneity in allelic effects on T2D has been reported at the CDKAL1 locus. Weaker evidence of ancestry-correlated heterogeneity was observed at the CDKN2A-CDKN2B locus (PHET = 0.0022), where there were marked differences in the effects on GDM of the lead SNV between GWAS undertaken in different ancestry groups (Supplementary Material, Fig. S5, Table S4). In contrast, there was no evidence of heterogeneity due to phenotype definition for any lead SNV, suggesting that differences in allelic effects between GWAS are more likely due to factors related to genetic ancestry than the use of a blood-based test in all women to screen for or diagnose GDM.

8Of the five GDM loci identified at genome-wide significance in the trans-ancestry meta-regression, four have been previously implicated in T2D susceptibility: MTNR1B, TCF7L2, CDKAL1 and CDKN2A-CDKN2B. In fact, in previously reported trans-ancestry GWAS meta-analyses of 180 834 T2D cases and 1 159 055 controls from the DIAMANTE Consortium, the lead T2D SNV is the same as we report for GDM at MTNR1B, TCF7L2 and CDKAL1, and is in strong linkage disequilibrium (LD) at the CDKN2A-CDKN2B locus (rs10811661, r2 = 0.91 across diverse populations in the 1000 Genomes Project). To further investigate the genetic correlation between the two diseases, we extracted GDM association summary statistics from our trans-ancestry meta-analysis for lead SNVs at 222 previously reported loci for T2D from the DIAMANTE Consortium (Supplementary Material, Fig. 1, Table S5). We observed a strong positive correlation in log-ORs for the T2D risk allele between the two diseases: Pearson r = 0.573 (P < 2.2 × 10−16). There was also a highly significant enrichment of GDM associations at T2D loci (50 of 222 lead SNVs with P < 0.05 and same direction of effect, binomial test P < 2.2 × 10−16), indicating that they would be discovered at genome-wide significance with larger effective sample sizes. Indeed, after excluding the four overlapping GDM-T2D loci, a weighted genetic risk score of lead T2D SNVs was significantly associated with GDM (P = 9.7 × 10−123, pseudo-R2 = 2.86%). Extending our analyses, genome-wide, using LD-score regression, we observed strong genetic correlation between GDM and T2D: rG (95% CI) 0.744 (0.052, 1.437). Weaker genetic correlations between GDM and other glycaemic traits, BMI and birth weight were also observed (Table 3, Supplementary Material, Table S6). These results are consistent with sharing of genetic determinants of GDM and T2D, although we acknowledge that LD score regression has limited statistical power because of the relatively small GDM sample size, and we note that the correlation from LD-score regression is not bound by −1 to 1, particularly when power is low.

9The most obvious difference in allelic effect sizes between GDM and T2D was observed at the MTNR1B locus (Fig. 1). The lead SNV, rs10830963, is the same for both diseases, but the allelic effect on GDM is substantially greater than on T2D: OR (95% CI) for GDM is 1.41 (1.35–1.47) and for T2D is just 1.09 (1.08–1.10). The MTNR1B lead SNV is associated, at genome-wide significance, with fasting glycaemic traits in non-diabetic individuals from the Meta-Analysis of Glucose and Insulin-related traits Consortium (MAGIC) Investigators. The GDM risk allele at the lead SNV is also associated with higher fasting plasma glucose and 1-hour plasma glucose in pregnant women from the Hyperglycemia and Adverse Pregnancy Outcomes (HAPO) Study. This SNV also has the strongest association of the maternal glucose-raising allele with higher offspring birth weight in women from the Early Growth Genetics Consortium, in line with the known effects of maternal hyperglycaemia on fetal growth. In non-diabetic individuals from the MAGIC Investigators, the MTNR1B lead SNV has a much larger impact on fasting glucose than those at TCF7L2, CDKAL1 and CDKN2A-CDKN2B (Supplementary Material, Table S7). Therefore, the difference in allelic effect sizes between GDM and T2D at MTNR1B may reflect the fact that thresholds of fasting plasma glucose used to diagnose GDM are lower than those used to diagnose T2D, meaning that a larger proportion of GDM than T2D cases will have higher fasting glucose that is regulated within the normal range.

10To gain insight into the molecular processes and tissues through which GDM association signals are mediated, genome-wide, we then undertook fGWAS enrichment analyses within three categories of functional and regulatory annotations: (i) genic regions; (ii) chromatin immuno-precipitation sequence (ChIP-seq) binding sites for 165 transcription factors and (iii) 13 unique and recurrent chromatin states in four diabetes-relevant tissues (pancreatic islets, liver, adipose, and skeletal muscle). We observed significant joint enrichment (P < 0.05) for GDM associations mapping to protein coding exons, binding sites for FOXA2, NFE2 and TFAP2, and chromatin states in adipose tissue and skeletal muscle that mark enhancers and transcribed regions (Supplementary Material, Table S8). FOXA2 is a pioneer factor involved in pancreatic and hepatic development, and T2D association signals have been previously reported to be enriched for FOXA2-binding sites. Skeletal muscle is the most prominent site of insulin-mediated glucose uptake in humans, and enhancers in skeletal muscle have been reported to overlap association signals for metabolic disorders, including T2D, insulin resistance and obesity. These enrichment analyses highlight molecular processes and tissues that are broadly consistent with those important in mediating T2D association signals, although the involvement of pancreatic islets appears to be less prominent for GDM.

11In contrast to the other GDM loci reported in this investigation, the lead SNV at the HKDC1 locus (rs9663238) demonstrates only weak statistical evidence of T2D association in previously reported trans-ancestry GWAS meta-analyses from the DIAMANTE Consortium (P = 0.0083, compared with P < 10−65 at the other four loci). GDM risk alleles at the lead SNV, and/or at variants in strong LD (European ancestry r2 > 0.9) with it, have been previously associated, at genome-wide significance, with higher 2-h plasma glucose (2HPG) in pregnant women in the HAPO Study and two replication studies of European ancestry, as well as with higher birth weight of first child (likely via greater maternal glucose availability), higher own birth weight (fetal effect independent of the maternal effect on birth weight) and comparative height and body size at age 10 in UK Biobank (Supplementary Material, Table S9). The lead SNV is also associated, more strongly in women than men, with higher alanine aminotransferase (ALT) levels in UK Biobank. Elevated ALT levels in early pregnancy have been associated with the risk of subsequent development of GDM and genome-wide, we observed positive genetic correlation between the two traits: rG (95% CI) 0.149 (0.005, 0.292). In addition to demonstrating the association of the maternal SNVs at this locus with GDM in the current study, we observed that 99% credible set variants are lead SNVs for HKDC1 expression quantitative trait loci in a range of tissues in the GTEx Project, including visceral adipose, subcutaneous adipose and pancreas (Supplementary Material, Table S10). HKDC1 (Hexokinase Domain Containing 1) catalyses the phosphorylation of hexose to hexose 6-phosphate and is involved in glucose homeostasis and hepatic lipid accumulation. Haplotypes of variants associated with 2HPG in pregnancy disrupt regulatory element activity and reduce HKDC1 expression across diverse tissues (including metabolically relevant liver stellate cells and pancreatic islet beta cells), which has been demonstrated to reduce hexokinase activity in multiple cellular models. Knockout of hepatic HKDC1 in pregnant mice has also been demonstrated to significantly impair glucose tolerance, highlighting the importance of liver HKDC1 on glucose metabolism during pregnancy. Taken together, the evidence from our study and others suggests a more important role for HKDC1 in glucose metabolism during pregnancy than outside of pregnancy, in addition to independent maternal and offspring effects on early growth, and highlights that while GDM shares many similarities with T2D, there are differences in at least one underlying pathway.

12Finally, we used two-sample MR to investigate causal effects on GDM of 282 metabolic measures and risk factors available in the MR-Base GWAS catalogue (www.mrbase.org), including metabolites, anthropometric measures, hormones, immune system phenotypes, kidney traits and metals (Supplementary Material, Table S11). We did not consider glycaemic traits (including HbA1c) because they are used to define GDM status. For each metabolic measure, we selected independent SNVs attaining genome-wide significance with the trait as instrumental variables. For each SNV, we extracted association summary statistics for GDM from the European ancestry-specific meta-analysis because we assessed independence of genetic instruments using LD from European ancestry haplotypes from the 1000 Genomes Project. Of the 282 exposures considered, only BMI demonstrated significant evidence for a causal effect on GDM risk at a false discovery rate of 5% (Supplementary Material, Table S11). The estimated causal effect of higher BMI on higher GDM risk was directionally consistent across multiple MR models (Fig. 2). The causal relationship of BMI with GDM is consistent with its effect on T2D.

13We have conducted the largest and most ancestrally diverse GWAS meta-analysis for GDM, where we identified associations mapping to MTNR1B, TCF7L2, CDKAL1, CDKN2A-CDKN2B and HKDC1. Our results demonstrated strong correlation in the effects of previously reported associations for T2D and those observed for GDM, and highlighted overlapping molecular mechanisms and tissues that mediate associations for both diseases. In contrast, variation at the HKDC1 locus is not strongly associated with T2D, but instead plays a more important role in glucose metabolism during pregnancy than outside of pregnancy. The genetic diversity of GWAS contributing to our meta-analysis enabled identification of ancestry-correlated heterogeneity in allelic effects on GDM at two loci. Such heterogeneity could reflect variable impact of different pathophysiology driving glycaemic dysregulation in pregnancy between ancestries and emphasizes the need for increased sample sizes in under-represented population groups. In contrast, results were consistent between GWAS in which all women had a universal blood-based test and those that did not, suggesting little impact from misclassification due to selective use of diagnostic tests only in those deemed to be at high-risk. Finally, MR analyses revealed a significant causal effect of higher BMI on GDM risk, consistent with the causal association observed with T2D. Taken together, these results provide further support for the hypothesis that T2D and GDM are part of the same underlying pathology. However, they also highlight there are pathways to GDM that impact on glucose regulation only in pregnancy, and that additional GDM-specific associations will be revealed through GWAS in larger sample sizes.

14All human research was approved by the relevant institutional review boards and conducted according to the Declaration of Helsinki. All participants provided written informed consent.

15Individuals were assayed with a range of GWAS genotyping arrays, with sample and SNV quality control undertaken within each study (Supplementary Material, Tables S2 and S3). Samples were pre-phased and imputed up to reference panels from the 1000 Genomes Project (phase 1, March 2012 release; phase 3, October 2014 release), Haplotype Reference Consortium or population-specific whole-genome sequencing (Supplementary Material, Table S3). SNVs with poor imputation quality (r2 < 0.3 or info < 0.4) and/or minor allele count < 5 were excluded from downstream association analyses (Supplementary Material, Table S3). Association with GDM was evaluated in a regression framework, under an additive model in the dosage of the minor allele, with adjustment for principal components and other study-specific covariates to minimize the population stratification effects (Supplementary Material, Table S3). Phenotype definition and covariate adjustments were not harmonized between GWAS because of differences in individual study design and availability of non-genetic risk factor information. Analyses accounted for structure (population stratification and/or familial relationships) by: (i) excluding related samples and adjustment for principal components derived from a genetic relatedness matrix (GRM) as additional covariates in the regression model; or (ii) incorporating a random effect for the GRM in a mixed model (Supplementary Material, Table S3). Allelic effects and corresponding standard errors that were estimated from a linear (mixed) model were converted to the log-odds scale. Study-level association summary statistics (P-values and standard error of allelic effects) were corrected for residual structure by means of genomic control if the inflation factor was > 1 (Supplementary Material, Table S3).

16To account for the different reference panels used for imputation across GWAS, we restricted our analyses to autosomal bi-allelic SNVs from the 1000 Genomes Project reference panel (phase 3, October 2014 release) that are also present in the Haplotype Reference Consortium reference panel. We considered only those SNVs with MAF > 0.5% in haplotypes in at least one of the five ancestry groups represented in the 1000 Genomes Project (phase 3, October 2014 release).

17Our primary multi-ancestry analysis utilized meta-regression, implemented in the MR-MEGA software, which allows for allelic effect heterogeneity between GWAS that is correlated with ancestry. We first constructed a distance matrix of mean effect allele frequency differences between each pair of GWAS across a subset of SNVs reported in all studies. We implemented multi-dimensional scaling of the distance matrix to obtain three principal components that defined axes of genetic variation to separate GWAS from the five ancestry groups (Supplementary Material, Fig. S1). For each SNV, we then modelled allelic log-ORs across GWAS in a linear regression framework, weighted by the inverse of the variance of the effect estimates, incorporating the three axes of genetic variation as covariates. Under this model, we tested for association with GDM allowing for allelic effect heterogeneity between GWAS that is correlated with ancestry. We corrected the meta-regression association P-values for inflation due to residual structure between GWAS using genomic control adjustment. We considered only those SNVs reported ≥50% of the total effective sample size in downstream analyses.

18For each SNV, we also conducted fixed-effects meta-analysis across GWAS under an inverse-variance weighting of allelic log-ORs using GWAMA. We corrected standard errors of the resulting effect estimates for inflation due to residual structure between GWAS by genomic control adjustment.

19We identified lead SNVs attaining genome-wide significant evidence of association (P < 5 × 10−8) in the multi-ancestry meta-regression that were separated by at least 500 kb. Loci were defined by the genomic interval mapping 500 kb up- and downstream of each lead SNV.

20Each GWAS was first assigned to one of the ancestry groups (Supplementary Material, Table S2) represented in the 1000 Genomes Project reference panel (phase 3, October 2014 release). Haplotypes in the panel that were specific to that ancestry group were used as a reference for LD between SNVs across loci for the GWAS in approximate conditional analyses implemented in GCTA. For each locus, we applied GCTA to each GWAS to condition on the lead SNV at the locus, using the study-level association summary statistics and matched LD reference. Allelic log-ORs from the approximate conditional analyses across GWAS were modelled in the multi-ancestry meta-regression framework implemented in MR-MEGA, incorporating the three axes of genetic variation as covariates, and weighted by the inverse of the variance of the effect estimates. The meta-regression association P-values were corrected for inflation due to residual structure between GWAS by using the same genomic control adjustment as in the unconditional analysis. If no SNVs attained genome-wide significant (P < 5 × 10−8) evidence of residual GDM association in the meta-regression, we concluded that there was a single association signal at the locus.

21We aggregated association summary statistics across GWAS from the same ancestry group via fixed-effects meta-analysis based on inverse-variance weighting of allelic log-OR to obtain effect size estimates using GWAMA. We corrected association P-values and standard errors of allelic effects from each ancestry group for residual inflation due to structure between GWAS by genomic control adjustment if the inflation factor was > 1. We estimated the mean effect allele frequency across GWAS from each ancestry group, weighted by the effective sample size of the study.

22We extended the meta-regression model implemented in the MR-MEGA software to investigate the impact of ancestry and the use of a universal blood-based test to define GDM status on heterogeneity in allelic effects on GDM at lead SNVs. We modelled allelic log-ORs across GWAS in a linear regression framework, weighted by the inverse of the variance of the effect estimates, incorporating a covariate indicating whether GDM status was defined by a universal blood-based test (Supplementary Material, Table S1) in addition to the three axes of genetic variation.

23We considered lead SNVs at 237 previously reported loci for T2D from the DIAMANTE Consortium obtained from a multi-ancestry meta-analysis of 180 834 cases and 1 159 055 controls (48.9% non-European ancestry). For each of the 222 SNVs that were reported in our multi-ancestry meta-analysis, we compared association summary statistics (risk allele, other allele, log-OR and P-value) for GDM and those reported for T2D. We excluded lead SNVs for T2D that also attained genome-wide significance for GDM. For the remaining SNVs, we regressed the log-ORs for GDM (weighted by their corresponding variances) on the log-OR for T2D, as implemented in grs. summary function of the gtx package in R version 3.4.2. We estimated the percentage of GDM variance explained, as measured by pseudo R2.

24We used LD Hub to perform LD score regression of the European ancestry association summary statistics for GDM on other glycaemic traits. We included T2D, fasting glucose, fasting insulin, fasting proinsulin, glucose 2 h post oral glucose tolerance test (adjusted for BMI), HbA1c, HOMA-B, HOMA-IR, BMI, birth weight and alanine aminotransferase. European ancestry association summary statistics for GDM were filtered so that only SNVs with minor allele frequency >0.01 was included before performing the LD score regression. Genetic correlations between the different glycaemic traits were obtained from the LD Hub lookup centre. Visualization was performed using the R package ggplot2 in R version 3.6.1.

25We mapped each SNV across the genome to three categories of functional and regulatory annotations. First, we considered genic regions, as defined by the GENCODE Project, including protein-coding exons, and 3′ and 5′ UTRs as different annotations. Second, we considered chromatin immuno-precipitation sequence (ChIP-seq) binding sites for 165 transcription factors: 161 proteins from the ENCODE Project and four additional factors assayed in primary pancreatic islets. Third, we considered 13 unique and recurrent chromatin states, including promoter, enhancer, transcribed, and repressed regions, in four diabetes-relevant tissues: pancreatic islets, liver, adipose and skeletal muscle. This resulted in a total of 220 genomic annotations for enrichment analyses.

26We tested for genome-wide enrichment of GDM associations that map to genomic annotations using fGWAS. To do this, we approximated the Bayes’ factor in favour of GDM association for the th SNV by,where is the deviance across Kj contributing GWAS contributing to the multi-ancestry meta-regression. We first considered each annotation separately and identified those with significant enrichment (P < 0.05). We then used an iterative approach to identify a joint model of enriched annotations from this set. At each iteration, we dropped the annotation from the joint model that minimized the reduction in the penalized likelihood. We continued until no additional annotations worsened the fit of the joint model at nominal significance (P < 0.05). We next used the cross-validation likelihood because the significance of parameter estimates from the penalized likelihood cannot be assessed using standard statistical approaches. For the selected joint model, we identified the penalty that maximized the cross-validation likelihood. Finally, we dropped any annotations from the joint model that resulted in a decrease in the cross-validation likelihood.

27At the HKDC1 locus, we calculated the posterior probability of driving the GDM association for each SNV under an annotation-informed prior model derived from the globally enriched functional and regulatory annotations identified by fGWAS. Specifically, for the th SNV at the locus, the posterior probability, where is the Bayes’ factor in favour of GDM association from the meta-regression, derived above. In this expression, the relative annotation-informed prior for the SNV is given by,where the summation is over the enriched annotations, is the estimated log-fold enrichment of the th annotation from the final joint model and is an indicator variable taking the value 1 if the th SNV maps to the th annotation, and 0 otherwise. We derived a 99% credible set for the locus by: (i) ranking all SNVs according to their posterior probability; and (ii) including ranked SNVs until their cumulative posterior probability attained or exceeded 0.99.

28We conducted a look-up of 99% credible set variants at the HKDC1 locus for significant (q < 0.05) cis-expression quantitative trait loci (eQTLs) across tissues in the GTEx Project. We reported only those 99% credible variants that were the lead SNV for the eQTL signal.

29We systematically searched the MR-Base GWAS catalogue (https://www.mrbase.org) for metabolic measures. We selected all subcategories of metabolites, which included ‘amino acid’, ‘carbohydrate’, ‘cofactors and vitamins’, ‘energy’, ‘fatty acid’, ‘keto acid’, ‘lipid’, ‘metabolite salt’, ‘metabolites ratio’, ‘NA’, ‘nucleotide’, ‘peptide’, ‘protein’, ‘unknown metabolite’ and ‘xenobiotics’. We also selected the following subcategories of risk factors: ‘anthropometric’, ‘hormone’, ‘immune system’, ‘kidney’ and ‘metal’. We identified European ancestry GWAS in MR-Base for each selected metabolic trait. Where more than one GWAS was available for a trait, we gave preference to women-specific studies with the largest sample sizes and numbers of SNVs. Any GWAS undertaken only in men were excluded.

30For each metabolic trait with more than five genetic instruments, we conducted MR analyses using a ‘mixture of experts’ (MoE) machine learning approach. This approach maximizes statistical power whilst minimizing the impact of horizontal pleiotropy by combining four instrument selection approaches to 14 different MR models. The four approaches for selecting genetic instruments using MoE were: (i) ‘top hits’ corresponding to independent variants associated at genome-wide significance (P < 5 × 10−8, r2 < 0.001 using 1000G CEU as the reference population); (ii) ‘directional filtration’ that removed instruments from ‘top hits’ that are likely to be related to the outcome through reverse causation using Steiger filtering; (iii) ‘heterogeneity filtering’ that removed instruments from ‘top hits’ that make a substantial contribution to Cochran’s Q statistic (P < 0.05); and (iv) combined ‘directional filtration’ and ‘heterogeneity filtering’. The 14 MR models were: seven mean-based methods (inverse variance weighting with fixed effects, IVW random effects, MR-Egger fixed effects, MR-Egger random effects and the three Rucker estimates), three median-based methods (simple, weighted and penalized median estimator) and four mode-based methods (simple and weighted mode, each weighted with or without the assumption of no measurement error in the exposure estimates). The best combination of instrument selection-MR method was identified using a variable predicted by MoE, scaled between 0 and 1, where 1 indicates best performance.

31For metabolic traits with five or fewer genetic instruments, the MoE approach could not be applied because many of the MR models require larger numbers of SNVs. For these metabolic traits, we used either the Wald ratio estimate (one SNV) or the inverse-variance weighted estimate (between two and five SNVs).

32All analyses were conducted in R version 3.6 using the packages ‘TwoSampleMR’ (version 0.5.4) and ‘MRInstruments’ for the MR analyses and ‘EpiCircos’ (https://github.com/mattlee821/EpiCircos).

33Meta-analysis summary statistics can be downloaded from: https://tools.gi.ut.ee/tools/GENDIP_PervjakovaEtAl2022.txt.gz.