Maternal Dietary Glycemic Index and Glycemic Load in Pregnancy and Offspring Cord Blood DNA Methylation

Glycemic index and glycemic load are measures to classify carbohydrate-containing foods according to the blood glucose response they evoke after consumption (1). Glycemic index represents the postprandial rise in blood glucose in response to a specific carbohydrate- containing food, fed as 50 g of net carbohydrate and compared with 50 g of pure glucose, and expressed as a percentage of the postprandial response to glucose. A higher glycemic index corresponds to a higher blood glucose–raising property of the food. Glycemic load incorporates both the glucose-raising properties of a food and the consumed amount of that food (2). Thus, a certain food always has the same glycemic index, but the glycemic load differs based on the amount consumed.

Individuals with a low glycemic–index diet or a low glycemic–load diet have a decreased risk of type 2 diabetes (3) and improved metabolic health (2,4,5). Additionally, compared with a conventional or high glycemic–index diet, a low glycemic–index diet reduces glycated hemoglobin (HbA_1c) in individuals with diabetes (6).

Evidence on the association of prenatal dietary glycemic index or load with neonatal outcomes is mixed (7–9). Of three previous observational studies, one reported a positive association of maternal dietary glycemic index with offspring adiposity measures at birth (8), one found negative associations (7), and one found no associations (9). One study in obese pregnant women did not find associations of a diet and physical activity intervention with offspring large-for-gestational-age status, despite, among other factors, a reduced dietary glycemic load (10).

A potential biological mechanism underlying associations of maternal diet with offspring health is differential DNA methylation. Some large-scale meta-analyses of epigenome-wide association studies (EWAS) have previously linked adverse maternal nutrition-related exposures during pregnancy (e.g., BMI [11] and plasma folate concentrations [12]) with cord blood DNA methylation. However, not much is known about the association of maternal dietary glycemic index or load during pregnancy with offspring DNA methylation. Two previous studies provided some first evidence that maternal dietary glycemic index and load may be related to DNA methylation in placental tissue (13) and in cord blood (14), although sample sizes were rather limited. DNA methylation may have functional consequences via changes in gene expression that potentially lead to changes in health outcomes. We therefore conducted a meta-analysis to investigate the association of maternal dietary glycemic index and load during pregnancy with offspring cord blood DNA methylation. In this EWAS, we took a hypothesis-free approach, analyzing DNA methylation at all cytosine-phosphate-guanine sites (CpGs) measured on the array in relation to glycemic index and load. After running the EWAS, we compared the results with previous EWAS of cardiometabolic phenotypes to examine whether the CpGs identified in our study were also related to cardiometabolic phenotypes. Lastly, because women with overweight are more likely to be insulin resistant (15), we therefore additionally stratified analyses by maternal weight status.

An explanation of the cohorts used for meta-analysis and all following statistical analyses is presented in the flowchart in Fig. 1.

Three population-based birth cohorts participated in this individual participant meta-analysis: 658 mother-offspring pairs from the Avon Longitudinal Study of Parents and Children (ALSPAC) from the U.K. (16,17), 998 mother-offspring pairs from the Generation R Study (Generation R) from the Netherlands (18), and 347 mother-offspring pairs from the INfancia y Medio Ambiente (INMA) Project from Spain (19). We excluded mothers with gestational diabetes or prepregnancy type 1 or 2 diabetes, because glycemic response changes with insulin resistance and with antidiabetes medication. After cohort-specific quality control on the food frequency questionnaires (FFQs), described in the cohort-specific Supplementary Methods, we additionally excluded mothers with maternal glycemic index and load values outside ±5 SD from the cohort mean. We further excluded multiple pregnancies and if mothers had multiple children (i.e., nontwin siblings) in the cohort. We included only one child per mother, based on completeness of data and, if equal, randomly. Subjects with missing data on any of the covariates were excluded from analyses.

GI_i (%) = the glycemic index of food i, relative to glucose; CHO_i = the amount of food i consumed (g/day, from the raw FFQ data set) × carbohydrate content from food i (g/g, from the reference database); and n = number of foods eaten per day.

Offspring cord blood DNA was isolated from cord blood drawn directly after birth, and DNA methylation was measured using the Illumina 450K array. Each cohort performed their own laboratory analyses, quality control, and normalization (Supplementary Methods). We used untransformed β values, with values between 0 and 1, as the outcome measure and excluded extreme DNA methylation outliers according to the Tukey method: outliers outside the range of (25th percentile − 3∗interquartile range [IQR]) to (75th percentile + 3∗IQR) were excluded. We removed control probes, X and Y chromosome probes, and cross-reactive probes (22,23), resulting in a maximum of 428,328 remaining CpG sites. After meta-analysis, we flagged potentially polymorphic sites (22,23) or methylation quantitative trait loci (mQTLs) (24).

Each cohort ran EWAS using robust linear regression models for glycemic index and load with each CpG site individually, using lmFit() in the limma R package (25), according to a prespecified analysis plan and R script. The fully adjusted models had maternal glycemic index or load as the exposure and cord blood DNA methylation as the outcome and were adjusted for sex, maternal education, age, smoking during pregnancy, total energy intake, batch, and cell types. We additionally ran these models stratified on maternal pre-pregnancy or early pregnancy (≤16 weeks of gestation) normal body weight (BMI 18 to <25 kg/m²) and overweight or obesity (BMI ≥25 kg/m²). Mothers with underweight were removed from these stratified analyses because the sample size was too low.

Sex of the child was obtained from birth records. Maternal educational level was based on self-reported questionnaires during pregnancy and assigned to two (ALSPAC and Generation R) or three (INMA) levels, depending on cohort preference (Supplementary Methods). Maternal age was self-reported in the questionnaire during pregnancy. Maternal smoking during pregnancy was self-reported in the questionnaire during pregnancy and defined as no smoking during pregnancy, as stopped smoking before the second trimester, as sustained smoking in Generation R and INMA, and as sustained smoking versus no smoking or quit before second trimester in ALSPAC. Maternal total energy intake during pregnancy (kcal/day) was calculated from the self-reported FFQ during pregnancy that was also used to calculate the glycemic index and load. Each cohort included their preferred variable to adjust for batch: Generation R used plate number as reported by the laboratory, ALSPAC used surrogate variable analysis (26), and INMA used ComBat (27). We used a cord blood methylation-based reference set to estimate cell type composition (28) (CD8 T cells, CD4 T cells, natural killer cells, B cells, monocytes, granulocytes, nucleated red blood cells).

All cohort EWAS results passed quality control using the QCEWAS R package (29). We then meta-analyzed the cohort-specific EWAS results using the fixed-effects inverse variance weighted method in METAL (30) at Erasmus MC and consecutively shadowed this meta-analysis at ISGlobal using GWAMA (31). Results were confirmed. We used a Bonferroni-corrected P-value threshold of <1.17 × 10⁻⁷ to define CpGs to take forward for further analyses, and we also present the false discovery rate (FDR).

We tested whether the CpGs that reached the Bonferroni threshold in the full group meta-analyses or in the analyses stratified on maternal BMI showed persistent associations with maternal glycemic index or load during pregnancy when DNA methylation was measured at later ages. We performed lookups in 1) early childhood: meta-analysis of INMA (mean: 4.5 years) and Generation R (6.1 years), 2) late childhood: meta-analysis of ALSPAC (7.5 years), INMA (8.8 years), Generation R (9.8 years), and 3) adolescence: ALSPAC (17.1 years). These analyses were adjusted for the same covariates as the cord blood analyses plus child age at DNA methylation measurement. For cell types, we used the adult blood methylation-based reference (32,33) (CD8 T ells, CD4 T cells, natural killer, B cells, monocytes, and granulocytes). We corrected for multiple testing in these lookups using a Bonferroni P-value cutoff of <0.0012 (0.05/41 tests).

We examined potential functionality of the identified CpGs using several analyses. First, we looked up the top hits from the full group and stratified meta-analyses in recently published EWAS on related exposures: maternal late-pregnancy glycemic traits (34), maternal early-pregnancy blood glucose and insulin concentrations (35), and maternal early-pregnancy BMI (11); and in EWAS of cord blood DNA methylation and child health outcomes: birth weight (36) and BMI in childhood or adolescence (37). We corrected for multiple testing in this lookup using a Bonferroni P-value cutoff <0.0012 (0.05/41 tests).

Second, we used the cell type-adjusted Human Early Life Exposome (HELIX) catalog of 13.6 million autosomal expression quantitative trait methylation (cis-eQTMs) in children to check whether the top CpGs from our EWAS were associated with gene expression in childhood blood (38). We additionally tested associations of the top CpGs from our EWAS with RNA transcript levels in subcutaneous adipose tissue in 223 participants (0–21 years old) from the Leipzig Childhood Adipose Tissue Cohort (39). Participating children underwent elective surgery, mostly orthopedic, and were free of severe diseases (diabetes, generalized inflammation, malignant disease, and genetic syndromes, or were not permanently immobilized) and medication potentially influencing adipose tissue biology. We focused on CpGs ≤500 kilobase of the transcription start site of the transcript. We applied Bonferroni correction using the total number of CpGs tested (33, since 8 CpGs were discarded prior to the statistical analysis) and the number of tests per each unique CpG-transcript pair. Detailed methods are described in the Supplementary Methods.

Third, we used missMethyl (40) to run functional enrichment pathway analyses using gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG). For these pathway analyses we used several selections of CpGs: 1) all CpGs that reached Bonferroni threshold in either the full group or stratified analyses for both glycemic index and load combined, to explore common pathways; and 2) CpGs with P < 0.0001 in each of the six analyses of maternal glycemic index or load in the full group or the stratified analyses separately.

Fourth, we explored in the EWAS Catalog (41) if the top hits from the full group meta-analyses were previously identified to be associated with other traits.

Fifth, for those CpGs that reached the Bonferroni threshold in the full group and stratified meta-analyses and that were flagged as potentially polymorphic sites (22,23) or methylation quantitative trait loci (mQTLs) (24), we created distribution plots using the Generation R data to check bimodality of the distribution.

Sixth, using all CpGs associated with glycemic index and load in the full group and stratified meta-analyses, we examined enrichment for tissue-specific DNaseI hypersensitivity regions using eForge version 2.0 (42).

Lastly, we examined between-cohort heterogeneity of the top CpGs with a combination of the I² statistic, forest plots, and leave-one-out plots.

We included 2,003 mother-child pairs from the three cohorts. Descriptives per cohort are presented in Table 1, Supplementary Fig. 1A–C, and the Supplementary Methods. The mean maternal glycemic index ranged between 50.8 and 58.3, and the mean maternal glycemic load ranged between 107.6 and 154.6.

In the full-group analyses, we observed an association between maternal glycemic index and cord blood DNA methylation level of cg21301148 (0.099% increase in DNA methylation per 1-unit increase in glycemic index, SE = 0.019, P = 1.00 × 10⁻⁷ (Table 2 and Supplementary Table 1). Maternal glycemic load was associated with cord blood DNA methylation levels of cg09874107, cg21301148, and cg27528695 (Table 2 and Supplementary Table 2). Manhattan plots are presented in Fig. 2.

In 1,509 mothers with normal weight, maternal glycemic index was associated with two CpGs, cg02079551 and cg21301148, and maternal glycemic load was associated with two CpGs, cg02920421 and cg26729101 (Table 2 and Supplementary Tables 3 and 4). In 463 mothers with overweight or obesity, maternal glycemic index was associated with DNA methylation of 31 CpGs, and maternal glycemic load was associated with DNA methylation of 4 CpGs (Table 2). Full epigenome-wide results for these stratified analyses are presented in Supplementary Tables 5 and 6. We found no overlap between the hits from these stratified analyses, except for cg21301148, which was associated with the glycemic index in the full group as well as in mothers with normal weight. Of the total of 39 CpGs that were identified in at least one of the stratified analyses, 19 (48.7%) showed opposite directions of association between the mothers with normal weight and those with overweight or obesity. Manhattan plots are presented in Fig. 2. In total, glycemic index or load were associated with 41 CpGs in either the full group or the stratified analyses. Results for all (full group and stratified) meta-analyses for these 41 CpGs are shown in Supplementary Table 7.

We tested persistence of associations of maternal glycemic index or load during pregnancy when DNA methylation was measured in early childhood (n = 554), late childhood (n = 1,255), and adolescence (n = 711). None of the CpGs persisted at later ages in relation to the same exposure. Although 7 of the 41 CpGs reached Bonferroni-corrected P < 0.0012 in one or more of the models in early childhood, late childhood, or adolescence (bold and underlined in Supplementary Tables 8–10), these associations were observed in other models than the initial meta-analysis in cord blood (e.g., these were associated with glycemic load in the lookup rather than glycemic index in the cord blood analysis). Additionally, nine CpGs reached nominal significance in the same model as in the main cord blood analyses.

None of the 41 CpGs from the full group and stratified analyses was associated with maternal early and late pregnancy glycemic traits or maternal BMI in previous EWAS (all P values >0.0012 [0.05/41]) (Supplementary Table 11). Similarly, none of the 41 CpGs was previously associated with birth weight or childhood or adolescence BMI (Supplementary Table 11).

In the eQTM analyses, we found three CpG-transcript cluster associations in childhood blood for cg24458009, with the PC-esterase domain containing 1B (PCED1B) gene annotated to this transcript cluster. We further found one CpG-transcript cluster association for cg23347399, with the protocadherin γ gene family (PCDHGA, PCDHGB, PCDHGC) annotated to this transcript cluster (Table 3). These two CpGs were both found in the EWAS for associations with glycemic index in mothers with overweight or obesity. In the analyses of associations between DNA methylation and RNA transcript levels in adipose tissue, 4 of the 561 CpG-transcript pairs tested reached statistical significance. These all involved cg27193519, showing a direct association between DNA methylation and expression levels of transcription factor AP-4 (TFAP4), zinc finger protein 500 (ZNF500), periplakin (PPL), and ankyrin repeat and sterile α motif domain containing 3 (ANKS3). Correlation plots of these four CpG-transcript pairs are shown in Supplementary Fig. 2.

Functional enrichment analyses on all 41 CpGs from the full group or stratified meta-analyses did not result in FDR significant GO or KEGG pathways (Supplementary Table 12). We further found no functional enrichment for the 76 and 231 CpGs with P < 0.0001 from the full group meta-analyses on glycemic index and load, respectively. Similarly, there was no functional enrichment for the 116 or 1,176 CpGs with P < 0.0001 from the stratified analyses for glycemic index in mothers with normal weight or overweight or obesity, respectively, nor for the 222 or 265 CpGs with P < 0.0001 from the stratified analyses for glycemic load, respectively.

We looked up the three CpGs from the full-group meta-analyses in the EWAS Catalog and found that cg21301148 was previously associated with (gestational) age, clear cell renal carcinoma, fetal versus adult liver gene expression, and alcohol consumption (Supplementary Table 13). Further, cg27528695 was associated with age and HIV infection, but cg09874107 was not found in the EWAS Catalog.

Distribution plots for all of the 41 CpGs that were flagged as either potentially polymorphic or mQTL mostly showed unimodal distributions. On the basis of visual inspection, only four CpGs had a marginal indication of bimodality: cg00006032, cg26474288, cg02725014, and cg06591466 (Supplementary Figs. 3–5). Further interpretation of these four CpGs needs consideration of potentially polymorphic or mQTL effects.

Finally, using the set of 41 CpGs, we found no evidence of enrichment for tissue-specific DNaseI hypersensitivity regions.

The I² value was >50 for 20 of the 41 CpGs, which might indicate between-cohort heterogeneity. To examine this in further detail, we created forest plots and leave-one-out plots for all 41 CpGs (Supplementary Figs. 6–17). Based on the forest plots, we found no indication that one cohort consistently caused heterogeneity in the meta-analysis. Leaving out one cohort at a time, the direction of effect never changed, as can be seen in the leave-one-out plots. However, for a limited number of CpGs, a combined exploration of the plots and I² statistic indicated one cohort had relatively strongly influenced the meta-analysis result. A sensitivity analysis rerunning the functional enrichment analysis only based on the 21 CpGs that reached the Bonferroni threshold and had an I² <50 also did not result in FDR significant GO or KEGG pathways. Of the three CpGs that were associated with gene expression, two had an I² >50.

This is the first large-scale EWAS to study associations of maternal dietary glycemic index and load during pregnancy with offspring DNA methylation. Overall, 41 unique CpGs were associated with glycemic index or load. We identified 3 of these in the full group analyses, 4 were identified when the sample was restricted to mothers with a normal weight, and 35 were identified when the sample was restricted to mothers with overweight or obesity. We found no persistence of the associations into childhood or adolescence. There was no overlap between our top CpGs and CpGs associated with maternal glycemic traits or BMI, or with birth weight or childhood BMI. We found no common pathways in functional enrichment analyses. However, two CpGs, associated with glycemic index in mothers with overweight or obesity were associated with gene expression in child blood, and one CpG was associated with gene expression in adipose tissue. On the basis of exploration of I², forest plots, and leave-one-out plots, some CpGs need careful interpretation. For example, especially interesting for interpretation of our expression findings are cg24458009 and cg27193519; these might have been driven by one cohort.

In the full-group meta-analyses, maternal glycemic index or load were associated with three CpGs: cg21301148, cg09874107, and cg27528695. First, cg21301148 was associated with both glycemic index and glycemic load, and it is annotated to the transcription start site of the myosin heavy chain 6 (MYH6) gene. MYH6 has previously been associated with cardiomyopathy (43). Second, cg09874107, associated with glycemic load, is annotated to the transcription start site of the chromosome 3 open reading frame 70 (C3orf70) gene, which may play a role in neural and neurobehavioral development (44). DNA methylation of another CpG in C3orf70 has been associated with fat-free mass in European children (45). Third, cg27528695, associated with glycemic load, is located close to the small Cajal body-specific RNA 2 (SCARNA2) gene. To the best of our knowledge, SCARNA2 has not been associated with cardiometabolic disease.

Because women with overweight are more likely to be insulin resistant (15) and opposite directions of association with blood glucose and insulin were previously observed between mothers with normal weight versus mothers with overweight or obesity (35), we tested a potential modifying effect of maternal weight status on the association between glycemic index or load and cord blood DNA methylation by running analyses stratified on maternal BMI categories. Only one CpG overlapped between two of our EWAS models: cg21301148 was associated with glycemic index in mothers with normal weight as well as in the full group. We found most hits in mothers with overweight or obesity. In these mothers, maternal glycemic index and load were associated with DNA methylation levels at 31 and 4 CpGs, respectively. This may indicate a stronger effect of glycemic index or load on biological pathways in offspring of mothers with overweight or obesity, for example, through altered insulin sensitivity in these mothers. Indeed, for 20 of the 39 CpGs from the stratified analyses, we observed stronger effect sizes in mothers with overweight or obesity. For 19 of the 39 CpGs, the effect sizes were in opposite directions between strata, explaining the fact that there were no findings for these CpGs in the full group.

Interestingly, the three CpGs that were associated with gene expression were all found in the EWAS for associations with glycemic index in mothers with overweight or obesity. In child blood, cg24458009 was associated with expression of PCED1B, and cg23347399 was associated with expression of PCDHG. These genes do not have a known function in processes associated with cardiometabolic phenotypes (46,47). In subcutaneous adipose tissue, cg27193519 was associated with expression of TFAP4, ZNF500, PPL, and ANKS3. TFAP4 has been associated with lipid metabolism (48). ANKS3 is associated with congenital kidney disease (49). ZNF500 and PPL are not known to be associated with cardiometabolic diseases. However, as mentioned, cg24458009 and cg27193519 showed relatively high heterogeneity and thus need to be carefully interpreted.

We found no overlap between our findings and previously reported associations of maternal glucose, insulin, and the glucose area under the curve for a 2-h oral glucose tolerance test or BMI during pregnancy (11,34,35) with cord blood DNA methylation. This could be because we measured glycemic traits based on dietary intake, while these previous studies measured BMI or different glycemic traits as blood concentrations; moreover, glycemic index values for a food are average for everyone, whereas individual biological responses to specific foods are likely to vary (50). It could also be because we excluded mothers with gestational diabetes, while the two previous glycemic trait EWASs did not (34,41). We could not establish a potentially mediating role of cord blood DNA methylation in the association of maternal glycemic index or load with later child health, since none of the CpGs from our EWAS were previously found to be associated with birth weight (36) or child or adolescent BMI (37). However, power might have been limited, and future studies with larger sample sizes need to investigate this further.

We found no persistence of differential DNA methylation into childhood or adolescence. However, seven CpGs reached the Bonferroni threshold for associations with the other glycemic phenotype (i.e., glycemic index instead of load or in mothers with normal weight instead of mothers with overweight or obesity). This may indicate that DNA methylation at those CpGs is associated with general glycemic processes, potentially with varying strengths at different time points. Additionally, nine CpGs reached nominal significance in the same model as in the main cord blood analyses, which might indicate a lack of power to detect persistence due to lower sample sizes. Future studies might focus on increasing sample size for such a lookup in older ages. We previously observed a similar lack of persistence of differential DNA methylation in relation to birth weight (36). Exposure to maternal factors rapidly diminishes after delivery, and many postnatal exposures may affect DNA methylation, which might explain the lack of persistence. However, even without this persistence, differential DNA methylation in the critical period of organ development may have already initiated structural functional alterations that might have long-term consequences independent of persistent differential DNA methylation or via differential DNA methylation in other tissues.

Limitations of this study include the self-reported dietary data, the difference in number of food items in the FFQs between cohorts, ranging from 47 to 293, and the difference in glycemic index reference databases used in the cohorts. However, the calculation of glycemic index and load and the EWAS analyses were performed in all cohorts according to a prespecified harmonized analysis plan and R script, and with extensive quality control of dietary data. Further, cohorts had comparable distributions of glycemic index and load, and the ranking of mothers based on glycemic index and load will not have differed substantially. Additionally, maternal dietary intake was registered at ∼12 weeks of gestation in Generation R and INMA, with the reference time window being the first trimester, while the FFQ was registered at ∼32 weeks in ALSPAC, with the reference time window “nowadays.” We do not expect this to be a problem, since it has been shown that micronutrient intake and diet quality generally do not change during pregnancy (51). Dietary intake might be influenced by many genetic, environmental, and behavioral factors. We adjusted our analyses for many covariates to exclude potential confounding, but as in all observational studies, residual confounding could still be present. The Illumina 450K array covers <2% of all CpG sites in the human genome. We therefore cannot exclude that maternal glycemic index or load is associated with DNA methylation at unmeasured CpG sites. Further, DNA methylation is tissue specific. Therefore, DNA methylation in pancreatic, liver, or adipose tissue might be more relevant than cord blood, which is also why we did a lookup of associations with gene expression in adipose tissue. However, in the large-scale birth cohorts that we used, cord blood is easily accessible, as opposed to these other tissues. Unfortunately, gene expression data were not available in cord blood. We were also limited in sample size for the stratified analyses, and would therefore recommend future studies to focus on increasing sample size specifically for analyses stratified on maternal weight status.

To conclude, in this meta-analysis we identified associations of maternal glycemic index and load during pregnancy with offspring cord blood DNA methylation levels of 41 CpGs, mostly in mothers with overweight or obesity. These CpGs did not overlap with previously published CpGs associated with maternal glycemic traits or with child birth weight or BMI. Although we found no enrichment for specific biological pathways, three CpGs were associated with expression of PCED1B, PCDHG, TFAP4, ZNF500, PPL, and ANKS3 and might therefore have functional effects. Further study is required to explore functionality in more detail, uncover causality of these findings, and study pathways to offspring health.

Acknowledgments. Cohort-specific acknowledgments are stated in the Supplementary Material.

Funding. This work was funded by the Joint Programming Initiative – A Healthy Diet for a Healthy Life to the NutriPROGRAM consortium: ZonMW, the Netherlands (529051022), MRC, U.K. (MR/S036520/1), and Instituto de Salud Carlos III, Spain (AC18/00006), and the PREcisE consortium: ZonMW, the Netherlands (529051023), German Federal Ministry of Education and Research, Germany (FKZ 01EA1905). Information regarding funding for the contributing cohorts and individual authors can be found in the Supplementary Material.

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

Author Contributions. L.K.K. and S.F.-B. performed the meta-analysis. L.K.K., S.F.-B., G.M., L.J., R.O., J.V., M.C., K.L., E.W.T., A.K., R.G., J.H.M.d.V., V.W.V.J., M.V., G.C.S., and J.F.F. interpreted the data. L.K.K., S.F.B., G.M., L.J., R.O., J.V., M.C., K.L., E.W.T., A.K., R.G., J.H.M.d.V., V.W.V.J., M.V., G.C.S., and J.F.F. contributed to critical revision for important intellectual content. L.K.K., S.F.B., G.M., L.J., R.O., J.V., M.C., K.L., E.W.T., A.K., R.G., J.H.M.d.V., V.W.V.J., M.V., G.C.S., and J.F.F. gave final approval of the version to be published. L.K.K., S.F.B., G.M., L.J., J.V., K.L., and M.C. conducted cohort-specific analyses. L.K.K. and J.F.F. designed the research. L.K.K. and J.F.F. drafted the manuscript. A.K., R.G., V.W.V.J., M.V., G.C.S., and J.F.F. acquired data. L.K.K. and J.F.F. are the guarantors of this work and, as such, had full access to all the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis.