Improved Performance of Dynamic Measures of Insulin Response Over Surrogate Indices to Identify Genetic Contributors of Type 2 Diabetes: The GUARDIAN Consortium

Type 2 diabetes (T2D) is a heterogeneous disorder in which complex interactions of peripheral insulin resistance with concomitant β-cell dysfunction lead to clinical presentation of disease. The “gold standards” for assessment of insulin resistance and β-cell dysfunction are the euglycemic-hyperinsulinemic and hyperglycemic clamps (1), respectively. An alternative approach, the frequently sampled intravenous glucose tolerance test (FSIGT) with minimal model (MINMOD) analysis (2), has been widely used and provides dynamic measures of glucose and insulin utilization, similar to the clamps, with correlated results across a range of glucose tolerance states (3–6). However, the expertise, time, and expense required by these measures as well as demands placed on the participant makes it difficult to perform these tests in large epidemiological studies. Consequently, basal estimates calculated from fasting glucose and insulin values, e.g., HOMA of insulin resistance (HOMA_IR) and of β-cell function (HOMA_B), in addition to simple measures of glycemic control have been widely used.

More than 80 loci (7) have been robustly implicated in T2D risk through evaluation of common variation across the genome in studies comprising up to 110,452 subjects (8). Collectively, however, these variants explain <10% of disease risk (8,9). Complementary efforts have explored the genetics of quantitative intermediate phenotypes of glucose homeostasis in normoglycemic individuals (8,10,11). To date, these studies have focused on the identification of variants modulating disease risk through assessment of basal insulin resistance or β-cell function. The majority of these loci appear to mediate their effects through β-cell function, while few loci have been identified that influence insulin resistance, despite an extensive literature documenting insulin resistance as a major component of T2D (12–16).

Despite the wide use of basal measures of glucose homeostasis to dissect the mechanistic heterogeneity of T2D, contributions to the pathophysiology remain unclear for many loci. Dynamic measures have the potential to elucidate contributors more proximal to the causal gene product resulting in the overt phenotype of T2D with attendant increases in power for discovery. Further, such analyses may more clearly identify the physiological path through which T2D susceptibility is transmitted. We evaluated the performance of basal (HOMA_IR and HOMA_B) and dynamic (acute insulin response [AIR_g] and insulin sensitivity [S_I]) measures of glucose homeostasis in the Genetics Underlying Diabetes in Hispanics (GUARDIAN) Consortium. Through statistical genetic comparison, we evaluated and contrasted the genetic basis of basal and dynamic measures of glucose homeostasis and used previously identified T2D susceptibility variants to evaluate the advantage of dynamic indices.

The GUARDIAN Consortium was established to evaluate the genetic basis of factors that predispose to T2D, including insulin resistance, metabolic clearance rate of insulin, and insulin response, in Mexican Americans (17). Participating cohorts were ascertained for various conditions including diabetes, gestational diabetes mellitus, or large family size and included persons with and without T2D who self-reported Mexican ancestry. Specific to this report, data were used from 2,560 Mexican American study subjects without T2D from four cohorts that measured glucose homeostasis by the FSIGT: the Insulin Resistance Atherosclerosis study (IRAS), the IRAS Family Study (IRASFS), BetaGene, and Troglitazone in the Prevention of Diabetes (TRIPOD). All participants provided written informed consent, and institutional review boards at the clinical, laboratory, and coordinating centers approved the study.

Dynamic measures of glucose homeostasis traits were measured in all participants by FSIGT with two modifications: an injection of insulin was used (TRIPOD injected tolbutamide) to ensure adequate plasma insulin levels for the accurate computation of S_I across a broad range of glucose tolerance (18), and a reduced sampling protocol was used (19). AIR_g was calculated as the increase in insulin concentrations at 2–8 min above the basal (fasting) insulin level after a bolus glucose injection at 0–1 min. S_I and glucose effectiveness (S_G) were derived from the FSIGT by mathematical modeling using the MINMOD program (20). Disposition index (DI) was calculated as the product of S_I × AIR_g. HOMA_IR and HOMA_B were modeled from fasting glucose and insulin measures using the updated HOMA model (21). A comprehensive description of study variables has previously been described (17,22).

Single nucleotide polymorphisms (SNPs) were selected for analysis with a bias toward variants for T2D and glucose homeostasis traits (e.g., fasting glucose), which have exhibited relatively large effect sizes and which have been widely replicated. This resulted in the selection of 57 variants (23–27) for analysis. Based on the a priori evidence of association, this discovery set yields increased power as well as increased probability of detecting effects across ancestries. Genotyping and quality control have been described in detail (22). Briefly, samples were genotyped on the Illumina HumanOmniExpress array. Samples with call rates >0.98 and SNPs with call rates >0.99 and minor allele frequency >0.001 passed laboratory quality control following usual best practices (e.g., sufficient signal and cluster separation with no replicate errors) (28). For family-based studies, pedigree structures were confirmed using standard procedures (e.g., Kinship-based INference for Gwas [KING] [http://people.virginia.edu/∼wc9c/KING/index.html]), and SNPs were examined for Mendelian inconsistencies using PedCheck (http://watson.hgen.pitt.edu/register/docs/pedcheck.html).

A variance components approach as implemented in Sequential Oligogenic Linkage Analysis Routines (SOLAR) (29) was used to compute estimates of heritability (h²) for each trait in the two family-based cohorts (BetaGene and IRASFS). Because BetaGene was ascertained for gestational diabetes mellitus, putting subjects at a higher risk to develop T2D, an ascertainment correction was implemented in SOLAR. When necessary, winsorization or transformation was applied to best approximate the distributional assumptions of conditional normality and homogeneity of variance. For traits warranting transformation, the same transformation was applied across both cohorts and included natural logarithm of the trait plus a constant (S_I), natural logarithm (fasting glucose, fasting insulin, HOMA_IR, and HOMA_B), and square root (AIR_g and DI); S_G was not transformed. Residual phenotypic variance, after accounting for covariates (age and sex ± BMI), was partitioned into additive genetic and nongenetic (environmental) components and tested using maximum likelihood methods in SOLAR.

Variance component models as implemented in SOLAR (29) were used to test for association in family cohorts and linear regression models as implemented in QSNPGWA from the SNPlash suite (https://github.com/guyrt/WFUBMC/tree/master/snplash) in nonfamily cohorts. All models included age, sex, BMI, study site (in multicenter recruitment studies), and admixture proportions. Traits were conditioned to approximate a normal distribution as described above. Population substructure was estimated using ADMIXTURE (http://www.genetics.ucla.edu/software/admixture/). In all tests for association, admixture proportions were included as covariates in the model such that the covariates were not collinear and tests of association did not exhibit evidence of inflation.

The primary inference was derived from the additive genetic model. The inverse variance–weighted method with weighting based on sample size was used to combine the evidence of association across cohorts as implemented in METAL (http://www.sph.umich.edu/csg/abecasis/metal/). A P < 8.77 × 10⁻⁴ (Bonferroni correction for 57 loci) was considered statistically significant. For each SNP-trait combination, we calculated the Wald statistic for comparison of the phenotypes on a unitless scale. With use of a matched pairs analysis as implemented in SAS (SAS Institute, Cary, NC), the Wilcoxon signed rank test was used to assess enrichment of previously reported loci for association with T2D-related quantitative traits.

Power for the association analysis (accounting for the familial correlations, with simulation-based estimations, resulting in an effective sample size of 92% of the total [n = 2,344]) was estimated to be 80% to detect SNP–quantitative trait associations that explain 1% and 0.56% of the variation in the quantitative traits at α = 5 × 10⁻⁸ and α = 1 × 10⁻⁴, respectively.

The study was performed with data from 2,560 Mexican Americans without T2D from four cohorts (Table 1). On average, the study subjects were overweight (BMI ≥25 kg/m²) and the majority of participants were female. The four cohorts varied in mean S_I from moderately insulin resistant (IRAS: mean S_I of 1.33 ± 1.24 × 10⁻⁴ ⋅ min⁻¹ ⋅ µU⁻¹ mL) to average S_I (IRASFS: mean S_I of 2.14 ± 1.86 × 10⁻⁴ ⋅ min⁻¹ ⋅ µU⁻¹ mL) to relatively insulin sensitive for the younger, largely female BetaGene cohort (mean S_I of 3.03 ± 1.63 × 10⁻⁴ ⋅ min⁻¹ ⋅ µU⁻¹ mL). Correspondingly, insulin response (AIR_g) was higher among more insulin resistant cohorts (IRAS and IRASFS), resulting in comparable DI values across the cohorts. The trend for the measure of insulin resistance derived from basal estimates (HOMA_IR) mirrored that of the FSIGT, and similarly, estimates of β-cell function (HOMA_B) were compensatory.

Genetic and environmental correlations (ρ_G and ρ_E, respectively) among the T2D-related quantitative traits are presented in Tables 2 and 3 (with BMI adjustment) and in Supplementary Table 1 (without BMI adjustment). AIR_g was positively correlated with the basal estimate of HOMA_B (ρ_G, IRASFS = 0.73; ρ_G, BetaGene = 0.28). However, directionality among the composite measures differed between IRASFS fasting glucose (fasting glucose ρ_G, IRASFS = −0.38; fasting insulin ρ_G, IRASFS = 0.68) and BetaGene fasting glucose (fasting glucose ρ_G, BetaGene = 0.07; fasting insulin ρ_G, BetaGene = 0.48). Among insulin resistance measures, S_I was negatively correlated with HOMA_IR (ρ_G, IRASFS = −0.73; ρ_G, BetaGene = −0.83) as well as the composite phenotypes of fasting glucose (ρ_G, IRASFS = −0.23; ρ_G, BetaGene = −0.23) and fasting insulin (ρ_G, IRASFS = −0.71; ρ_G, BetaGene = −0.84).

Heritability estimates from IRASFS and BetaGene for dynamic and basal measures of glucose homeostasis are presented in Tables 2 and 3 (with BMI adjustment) and in Supplementary Table 1 (without BMI adjustment). The dynamic measure of β-cell function, assessed by AIR_g, was the most consistent and highly heritable (h² = 0.47–0.56) measure assessed. In comparison, basal measures of β-cell function were lower in IRASFS (h² = 0.34) but comparable in BetaGene (h² = 0.55), which could be attributed to sample ascertainment differences between population-based and a family history of T2D, respectively. Measures of insulin resistance were also heritable (h² = 0.33–0.34). In contrast to HOMA_B, basal heritability estimates of insulin resistance were higher in BetaGene (h² = 0.48) but again comparable in IRASFS (h² = 0.31). Furthermore, an examination of the heritability estimates and associated SEs for the dynamic and basal measures revealed that S_I had the most consistent heritability estimates between the two studies, while the basal measures (HOMA_B, HOMA_IR, fasting glucose, and fasting insulin) had nonoverlapping point estimates within the SE-defined CIs.

Significant genetic association results observed among the 57 SNPs with dynamic (AIR_g, S_I, DI and S_G) and basal (HOMA_B, HOMA_IR, fasting glucose, and fasting insulin) glucose homeostasis traits are summarized in Table 4 (full results for SNP-trait combinations are presented in Supplementary Table 2). The most significant associations observed were among two modestly correlated SNPs, rs10830963 and rs1387153 (r² = 0.69), at the melatonin receptor 1B gene (MTNR1B) with AIR_g that reached genome-wide significance (P < 5.00 × 10⁻⁸). The MTNR1B locus was initially identified as a locus for fasting glucose and subsequently evaluated for association with T2D (30). Among the GUARDIAN cohorts, evidence of association between MTNR1B SNPs and other traits, such as fasting glucose, was comparatively modest (P > 1.50 × 10⁻⁶), and no evidence of association was observed with a basal measure of β-cell function (HOMA_B, P > 0.020).

The majority of significant associations (Bonferroni-corrected P < 8.77 × 10⁻⁴) among the 57 previously reported T2D-associated SNPs were with AIR_g (n = 10), with eight SNPs showing consistent association of the T2D-associated allele with decreased β-cell function. Further, 13 additional T2D variants were nominally associated with AIR_g (P < 0.05), with eight SNPs showing a consistent direction of effect. Comparison of the effect sizes for measures of β-cell function using the Wilcoxon signed rank test revealed a significant nonzero shift toward AIR_g (P = 7.0 × 10⁻⁴). The proportion of variants associated with AIR_g (17.5%) was more than expected by chance (P = 1.12 × 10⁻²⁰), consistent with the previous observation of enrichment of T2D-associated loci that contribute to β-cell function.

In addition to MTNR1B, there was a consistent pattern of variants showing association with AIR_g but little or no association with the basal measure of insulin response, i.e., HOMA_B, including KCNQ1 and TCF7L2, which is consistent with prior literature (31–33). Strikingly, only a single variant at the glucose-6-phosphatase catalytic subunit gene (G6PC2) (rs560887, P = 1.25 × 10⁻⁴) showed nominal evidence of association with basal estimates of β-cell function (HOMA_B).

Measures of insulin resistance, both dynamic (S_I) and basal (HOMA_IR), failed to show evidence of association among the 57 T2D-associated SNPs. More nominal evidence of association (P < 0.05) was observed among 11 SNPs with S_I (P = 0.48–0.0013) and 10 SNPs with HOMA_IR (P = 0.048–0.0023), only six of which overlapped between traits. Notably, variants in the potassium inwardly-rectifying channel gene (KCNJ11) were nominally associated with the dynamic measure of insulin resistance (S_I) (rs5219, P = 0.032). Effect size comparisons among insulin resistance loci using the Wilcoxon signed rank test were nonsignificant (P = 0.60).

Among additional phenotypes obtained from the FSIGT, three variants at two loci were significantly associated with S_G, which captures the ability of glucose to enhance its own disposal (34). SNP rs780094 (P = 5.38 × 10⁻⁶) is located in the glucokinase regulator gene (GCKR), and SNPs rs10830963 (P = 1.09 × 10⁻⁴) and rs1387153 (P = 6.85 × 10⁻⁴) are located near the MTNR1B locus.

Seven variants were significantly associated with DI, which is thought to be a good predictor of diabetes onset (35) taking into account the contributions of both insulin sensitivity and response, S_I and AIR_g, respectively. In each case the association with DI was one to two orders of magnitude more significant than that observed with AIR_g. Notably, two variants at the MTNR1B locus were more strongly associated with DI (rs10830963, P = 1.11 × 10⁻¹⁴ and rs1387153, P = 1.40 × 10⁻¹¹) than AIR_g despite the lack of contribution from S_I (P > 0.089). Two additional variants that were previously identified through association with fasting glucose, rs11708067 in the adenylate cyclase 5 gene (ADCY5) (36) and rs560887 in the glucose-6-phosphatase 2 gene (G6PC2) (37), were also associated with DI. Notably, neither SNP was significantly associated with fasting glucose (P = 0.18 and 0.033, respectively). Variants in KCNQ1 and insulin-like growth factor 2 mRNA–binding protein 2 gene (IGF2BP2) were also significantly associated with DI (P < 4.93 × 10⁻⁴).

Genome-wide association studies (GWAS) of T2D have identified >80 susceptibility loci to date (8,9,26,38–41). However, there is a diminishing return on investment, as contemporary studies require increasingly large sample sizes; e.g., the most recent analysis by Morris et al. (9) analyzed 34,840 case and 114,981 control subjects. In an effort to reduce the phenotypic heterogeneity, quantitative intermediate phenotypes of glucose homeostasis have been assessed. These studies have focused most frequently on basal measures of glucose homeostasis derived from easily obtained fasting measures. These genetic studies have confirmed the contribution of 53 loci, 33 of which impact T2D risk, but required similarly large sample sizes; e.g., Scott et al. (10) analyzed 133,010 individuals. Even when such large samples are used, much of the genetic component of T2D and its underlying glucometabolic phenotypes remain unknown. This challenge is amplified in populations such as African Americans and Mexican Americans where available samples sizes are appreciably smaller and genetic admixture complicates analysis. The current study demonstrates the potential value of further refining quantitative intermediate phenotypes of T2D through analysis of dynamic quantitative measures of insulin sensitivity and β-cell function.

Basal and dynamic measures of glucose homeostasis exhibit a differential genetic basis. Basal metabolic measures are derived from fasting measures of glucose and insulin via the HOMA approach, while dynamic phenotypes characterize an elicited response. For example, in this study AIR_g and S_I were measured in response to an intravenous glucose load using the minimal modeling approach. Extending our previous work (42) with an increased sample size and inclusion of contemporary genetic data, we observed that S_I was significantly correlated with HOMA_IR (ρ_G, IRASFS = −0.73; ρ_G, BetaGene = −0.83) with a stronger contribution from fasting insulin (ρ_G, IRASFS = −0.71; ρ_G, BetaGene = −0.84) compared with fasting glucose (ρ_G, IRASFS = −0.23; ρ_G, BetaGene = −0.23). In IRASFS, heritability of S_I (h² = 0.34) was modestly greater than that for HOMA_IR (h² = 0.31) with genetic background accounting for a greater proportion of S_I, while environmental factors made a stronger contribution to HOMA_IR. Notably, a direct assessment of the heritability and 95% CIs for basal and dynamic measures revealed that S_I was more similar between studies than HOMA_IR. Although not assessed herein, basal measures of insulin sensitivity have failed to adequately capture longitudinal change despite good correlation in the cross-sectional setting (43). As expected, the measures of β-cell function were positively correlated (AIR_g-HOMA_B ρ_G, IRASFS = 0.73 and BetaGene = 0.28) with a relatively lower correlation observed among the component fasting measures, particularly fasting glucose (AIR_g–fasting glucose ρ_G, IRASFS = −0.38 and BetaGene = 0.07), suggesting that the result is driven by the contribution of fasting insulin (AIR_g–fasting insulin ρ_G, IRASFS = 0.68 and BetaGene = 0.48).

Using previously identified T2D-susceptibility variants, the current study demonstrates the value of high-quality dynamic measures of glucose homeostasis. Most notable among these observations is the association of MTNR1B with AIR_g. This locus was first identified in a fasting glucose GWAS of 36,610 individuals of European descent (rs10830963, minor allele frequency =0.30, P = 3.2 × 10⁻⁵⁰) (30) and was only subsequently attributed to association with T2D in analyses testing up to 40,655 case and 87,022 control subjects (rs10830963, P = 8.0 × 10⁻¹³) (36). Comparatively, by targeting a precise measure of β-cell function, i.e., first-phase insulin response (AIR_g), we identified genome-wide significant association with this locus (rs10830963, P = 9.46 × 10⁻¹²) in a sample size of just 2,548 subjects, despite more nominal associations with fasting glucose (P = 1.50 × 10⁻⁶) (22). This association is consistent with the reported biology, i.e., the colocalization of MTNR1B with insulin in human islets (44). This pattern of superior performance of AIR_g was repeated for additional well-established T2D loci KCNQ1 and TCF7L2, which attained study-wise levels of statistical significance with much more nominal evidence of association with basal measures. Impairment of insulin response is believed to be the mechanism of action for both KCNQ1 (32,33) and TCF7L2 (31), although these results suggest a direct role in first-phase insulin response as opposed to significant contributions from changes in incretin secretion, as has been suggested for KCNQ1 (32). To compliment these dynamic measures of glucose homeostasis, we also evaluated DI, which is the product of S_I × AIR_g. Notably, this measure outperformed component phenotypes with comparatively more significant associations observed at the IGF2BP2 locus (AIR_g P < 0.023 vs. DI P < 4.93 × 10⁻⁴) and may indicate a more direct involvement in physiological cross talk mechanisms used to maintain glucose homeostasis.

Huyghe et al. (45) have genetically assessed a battery of quantitative intermediate phenotypes characterizing insulin processing, secretion, and glycemic traits in the Metabolic Syndrome in Men (METSIM) study. While the approach herein used metabolic phenotypes derived from the FSIGT, clinical testing in the METSIM study used the oral glucose tolerance test (OGTT). This study identified associations with fasting proinsulin levels at previously reported GWAS loci as well as novel genes associated with fasting proinsulin and the insulinogenic index. It is noteworthy that association with MTNR1B was not reported. In contrast, Prokopenko et al. (46) identified significant association at MTNR1B with decreased insulin secretion (corrected insulin response [CIR], P = 6.71 × 10⁻²⁸) obtained from the OGTT in a comparable sample size from the Meta-Analysis of Glucose- and Insulin-related traits Consortium (MAGIC). The OGTT-derived measure of insulin secretion represents stimulated response to oral glucose administration and may highlight additional component pathways toward development of the overt phenotype of T2D with variable contribution by MTNR1B.

Among additional novel dynamic phenotypes obtained from the FSIGT, the ability of glucose to enhance its own disposal is captured in the form of S_G. Among significant results, SNP rs780094 located in the glucokinase regulator gene (GCKR) was associated with S_G (P = 5.38 × 10⁻⁶). This locus is supported biologically by glucokinase, which catalyzes the ATP-dependent phosphorylation of glucose, the first and rate-limiting step in liver glucose metabolism (47). Extending upon the current literature, we observed association of MTNR1B rs10830963 with S_G (P = 1.09 × 10⁻⁴), which warrants additional follow-up studies for a role in glucose tolerance, as has been previously suggested (48). However, among the candidates evaluated, G6PC2 has been implicated in the alteration of hepatic glucose production (49,50). Although no association was observed with S_G (P = 0.15) that would represent the most proximal phenotype in GUARDIAN, further work is needed to accurately measure this metabolic pathway.

The observations described here suggest that basal and dynamic measures, the latter resulting from either oral (OGTT) or intravenous (FSIGT) stimulation, provide different estimates with differing levels of sensitivity to discrete elements of glucose homeostasis. Thus, T2D risk polymorphisms may selectively be associated with distinct measures. While these measures are correlated, the association results with AIR_g are striking, with rs10830963 having association P values 4 orders of magnitude stronger than fasting glucose and 10 orders of magnitude stronger than HOMA_B. Multiple other T2D variants (e.g., KCNQ1, TCF7L2) showed similar if less dramatic differences. These results suggest that there is discrete involvement of these genes in first-phase insulin response. In a similar vein, Huyghe et al. (45) identified a variant associated with insulin processing (C-peptide), a phenotype that is not available in the GUARDIAN cohorts. The utility of this type of approach for understanding the genetic contribution to disease progression is the proximity of the phenotypes to the underlying genetic variation. Thus, increased power is observed by reducing phenotypic heterogeneity; e.g., insulin resistance precedes development of impaired insulin response, yet few studies of T2D as a qualitative trait assess insulin resistance among control subjects. Together these studies are consistent with multiple metabolic contributions to T2D that are revealed and available for investigation only when detailed physiological phenotyping has been performed.

Further studies of precise metabolic phenotyping are needed to identify informative intermediate phenotypes of glucose homeostasis, and complimentary studies built upon this knowledge are needed, particularly in ethnic minority populations who are disproportionately burdened by T2D. Although more costly to attain, these measures, when analyzed in a comparatively small study population, yielded significant results that could be considered physiologically “closer” to the causal pathway. More broadly, the results presented here argue for detailed metabolic phenotyping in the further search for diabetogenic loci and as a way to gain insight into the discrete mechanisms of action.

Acknowledgments. The authors thank the investigators, staff, and participants of the studies for their valuable contributions.

Funding. This research was supported by the GUARDIAN study through funding provided by the National Institute of Diabetes and Digestive and Kidney Diseases (NIDDK) (DK-085175) and from National Heart, Lung, and Blood Institute grants HL047887, HL047889, HL047890, and HL047902 for IRAS and HL060944, HL061019, and HL060919 for IRASFS; NIDDK grant DK-061628; and the American Diabetes Association Distinguished Clinical Scientist Award for BetaGene. The provision of genotyping data was supported in part by National Center for Advancing Translational Science grant UL1-TR-000124 (Clinical and Translational Science Institute) and NIDDK grant DK-063491 (Diabetes Research Center). Computing resources were provided in part by the Wake Forest School of Medicine Center for Public Health Genomics.

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

Author Contributions. N.D.P. and D.W.B. contributed to the writing of the manuscript. N.D.P., L.E.W., J.I.R., R.M.W., and D.W.B. contributed to the study design, management, and coordination of the project. C.D.L., N.W., and A.H.X. contributed to the initial analysis in GUARDIAN. C.D.L. contributed to meta-analysis of data. All authors gave final approval of the manuscript. T.A.B., H.A., R.N.B., L.J.R., Y.-D.I.C., M.O.G., and T.F. contributed to the phenotyping. T.H. and K.D.T. contributed to the genotyping in GUARDIAN. N.D.P. and D.W.B. 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.

Prior Presentation. Parts of this study were presented in abstract form at the 73rd Scientific Sessions of the American Diabetes Association, Chicago, IL, 21–25 June 2013.