A Mendelian Randomization Study of Metabolite Profiles, Fasting Glucose, and Type 2 Diabetes

Type 2 diabetes is a progressive metabolic disease characterized by hyperglycemia, initially as a result of insulin resistance and in later stages also as a result of insulin insufficiency. Type 2 diabetes is also associated with dyslipidemia, including higher circulating concentrations of triglycerides and lower concentrations of HDL cholesterol. In addition, several circulating molecules have previously been shown to be dysregulated in type 2 diabetes, including phospholipids, branched-chain amino acids, keto-acid metabolites, and other metabolites such as acyl-carnitines (1–3). However, the causal paths between these metabolites and glucose/type 2 diabetes in humans remain unclear from observational studies and require randomized controlled trials that are difficult to conduct.

As an alternative, Mendelian randomization (MR) is an instrumental variable method that has gained popularity over the last decade to investigate causal effects of traits using genetic predictors. MR uses the principle that the allocation of genetic variants that affect a specific trait from parents to offspring is random and unrelated to factors other than the trait (4). Furthermore, associations between the genotype and the outcome will not be affected by reverse causation because disease will occur after the meiosis. MR has previously been used to determine whether metabolic markers such as classic blood lipids are causally involved in type 2 diabetes (5–11) and has yielded contradicting results. One reason for this could be that these studies are affected by the heterogeneous nature of the metabolic markers chosen, such as in the example of total HDL cholesterol, which in reality consists of a collection of different-sized HDL particles possibly with different functions. This may dilute the causal effects of single nucleotide polymorphisms (SNPs) when only combined (total) HDL is considered. However, false signals may also be due to pleiotropic effects of the chosen genetic variants leading to possibly invalid instrumental variables. As high-throughput analyses techniques improve, the quantification of circulating molecules is becoming ever more detailed and precise. For instance, instead of LDL cholesterol, HDL cholesterol, and total triglycerides determined by routine clinical biochemistry, lipoprotein particle size distribution and content as well as tens of biochemical components can now be measured using nuclear magnetic resonance (NMR) spectroscopy– and mass spectrometry (MS)-based approaches (12,13). These additional measures offer an opportunity to gain novel insight into the pathogenesis of diseases like type 2 diabetes. With the knowledge of genetic determinants of metabolites gained from genome-wide association studies (GWAS) (14–16), one can use MR for causal inference given the specific conditions encoded in Fig. 1. In the current study, with the aim of unraveling potentially causal metabolic paths that underlie the observed associations, we used genetic predictors from published metabolite GWAS, guided by pathway-based evidence to select instrumental variables, and performed MR between selected metabolic markers and glucose/type 2 diabetes.

The observational associations between metabolites and fasting glucose/type 2 diabetes were tested in the Erasmus Rucphen Family (ERF) study, which is a prospective family-based study with 3,465 individuals in the southwest of the Netherlands. The study protocol for ERF was approved by the medical ethics board of the Erasmus Medical Center, Rotterdam, the Netherlands (17). The baseline demographic data and measurements of the ERF participants were collected between 2002 and 2006. Venous blood samples were collected after at least 8 h fasting. The detailed description of the ERF study and related measurements were reported previously (17). Type 2 diabetes was defined according to a fasting plasma glucose ≥7.0 mmol/L and/or antidiabetes treatment. The analytical sample included 2,564 participants without diabetes and 212 participants with diabetes.

Metabolic markers were measured by five different metabolomics platforms using the methods that have been described in earlier publications (15,16,18,19). In total, 562 metabolic markers including subfractions of lipoproteins, triglycerides, phospholipids, ceramides, amino acids, acyl-carnitines, and small intermediate compounds, which throughout this article will be referred as “metabolites,” were measured either by NMR spectrometry or by MS. The platforms used in this research are the following: 1) liquid chromatography-MS (LC-MS) (116 positively charged lipids comprising 39 triglycerides, 47 phosphatidylcholines, 8 phosphatidylethanolamines, 20 sphingolipids, and 2 ceramides available in up to 2,638 participants) measured in the Netherlands Metabolomics Centre, Leiden, using the method described previously (18); 2) electrospray-ionization MS (ESI-MS) (in total, 148 phospholipids and sphingolipids comprising 16 plasmologens, 72 phosphatidylcholines, 27 phosphatidylethanolamines, 24 sphingolipids, and 9 ceramides available in up to 878 participants) measured in the Institute for Clinical Chemistry and Laboratory Medicine, University Hospital Regensburg, Regensburg, Germany, using the method described previously (14); 3) small molecular compounds window-based NMR spectroscopy (NMR-COMP) (41 molecules comprising 29 low–molecular weight molecules and 12 amino acids available in up to 2,639 participants) measured in the Center for Proteomics and Metabolomics, Leiden University Medical Center (19,20); 4) lipoprotein window-based NMR spectroscopy (NMR-LIPO) (104 lipoprotein particle subfractions comprising 28 VLDL components, 30 HDL components, 35 LDL components, 5 IDL components, and 6 plasma totals available in up to 2,609 participants) measured in the Center for Proteomics and Metabolomics, Leiden University Medical Center (lipoprotein subfraction concentrations were determined by the Bruker algorithm [Bruker BioSpin GmbH, Germany] [16]); and 5) AbsoluteIDQ p150 Kit of Biocrates Life Sciences AG (153 molecules comprising 14 amino acids, 91 phospholipids, 14 sphingolipids, 33 acyl-carnitines, and hexose available in up to 989 participants) measured with the experiments carried out at the Metabolomics Platform of the Genome Analysis Center at the Helmholtz Zentrum München, Germany, per manufacturer instructions (15). The laboratories had no access to phenotype information.

We assessed the pairwise partial correlation between each metabolite and each glycemic trait (i.e., fasting glucose, fasting insulin, HOMA of insulin resistance [HOMA-IR], BMI, and waist-to-hip ratio [WHR]) in the group of participants without diabetes. We included age, sex, and lipid-lowering medication as covariates. Bonferroni correction was applied based on the number of independent vectors in the data. By the matrix spectral decomposition method (21), we estimated 191 independent vectors using the pairwise bivariate correlation matrix of the 562 metabolites. A P value < 5.24 × 10⁻⁵ (0.05/191/5) adjusted by number of independent vectors and number of outcomes was used as the threshold for metabolome-wide significance. The metabolites associated with glucose in the ERF study were taken forward (n = 124) as candidates for MR. In this set of 124 metabolites, we also tested the association with type 2 diabetes using logistic regression.

For each metabolite associated with glucose, we performed two-sample bidirectional MR. The same method on two-sample MR has been performed in the previous MR studies on type 2 diabetes (6,9,22). We tested whether genetically varying levels of a particular metabolite affect the risk for elevated glucose and type 2 diabetes (we call this the forward approach) and whether genetically increased risk of type 2 diabetes or elevated glucose is associated with circulating levels of a particular metabolite (we call this the backward approach). The associations between the instrumental variables and the exposure and the outcome are estimated from different studies, either the metabolite GWAS (14–16,19) or fasting glucose/type 2 diabetes GWAS published by Meta-Analyses of Glucose and Insulin-related traits Consortium (MAGIC) and DIAbetes Genetics Replication and Meta-analysis (DIAGRAM) Consortium (23,24), using the genetic risk score method. The effect of the genetic risk score was constructed by summing up the weighted effects of genome-wide significant SNPs on the exposure variable in relation to their effects on the outcome, as detailed in a previous publication (6). This was performed using summary statistics level data utilizing the method described by Dastani et al. (25) and implemented in the R package gtx. Figure 2 shows the overview of the instrumental variable construction. All SNPs were mapped to human genome build hg19. Given that MR assumes no pleiotropic effect beyond that on the risk factor of interest (i.e., exposure), we excluded the top SNPs from previously published BMI and WHR GWAS (26,27), as well as any SNPs within a 1-Mbp window distance of these, from the genetic score. We additionally excluded the genetic loci (1-Mbp window) of glucose, type 2 diabetes, insulin, and HOMA-IR extracted from previous publications (23,24,28) in the forward MR and the genetic loci (1-Mbp window) of the particular metabolite of interest, using the published GWAS information in the backward MR (see Supplementary Table 1 for a list of genetic loci excluded at this stage). We restricted the SNP lists to a set of independent SNPs in low linkage disequilibrium (pairwise R² <0.05) for each test (29) based on the genotype data in ERF. SNPs with disproportionate effects in the risk score were excluded to reduce pleiotropy (see Supplementary Table 2 for a list of SNPs excluded at this stage). Genetic risk scores comprising >5 SNPs that explain >1% of variance in exposure were taken forward. This effort yielded 20 metabolite–glucose/type 2 diabetes sets in the forward MR and 76 glucose–metabolite sets and 79 type 2 diabetes–metabolite sets in the backward MR. A false discovery rate (FDR) of 0.05 was used as the significance threshold for the four series (i.e., metabolite–glucose, metabolite–type 2 diabetes, glucose–metabolite, and type 2 diabetes–metabolite series).

Although we applied several restrictions on the SNPs in the genetic risk scores as explained above, the instrumental variable assumption that the locus is associated with the outcome only via the association with the exposure (Fig. 1) is still hard to justify in practice. We harnessed the extensive background biological knowledge available to make the additional semiparametric assumption to get the MR estimates of the causal effect. That is, for each set, we evaluated whether we could identify the gene in proximity to the locus that could explain the change in exposure levels. If a gene codes for an enzyme that catalyzes the exposure or a related compound or if it is present in a signaling cascade that affects the exposure, we assumed that the link between the instrumental variable and the exposure was direct and not mediated by the outcome. For the forward approach, we checked the biological link between the locus and the target metabolite, and for the backward approach, we checked the link with glucose. As the pathway in the disease type 2 diabetes is complex, we did not check the biological link with type 2 diabetes in the backward approach. If the gene directly links to the exposure, the related SNPs are taken forward to calculate the genetic risk score. Next, MR is performed for any genetic risk score (comprising >5 SNPs) that explains >1% of variance in exposure. To explore potential mechanistic links between the locus and the exposure, we used an automated workflow that was developed in house to gather gene-specific knowledge of all genes in proximity to each locus. In detail, we downloaded a number of online databases from the respective FTP servers and integrated them offline in MATLAB. Subsequently, for each SNP we selected genes within a window of 100 kbp, with coordinates based on the dbSNP (30) and National Center for Biotechnology Information Gene (http://www.ncbi.nlm.nih.gov/gene; GRCh37), as well as genes whose expression is affected by the locus (GTEx-eQTL). Next, for each gene we gathered protein-related knowledge from UniProtKB (31) and affected pathways from ConsensusPathDB (32). Finally, for each protein we investigated metabolic activity by checking whether it concerned a transporter protein in TCDB (33) or enzyme in ExPASy (34) and, if so, checking what the catalyzed metabolic reaction was in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (35). For the KEGG database, the last freely available version was used. The database integration pipeline generated one HTML file for each locus, containing gene-specific knowledge and hyperlinks to the original database entries, which was then inspected to find a mechanistic link with the exposure. The strength of this approach is that it identifies loci for which the instrumental variable assumption can be validated using genetic and biochemical evidence from online databases. We have successfully applied this workflow in earlier studies (15,19,36–37). Heemskerk et al. (37) give the best example of the power of our method, where we reanalyzed published results of a GWAS on metabolite levels (38) and confirmed the annotation by an in vitro experiment.

Characteristics of the current study population are given in Table 1. Participants with type 2 diabetes were older and tended to be more often male and more likely to be on lipid-lowering medication. They had higher BMI, WHR, systolic blood pressure, triglycerides, fasting glucose, insulin, and HOMA-IR and lower HDL cholesterol, adiponectin, and LDL cholesterol.

We identified 124 metabolites that observationally associate with fasting glucose in the control population with P value < 5.24 × 10⁻⁵ (Fig. 3). These consisted of 36 phospholipids (Fig. 3A), 20 triglycerides (Fig. 3B), 24 small molecular compounds (Fig. 3C), and 44 lipoprotein particle subfractions (Fig. 3D). Correlation coefficients and P values as well as the overlap with previous research for all 124 metabolites are given in Supplementary Table 3. A clustered heatmap of correlation structure in between the 124 selected metabolites is shown in Supplementary Fig. 1. Among the 124, 112 of them also associated with type 2 diabetes (P < 0.05) and their associations with type 2 diabetes and glucose were in the same direction. In addition, their associations with BMI, WHR, fasting insulin, and HOMA-IR were in line with the direction of their associations with glucose. Out of the 124 metabolites, 90 correlated positively and 34 correlated negatively with fasting glucose. We observed negative correlation between glucose and alkyl-acyl and diacyl phosphatidylcholines (mostly of the polyunsaturated type), lysophosphatidylcholines (mostly of the saturated type), and parts of the lipoprotein subfractions from LDL and HDL. These lipoprotein subfractions particularly consisted of lipid components of extra large (XL) and large (L) LDL particles and XL-HDL and L-HDL particles, as well as total HDL measurements. The second cluster of metabolites that we observed to correlate positively with glucose included several phospholipids, phosphotidylcholines, phosphatidylethanolamines, and lysophosphatidylcholines. Amino acids and low–molecular weight compounds also correlated positively with glucose, in addition to lipid side groups and triglycerides. Finally, from the lipoprotein subfractions, small (S), extra small (XS), medium (M), and XL-VLDL particles, as well as the total VLDL components, followed by IDL and LDL triglycerides, XS-LDL to M-LDL particle components, and the ApoA1 and triglyceride content of S-HDL particles were correlated positively with fasting glucose in the population without diabetes.

Table 2 shows the significant results from the association of the relevant metabolites with fasting glucose using MR. Among the 20 eligible metabolite–glucose/type 2 diabetes sets, genetically decreased levels of 8 metabolites associated significantly with fasting glucose (FDR <0.05). These include XL-HDL cholesterol (FDR = 0.03), XL-HDL phospholipids (FDR = 2.76 × 10⁻³), XS-VLDL phospholipids (FDR = 0.04), XL-HDL–free cholesterol (FDR = 0.01), L-HDL cholesterol (FDR = 0.01), L-HDL–free cholesterol (FDR = 2.76 × 10⁻³), HDL cholesterol (FDR = 0.04), and IDL phospholipids (FDR = 0.04). After the pathway-based subset analysis, a causal role for IDL phospholipids was not supported (FDR = 0.17). At the same time, pathway-based sensitivity analysis revealed possibly causal roles for three additional metabolic markers, including S-VLDL triglycerides (FDR = 0.04), S-HDL triglycerides (FDR = 0.04), and plasma triglycerides (FDR = 0.04). Table 3 shows the suggested causal effects of metabolites on type 2 diabetes, i.e., XS-VLDL phospholipids, IDL phospholipids, and plasma triglycerides. Interestingly, the statistical significance for both XS-VLDL phospholipids and IDL phospholipids in the initial results is filtered out after the sensitivity analysis (FDR: XS-VLDL phospholipids 0.03 vs. 0.31; IDL phospholipids 0.01 vs. 0.24), while plasma triglycerides shift to being borderline significant (FDR = 0.07 vs. 0.046). The results from the full lists of performed forward MR tests are given in Supplementary Table 4, and the SNPs included in all the genetic risk scores are given in Supplementary Table 5.

The significant results of the backward MR are shown in Table 4. We found that genetic predisposition to type 2 diabetes is associated with lower levels of phosphatidylcholine alkyl-acyl 42:5 (FDR = 0.02) and phosphatidylcholine alkyl-acyl 44:4 (FDR = 0.02) and higher levels of alanine (FDR = 0.02). The details of all the tested SNP sets are shown in Supplementary Table 6 and Supplementary Table 7. No possible causal role for glucose was supported. As the genetic risk scores of the glucose explained <1% of variance, the backward MR with pathway analysis was not performed. Figure 4 displays the suggested paths discovered by the MR approach after the pathway-based sensitivity analysis. Overall, the associations estimated by MR were in the consistent direction with the observed associations in the ERF study.

We selected 124 metabolites that are correlated with glucose in the population without diabetes, and using MR, we tested whether this metabolic profile points to any causal paths involved in glucose level or type 2 diabetes. Combining metabolomics and MR, we detected 14 candidate causal associations: 10 metabolites influencing fasting glucose, 1 influencing type 2 diabetes, and 3 influenced by type 2 diabetes.

Our initial observational association tests yielded correlation estimates within the expected range of power calculations for the 124 glucose-associated metabolites. To our knowledge, 35 of these metabolites were previously shown to be associated with glucose or type 2 diabetes, including 31 concordant and 4 discordant results (Supplementary Table 3) in studies with very limited sample size (39,40). Our significant results on subfractions of lipoproteins yielded resolution on the established association of dyslipidemia, especially for the HDL subfractions.

One of our main findings is that genetically increased cholesterol and free cholesterol contents of circulating XL-HDL and L-HDL particles, XL-HDL phospholipids, and HDL cholesterol, together with XS-VLDL phospholipids, associate with decreased glucose levels. Our second finding is that triglyceride contents of S-HDL and S-VLDL particles, as well as total plasma triglycerides, seem to have a glucose-increasing causal effect and, considering the total triglycerides, this effect has been extended to the outcome type 2 diabetes. Finally, we showed that genetic predisposition to type 2 diabetes associates with lower levels of two alkyl-acyl phosphatidylcholines and a higher level of alanine. Our report is the first using higher-resolution (metabolomics-driven), lipoprotein-based exposure variables. Hence, no other study exists for comparison, although HDL cholesterol, LDL cholesterol, and total triglycerides (from routine biochemistry) have been previously studied as exposure variables for MR in order to understand their causal effects on type 2 diabetes and glucose (an overview is given in Table 5). Our method is similar to the method of White et al. (9) and Fall et al. (6) in terms of the application of the genetic risk score function utilizing the DIAGRAM/MAGIC data sets. White et al. (9) showed that high levels of all three blood lipids (HDL cholesterol, LDL cholesterol, and plasma triglycerides) were genetically associated with a lower risk of diabetes, although the results for triglycerides were inconsistent. However, the study did not consider the genetic variants that might be involved in the confounding phenotypes such as BMI or WHR, nor did it exclude the SNPs that are involved in type 2 diabetes directly. Fall et al. (6) showed that the association between total HDL cholesterol risk score and low fasting glucose was attenuated when adjusted for the effects of SNPs on LDL cholesterol, triglycerides, and surrogates of adiposity. Different from the two studies mentioned above, the MR in our study was done in a broad spectrum of metabolites and included a detailed subclassification of lipoproteins that have not been tested before. Using such high-resolution phenotypes, we demonstrate that the decreasing effect of HDL cholesterol on fasting glucose is more specific to the L-HDL or XL-HDL subclasses, whereas for S-HDL triglycerides, an increasing effect exists. These results advocate that a higher resolution of high density in lipoproteins may reveal the observed epidemiological associations or biological functions of HDL cholesterol more accurately and will uncover the mystery of complex lipids such as HDL. Certain HDL subfractions and characteristics of these subfractions may have independent associations with glucose, particularly for the small versus large size particles. Such a different role for HDL triglycerides and HDL large fractions occurred upon sleeve gastrectomy of obese patients and was associated with reduced insulin resistance and HDL remodeling (41). In addition, as experimentally shown, HDL indeed may mediate glucose regulation in the pathophysiology of type 2 diabetes (42). Suggested mechanisms include the following (43): 1) insulin secretion from pancreatic β-cells combating cellular lipid accumulation and lipotoxicity (44) and endoplasmic reticulum stress and apoptosis (45,46), 2) insulin-independent glucose uptake by muscle via the AMP-activated protein kinase (47) and calcium/calmodulin activated protein kinase (48), and 3) insulin sensitivity (49). The Investigation of Lipid Level Management to Understand its Impact in Atherosclerotic Events (ILLUMINATE) trial (50) demonstrated that in a subgroup of participants with diabetes, statin treatment led to increased glucose levels, while this effect was not observed in participants treated with combination of statin and CETP inhibitor torcetrapib, suggesting that CETP inhibition and consequent HDL cholesterol elevation may improve glycemic control in patients with diabetes. It is of note that the CETP gene is a major determinant of XL-HDL and was included in our MR experiment.

We have detected three associations potentially pointing out an influence of type 2 diabetes over the metabolome. The first two are long-chain polyunsaturated alkyl-acyl phosphatidylcholines that are decreased in type 2 diabetes. This is interesting considering our previous report, which showed that three shorter-chain polyunsaturated alkyl-acyl phosphatidylcholines are increased in patients with type 2 diabetes and decreased in patients using the glucose-lowering drug metformin (51). The other molecule affected by diabetes was alanine, which is a nonessential amino acid and can be synthesized in the body from pyruvate and branched-chain amino acids such as valine, leucine, and isoleucine. Alanine has been previously implicated in glucose response (52). The enzyme alanine aminotransferase (ALT) catalyzes the conversion of alanine to pyruvate and glutamate, and high levels of ALT indicate liver damage.

Our study differs from previous reports in three ways. First, we used a bidirectional approach and included a wide range of molecular markers to be tested, using high-resolution phenotypes, measured by MS or NMR. Second, we exploited pathway knowledge that was gathered through an automated workflow to perform subset analysis in MR. Statistical methods that deal with pleiotropy in MR analyses, such as the Egger regression method (9,53,54), exist but are still being refined. They all rely on additional strong assumptions about the unobserved pleiotropy, such as the InSIDE assumption, and are sensitive to violations of these assumptions. They also suffer from a lack of power. However, one can harness the available genetic and biological knowledge in online databases in order to maximize the uniqueness of the genetic risk score for the exposure variable for this purpose and to validate the chosen instruments. It has to be mentioned that although powerful for most metabolites, our approach with the genetic and biological knowledge is also conservative because it ultimately relies on the comprehensiveness of the content of the databases that are included. As a consequence, all loci for which no strong evidence is present that a nearby gene directly affects the exposure, e.g., because the involved gene is affected through a yet-unknown regulatory mechanism, are excluded. Considering glucose for which the instrumental strength was initially lower compared with the others, the pathway approach yielded lower explained variance in exposure (R² <1%). While one can argue that this would lead to lack of power, it may also reflect the fact that polygenic traits such as glucose may not be the most suitable exposure variables for an MR analysis. To limit this, we utilized the large population–based GWAS of broad-spectrum metabolites and fasting glucose/type 2 diabetes with the combined-instrument MR approach (25). We want to point out that although we controlled the pleiotropic effects between the outcome and exposure by 1) excluding the known predictors, 2) heterogeneity tests, and finally 3) pathway analysis, we cannot exclude a correlation between the genetic instruments tested, especially for the HDL subfractions, for which the coding genes overlap. While effect alleles were weighted by their original effects estimates from each GWAS (of exposure variables), there was strong overlap in the SNPs used for different lipid subfractions, meaning the genetic instruments were not highly specific to these subfractions.

In conclusion, using MR, the current study provides evidence for potentially causal metabolic paths of glucose homeostasis and type 2 diabetes. Our results indicate that an increase of large HDL particles might have a decreasing effect on glucose, while an increase of small HDL particles might have an increasing effect, refining earlier MR findings suggesting a possible causal effect of HDL on glucose levels as well as pointing out these particles as targets for glucose management. We further found evidence that type 2 diabetes may alter levels of alkyl-acyl phosphatidylcholines and alanine, which also here can be translated into prevention of disease complications and prognosis.

Acknowledgments. The authors are grateful to all study participants and their relatives, general practitioners, and neurologists for their contributions to the ERF study; P. Veraart (Department of Epidemiology, Erasmus Medical Center, Rotterdam, the Netherlands) for her help in genealogy; J. Vergeer (Department of Epidemiology, Erasmus Medical Center, Rotterdam, the Netherlands) for the supervision of the laboratory work; and P. Snijders (Department of Epidemiology, Erasmus Medical Center, Rotterdam, the Netherlands) for his help in data collection.

Funding. The ERF study was supported by the Netherlands Consortium for Systems Biology, within the framework of the Netherlands Genomics Initiative/Netherlands Organisation for Scientific Research. The ERF study as a part of the European Special Populations Research Network was supported by European Commission Sixth Framework Programme STRP grant no. 018947 (LSHG-CT-2006-01947) and also received funding from the European Community Seventh Framework Programme (FP7/2007-2013)/grant agreement HEALTH-F4-2007-201413 by the European Commission under the program “Quality of Life and Management of the Living Resources” of Fifth Framework Programme (no. QLG2-CT-2002-01254), as well as FP7 project EUROHEADPAIN (no. 602633). High-throughput genetic analysis of the ERF data was supported by a joint grant from the Netherlands Organisation for Scientific Research and the Russian Foundation for Basic Research (NWO-RFBR 047.017.043). High-throughput metabolomics measurements of the ERF study have been supported by Biobanking and BioMolecular resources Research Infrastructure the Netherlands (BBMRI-NL). Biocrates platform measurements were supported by the European Community's Seventh Framework Programme (FP7/2007–2013), ENGAGE Consortium, grant agreement HEALTH-F4-2007-201413. Lipidomics analysis was supported by the European Commission FP7 grant LipidomicNet (2007-202272). J.L., C.M.v.D., and A.D. have used exchange grants from the Personalized pREvention of Chronic DIseases consortium (PRECeDI). S.S. has been awarded the Erasmus Mundus–Western Balkans (ERAWEB) mobility program academic scholarship. A.D. is supported by a Veni grant (2015) from ZonMw.

The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Duality of Interest. No potential conflicts of interest relevant to this article were reported.

Author Contributions. J.L. and J.B.v.K. analyzed the data. J.L., J.B.v.K., S.S., K.W.v.D., A.V., T.H., A.C.H., E.S., N.A.S., C.M.v.D., and A.D. explained the results and wrote the manuscript. A.V., T.H., and A.C.H. generated the metabolomics data. C.M.v.D. and A.D. designed the study. J.L. 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.