Genome-wide association study of 1,391 plasma metabolites in 6,136 Finnish men identifies 303 novel signals and provides biological insights into human diseases

GWAS on 1,391 metabolites

We assayed 1,544 plasma metabolites using the Metabolon DiscoveryHD4 mass spectrometry platform (Supplementary Tables 1-2) in 6,136 randomly-selected METSIM participants who were non-diabetic at baseline and passed quality control (QC) (Supplementary Table 3; Supplementary Fig. 1). 1,391 metabolites were successfully quantified in ≥500 of these 6,136 participants. We created a METSIM imputation reference panel of >26M genetic variants by integrating genome and exome sequence and array genotypes in 2,922 METSIM participants (Methods; Supplementary Table 4). We used this reference panel to impute genotypes in all METSIM participants.

We carried out GWAS across >16M variants with imputation r²≥0.3 and minor allele count (MAC)≥5 in the 6,136 METSIM participants for the 1,391 (correlated) metabolites (Methods; Supplementary Table 5; Supplementary Fig. 1). Single-variant association tests identified 305,555 associations at 109,368 variants for 803 metabolites at P<7.2×10⁻¹¹=5.0×10⁻⁸/692 (Bonferroni correction for 692 principal components that together explained 95%¹⁵ of phenotypic variance for the 1,391 correlated metabolites; Methods). The GWAS p-values for each metabolite were well calibrated (genomic control inflation factor median=1.00, range=0.92-1.07; Supplementary Fig. 2). We built a multi-phenotype GWAS browser (PheWeb) (https://pheweb.org/metsim-metab/) to visualize and make publicly available our results for all 1,391 GWAS (Fig. 1; see Discussion).

• Download figure
• Open in new tab

Figure 1:

Our METSIM Metabolomics PheWeb facilitates the characterization of genetic associations and gene activities. a) Manhattan plot for N-acetylkynurenine highlighting the roles of the associated genes (https://pheweb.org/metsim-metab/pheno/C100006378). Chemical structure for N-acetylkynurenine and activities for the associated genes are added manually on top of the Manhattan plot. b) Stacked PheWeb plots show significant associations between rs6705977 (NAT8, https://pheweb.org/metsim-metab/variant/2:73622043-C-G) and fifteen N-acetylated molecules, and the more restricted set of associations between rs948445 (ACY3, https://pheweb.org/metsim-metab/variant/11:67647021-C-T) and four N-acetylated aromatic amino acids. LI: lipid; XE: xenobiotics; AA: amino acid; CA: carbohydrate; NU: nucleotide; PE: peptide; CV: cofactor and vitamin; EN: energy; PC: partially characterized; UN: unnamed.

Since body mass index (BMI) influences levels of many metabolites¹⁸, we repeated all 1,391 GWAS with BMI as an additional covariate. Results with and without BMI adjustment were generally very similar, with Pearson correlation coefficient r=0.999 for effect size estimates and -log₁₀p-values for variant-metabolite pairs with P<7.2×10⁻¹¹ in either of the two analyses (Supplementary Fig. 3). Supplementary Table 6 lists the 83 associations with substantially different effect sizes (ratio≥1.20) with and without BMI adjustment. In what follows, we present results for analyses without BMI adjustment.

Detecting independent association signals

To identify (nearly) independent association signals, we carried out chromosome-wide stepwise conditional analysis for each chromosome-metabolite pair with ≥1 association at P<5.0×10⁻⁸. Conditional analysis identified 2,030 association signals at 1,143 index variants for 803 metabolites at P<7.2×10⁻¹¹ (Table 1; Supplementary Table 7; Supplementary Figs. 4-5). The 1,143 index variants were of high imputation quality (r² median=0.99, range=0.63-1.00). 311 (27.2%) of the 1,143 index variants were associated (P<7.2×10⁻¹¹) with ≥2 metabolites, suggesting widespread pleiotropy (Supplementary Fig. 6). Among the 1,143 index variants, 121 (for 125 metabolites) are rare in METSIM and 99 (for 148 metabolites) have MAF>10-fold greater in METSIM than in NFE (gnomAD v3.1); 58 of these variants are both rare and enriched in Finns (Fig. 2a; Supplementary Table 7).

• Download figure
• Open in new tab

Figure 2:

Characterization of the 2,030 significant metabolite genetic association signals in the 6,136 METSIM participants with Metabolon metabolomics data. Comparison of MAFs for the 1,143 index variants between METSIM and non-Finnish Europeans in gnomAD v3.1; index variants are colored a) purple if MAF>10-fold greater in METSIM than in non-Finnish Europeans; or b) blue they represent novel association signals. The dashed line is of slope one through the origin. c) Overlaid Manhattan plots of the 1,391 metabolite GWAS. The red dashed line depicts P=7.2×10⁻¹¹. The associations at 40 novel putative causal genes within novel regions (blue) and 18 novel putative causal genes within previously reported regions (maize) are highlighted. The seven novel putative causal genes implicated only by fine-mapping analysis are starred. HADHA/B represents the HADHA and HADHB genes and ARSD/L the ARSD and ARSL genes.

Index variants explained from 0.7% to 62.0% (median=1.4%; Supplementary Fig. 7) of the phenotypic variance of the corresponding metabolite; 99 index variants explained ≥10% of the variance (Supplementary Table 8), including three missense variants with >10-fold greater frequency in METSIM than in NFE. For example, the putatively-deleterious AFMID missense variant p.Ala41Pro (rs77585764; MAF=5.4% in METSIM vs. 0.38% in NFE) explained 15.5% of the variance in N-formylanthranilic acid. AFMID encodes arylformamidase, an enzyme which catalyzes N-formylanthranilate to produce anthranilate and formate¹⁹.

Fine mapping

To fine map the causal variants for the 2,030 association signals, we created 2 Mb regions centered on each index variant and merged overlapping regions associated with the same metabolite, resulting in 1,501 regions. We used Bayesian fine-mapping²⁰ with a uniform prior to calculate the variant posterior inclusion probability (VPIP) that each variant is causal and the signal posterior inclusion probability (SPIP), the sum of the VPIPs for the variants in a region (Methods). This method can identify multiple independent signals in a region. In the 1,501 regions, we identified 2,435 signals with SPIP≥0.95, 1,952 of which are among the 2,030 association signals identified in conditional analysis. We built 95% credible sets²¹ of potential causal variants for the 1,952 signals (Fig. 3a; Supplementary Table 9). These credible sets included 1 to 544 variants (median=6); notably, 334 credible sets included only one variant (168 distinct variants).

• Download figure
• Open in new tab

Figure 3:

The 1,952 of the 2,030 metabolite genetic association signals identified in stepwise conditional tests with SPIP_≥0.95 in DAP-g Bayesian fine mapping. a) Numbers of variants in the 95% credible sets and distribution of variant posterior inclusion probabilities (VPIPs) for the most likely causal variants within the 95% credible sets. b) Density plot of largest VPIPs highlights the variants with >10-fold greater frequency in METSIM than non-Finnish Europeans (gnomAD v3.1; blue) have larger VPIPs than all other variants (maize).

In each of the 1,952 credible sets, we identified the variant with the largest VPIP. This list comprised 1,119 distinct variants, 100 with MAF>10-fold greater in METSIM than in NFE. VPIPs for these 100 variants were greater than those for the remaining 1,019 (VPIP mean=0.73 vs. 0.47; t-test P=1.7×10⁻²¹; Fig. 3b).

Of the 1,119 variants, 263 had VPIP≥0.8 in 547 credible sets. Among these 263 variants, 46 are rare in METSIM and 47 have MAF>10-fold greater in METSIM than in NFE; 28 of these variants are both rare and enriched in Finns (Supplementary Table 9). The 263 variants include 11 protein-truncating (PTV) and 69 missense variants across 66 genes, and 183 other (mostly non-coding) variants (Supplementary Table 9). Given their likely impact on gene function, we focused on the 80=11+69 PTV and missense variants, which suggested causal roles for the corresponding 66 genes. These 80 variants had VPIP≥0.8 in credible sets for 208 signals with 173 metabolites. Among the 80 variants, 26 (5 PTV and 21 missense) are rare and 30 (6 PTV and 24 missense) have MAF>10-fold greater in METSIM than in NFE; 16 of these variants are both rare and enriched in Finns.

Identifying novel association signals at rare and Finnish enriched variants

To determine which of the 2,030 association signals are distinct from previous metabolite GWAS findings, we repeated metabolite association analysis conditioning on all variants that were (a) ≤1 Mb of the index variant and (b) previously reported as associated with any metabolite in a curated list of 381 publications (Methods; Supplementary Table 10). 303 association signals at 229 index variants remained significant for 201 metabolites (P_condition<7.2×10⁻¹¹; Fig. 2b; Supplementary Table 7). The 303 novel signals included 64 signals (for 58 metabolites) at 51 rare variants and 79 signals (for 71 metabolites) at 47 variants >10-fold more frequent in METSIM than in NFE (Table 2); 33 of these signals are at variants both rare and enriched in Finns (Supplementary Table 7). In addition, 17 signals for 16 metabolites are on the X chromosome, and 78 signals are for 44 metabolites identified since 2015 on the Metabolon DiscoveryHD4 platform.

Multiple novel associations arose at the same index variants. For example, we identified novel association signals with 19 metabolites at the putatively-deleterious SLC23A3 missense variant p.Asn336Lys (rs192756070; Supplementary Fig. 8). p.Asn336Lys has 107-fold greater frequency in METSIM than in NFE (MAF=2.3% vs. 0.022%) and is likely the causal variant for most or all 19 metabolite associations (VPIP median=0.98, range=0.57-1.00). SLC23A3 encodes an SLC23 ascorbic acid transporter without demonstrated nucleobase transport²². These novel associations suggest a wide range of transport functions for SLC23A3.

Among the novel association signals at index variants enriched in Finns, we identified an association with 3-amino-2-piperidone at the putatively-deleterious OAT missense variant p.Leu402Pro (rs121965043, β=1.91, P=3.7×10⁻³⁵). p.Leu402Pro has 100-fold greater frequency in METSIM than in NFE (MAF=0.35% vs. 0.0031%) and is the likely causal variant for this association (VPIP=0.997). OAT encodes the key mitochondrial enzyme ornithine aminotransferase which converts arginine and ornithine into glutamate and gamma aminobutyric acid²³. OAT has not previously been implicated in metabolite GWAS, but inactivation of OAT is responsible for the Finnish heritage disease gyrate atrophy characterized by hyperornithinemia²⁴. Previous studies have found increased 3-amino-2-piperidone levels in the urine of individuals with gyrate atrophy²⁵.

Among the novel association signals on the X chromosome, we identified an association for tiglylcarnitine at the putatively-deleterious HSD17B10 missense variant p.Ala95Thr (rs201378370, β=0.94, P=5.2×10⁻¹²²). p.Ala95Thr has 76-fold greater frequency in METSIM than in NFE (MAF=2.6% vs. 0.034%) and is the likely causal variant for this association (VPIP=0.997). HSD17B10 encodes 17-β-hydroxysteroid dehydrogenase X, a mitochondrial enzyme which catalyzes oxidation of neuroactive steroids and degradation of isoleucine²⁶. Mutations in HSD17B10 that abolish enzyme activity lead to HSD10 deficiency, an infantile neurodegenerative disorder in which tiglylcarnitine level is elevated²⁷. In contrast, HSD17B10 is overexpressed in brains of individuals with Alzheimer’s disease²⁸, in which tiglylcarnitine level is decreased²⁹.

A previous study reported an association between the GNPTAB intronic variant rs7964859 and aspartate¹⁵. We identified an independent aspartate association signal with the GNPTAB frameshift variant p.Cys528ValfsTer19 (rs1209353188) (LD r²=0.01; β=0.91, P_condition=5.2×10⁻¹⁵ conditioning on rs7964859). p.Cys528ValfsTer19 is rare in METSIM (MAF=0.58%), absent in gnomAD NFE, and is the likely causal variant for the METSIM aspartate association (VPIP=0.996). GNPTAB encodes the alpha- and beta-subunits of N-acetylglucosamine-1-phosphotransferase which catalyzes the N-linked glycosylation of asparagine residues with mannose-6-phosphate³⁰.

Nominating putative causal genes

To nominate putative causal genes for metabolite association signals, we applied two approaches. First, we nominated 66 putative causal genes for the 208 association signals for which fine-mapping analysis identified PTV or missense variants at VPIP≥0.8 (see Fine mapping). Second, we implemented a knowledge-based approach to integrate biological information about the metabolite and the 20 protein-coding genes nearest the corresponding index variant for the 1,666 of the 2,030 association signals with named metabolites (Methods). The knowledge-based approach nominated 215 single genes for 1,033 association signals with 480 metabolites (Supplementary Fig. 9) and 19 sets of 2 to 7 genes with similar biochemical activity (62 additional genes) for 324 association signals with 208 metabolites (Supplementary Tables 7, 11).

We compared gene nominations for the 138 metabolite association signals for which both approaches nominated causal genes. Compared to the fine-mapping analysis, the knowledge-based approach nominated the same gene for 119 signals, multiple paralogs including the same gene for 18, and a different gene for 1, for an overall consistency >99% (Supplementary Table 12).

The 277=215+62 genes identified by the knowledge-based approach and the 66 identified by fine mapping together comprised 290 genes. 204 (70%) of the 290 genes are the closest genes to the index variants, including 188 (68%) of the 277 genes identified by the knowledge-based approach. 58 of the 290 genes have not previously been implicated in metabolite GWAS (Fig. 2c). Of the 58 novel genes, 51 were identified by the knowledge-based approach (Box 1), 21 by fine-mapping analysis (Supplementary Table 9), and 14 by both. Of the 58 novel genes, 40 represent novel loci and 18 are within loci in which metabolite associations have previously been identified but the genes we nominated have not previously been implicated.

Novel genes nominated based on association signals with amino acid levels provide insight into how the encoded enzymes or transporters contribute to modifications of amino acid derivatives. As a first example, we identified a novel association at the HDAC6 missense variant p.Arg832His (rs61735967) with N6-acetyllysine (MAF=2.9%, β=0.71, P=3.6×10⁻⁸⁰) and suggested p.Arg832His is the likely causal variant (VPIP=0.998). Both the fine-mapping and knowledge-based approaches nominated HDAC6 as the putative causal gene. HDAC6 encodes a lysine deacetylase that removes the acetyl group from acetyllysine in histones. Increased HDAC6 expression has been found in brains of individuals with Alzheimer’s disease³¹. Elevated levels of N6-acetyllysine were recently found in an Alzheimer’s disease mouse model³².

As a second example, we identified a novel association between the QPCT intronic variant rs77684493 and pyroglutamylglutamine (MAF=6.2%, β=-0.55, P=9.0×10⁻³¹). rs77684493 is in near-perfect LD (r²=0.996) with the putatively-deleterious QPCT missense variant p.Arg54Trp (rs2255991), which was also associated with pyroglutamylglutamine (β=-0.54, P=1.6×10⁻³⁰) and has >7-fold greater frequency in METSIM than in NFE (MAF=6.3% vs. 0.89%). Our knowledge-based approach nominated QPCT as the putative causal gene for this association. QPCT encodes the enzyme glutaminyl-peptide cyclotransferase, which performs cyclization of the N-terminal glutamine residues and results in the pyroglutamine residue³³. QPCT has been implicated in a schizophrenia GWAS³⁴ and suggested as a druggable target for Huntington’s disease³⁵.

GWAS of metabolites recently identified on the Metabolon platform helped nominate novel putative causal genes with high biochemical relevance in known metabolite-associated regions. For example, a previous study in a Japanese sample identified associations for blood creatinine and uracil levels at the LRIG1 missense variant p.Thr792Met (rs202007714)³⁶, which is monomorphic in METSIM and gnomAD Finns. We identified an association at the nearby (13 kb) SLC25A26 missense variant p.Thr208Met (rs13874) with 2,3-dihydroxy-5-methylthio-4-pentenoate (DMTPA) (MAF=48.3%, β=0.17, P=2.3×10⁻²¹). We suggest SLC25A26 as the causal gene for the DMTPA association. DMTPA, an S-adenosylmethionine, was recently identified on the Metabolon DiscoveryHD4 platform. SLC25A26 is the only known mitochondrial S-adenosylmethionine transporter.

Seven of the 58 novel genes were identified only by fine mapping. Among them, we identified a novel association for glycocholenate sulfate at the rare ADCK5 missense variant p.Ala508Thr (rs552968665; β=1.31, P=3.6×10⁻¹²), which is >79-fold more frequent in METSIM than in NFE (MAF=0.25% vs. 0.0031%). Fine-mapping analysis suggested p.Ala508Thr is the likely causal variant (VPIP=0.89), implicating a causal role for ADCK5. ADCK5 encodes the aarF domain containing kinase 5. These results suggest ADCK5 plays a role in human bile acid metabolism.

Colocalization of metabolites with human diseases

Integrating metabolite and disease genetic associations can improve fine-mapping resolution³⁷ and clarify the potentially causal variants and disease genes. We performed Bayesian colocalization analysis^37,38 based on the probabilistic fine-mapping results of METSIM metabolites and of 980 disease and disease-related dichotomous traits (henceforth, disease traits) in 176,899 Finns in FinnGen release 4 (Methods; Supplementary Table 13). We calculated the regional colocalization probability (RCP) of a putative causal variant shared between a METSIM metabolite and a FinnGen disease trait (Methods). We identified 946 colocalizations involving 248 metabolites and 105 interrelated disease traits (RCP≥0.5; Supplementary Table 14).

Integrating metabolite associations substantially increased the fine-mapping confidence of FinnGen disease trait association signals (SPIP median=0.91 vs. 0.73; paired t-test P=2.5×10⁻⁷⁰) and the probability assigned to the most likely disease variant (maximum VPIP median=0.54 vs. 0.11; paired t-test P=9.1×10⁻¹²⁸; Supplementary Fig. 10). For example, the putatively-deleterious SERPINA1 missense variant p.Glu366Lys (rs28929474) was associated with N-acetylglucosaminylasparagine level in METSIM (MAF=2.3%, β=0.56, P=4.9×10⁻¹⁶) and risk of cholestasis of pregnancy in FinnGen (odds ratio (OR)=6.23, P=8.1×10⁻¹⁷). Our fine-mapping analysis suggested a causal role of SERPINA1 p.Glu366Lys for N-acetylglucosaminylasparagine (VPIP=0.81). We detected colocalization in this region between signals for N-acetylglucosaminylasparagine and cholestasis of pregnancy (RCP=0.99). Colocalizing these association signals increased the SPIP for cholestasis of pregnancy from 0.69 to 0.99, and the VPIP of SERPINA1 p.Glu366Lys from 0.37 to 0.80. SERPINA1 encodes a serine protease inhibitor produced mainly in the liver. SERPINA1 mutations have been associated with familial intrahepatic cholestasis³⁹, and SERPINA1 p.Glu366Lys with liver diseases⁴⁰ and circulating liver enzymes⁴¹.

Campesterol and gallstones: potential causal link

Gallstones affect 10-20% of adults worldwide⁴². Aberrant cholesterol homeostasis, particularly the physical-chemical imbalance of cholesterol solubility in bile, induces gallstones⁴³. Blood campesterol levels have been associated with gallstones⁴⁴, but it is uncertain whether the relationship is causal. We identified associations at the ABCG8 intronic variant rs6544713 with lower campesterol level in METSIM (MAF=20.1%, β=-0.33, P=2.7×10⁻³⁷) and higher gallstone risk in FinnGen (OR=1.32, P=8.0×10⁻⁶⁵). In 4,689 METSIM participants with observed campesterol levels, 199 with gallstones, plasma campesterol level was inversely associated with gallstone risk (β=-0.52, P=3.7×10⁻⁵).

Colocalization analysis suggested campesterol and gallstones share the same causal variant in this region (RCP=0.65; Fig. 4a) and nominated rs6544713 (SCP=0.45) as the most likely causal variant. rs6544713 resides in active regulatory units in intestinal tissue⁴⁵ and its campesterol-decreasing allele is associated with higher ABCG8 expression in colon tissue⁴⁶ (β=0.37, P=7.2×10⁻¹⁶), but not in the liver⁴⁶. ABCG8 and its nearby paralog ABCG5 have previously been suggested as candidate genes for gallstones⁴⁷. Mutations in ABCG5 and ABCG8 cause sitosterolemia, characterized by elevated campesterol⁴⁸.

• Download figure
• Open in new tab

Figure 4:

Colocalization and causal relationship between campesterol and gallstones. a) Stacked regional association plots for campesterol and gallstones (cholelithiasis, K11_CHOLELITH in FinnGen release 4) in the ABCG5/ABCG8 region. The index variants identified in stepwise conditional analysis (campesterol) and approximate conditional analysis (gallstones) are labeled and variants colored by their linkage disequilibrium (LD) to the index variant with which they are in strongest LD in METSIM. The campesterol signal (index variant rs6544713) is colocalized with the gallstone signal (rs4299376, pairwise LD r²=0.993, RCP=0.65) shown in the gray box. In contrast, no colocalization was detected between the signals indexed by rs4614977 and rs11887534. No coding variants within 1 Mb have LD r²>0.2 with rs6544713 in METSIM. b) Comparison of effect sizes for the 15 instrumental variables genome-wide without significant heterogeneity (P>0.05) used in Mendelian randomization analysis between campesterol and gallstones. rs6544713 is in blue. The slope of the blue dashed line depicts the estimated causal effect size of campesterol on gallstones. The Egger regression intercept is deemed not significant (P=0.15). c) Negative relationship between instrumental variable and risk of gallstones. OR: odds ratio.

Mendelian randomization analysis using 15 independent variants for campesterol suggested a causal effect of lower plasma campesterol level on higher gallstone risk (Methods; β=-0.70, P=7.2×10⁻⁸; Fig. 4b-c). Interestingly, we also identified a causal effect of gallstones on lower plasma campesterol level (β=-0.49, P=1.3×10⁻³⁹) in Mendelian randomization analysis using 151 independent variants. ABCG5 and ABCG8 together encode a heterodimeric ATP-binding cassette transporter that facilitates secretion of cholesterol and non-cholesterol sterols in the intestine and bile. High plasma campesterol levels might compete with cholesterol for ABCG5/ABCG8 transporters during biliary cholesterol secretion, resulting in decreased biliary cholesterol levels and reduced risk of gallstones^49,50 (Supplementary Fig. 11).

DBH influence on vanillylmandelate and hypertension: distinct pathways

Stepwise conditional analysis identified the putatively-deleterious DBH missense variant p.Arg79Trp (rs77273740) as associated with lower vanillylmandelate (β=-0.38, P=1.8×10⁻¹⁵); p.Arg79Trp is >10-fold more frequent in METSIM than in NFE (MAF=4.5% vs. 0.39%). Both fine-mapping and the knowledge-based approach suggested a causal role for DBH. The knowledge-based approach suggested DBH could exhibit an effect on vanillylmandelate in two ways (Supplementary Fig. 11): by converting dopamine to norepinephrine, a vanillylmandelate precursor⁵¹, or by transforming homovanillate acid to vanillylmandelate through hydroxylation⁵².

In FinnGen, the DBH p.Arg79Trp vanillylmandelate-decreasing allele was significantly associated with lower hypertension risk (OR=0.84, P=5.2×10⁻¹³), consistent with previous associations with systolic and diastolic blood pressures^53,54. No other variants were associated with FinnGen hypertension or METSIM vanillylmandelate in this region (P>10⁻⁶; Supplementary Fig. 12). In 5,173 METSIM participants with observed vanillylmandelate levels, 1,233 with hypertension, plasma vanillylmandelate level was associated with hypertension (β=0.52, P=3.2×10⁻⁹).

Colocalization analysis suggested that hypertension colocalized with vanillylmandelate (RCP=0.996). However, Mendelian randomization found no evidence for a causal effect of vanillylmandelate on hypertension (P=0.15; 10 independent variants; Supplementary Fig. 12) or hypertension on vanillylmandelate (P=0.17; 157 independent variants), suggesting this signal conferred effects on hypertension risk and vanillylmandelate through different pathways, consistent with the two possible DBH effects identified by our knowledge-based approach (Supplementary Fig. 11).

METabolic Syndrome In Men (METSIM) study

METSIM is a study of 10,197 Finnish men from Kuopio in the late-settlement region of northeast Finland designed to investigate factors associated with type 2 diabetes and cardiovascular diseases³ (Supplementary Table 3). Participants were aged 45 to 74 (median=58) years during baseline visits from 2005-2010. Participants provided demographic, diet, exercise, disease, and medication history information, and underwent laboratory measurements, including oral glucose tolerance test, after ≥10-hour overnight fast. Morbidity, mortality, and drug treatment information was periodically updated for participants who consented to use of their hospital admission, drug reimbursement, and prescription records in Finnish national registries. Due to funding constraints, we randomly selected 6,490 of the 8,777 METSIM participants who at baseline were neither diagnosed with diabetes nor taking diabetes medications that might broadly impact metabolomics levels for the Metabolon metabolomics assay. After exclusion of participants who subsequently developed diabetes (n=264), lacked array genotypes (n=65) or body mass index (BMI) measurement (n=1), had sex mismatch (n=3), and/or were non-Finnish (n=21), our analysis set comprised 6,136 participants (Supplementary Table 3; Supplementary Fig. 1). This study was approved by the Ethics Committee at the University of Eastern Finland and the Institutional Review Board at the University of Michigan. All participants provided written informed consent.

Metabolomics profiling and data processing

Non-targeted metabolomics profiling was performed at Metabolon, Inc. (Durham, North Carolina, USA)⁶¹ on EDTA-plasma samples obtained after ≥10-hour overnight fast during METSIM baseline visits. Briefly, methanol extraction of biochemicals followed by a non-targeted relative quantitative liquid chromatography–tandem mass spectrometry (LC-MS/MS) Metabolon DiscoveryHD4 platform was applied to assay named (n=1,240) and unnamed (n=304) metabolites (Supplementary Tables 1-2). Samples were randomized across batches.

Batches contained ∼144 METSIM samples and 20 well-characterized human-EDTA plasma samples for quality control (QC). All 6,490 samples were processed together for peak quantification and data scaling. We quantified raw mass spectrometry peaks for each metabolite using the area under the curve. We evaluated overall process variability by the median relative standard deviation for endogenous metabolites present in all 20 technical replicates in each batch. We adjusted for variation caused by day-to-day instrument tuning differences and columns used for biochemical extraction by scaling the raw peak quantification to the median for each metabolite by Metabolon batch.

Array genotyping and exome sequencing

All METSIM participants were array genotyped on the Human OmniExpress-12v1_C BeadChip (OmniExpress) and Infinium HumanExome-12 v1.0 BeadChip (exome array) platforms⁶². We excluded individuals for sex or relationship mismatch, apparent sample duplication, or ancestry outliers based on genetic principal component analysis (PCA). We removed variants with genotype call rate <95% (OmniExpress) or <98% (exome array), or Hardy-Weinberg equilibrium (HWE) P<10⁻⁶ (either array)⁶².

We captured exomes for all METSIM participants by SeqCap EZ HGSC VCRome kit (Roche) and sequenced them by HiSeq2000 (Illumina) (average depth 45x)². For exome sequences, we excluded samples with estimated contamination >3% or sample swaps compared to the array genotype data⁶² and required single-nucleotide variant (SNV) array genotype concordance >90% if array data were available. We filtered variants with genotype call rate <98%, HWE P<10⁻⁶, or overall low allele balance (alternate allele count/sum of total allele count <30%)². The resulting array-genotype dataset consisted of n=10,066 METSIM participants with 679,866 SNVs. The exome-sequence dataset consisted of n=9,957 participants with 583,947 SNVs and 40,270 small insertions/deletions (indels).

Genome sequencing

We whole genome sequenced METSIM participants in two waves. In wave 1, we genome sequenced 3,074 METSIM participants (average depth 23x)⁶³. Genomic DNA was fragmented on a Covaris LE220 instrument and size-selected to ensure an average insert size of 350-375 base pairs (bp). Libraries were constructed with the Illumina TruSeq or KAPA Hyper PCR-free library prep kit. qPCR was used to determine concentration of each library. Libraries were subsequently pooled and sequenced with 2×150 bp paired-end reads using HiSeq X (Illumina). We filtered read alignments with mismatch rate ≥5%, inter-chromosomal rate ≥5%, discordance rate of paired reads ≥5%, or haploid coverage <19.5x. We generated QC statistics in Picard v2.4.1 (http://broadinstitute.github.io/picard/), Samtools v1.3.1⁶⁴, and VerifyBamID v1.1.3⁶⁵. We called SNVs and small indels and performed base quality score recalibration in GATK v3.5 (https://gatk.broadinstitute.org/). We excluded variants with missingness >2%, HWE P<10⁻⁶ in unrelated individuals, or allele imbalance <30%. The resulting genome sequence consisted of n=3,074 participants genotyped for 23,849,428 SNVs and 2,914,167 indels. We used wave 1 as part of our imputation reference panel (see METSIM integrative panel and genotype imputation).

In wave 2, we sequenced 2,875 additional METSIM participants using the same methods used for wave 1. We generated a combined wave 1+2 call set of n=5,949 using the same methods, resulting in calls for 55,648,111 SNVs and 12,850,837 indels. Wave 2 data became available only after the main analysis for this paper was complete; we used wave 1+2 combined data to determine linkage disequilibrium (LD) proxies for previously-identified metabolite associated variants that were missing in wave 1 but present in wave 1+2 combined data (see Identification of novel associations).

METSIM integrative panel and genotype imputation

Using the 3,074 METSIM participants with wave 1 genome sequence data, we generated an integrated list of genetic variant sites by merging site lists from the genome and exome sequence data, and the OmniExpress and exome array data. Of the 3,074 participants, 3,055 had OmniExpress and exome array data, and 3,037 had exome sequence data. We calculated genotype likelihoods for each individual at each site as the product of genotype likelihoods assuming independent data across platforms⁶⁶. For OmniExpress and exome array genotypes, we converted genotype calls to genotype likelihoods assuming a genotype error rate of 10⁻⁶. We then phased genotypes using integrated genotype likelihoods in Beagle v4.1 with 50,000 markers per chunk and 3,000 overlapping genetic markers between consecutive chunks⁶⁷. We subsequently excluded 1 individual who self-identified as non-Finnish, 2 individuals identified as population outliers in genetic PCA, and 149 close relatives (estimated kinship ≥0.125 in KING v2.2.1⁶⁸). The resulting integrative panel comprised 2,922 individuals genotyped for 23,294,337 SNVs and 2,851,848 indels (Supplementary Table 4). 2,670 (91.4%) of the 2,922 individuals had Metabolon metabolomics data.

We imputed genotypes for the 6,490 study participants on the framework of their OmniExpress genotypes using the METSIM integrative panel with Minimac v4⁶⁹. We excluded imputed variants with imputation r²<0.3, leaving 19,182,997 SNVs and 2,404,717 indels for downstream analysis (Supplementary Table 5).

Variant functional annotation

We annotated all variants using the Ensembl Variant Effect Predictor (VEP) version 99⁷⁰. We used the “-pick_order” option to annotate each variant using a single transcript, with transcripts prioritized in the following order: transcript support level (i.e. well-supported and poorly-supported transcript models based on the type and quality of the alignments used to annotate the transcript), transcript biotype (protein coding preferred), APPRIS isoform annotation (i.e. annotation based on a range of computational methods to identify the most functionally important transcripts from cross-species conservation), deleteriousness of annotation as estimated by Ensembl, transcript CCDS status (i.e. amount and type of evidence that supports the existence of a variant), canonical status of transcript (https://m.ensembl.org/Help/Glossary), and transcript length⁷¹. We used the dbNSFP (version 4.0)⁷² plugin to generate additional predictions of variant deleteriousness from five in silico algorithms: Polyphen2 HDIV⁷³, Polyphen2 HVAR⁷³, SIFT4G⁷⁴, MutationTaster⁷⁵, and the Likelihood Ratio Test (LRT)⁷⁶.

Trait transformatio

For each metabolite, we inverse normalized the Metabolon-reported metabolite level, regressed on covariates (age at sampling, Metabolon batch, and lipid-lowering medication use status for lipid traits only), and inverse normalized the residuals. In the single-variant association analyses with BMI adjustment, we also included BMI among the covariates.

Single-variant association analysis

We carried out single-variant association tests using a linear mixed model in EPACTS (v3.2.6) (https://github.com/statgen/EPACTS) on the normalized residual metabolite values. We limited analysis to the 1,391 metabolites successfully measured on ≥500 METSIM participants and to the genetic variants with minor allele count (MAC)≥5. This resulted in 10,914,098 to 16,172,108 variants (median=16,042,879) tested across the 1,391 metabolites, since the number of variants with MAC≥5 varied with the set of individuals successfully measured for each metabolite.

To choose a study-wise significance threshold for the 1,391 parallel metabolite GWAS, we carried out PCA across the metabolites to determine the number of principal components required to explain metabolite variation. To account for missing data (Supplementary Fig. 13), we first imputed missing metabolite values using the K-nearest neighbors approach⁷⁷ with K=5. PCA of the imputed data showed that 692 principal components explained 95% of phenotypic variation for the 1,391 metabolites. We therefore used a study-wise significance threshold of P<5.0×10⁻⁸/692=7.2×10⁻¹¹ for our single-variant analyses. A metabolite quantified in n=500 participants provided >80% power at P<7.2×10⁻¹¹ to detect variants that explained phenotypic variance≥11% and had MAC≥5.

PheWeb browser

We built a PheWeb browser⁵⁹ of the 1,391 metabolite GWAS to support interactive visualization, exploration, and download of these results. This PheWeb (https://pheweb.org/metsim-metab/) includes Manhattan and quantile-quantile plots, summary statistics, and links to biochemical characteristics and functions in the Human Metabolome Database (HMDB)⁶⁰ for all 1,391 metabolites.

Stepwise conditional analysis

We carried out stepwise conditional analysis in EPACTS (v3.2.6) (https://github.com/statgen/EPACTS) to identify near-independent association signals. For each metabolite-chromosome pair with at least one single-trait genome-wide significant association (P<5.0×10⁻⁸), we first conditioned on the most significant associated variant and continued conditioning on the most significant remaining variant until no variant attained P<5.0×10⁻⁸.

Fine mapping and credible sets

For each of the 2,030 nearly-independent association signals, we built genomic regions of 1 Mb on either side of the index variant, less near chromosome ends. We merged overlapping regions for the same metabolite, resulting in 1,501 genomic regions of 1.2 to 3.1 Mb. To identify potential causal variants within each region, we performed fine-mapping analysis using the Deterministic Approximation of Posteriors (DAP-g) method²⁰ (https://github.com/xqwen/dap), assigning equal priors to all candidate variants. DAP-g uses individual-level metabolite, genotype, and covariate data to produce fine-mapping results. Since DAP-g does not allow for related participants, we corrected for relatedness approximately by including the first ten genetic principal components as covariates; repeating the DAP-g analysis with 0, 20, or 100 principal components yielded similar results.

DAP-g allows for multiple independent association signals within each region. For each identified signal, DAP-g computes (1) a signal posterior inclusion probability (SPIP) that there is at least one causal variant in the signal; and (2) a posterior inclusion probability for each variant (VPIP) that the variant is causal for the signal. For each of the 1,952 signals identified in stepwise conditional tests that had SPIP≥0.95 in DAP-g, we constructed a 95% credible set of potential causal variants by ranking the variants in descending VPIP and including variants until their summed VPIP was ≥0.95.

Identification of novel associations

To assess which of our metabolite associations were novel, we compiled a list of 381 published metabolite GWAS papers (Supplementary Table 10): 354 from the NHGRI-EBI GWAS catalog⁷⁸ (https://www.ebi.ac.uk/gwas/; release date December 1, 2020); and 27 others from the list curated by Kastenmüller et al.¹¹ (accessed April 1, 2021). From these papers, we identified 8,502 variants with metabolite associations at P<5.0×10⁻⁸ or at the significance threshold used in the paper, whichever was more stringent (Supplementary Table 10). Among these 8,502 published variants, 7,807 were present in the METSIM imputed genotype data. For 194 of the published variants not present in the METSIM imputed genotype data, we identified proxies (LD r²≥0.8 and ≤500 kb) using the wave 1+2 genome sequence dataset of 5,949 METSIM participants. The 7,807 variants present in the METSIM imputed genotype data and the 194 LD proxies for missing variants together comprised 8,000 unique variants. To avoid problems with multicollinearity, we pruned these 8,000 variants at METSIM LD r²>0.99 and ≤1 Mb, yielding 6,501 LD-pruned variants. Then, for each of the 2,030 association signals, we repeated the conditional association analysis including the subset of these LD-pruned variants within ≤1 Mb of the corresponding index variant as covariates. We considered as novel signals those index variants with conditional P<7.2×10⁻¹¹ and location >500 kb from any of the 8,502-7,807-194=501 published variants, which were neither present nor with proxies in the METSIM imputed genotype data. Among the 501 variants, 380 were monomorphic in gnomAD v3.1 Finns (n=5,316).

Knowledge-based approach to gene nomination

To nominate putative causal genes for the 1,666 of 2,030 signals associated with named metabolites, we employed a two-stage knowledge-based approach¹⁵. In stage 1, for each variant, we identified the 20 closest protein-coding genes using the minimum distance from the index variant to the refSeq genes’ transcription start or end sites. We employed an algorithm to look for lexical overlaps between the associated metabolite and each of the 20 genes.

Specifically, we searched automatically for matching strings using customized scripts between: (1) the HMDB⁶⁰ metabolite name and synonyms and Entrez gene names⁷⁹; (2) the metabolite and Entrez gene names listed in HMDB as interacting with the metabolite; (3) the metabolite name and Uniprot protein names⁷⁵ and their synonyms; (4) the metabolite and its parent classes as defined in HMDB and the Uniprot protein names and their synonyms; (5) the metabolite name and rare disease names linked to each gene in OMIM (Online Mendelian Inheritance in Man, https://omim.org/, accessed January 1, 2021) after removing the non-specific substrings uria, emia, deficiency, disease, transient, neonatal, hyper, hypo, defect, syndrome, familial, autosomal, dominant, recessive, benign, infantile, hereditary, congenital, early-onset, idiopathic; (6) the metabolite and its parent classes and Gene Ontology (GO) biological process names⁸⁰ associated with each gene after removing the non-specific substrings metabolic process, metabolism, catabolic process, response to, positive regulation of, negative regulation of, regulation of (we only considered gene sets of <500 genes); and (7) Kyoto Encyclopedia of Genes and Genomes (KEGG)⁸¹ maps (https://www.kegg.jp/) containing the metabolite (as defined in HMDB) and KEGG maps containing each gene (as defined in KEGG) omitting the large “metabolic process map”. For each of these pairs of terms, we calculated a Pair Distance score ranging from 0 to 1 using the Ruby gem “fuzzy_match” (https://github.com/seamusabshere/fuzzy_match), and considered a score >0.5 as a match.

In stage 2, we manually reviewed the evidence collected at stage 1. We selected the biologically most plausible causal gene if we identified experimental evidence linking the gene to the metabolite. >1 putative causal genes could be nominated if >1 gene was suggested in stage 1 and/or 2; this happened most often when a locus contains multiple paralogs with similar biochemical activity. If no clear experimental evidence existed for any of the 20 genes, no causal gene was selected.

Colocalization of FinnGen disease traits and METSIM metabolites

To identify shared causal variants between METSIM metabolites and FinnGen disease traits, we carried out Bayesian pairwise colocalization analysis using fastENLOC^37,38 (https://github.com/xqwen/fastenloc). We downloaded FinnGen release 4 (https://www.finngen.fi/en/access_results) FINEMAP⁸²-based fine-mapping results for 980 disease traits with at least one association at P<5.0×10⁻⁸. fastENLOC used these FinnGen fine-mapping results and our DAP-g-based fine-mapping results for METSIM metabolites to carry out colocalization analysis assuming a single causal variant. For each FinnGen disease trait, we estimated its degree of enrichment for genome-wide associations in metabolite GWAS using TORUS⁸³ (https://github.com/xqwen/torus) and used this enrichment estimate as the prior for Bayesian analysis in fastENLOC. fastENLOC computes two probabilities. The regional colocalization posterior probability (RCP) is the probability of the same causal variant within a region for both the metabolite and the FinnGen disease trait. The variant colocalization posterior probability (SCP) is the probability a specific variant is causal for both traits. We limited colocalization analysis to the 1,952 metabolite stepwise association signals with SPIP≥0.95 for 792 metabolites and present colocalizations for metabolite-FinnGen disease trait pairs with RCP≥0.5 (Supplementary Table 14).

Associations of campesterol with gallstones and vanillylmandelate with hypertension in METSIM

Among the 4,698 METSIM participants with measured plasma campesterol level at baseline, we identified 199 with gallstones in METSIM (December 2020). To test for association between plasma campesterol level and presence of gallstones, we used logistic regression with covariates baseline study age, Metabolon batch, and lipid medication use. Among the 5,173 METSIM participants with measured plasma vanillylmandelate level at baseline, we identified 1,233 individuals with hypertension, and used logistic regression with covariates baseline study age, Metabolon batch, and hypertension medication use to test for association between plasma vanillylmandelate level and hypertension status.

Causal effects between metabolites and FinnGen disease traits

To infer the potential causal effects of plasma campesterol on FinnGen gallstones (phenocode: K11_CHOLELITH) and plasma vanillylmandelate on hypertension (phenocode: I9_HYPTENS), we applied four two-sample Mendelian randomization methods: inverse variance weighted⁸⁴, weighted median⁸⁵, MR-PRESSO⁸⁶, and MR-Egger⁸⁷. These methods make different assumptions and use different strategies to account for horizontal pleiotropy, which can result in false positive inference of causality. For each metabolite, we identified nearly-independent genetic instrument variables (LD r²<0.1, distance≥500 kb) with unconditional single-variant association P<10⁻⁶. We also ran Mendelian randomization analyses to infer the causal effect of FinnGen gallstones (phenocode: K11_CHOLELITH) on campesterol and FinnGen hypertension (phenocode: I9_HYPTENS) on vanillylmandelate in a similar way. We considered findings significant if they had the same effect direction and P<0.05 for all four Mendelian randomization methods. We present MR-PRESSO effect estimate and p-values in the main text.