Genetic Susceptibility Determines β-Cell Function and Fasting Glycemia Trajectories Throughout Childhood: A 12-Year Cohort Study (EarlyBird 76)

Diabetes is now one of the most common noncommunicable diseases in the world. The World Health Organization (WHO) reports that diabetes affects >422 million people or 8.5% of the world’s adult population (1). It has been projected that one in every three individuals born in the U.S. in the year 2000 will develop diabetes during their lifetime (2). However, the increasing burden of diabetes will be felt, according to the WHO, most acutely in low- and middle-income countries, with catastrophic economic consequences (1). Therefore, there is a pressing need to improve the prediction and early prevention of diabetes.

Diabetes results from impaired insulin secretion, resistance to the action of insulin, or a combination of both factors. Failure of insulin secretion predominates in type 1 diabetes, whereas resistance to insulin (IR) and relative insulin deficiency characterize type 2 diabetes (T2D). Because of the rising prevalence of obesity, T2D has become increasingly common in children and adolescents (3). Adolescence is a period of high vulnerability to the onset of diabetes in children because of rapid endocrine changes, increasingly accompanied by obesity, and on a variable background of genetic risk (4). Yet, very little is known about the molecular pathways linking genetic risk factors, weight gain, and puberty to the risk of developing diabetes in adolescents. Deciphering these pathways is of fundamental importance to understanding disease progression in children and how prediabetes could be detected at an earlier stage to enable preventive measures to be taken to tackle the worldwide epidemic of diabetes (5).

The EarlyBird study is a landmark prospective cohort study investigating the origins of T2D in children. This cohort of healthy children has been followed from 5 to 16 years of age with annual clinical, anthropometric, and physiological measurements (6). We were among the first to report the occurrence of an early defect in β-cell function among children who go on to develop prediabetes (6). Some 17% of the nominally healthy children in the EarlyBird cohort showed impaired fasting glycemia (IFG) by 15 years of age. Furthermore, children who developed IFG already exhibited higher fasting blood glucose levels at 5 years of age compared with those who did not subsequently develop IFG, and this effect was independent of BMI (6). However, prediabetes did not appear until puberty, when IR was at its highest. This is consistent with the U.S. National Health and Nutrition Examination Survey 2005 to 2006, which also included children of normal weight and suggested that the prevalence of prediabetes in adolescence was strongly influenced by IR (7).

Information from longitudinal studies of healthy-weight children is required to determine how prediabetes results from the interaction of genetic risk factors, weight gain, changing levels of IR and other endocrine parameters during puberty, and other nongenetic risk factors (8). Therefore, the present analysis of the EarlyBird cohort was designed to examine how genetic variants, puberty, and weight gain interact to influence insulin action and blood glucose levels at an early age of high vulnerability to diabetes.

We hypothesized that genes and single nucleotide polymorphisms (SNPs), known to be associated with increased risks of T2D, would be associated with trajectories of fasting glucose, insulin, and HOMA indices of insulin resistance (HOMA-IR) and β-cell function (HOMA-B) in children. We also explored if insulin action may be influenced by genetic variations in the growth hormone–IGF-1 axis and adrenal and sex steroid activity. The aim of this study was to identify specific genes and SNPs associated with these glycemic traits independently of body weight that could identify young people at high risk of diabetes. The long-term goal of this work is to develop risk-modifying interventions before adolescence.

The study was conducted in accordance with the principles of the Declaration of Helsinki II. Ethical approval was granted by the Plymouth Local Research Ethics Committee (1999; Plymouth, U.K.), and parents gave written consent and children verbal assent. The EarlyBird Diabetes Study incorporates a 1995/1996 birth cohort recruited in 2000/2001 when the children were 5 years old (307 children [170 boys]) (9). Most of the children were white Caucasian, and five children were of mixed race, reflecting the racial mix of the city of Plymouth. According to the pediatric thresholds for overweight and obesity proposed by the International Obesity Task Force, 13% of the EarlyBird boys and 26% of the girls were overweight at baseline, which included 4% and 5%, respectively, who were obese. The thresholds approximate to the 91st and 98th centiles of the 1990 BMI reference curves for the U.K. and are deemed to correspond to equivalent thresholds in adulthood. At baseline, the number of children with family history of T2D were as follows: mother (N = 0), father (N = 2), maternal grandmother (N = 14), maternal grandfather (N = 19), paternal grandmother (N = 18), and paternal grandfather (N = 18). The collection of data from the EarlyBird cohort is composed of clinical and anthropometric variables measured on an annual basis from 5 to 16 years of age.

BMI was derived from direct measurement of height (Leicester Height Measure; Child Growth Foundation, London, U.K.) and weight (Solar 1632 electronic scales; Tanita Europe BV), performed in duplicate and averaged. BMI SD scores were calculated from the British 1990 standards (10).

The children were fasted overnight for 10 h before venesection. HOMA2-IR and HOMA2-B were determined each year from fasting glucose (Cobas Integra 700 analyzer; Roche Diagnostics) and insulin (DPC Immulite) (cross-reactivity with proinsulin, 1%) using the HOMA program, which has been validated in children (11). Peripheral whole blood and serum were collected annually and stored at −80°C for analysis. Dehydroepiandrosterone sulfate (DHEAS), androstenedione, and testosterone were measured in serum by liquid chromatography–tandem mass spectrometry using the Waters Acquity UPLC system and Quattro Premier Tandem Quadrupole mass spectrometry (Waters Corporation, Milford, MA). Free testosterone (FT) was calculated using the formula of Vermeulen et al. (12). Serum sex hormone–binding globulins (SHBGs) were assayed using the Roche Cobas method on the E170 Modular Analytics system, and IGF-1 was measured by chemiluminescence immunoassay (Quest Diagnostics Nichols Institute Diagnostics, San Juan Capistrano, CA) using standards referenced in the WHO First International Reference Reagent 1988 (IGF‐1 87/518).

A total of 1,793 associations with T2D are reported in the GWAS Catalog (https://www.ebi.ac.uk/gwas/) (13). Results from small studies (n < 1,000 subjects) were excluded. Only SNPs with genome-wide significance for T2D itself, HOMA-IR, and response to metformin (but not T2D-associated comorbidities) and reported in European populations were selected. Excluding overlapping markers among those phenotypes, the GWAS Catalog yields 116 associations, 30 associations, and 1 association, respectively (147 SNPs). We were able to design genotyping assays for 136 of these SNPs that were included in the analysis. Descriptive information on selected SNPs is reported in Supplementary Table 1. Genomic DNA was extracted from blood using the QIAsymphony DSP DNA Midi Kit (96) on a QIAsymphony automation platform (Qiagen). The DNA concentration was measured with a fluorimetric method (PicoGreen; Thermo Fisher Scientific). Genotyping was performed using the SNP Type assay (Fluidigm Corporation), which relies on allele-specific PCR reactions that use three primers and two universal probes to distinguish between two alleles. Genotyping was performed using microfluidic chips requiring 2 chips per batch of 96 samples. Each batch was actually made of 86 samples of interest, 3 nontemplate controls with preamplification, 3 nontemplate controls without preamplification, and 2 control samples of known genotype run in duplicate (HapMap NA12891 and NA12892, obtained from the National Institute of General Medical Sciences Human Genetic Cell Repository at the Coriell Institute for Medical Research). DNA was normalized to 10 ng/μL and randomized, and 25 ng DNA was preamplified for 14 cycles with a pool of all 136 assay primer pairs using the Qiagen Multiplex PCR Kit, following Fluidigm Corporation’s recommendations. The amplification mixes were subsequently prepared from diluted preamplified products (1/100 in low Tris-EDTA buffer) and loaded together with the SNP Type assay mixes in 96.96 Dynamic Arrays IFC (Fluidigm Corporation). The final amplification and data collection were performed on a Biomark HD (Fluidigm Corporation). The analysis was performed with SNP Genotyping Analysis software (version 4.1.2; Fluidigm Corporation), which uses a cluster analysis method to automatically call genotypes (confidence threshold set at 65, SNP Type normalization, and K-means clustering method). Each SNP Type assay was manually quality controlled on each individual chip by evaluating the background fluorescence levels from nontemplate controls and the actual signal accuracy from the two positive controls. All genotypes were validated after manual inspection of each cluster, resulting in genetic data available for 121 SNPs. Quality control consisted in the exclusion of candidate SNPs with >5% missing values and subjects with >10% missing values. The sample consisted of 318 children with clinical and genetic information for 121 SNPs after this quality control. Pairwise linkage disequilibrium of SNPs was also calculated and reported in Supplementary Table 2. Comparison of minor allele frequencies using EarlyBird genetic data, as well as 1000 Genome British in England and Scotland (GBR) population data, was conducted and reported in Supplementary Table 3.

For insulin, glucose, and HOMA-IR and HOMA-B parameters collected over the 12-year period (subjects 5–16 years old), outlier detection filtering was based on interquartile range (IQR). Values below the lower bound (Q1 − Q4*IQR) or above the upper bound (Q3 + Q4*IQR) were set as missing for each period of age, in which Q1 and Q2 correspond to the value for the first and third quartile, respectively. Assuming subsequent repeated measurements analysis, subjects with missing data for at least 4 visits were excluded, resulting in the selection of 224 children.

Fasting insulin, HOMA-IR, HOMA-B, and fasting glucose were the tested outcome variables. Sex and BMI were considered as confounders. Fasting insulin and HOMA-IR distributions were skewed because of numerous values below the limit of detection. This proportion was highly variable, from 5% at 12 years of age to 39% at 5 years of age. Under such conditions, means and SEs become unreliable, whereas the median remains an appropriate measure of the distribution of the variable. Therefore, data analysis based on quantile regression (QR) was used. QR is a nonparametric method based on linear regression that makes no assumptions on the underlying distribution and is better suited to deal with data skewness (14). Assuming that the number of values below the limit of detection did not exceed 40% at any single visit, the median (i.e., 50th percentile) was considered as an analog to the mean of each trait. However, in genetic studies of metabolic traits, genetic variants may display quantile-specific effects (15–17). QR was thus performed on the 75th percentile for each trait to identify common genetic variants associated with high values of glycemic and insulin traits. Longitudinal association analysis testing the association between candidate SNPs and traits trajectory was conducted using mixed-effect QR as implemented in the lqmm R package (18). SNPs were encoded assuming an additive effect according to the risk allele definition available from PhenoScanner. In addition to the SNP main effect, the SNP × age interaction was also tested. For each trait, multiple testing was controlled by computing the false discovery rate (FDR). Because our analysis was based on candidate SNPs already identified for related glucose/insulin traits, we applied a relaxed 20% FDR significance cutoff. Additional genetic risk scores were computed and are reported in the Supplementary Data and Supplementary Tables 4–9. Additional longitudinal analyses were performed using the same methodology and criteria to test associations between SNPs and DHEAS, androstenedione, FT, 17-hydroxyprogesterone (17-OHP), SHBG, and IGF-1. Association analysis with 50th and 75th percentiles of the hormone distributions was performed only for SNPs significantly associated with 50th and 75th percentile insulin/glucose-related traits, respectively.

Data may be available upon request from J.P. and F.-P.M., subject in particular to ethical and privacy considerations.

The population characteristics are reported for the 12-year period in Fig. 1 and Supplementary Table 10. The glycemic and insulin traits followed the same trajectories in males and females. However, females had higher HOMA-B levels (P = 0.02 for a sex main effect) (Fig. 1). Fasting glucose increased from 5 years of age to reach a plateau between 13 and 15 years of age and then tended to decrease. Trajectories were very similar for insulin and HOMA indices, with decreases from 5 to 7 years of age and increases until 12 years in females and 14 years in males, after which a plateau was reached, followed by a fall.

To examine how genetic variants, puberty, and weight gain interact to influence insulin action and blood glucose levels, the primary analysis studied the genetic associations with the 50th percentile traits for fasting glucose, fasting insulin, HOMA-IR, and HOMA-B across the 12-year period (Table 1 and Supplementary Tables 11–14). Assuming a 20% FDR cutoff, eight SNPs were associated with one or more traits (Table 1), and two of them interacted with age (Supplementary Table 15).

The SNPs rs780094 and rs4457053, two intronic genetic variants of GCKR and ZBED3, respectively, showed a positive main effect on fasting glucose. The first variant showed an age-dependent interaction with the glucose trait (interaction test, P = 9.6e−4; adjusted, P = 0.11).

HOMA-B carried most of the replicated genetic signal, with six negatively associated common variants: rs11257655 located within DNase1-hypersensitive sites in the regulatory region of CDC12; rs12779790 within the CDC123–CAMK1D region; rs1111875 within the HHEX–IDE gene regions; rs7178572, an intronic SNP of HMG20A; rs9787485 within the promoter region of NRG3; and the intronic SNP rs1535500 of KCNK16. Three of these SNPs—rs9787485, rs11257655, and rs12779790—were also associated with fasting insulin, and the rs9787485 SNP was associated with HOMA-IR as well. For the rs9787485 SNP, homozygous carriers of the T allele tended to present a different time course for insulin and HOMA indices when compared with allele C carriers (Supplementary Fig. 1).

Because elevated fasting glucose and insulin are associated with impaired fasting glucose or T2D, a secondary analysis was performed using QR on the 75th percentile for each trait to identify common genetic variants associated with high values of glycemic and insulin traits (Table 2 and Supplementary Tables 16–19). All but the rs1535500 SNP replicated an association with one or more 50th percentiles for insulin or glucose traits (Table 2 and Supplementary Table 20). Altogether, association signal was extracted for 50 SNPs in or close to 36 genes. The analysis reveals an additional set of SNPs targeting 15 genes likely involved in regulation of high-level fasting glucose. For insulin-associated traits, no new SNPs were identified.

Puberty is a time during which rapid and dynamic changes occur in hormonal regulations due to an increase in sex hormones that induce transient changes in insulin sensitivity and insulin secretion. As a secondary exploratory analysis, we explored if insulin action may be influenced by genetic variations in the growth hormone–IGF-1 axis and adrenal and sex steroid activity. Genetic associations with 17-OHP, DHEAS, androstenedione, FT, SHBG, and IGF-1 were investigated over the 12-year period only for SNPs previously associated with insulin/glucose-related traits. The analysis was first conducted using the 50th percentile traits for each endocrine parameter (Supplementary Tables 21–26). As seen previously, the analysis was also performed using the 75th percentile endocrine traits (Supplementary Tables 27–32). Of note is the association of the genetic variant rs1111875 within the HHEX–IDE gene region with concentrations of SHBG (main effect test, P = 0.0059; FDR = 0.047), IGF-1 (age interaction test, P = 0.013; FDR = 0.106), and 17-OHP (main effect test, P = 0.0085 and FDR = 0.068; age interaction test, P = 0.0081 and FDR = 0.065), suggesting the genotype-related HOMA-B trajectory is associated with more profound physiological changes through pubertal development.

This novel work on the EarlyBird cohort has demonstrated that some genetic risk markers associated with increased risks of T2D also influence the trajectory of insulin and glucose in children independently of BMI. We report how longitudinal interaction analysis of SNPs and metabolic phenotypes of children have the potential to shed further light on the molecular disturbances associated with IR in childhood. Some of these genetic risk markers may be particularly relevant for the early identification of children at risk for prediabetes, independently of obesity.

Genome-wide association studies (GWAS) have identified SNPs that are associated with tissue-specific IR, β-cell dysfunction, or both, making individuals with these variants more prone to the adverse metabolic effects of obesity and T2D (19,20). Yet, due to the difficulty in recruiting large children’s cohorts, most GWAS have been undertaken in adults. Information about the genetics of prediabetes in children and adolescents remains scarce and mainly limited to cross-sectional studies on children who were already obese (21,22).

The present analysis of children from the EarlyBird cohort is a significant advance because it sheds light upon changes in insulin secretion, insulin action, and glycemia over time. We found that SNPs associated with adult diabetes susceptibility not only replicated their association with insulin and glycemic traits in these children independently of BMI but also showed age-specific interactions with some of these traits (50th percentile traits analysis), including rs780094 (GCKR), rs4457053 (ZBED3), rs11257655 (CDC123), rs12779790 (CDC123 and CAMK1D), rs1111875 (HHEX), rs7178572 (HMG20A), rs9787485 (NRG3), and rs1535500 (KCNK16). A greater number of common genetic variants were associated with variation of high levels of fasting glucose distribution (75th percentile). This observation is consistent with the understanding that elevated fasting plasma glucose is a risk factor for T2D and previous evidence that selected SNPs were associated with T2D in adult populations (23).

Two intronic genetic variants of GCKR and ZBED3 were associated with different trajectories of fasting glucose from an early age throughout childhood for both the median and higher (75th percentile) trait distributions. The SNP in ZEBD3 may contribute to the risk of T2D through elevated WNT activity (24). The WNT pathway involves pancreas β-cell genesis, glucagon-like peptide 1–mediated proliferation, and synthesis of glucagon-like peptide 1 (24). Moreover, the activity of this pathway has been shown to be modulated by short-chain fatty acids produced during digestion of dietary fibers and thus partially mediated their antidiabetogenic effects (24). Because the SNP in ZBED3 was shown previously to influence the effect of fiber intake on the incidence of T2D (24), it may be a relevant genetic marker for children who would benefit from optimized dietary fiber intake. GCKR regulates the activity of glucokinase in the liver and the pancreatic β-cells, with a pivotal role in glucose-stimulated insulin release and systemic glucose homeostasis (25). The intragenic SNP of GCKR contributes to the risk of T2D and dyslipidemia in different populations and may be associated with lower fasting blood glucose levels in adults (25). We found that the GCKR genotype was negatively associated with an age-dependent course of glucose, yet showed a positive association with blood glucose at 5 years of age. Glucose metabolism of β-cells also regulates the adaptive response to fuel loads (26), and it has been shown that children exhibit increased carbohydrate oxidation during pubertal development (27). Therefore, understanding of the role of β-cells in regulating fuel homeostasis during growth and development in healthy children is likely to be relevant for diabetes risk management (26).

Obese children and adolescents developing (pre)diabetes generally have a higher genetic predisposition related to gene variants modulating the early, dynamic phase of insulin secretion (20). In contrast with previous studies, we identified that relative defects in β-cell function distinguished children who developed prediabetes in the EarlyBird cohort, independently of BMI (6). Our analysis has found that among common variants associated with diabetes susceptibility, HOMA-B carried most of the replicated genetic associations, independently of BMI. Some of the replicated genetic variants, namely, SNPs in the CDC123, HHEX–IDE, and KCNK16 loci, showed a negative association with HOMA-B in the EarlyBird cohort. Our observations are in agreement with reports in various populations of their associations with reduced β-cell function and insulin secretion (20,28,29) and the proposition that SNPs in the CDC123 loci may represent a proxy for β-cell mass (28). Therefore, these genetic variants, previously identified as risk markers, also influence the normal course of β-cell function during growth and development of healthy children.

We explored whether genotypes also influence endocrine traits and therefore the maturation of other biological processes in childhood. We found some evidence for this in the interactions of the SNP in HHEX–IDE loci with the trajectories of SHBG, IGF-1, and 17-OHP. These findings support the biologically plausible idea that insulin action is also influenced by genetic variations in the growth hormone–IGF-1 axis and adrenal and sex steroid activity. Of particular note is the negative association with both HOMA-B and SHBG. SHBG is a glycoprotein that transports sex steroids in the circulation and regulates their access to target cells and for which low levels have been linked to diabetes and early puberty (30). Because the HHEX–IDE SNP relates to insulin secretion capacity in early life (31), our observations suggest effects beyond fetal development into pubertal growth and development.

This study has several limitations, including limited sample size and ethnic homogeneity of the EarlyBird cohort. However, the inclusion of healthy children who were of normal weight and insulin sensitivity reduces the variance of insulin and HOMA indices in the study population. A significant strength of the EarlyBird study is that it is a truly longitudinal study of a cohort of healthy children, with detailed annual phenotypic measurements from 5 to 16 years of age. Importantly, the majority of the children were of normal weight, and the IR observed was within the physiological range. Therefore, the EarlyBird study is not confounded by the effects of pre-existing obesity and IR; rather, it provides insights into the beginnings of these conditions.

In conclusion, we have demonstrated that SNPs previously associated with diabetes in adults also influence the course of glycemic and insulin traits during childhood independently of BMI. Our study has potential clinical applications, since β-cell dysfunction is an early event in the pathogenesis of diabetes. Further weight gain and demand for insulin will aggravate the progression from prediabetes to overt diabetes (32). Eight of these genetic risk markers may be particularly relevant for future studies aiming at early identification of children at risk for prediabetes, before puberty and prior to the development of obesity. Individuals at risk potentially could be offered an early lifestyle modification program to reduce the risks of progression toward diabetes. An important aim for future research will be to understand mechanisms of how genetic factors, puberty, and weight gain interact in the pathogenesis of prediabetes. This knowledge could allow the development of disease-prevention strategies in children, using personalized lifestyle and dietary interventions.

Acknowledgments. The authors acknowledge the life and work of our former colleague Terence Wilkin (1945–2017), Professor of Endocrinology and Metabolism, Peninsula Medical and Dental Schools at the University of Plymouth, whose vision and original thinking led to the creation of the EarlyBird Study and the establishment of the collaboration that made possible the genetic study reported in this article. The authors thank Ondine Walter (Nestlé Research, Lausanne, Switzerland) for biobanking, sample handling, and preparation at Nestlé and for support for compliance with the Human Research Act and Kei Sakamoto (Nestlé Research) and Ralf Heine (Nestlé Health Science, Lausanne, Switzerland) for critical review of the manuscript.

Funding. The EarlyBird study was supported by the Bright Future Trust, Kirby Laing Foundation, Peninsula Medical Foundation, Diabetes UK, the U.K. National Institute for Health Research, and the EarlyBird Diabetes Trust. The genotyping and genetic analyses, data interpretation, and writing of the manuscript were funded and performed in collaboration with Nestlé Research.

This report is independent research supported by the U.K. National Institute for Health Research Applied Research Collaboration South West Peninsula. The views expressed in this publication are those of the authors and not necessarily those of the National Institute for Health Research or the Department of Health and Social Care.

Duality of Interest. J.C., J.M., A.C., S.M., J.Ha., and F.-P.M. are or were employees of the Nestlé Group at the time the work was performed. J.Ho., J.P., and A.J. are employees of the Peninsula Medical and Dental Schools at the University of Plymouth. J.Ho. and A.J. have received funding from the Nestlé Group. No other potential conflicts of interest relevant to this article were reported.

Author Contributions. J.Ha. and F.-P.M. designed the study. J.M., A.C., S.M., and A.J. were involved in the acquisition of the data. J.C., J.Ho., J.P., J.M., J.Ha., and F.-P.M. were involved in the analysis and interpretation of the data. J.C., J.Ho., J.P., J.Ha., and F.-P.M. drafted the manuscript, and all co-authors approved the final version. J.P. is the guarantor of this work and, as such, had full access to all of the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.