Genetic Studies of Metabolomics Change After a Liquid Meal Illuminate Novel Pathways for Glucose and Lipid Metabolism

Metabolites, the products of cellular activity and the environment, can be considered as intermediates between genes and clinical phenotypes (1). Recent developments in high-throughput metabolomic profiling based on mass spectrometry (1–4) and nuclear magnetic resonance (NMR) spectroscopy (2,5–9) platforms have opened new avenues to explore gene-metabolite associations by genome-wide association studies (GWAS). In 2014, an atlas describing 145 genetic loci associated with a broad spectrum of blood measures covering amino acids, carbohydrates, lipids, and peptides was generated (10). In a subsequent 2016 study evaluating circulating fasting metabolites (mainly lipoprotein subclasses), >60 genetic loci associated with at least one measure were observed (11). These and other studies have provided considerable mechanistic insight into physiological pathways of diseases. However, the predominant focus of most studies has been on the genetics of metabolomic measures under fasting conditions.

Because of frequent food intake, humans spend the majority of their waking hours in a nonfasting state. However, insight into the genes that affect plasma metabolite concentrations in response to food intake is relatively limited. It seems more than likely that cumulative and prolonged exposure to multiple plasma metabolites in response to food intake may have pathological consequences, as has been well documented for certain lipids and lipoproteins (12–16). Furthermore, it has long been observed that while there may be accepted relationships between consumption and metabolite abundance, there is also natural variation in this response (17).

In clinical practice, an oral glucose tolerance test (OGTT) is commonly used for the screening of suspected diabetes and is performed by determining glucose levels 2 h after ingestion of a fixed dose of glucose. Previous GWAS have identified genetic loci that are associated with glucose and insulin responses from OGTT (18), which only partly overlapped with GWAS findings on fasting glucose and insulin measures. The measurements from OGTT have expanded the understanding of the genetics and pathophysiology of diabetes. However, human food intake is more than sugar alone, and similar to OGTT, meal responses reflected by metabolite profiles are likely to exhibit a large amount of variability, part of which will be heritable. A recent candidate gene study showed that the genetics of fasting and postprandial metabolite concentrations overlap (19). A recent twin study indicated that up to 48% of the variation in the 2-h glucose response to food can be attributed to genetics, whereas the genetic contribution to variation in the triglyceride (TG) and insulin responses was found to be nil to low (0% and 9%, respectively) (12).

In the current study, we aimed to further explore the genetics and heritability of plasma metabolomic measures as assessed by NMR spectroscopy before and after a liquid mixed meal. For this, we performed GWAS on fasting and postprandial metabolomic measure concentrations. We also performed GWAS on the response of these measures using various methods to assess these responses (the change from fasting state to postprandial state after a liquid mixed meal) in a large (n = 5,705) population-based cohort study, the Netherlands Epidemiology of Obesity (NEO) study.

This study was performed in a population-based prospective cohort, the NEO study (20). All 6,671 participants gave written informed consent, and the Medical Ethical Committee of the Leiden University Medical Center approved the study design. Initiated in 2008, the NEO study was designed to study pathways that lead to obesity-related diseases. Detailed information about the study design and data collection has been described elsewhere (20). Briefly, men and women aged between 45 and 65 years with a self-reported BMI of ≥27 kg/m² living in the greater area of Leiden (in the west of the Netherlands) were eligible to participate in the NEO study. In addition, all inhabitants aged between 45 and 65 years from one municipality (Leiderdorp) were invited, irrespective of their BMI.

Participants were invited for a baseline visit at the NEO study center in the Leiden University Medial Center after an overnight fast. Before their visits, participants completed a questionnaire at home, providing demographic, lifestyle, and clinical data. At the baseline visit, fasting blood samples were drawn. Within the next 5 min after the fasting blood draw, a liquid mixed meal (400 mL of 600 kcal, with 16% of energy derived from protein, 50% of energy from carbohydrates, and 34% of energy from fat) was consumed, and subsequent blood samples were drawn 30 and 150 min after the liquid mixed meal. Individuals were excluded from the analyses (Fig. 1) when 1) taking any lipid-lowering medication, 2) violating overnight fasting, or 3) violating liquid mixed meal challenge protocol (i.e., did not finish the liquid mixed meal completely).

DNA was extracted from 6,671 venous blood samples obtained from the antecubital vein. Genotyping was performed in the Centre National de Génotypage (Evry Cedex, France), using HumanCoreExome-24 BeadChip (Illumina, San Diego, CA). The detailed quality control process has been described previously (21). Genotypes were further imputed to the Haplotype Reference Consortium release 1.1 (22). All genetic variants with an imputation quality <0.4 or a minor allele frequency (MAF) <0.01 were not considered for the analyses in the current study. As such, a total of 5,705 individuals with genotype data for 7,381,632 variants were used in our association analysis.

Metabolomic measurements were performed in both fasting and postprandial (t = 150 min after the liquid mixed meal) plasma samples using the Nightingale high-throughput NMR metabolomics platform (23). Each box of 94 samples (including both fasting and postprandial) was run in sample order, which was based on the sampling date and guaranteed that fasting and postprandial samples from the same participant were run in the same batch. Due to budget constraints and the fact that metabolite levels at 150 min change more significantly than 30 min after the meal (24), there were no metabolomic measurements at 30 min. We note that metabolites are commonly defined as biological molecules <1.5 kDa in size and that many of the molecules (lipids) assayed here are larger than this threshold. Nevertheless, for simplicity, we refer to all of them here as metabolites. The metabolomics platform provides measurements for 148 metabolites (Supplementary Table 1) from 11 substance classes: lipoprotein subclasses (n = 98), lipoprotein particle sizes (n = 3), apolipoproteins (n = 2), fatty acids and saturation (n = 11), cholesterol (n = 9), glycerides and phospholipids (n = 9), amino acids (n = 8), ketone bodies (n = 2), inflammation (n = 1), glycolysis-related metabolites (n = 3), and fluid balance (n = 2). The NMR-based metabolomics platform and the experimental procedure have been described in detail previously (25).

In addition to metabolomics measurements, fasting glucose and insulin as well as postprandial glucose and insulin at 30 min and 150 min after a liquid meal were measured by clinical chemistry approaches (26). Fasting hemoglobin A_1c (HbA_1c) concentrations were expressed as percentages. By using both fasting and postprandial measurements at 30 min and 150 min after the liquid mixed meal, insulin disposition index (IDI), insulinogenic index (IGI), Matsuda insulin sensitivity index (ISI), HOMA for insulin resistance, and HOMA for β-cell function (HOMA-β) indices were derived on the basis of standard formulas.

To remove samples of low quality and with measurement errors, individuals were excluded when metabolite concentrations 1) deviated >10 SDs of the mean from the entire NEO population and 2) had >30% missingness across all 148 metabolites in either fasting or postprandial states (Fig. 1). Metabolite concentrations were inverse rank normal transformed (27) using an edited version of the rntransform() function from the GenABEL package in R (28), which randomly ranks tied values. When analyzing metabolites in the fasting or postprandial states alone, each state was transformed independently of each other. However, when analyzing or deriving the response phenotype, the data from the two states were merged before data transformation. As previously observed that sampling date (a compound variable composed of blood sampling year and month) has an appreciable effect on metabolite concentration in this data set, we included this variable as a covariate before/during linear regression (29).

Metabolite response, or the change in concentration of a metabolite between the fasting and postprandial state, was analyzed in two manners. First and foremost, for each metabolite, we derived a response phenotype defined as the residuals of an orthogonal nonlinear least squares regression (OrNLSr), where the postprandial state was set as the dependent variable and the fasting state as the independent variable in a univariate analysis. The technical details regarding the OrNLSr method and comparison with other methods are discussed in the Supplementary Material. Response could have been defined as a simple estimation of change or Δ between postprandial and fasting states. However, we observed that for 30 of the 148 metabolites analyzed, the data were best explained by a nonlinear curve as opposed to a linear one possibly because of a physiological plateauing effect for some metabolites (Supplementary Fig. 1). Simple Δ estimates of change or residuals derived from linear regressions between the two states would not accurately capture the variation of response or metabolite change.

We also evaluated response via linear mixed models (LMMs) that included an interaction term between dietary state (either fasting or postprandial state) and genotype alongside individual random effects. Under this framework, we simultaneously measured the effect of dietary state, genotype, and the interaction of state and genotype on metabolites. Here, the interaction term provides a measurement of genotype on response or specifically, a differential effect across the two states. Given the computational expense of the LMM, we chose to remove all genetic variants with an imputation quality <0.4 or MAF <0.01 in the double-sized genotype data (n = 11,410) before genome-wide association analysis for a specific metabolite, thereby also reducing our testing burden.

Four unique GWAS (fasting, postprandial, OrNLSr response, and LMM response) were performed on 148 metabolomic measures across 4,734, 4,348, 4,292, and 11,410 individuals, respectively (Fig. 1), and 7,701,709 (7,568,622 for the LMM) genetic variants across the 22 autosomal chromosomes. For fasting, postprandial, and OrNLSr response, linear regression analyses, assuming an additive genetic model, were performed with SNPTEST v2 (30) on the residuals of inverse rank normal transformed metabolomic measures after adjusting for age, sex, the first 10 principal components, and the batch effect variable (a compound variable composed of blood sampling year and month). By considering that patients with diabetes are known to have a dramatically altered metabolism of macronutrients, we conducted sensitivity analyses by adding diabetes status as a covariate in the model and reran association analyses for top signals. The explained variance was estimated as the partial R² from the linear regression model, with the single nucleotide polymorphism (SNP) as independent variable and metabolite response represented by the OrNLS residual (with inverse normal transformation) as dependent variable. For illustration purposes, we also reported the response to a liquid meal (i.e., OrNLSr response GWAS results) expressed as percent change from the fasting state. To obtain this percent change, we reran the association analysis on the traits of interest using the raw OrNLS residuals (i.e., without inverse rank normal transformation) divided by the baseline fasting measurements times 100, which gives the percent change per allele. The P value from this analysis is distorted because of a nonnormal distribution of phenotypes, but the effect size obtained corresponds to the value of a percent change from the fasting state.

The LMM framework essentially tested for an interaction between genotype and dietary state. To reduce the computational burden, the LMM was implemented in two steps. In the first step, fasting and postprandial metabolomic measures were concatenated into a long-format phenotype file, and an edited version of inverse rank normal transformation (mentioned above) was applied to the raw measurements. After transformation, the LMM with individual random effects was applied using the function lmer in the R package lme4 (31) to the inverse rank normal transformed metabolite measures, further adjusting for age, sex, the first 10 principal components, and the batch effect variable sample date. Conditional residuals derived from the first step were subsequently used as dependent variables for linear regression models against genotype and dietary state interaction in the second step performed in ProbABEL, with sandwich SE being estimated (32).

To identify independent loci in each metabolite-by-dietary state GWAS, we first identified linkage disequilibrium (LD)–independent blocks by clumping all variants with a standard GWAS significance level (α) of P < 5 × 10⁻⁸ in PLINK (33). All genetic variants with a P value below the α threshold (–clump-p1) were set as “index” SNPs, and all the other SNPs were clustered into different clumps or LD blocks on the basis of their LD and physical proximity to the index SNPs controlled by clump-r2 and clump-kb separately in the command (4). In the current analyses, the following parameters were adopted: –clump-p1 5.0 × 10⁻⁸, –clump-r2 0.5, –clump-kb 1,000. LD patterns were based on the 1000 Genomes v3 20101123 reference set of Utah Residents with Northern and Western European Ancestry population (34).

In addition, primary and secondary signals were subsequently identified through stepwise conditional analyses using the genome-wide complex trait analysis (GCTA) tool v1.24.4 (35), with a parameter of MAF >0.01. Conditional P < 5 × 10⁻⁸ was considered to be genome-wide significant. These primary and secondary signals identify tagging SNPs for an association and represent statistically independent (at P < 5 × 10⁻⁸) associations.

For each dietary state, we also defined “genomic regions” of association by collapsing all associated genetic variants that were within 500 kb of one another into a single region. Finally, we compared the genomic regions defined for each dietary state with one another to define “study genomic regions.” If multiple genomic regions across dietary states overlapped in their genomic coordinates, they were merged into one region.

To account for the strong intercorrelations of lipoprotein subclasses (Supplementary Fig. 2) in false discovery rate–based multiple test corrections, we applied the variance decomposition method proposed by Li and Ji (36). This resulted in 39, 38, and 44 independent components underlying 148 metabolomic measures in the fasting state, postprandial state, and response as determined by OrNLSr, respectively. Accordingly, metabolome-wide significance was set to 1.28 × 10⁻⁹, 1.32 × 10⁻⁹, and 1.14 × 10⁻⁹ for fasting, postprandial, and response GWAS (OrNLSr and LMM GWAS), respectively. Genome-wide significance was set to 5 × 10⁻⁸, which was used to derive the metabolome-wide values by variance decomposition (36). Our study-wide α-threshold is 3.0 × 10⁻¹⁰ (5 × 10⁻⁸ / [39 + 38 + 44 + 44]). The analysis workflow is shown in Fig. 1. Associations that surpassed the genome-wide α-threshold were used to identify primary and secondary associations as well as LD-independent blocks and genomic regions of association.

The independent top variants were entered in GTEx Portal (37) to identify the tissue-specific expression patterns (expression quantitative trait loci [eQTLs]) of mRNAs associated with these variants. The phenome-wide association studies (Phewas) of independent top variants was performed by GeneAtlas on the basis of 118 nonbinary and 599 binary traits of 408,455 related and unrelated UK Biobank participants (38).

Narrow-sense SNP-based heritability for each metabolomic measure under fasting, postprandial, and OrNLSr-derived response states was estimated by restricted maximum likelihood under the framework of GCTA. Genetic variants (both genotyped and imputed genotype data) with MAF >1% were retained for the analysis (32).

The individual-level data sets that are used for current study are not publicly available because of a privacy issue, but all the summary statistics results from the GWAS are available from the corresponding author upon reasonable request. No applicable resources were generated or analyzed during the current study.

Postprandial metabolite concentrations were correlated with their fasting levels (median absolute Pearson r = 0.29 [interquartile range 0.12, 0.59]). Notably, for 30 metabolomic measures, these associations were nonlinear (Supplementary Fig. 1). Since the effects of fasting state levels on the response was regressed out by OrNLSrs, the metabolite response residuals as determined by OrNLSr showed low correlations to either fasting (median absolute Pearson r = 0.088 [interquartile range 0.043, 0.14]) or postprandial (median absolute Pearson r = 0.11 [interquartile range 0.047, 0.19]) state measures.

In total, 32,212 SNP-fasting metabolite associations were discovered at P < 5 × 10⁻⁸ (Supplementary Table 2). These associations include 2,249 unique genetic variants and 144 of the 148 tested metabolites (Fig. 2). A total of 743 independent, primary associations were identified, which involve 119 unique SNPs and mapped to 39 genomic regions (Supplementary Table 2 and Supplementary Fig. 3). In addition, 16 secondary associations were identified involving six unique SNPs and 15 metabolites (Supplementary Table 3). The 16 secondary signals only surpassed the GWAS α-threshold in the conditional analysis and were not among the initial 32,212 associations. In total, 512 fasting primary associations surpassed our study-wide α-threshold of 3.0 × 10⁻¹⁰ (Supplementary Table 2).

Among the 62 associations previously reported by Kettunen et al. (11), 46 could be tested in the current fasting GWAS; the other 16 associations could not be evaluated because of our quality control: SNP MAF <0.01 or imputation quality <0.4. Among these 46 identified associations, 36 were successfully replicated (P < 0.05) in the NEO study (Supplementary Table 4).

At P < 5 × 10⁻⁸, the postprandial metabolite GWAS identified 30,747 SNP-metabolite associations (Supplementary Table 5). These associations include 1,817 unique genetic variants and 140 of the 148 tested metabolites (Fig. 2). A total of 675 independent, primary associations were identified using joint conditional analysis, which involve 94 unique SNPs and mapped to 28 genomic regions (Supplementary Table 5 and Supplementary Fig. 4). In addition, 14 secondary associations were identified involving eight tagging SNPs and 13 metabolites (Supplementary Table 3). In total, 475 postprandial primary associations surpassed our study-wide α-threshold of 3.0 × 10⁻¹⁰ (Supplementary Table 5).

In the GWAS of residuals derived by OrNLSr, we identified 234 SNPs that contributed to variation in metabolite response to a meal (P < 5 × 10⁻⁸) (Supplementary Table 6). These associations included 74 unique genetic variants and 23 of the 148 tested metabolites. A total of 23 primary associations were identified using joint conditional analysis, which involved 14 unique SNPs and mapped to 12 genomic regions (Supplementary Table 6 and Supplementary Fig. 5). No secondary associations were observed.

Three primary associations surpassed the metabolome-wide α in the OrNLSr GWASs and mapped to two genomic regions (Table 1). The strongest association, which also surpasses our study-wide α, is on chromosome 11, is tagged by rs10830963:G (MAF 0.26, β [SE] −0.23 [0.03], P = 2.15 × 10⁻¹⁹) at the MTNR1B locus (Fig. 3A and Table 1) and is associated with glucose response. rs10830963 explained 1.9% of total variance in glucose response represented by OrNLS residuals. In addition, rs458741:C located on chromosome 5 (ANKRD55 locus) was associated with extremely large VLDL total cholesterol (XXLVLDLC) levels (MAF 0.23, β [SE] 0.17 [0.03], P = 5.76 × 10⁻¹⁰) and XXLVLDL cholesterol ester (XXLVLDLCE) levels (MAF 0.23, β [SE] 0.17 [0.03], P = 9.74 × 10⁻¹⁰). Although not reaching study-wise significance, signals led by rs467022:T and rs173964:G at the ANKRD55 locus also showed suggestive associations with extremely large VLDL total lipids (XXLVLDLL), extremely large VLDL particle concentration (XXLVLDLP), extremely large VLDL triglycerides (XXLVLDLTG), extremely large VLDL free cholesterol (XXLVLDLFC), and extremely large VLDL phospholipids (XXLVLDLPL) (Table 1). The three SNPs (i.e., rs458741, rs467022, rs173964) located in ANKRD55 locus in total explained 0.85–0.98% of total variance of seven XXLVLDL-related metabolite response.

To further assess whether genetic variants associated with fasting or postprandial metabolomic measure concentrations affected response, we performed a genome-wide LMM analysis for each metabolite in an alternative assessment of meal response. We identified a total of 37 genetic variants that contributed to variation in metabolite response to a meal in this mixed model analysis (P < 5 × 10⁻⁸) (Supplementary Table 7). These associations include 15 metabolites and 31 SNPs (15 unique SNPs) mapping to 12 genomic regions (Table 2). The strongest effects were observed for very large HDL particles total lipids (XLHDLL) and citrate, mapping to rs116041093 and rs632200, respectively (Table 2). Both these genetic variants had a low-effect allele frequency (0.03 and 0.01 separately) and did not reach genome-wide significance in the OrNLSr response GWAS. The third strongest effect was the glucose response and mapped to rs10830963:G in the intron of MTNR1B gene (genotype-by-dietary state interaction: β [SE] −0.16 [0.02], P = 4.66 × 10⁻⁹), which is consistent with findings in the OrNLSr response GWAS (Supplementary Fig. 6). The sensitivity analysis by adding diabetes status as a covariate showed very similar effect size estimations (data not shown).

The fasting and postprandial GWASs shared an abundance of overlap. In total, 86% (26,419 of 30,747) of the postprandial associations were also observed among the fasting associations. Moreover, 397 of the 675 independent primary SNP-metabolite postprandial associations were also primary associations in the fasting data set. This overlap includes 131 metabolites and 50 SNPs, the most common of which are SNP rs964184 at band 11q23.3 near the APOA5 locus (48 metabolites), previously associated with hypertriglyceridemia (39) and the APOE exonic missense SNP rs429358 at band 19q13.32 (59 metabolites) previously associated with hyperlipoproteinemia (40) and Alzheimer disease (41,42). In contrast, 278 primary associations were unique to the postprandial state, while 346 were unique to the fasting dietary state.

However, if we look not at overlap in tagging SNPs but at the genomic loci they map to, we find that there are 22 genomic regions involving 139 metabolites that carry independent, primary associations (546 genomic region-metabolite pairs) in both the fasting and the postprandial state (Supplementary Fig. 7 and Supplementary Table 8). In contrast, 15, 5, and 10 genomic regions harbor SNP-metabolite primary associations in just one dietary state: fasting, postprandial, and response, respectively. This includes glucose at 1) the G6PC2 locus (rs560887, band 2q31.1, Study-GR8), 2) the AGMO locus (rs10231021, band 7p21.2, Study-GR20), 3) the GCK locus (rs2971670, band 7p13, Study-GR21), and 4) the MTNR1B locus (rs10830963, band 11q14.3, Study-GR31), which each exhibit an association with glucose in the fasting state (P < 2.20 × 10⁻⁸) but not in the postprandial state (P > 0.0037).

A previous study showed that rs10830963 is an eQTL in human islets conferring increased MTNR1B mRNA expression. The G-allele carriers exhibited a stronger inhibition of insulin secretion after melatonin treatment (43). Given the consistently observed association between rs10830963:G and glucose response from our metabolite response GWAS, glucose as measured by a clinical chemistry laboratory at fasting state, 30 min, and 150 min after a liquid mixed meal were examined in the NEO cohort stratified by rs10830963:G genotype. The postprandial glucose excursions, defined as the change in glucose concentrations from before to after a meal, showed evidence of difference (fasting status and genotype interaction P = 7.6 × 10⁻⁴) (Fig. 3B). rs10830963:G has been linked to several glycemic traits in a previous large-scale GWAS meta-analysis, with the strongest signals on fasting glucose (Fig. 3C) followed by HOMA-β.

The three SNPs (rs458741, rs467022, and rs173964) are in high LD with one another and belong to the same ANKRD55 locus. By a Phewas in the UK Biobank, the ANKRD55 locus is significantly associated with several traits reflecting body composition (e.g., impedance of arms and legs, waist-to-hip ratio, trunk fat percentage, body fat percentage), blood traits (e.g., red blood cell distribution width, platelet distribution width, hemoglobin concentration), and diabetes, with P < 5 × 10⁻⁸ (Fig. 5A). Except for rs458741without eQTL information, rs467022 and rs173964 both are eQTLs in multiple tissues conferring increased mRNA expression levels, with the strongest expression levels in kidney (for rs467022) and esophagus as well as for whole blood (for rs173964) separately (Fig. 4A). A previous study performed by Scott et al. (44) reported a type 2 diabetes (T2D) association signal led by rs173964 at the ANKRD55 locus, with an odds ratio of 1.08; however, we did not replicate this association in the NEO study (participants with and without T2D 254 and 4,215, respectively), which is likely due to the lack of power. To further understand the roles of chylomicrons and XXLVLDL particle response in the risk of T2D, we investigated the associations between responses of XXLVLDL particles (rXXLVLDL) and T2D risks as well as other glycemic traits in the NEO study. Although there was no significant association to T2D, we identified strong links between the responses of chylomicrons and rXXLVLDL particle to IGI, Matsuda ISI, HOMA-β, and HbA_1c (Fig. 4B and Supplementary Table 10), after adjusting for age, sex, and BMI.

On average, the SNP heritability was 31% for fasting metabolomic measure concentrations, which was higher than the average SNP heritability for postprandial metabolomic concentration and response measures (27% and 12%, respectively) (Supplementary Table 9). SNP heritability for fasting-state metabolomic measures was higher than the counterpart postprandial measures (paired one-sided Wilcoxon rank sum test P = 2.2 × 10⁻¹¹). Overall, the heritability of metabolomic measure responses was much lower than either fasting- or postprandial-state measures. However, the response measures for glucose, amino acids (except for histidine), XXLVLDL, and XLHDL showed a heritability estimate of ∼25% (Fig. 5). In fact, heritability point estimates for glucose (0.27, SE 0.07), alanine (0.29, SE 0.08), phenylalanine (0.22, SE 0.07), valine (0.31, SE 0.07), and XXLVLDL traits (mean 0.28) were larger in the response state than they were in fasting and postprandial states.

By performing GWAS on measures before and after a liquid mixed meal and on parameters describing meal response. We observed highly overlapping genetic association signals between fasting and postprandial metabolites. By using baseline-adjusted nonlinear residuals to determine the meal response, rs10830963:G, located in the intron of gene MTNR1B, was observed to be associated with glucose response, which was also found using an LMM-based approach. Importantly, we identified a previously unreported distinct association between ANKRD55 locus led by rs458741, rs467022, and rs173964 and rXXLVLDL particle (all seven XXLVLDL particles). The heritability of the meal response measured by OrNLSr was much lower, on average, than either fasting or postprandial state measures. This suggests that variation in the meal response is to a larger extent more attributable (or susceptible) to nongenetic determinants, such as environmental exposures like the gut microbiome (45).

MTNR1B plays a role in central circadian clock regulation to accommodate diurnal rhythms (46). rs10830963:G is located in the middle of the single intron of the MTNR1B gene and was found to be associated with fasting glucose levels and T2D risk in previous GWAS studies (47,48). Concurrently, this variant was shown to be associated with decreased early-phase insulin secretion (49,50), which is normally considered as the earliest detectable abnormality in individuals who are prone to develop T2D (51). In line with this, G-allele carriers were observed to have a 20% increased risk to develop prediabetes (hazard ratio 1.20 [95% CI 1.15, 1.27]), whereas no additional risk was observed for the progression to T2D from an impaired fasting glucose state (hazard ratio 0.98 [95% CI 0.89, 1.07]) (52). In the current study, glucose response was induced by a liquid mixed meal with 50% of energy from carbohydrates, which more closely resembles regular meal intake. We observed a strong effect of the rs10830963:G risk allele on plasma glucose levels at 30 min after a meal (Fig. 3B). Theoretically, this translates to a lifelong accumulative effect from this genotype on glucose exposure after meals and likely plays an important role in the development of glucose intolerance and insulin resistance. A recent study showed that glycemic response to a meal for the same individuals was on average twofold higher at lunch than at breakfast (45). Taking the fact of meal timing into account, the effect of rs10830963:G risk allele on glucose response identified in the current study would be even more striking, as the meal challenge intervention was performed between 8:00 a.m. and 12:00 p.m. (i.e., before lunch) in the NEO study.

The liquid mixed meal used in the current study provided 34% of energy from fat, which confers an opportunity to investigate the exogenous pathway for lipid metabolism. Dietary cholesterol and fatty acids are absorbed into small intestine, where micelles are formed and subsequently transported to enterocytes. TGs are formed from free fatty acids and glycerol, and cholesterol is esterified. Inside the Golgi body, TGs combine with proteins (e.g., ApoB-48) synthesized in rough endoplasmic reticulum to form chylomicrons (Fig. 5B), which enter the circulation and travel to peripheral sites. ANKRD55 locus led by rs458741, rs467022, and rs173964 showed strong associations with the responses of chylomicrons and XXLVLDL particles in the blood, implying its role in chylomicron synthesis and transportation. Interestingly, ANKRD55 locus is associated with both body composition and diabetes. The association between rs173964 and the risk of T2D has been reported previously by Scott et al. (44); however, the function of rs173964 was not elucidated in their article. By deriving endophenotypes of chylomicrons and XXLVLDL particle responses to a liquid meal, we observed that the responses of chylomicrons and XXLVLDL particles are strongly associated with glycemic indicators for β-cell function (i.e., IGI, HOMA-β), insulin resistance (i.e., Matsuda ISI), and a surrogate for the average blood glucose levels in the previous 3 months (i.e., HbA_1c). In summary, ANKRD55 locus could potentially affect chylomicron synthesis and transportation after fat intake and consequently influence β-cell function and ultimately lead to a higher risk of T2D.

Personalized nutrition aims to exploit the interactions between genetics and modifiable environmental factors to determine food choices that lead to favorable glycemic and lipidemic responses and presumably decreased disease risk. Predictors from the gut microbiome, genetic variants, habitual diet, and meal context have been explored in prediction models of glycemic responses (45,53) and lipidemic responses (i.e., TG response) to different types of food (45). It was recently shown that the genetic component explained nearly one-half of the variation in the glucose response; however, the influence on TG and insulin response (reflected by C-peptide) was minimal (0% and 9%, separately). In the current study, we found that the SNP-based heritability of the glucose response to a meal was 27% (SD 7.2%), whereas only 7% (SD 7.4%) of variation in the serum TG response was due to genetics. Our findings are in line with the previous study (45), and in addition, we extended the heritability estimations of a meal response to a spectrum of 148 metabolomic measures. In general, we found that genetics was not a predominant determinant of meal response for most metabolites (on average 12%), which clearly challenges the application of genetics for personalized nutrition advice.

Several methodological aspects should be considered. The main strength of this study is the liquid mixed meal that was provided to all the NEO participants, which more closely resembles normal meal consumption than a glucose tolerance test to assess glucose metabolism after a meal. We generated novel insight into the genetic basis for fasting and postprandial metabolite concentrations in a general population. Moreover, to assess metabolite responses, we used two different methods to account for the potential bias introduced by baseline adjustment. When genetics contributions to metabolite responses after a meal are assessed, it is important to realize that these responses are affected by their fasting baseline abundance. In the original analysis plan, we adopted plain Δ (i.e., the difference between postprandial and fasting measures) to represent metabolomic responses, as had several other GWAS efforts on response measures (26,54). However, once the theoretical framework of the OrNLS method had been fully developed, we were convinced that we also needed to use it in our present study because OrNLS corrects the metabolite response measure much better for the baseline values than plain Δ and an LMM given that it takes nonlinear response profiles into account. Since LMM is widely considered as the de facto statistical approach to analyze repeated-measures data, even though we believe that it has certain shortcomings in the setting of our present study (e.g., LMM cannot address the nonlinear relationships between fasting and postprandial measurements that we observed in our data), we decided to still include the results from the LMM as a reference. Nonetheless, the sample size of the GWAS was relatively small to identify low-frequency and rare genetic variants. In addition, it is not known whether the association between rs10830963:G genotype and glucose response is generalizable to different periods of the day. Finally, because of the uniqueness of the current study design administering all participants a liquid mixed meal, we could not replicate our findings in other independent cohorts, which calls for further efforts to investigate our findings.

The majority of metabolomic response variation does not have a large genetic component, which raises questions about the applicability of genetic variants to predict postprandial responses to food. We observed that the rs10830963:G variant in the MTNR1B gene is a genetic determinant of both fasting glucose levels and fasting-independent postprandial glucose response to a liquid mixed meal. ANKRD55 locus led by rs458741, rs467022, and rs173964 showed strong associations with XXLVLDL particle response that were linked to β-cell function and insulin resistance. Further studies are warranted to corroborate the biological pathways of the ANKRD55 locus underlying diabetes.

Acknowledgments. The authors thank all individuals who participated in the NEO study and all participating general practitioners for inviting eligible participants. The authors also thank the NEO study group, Pat van Beelen, Department of Clinical Epidemiology, Leiden University Medical Center, Leiden, the Netherlands, and all research nurses for data collection; Petra Noordijk, Department of Clinical Epidemiology, Leiden University Medical Center, Leiden, the Netherlands, and her team for laboratory management; and Ingeborg de Jonge, Department of Clinical Epidemiology, Leiden University Medical Center, Leiden, the Netherlands, for data management of the NEO study.

Funding. The authors acknowledge support from the Netherlands Cardiovascular Research Initiative, an initiative supported by the Dutch Heart Foundation (CVON2014-02 ENERGISE), and Nutricia Research (Utrecht, the Netherlands) for providing the mixed meal. N.J.T. is a Wellcome Trust Investigator (202802/Z/16/Z), is the principal investigator of the Medical Research Council (MRC) Avon Longitudinal Study of Parents and Children (WT 102215/2/13/2), is supported by the University of Bristol National Institute for Health Research Biomedical Research Centre (BRC-1215-20011), the MRC Integrative Epidemiology Unit (MC_UU_12013/3), and works within the CRUK Integrative Cancer Epidemiology Programme (C18281/A19169). D.A.H. is supported by the Wellcome Trust Investigator Award (202802/Z/16/Z). D.O.M.-K. is supported by the Dutch Science Organization (ZonMW-VENI grant 916.14.023).

Duality of Interest. R.L.-G. is a part-time consultant for Metabolon, Inc. No other potential conflicts of interest relevant to this article were reported.

Author Contributions. R.L.-G. and D.A.H. performed the analysis and wrote and edited the manuscript. J.B.v.K. derived the OrNLS response measurements. R.d.M. and F.R.R. contributed to the study design and manuscript review. D.O.M.-K., N.J.T., and K.W.v.D. contributed to the study design, conceived the idea of the current study, interpreted the results, and reviewed the manuscript. R.L.-G. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.

Prior Presentation. Parts of this study were presented in abstract form at the American Society of Human Genetics Conference, Houston, TX, 15–19 October 2019, and as a poster at the 13th Annual Meeting of the International Conference on Genomics (ICG-13), Shenzhen, China.