Genetic Variants Associated With Quantitative Glucose Homeostasis Traits Translate to Type 2 Diabetes in Mexican Americans: The GUARDIAN (Genetics Underlying Diabetes in Hispanics) Consortium.

Insulin sensitivity, insulin secretion, insulin clearance, and glucose effectiveness exhibit strong genetic components, although few studies have examined their genetic architecture or influence on type 2 diabetes (T2D) risk. We hypothesized that loci affecting variation in these quantitative traits influence T2D. We completed a multicohort genome-wide association study to search for loci influencing T2D-related quantitative traits in 4,176 Mexican Americans. Quantitative traits were measured by the frequently sampled intravenous glucose tolerance test (four cohorts) or euglycemic clamp (three cohorts), and random-effects models were used to test the association between loci and quantitative traits, adjusting for age, sex, and admixture proportions (Discovery). Analysis revealed a significant (P < 5.00 × 10(-8)) association at 11q14.3 (MTNR1B) with acute insulin response. Loci with P < 0.0001 among the quantitative traits were examined for translation to T2D risk in 6,463 T2D case and 9,232 control subjects of Mexican ancestry (Translation). Nonparametric meta-analysis of the Discovery and Translation cohorts identified significant associations at 6p24 (SLC35B3/TFAP2A) with glucose effectiveness/T2D, 11p15 (KCNQ1) with disposition index/T2D, and 6p22 (CDKAL1) and 11q14 (MTNR1B) with acute insulin response/T2D. These results suggest that T2D and insulin secretion and sensitivity have both shared and distinct genetic factors, potentially delineating genomic components of these quantitative traits that drive the risk for T2D.

Introduction

The pathophysiologic basis of type 2 diabetes (T2D) reflects derangements in both insulin sensitivity and β-cell function (1). Alterations in insulin clearance and glucose effectiveness may also contribute to the development of T2D (2). Genome-wide association studies (GWAS) of T2D have focused almost entirely on clinical presentation of disease and not on these underlying pathophysiologic traits. Expanding the focus to include the genetic basis of insulin sensitivity and β-cell function could expand our knowledge of the pathophysiologic pathways underlying T2D. To date, GWAS of T2D and related traits have been conducted primarily in populations of European origin (3). However, the prevalence of T2D and related traits varies by ethnicity, suggesting that differential genetic architecture will provide important insight into T2D diathesis.

GWAS in case/control samples of T2D have had a substantial impact on the current understanding of genetic susceptibility to disease, implicating variants in at least 70 genes/regions, each of which has relatively small individual effects but is common in the general population (4). Most identified T2D genes appear to mediate their influence through the β-cell and not through insulin resistance. These data contrast with other evidence and the widely accepted belief that insulin resistance is a major (5,6) heritable (7–10) component of T2D susceptibility. This suggests that insulin resistance is a part of the necessary milieu but is insufficient to cause frank T2D in isolation.

GWAS of the underlying pathophysiologic traits of insulin sensitivity and β-cell function have relied almost entirely on surrogate measures, such as HOMA parameters (11). Although these fasting measures do not reflect the dynamic processes of glucose homeostasis, new T2D loci have been identified through GWAS of basic T2D-related traits, such as fasting glucose (11,12). We recently documented substantial heritability of direct measures of insulin resistance and insulin clearance in Mexican Americans (13), suggesting that genetic factors underlying these traits should be investigated to identify new loci underlying disease susceptibility. In addition, strong genetic correlation was observed between these traits, raising the possibility of shared genetic determinants (13).

Only two GWAS of T2D (14,15) have been conducted in Mexican-origin populations, whose disease risk is nearly two times greater than that of European-origin populations (16). The recent SIGMA (Slim Initiative in Genomic Medicine for the Americas) T2D Consortium identified a novel risk variant in SLC16A11, which is rare in European and African individuals, suggesting a possible role for triacylglycerol metabolism in T2D (15). Thus, the study of detailed physiologic traits in individuals of Mexican ancestry could uniquely expand our understanding of T2D.

The Genetics Underlying Diabetes in Hispanics (GUARDIAN) Consortium was designed to overcome numerous gaps in the field of T2D genetics. GUARDIAN conducted a GWAS in multiple Mexican ancestry cohorts with highly detailed glucose homeostasis measures. In the Discovery phase, measures were obtained through gold-standard protocols (i.e., euglycemic clamp, frequently sampled intravenous glucose tolerance test [FSIGT]). Genomic regions associated with these quantitative traits were carried forward in a Translation phase that evaluated association with the clinical outcome of T2D. Using this approach, GUARDIAN has found novel and known risk variants in a Mexican ancestry population that specifically translate to T2D. The study identified possible new pathways of disease etiology and discovered risk variants for glucose homeostasis traits that do not associate with overt T2D, providing unique opportunities to understand physiologic regulation of glucose homeostasis traits within the normal range.

Research Design and Methods

Study Populations

Discovery Cohorts

Seven cohorts were included in the Discovery phase: five family-based studies [Insulin Resistance Atherosclerosis Family Study (IRAS-FS) (17), BetaGene (18), Hypertension-Insulin Resistance Family (HTN-IR) study (10), Mexican-American Coronary Artery Disease (MACAD) study (19), and NIDDM-Atherosclerosis Study (NIDDM-Athero) (20)] (n = 3,925) and two non–family-based studies [IRAS (1) and Troglitazone in the Prevention of Diabetes (TRIPOD) study (21)] (n = 411). Cohorts were ascertained based on various conditions, including diabetes, gestational diabetes mellitus, hypertension, and atherosclerosis (Supplementary Data). Cohorts included persons without T2D who self-reported Mexican ancestry. Four studies measured glucose homeostasis traits by FSIGT (22) (IRAS, IRAS-FS, BetaGene, and TRIPOD) and three by euglycemic clamp (23) (MACAD, HTN-IR, and NIDDM-Athero). The primary traits of interest were insulin sensitivity (SI from FSIGT or glucose infusion rate [M] from clamp as well as a meta-analysis combining these, denoted as SI + M), metabolic clearance rate of insulin (MCRI), acute insulin response (AIRg), disposition index (DI), and glucose effectiveness (SG). All participants provided written informed consent, and institutional review boards at the clinical, laboratory, and coordinating centers approved the study.

Phenotyping

Glucose homeostasis traits were measured by hyperinsulinemic-euglycemic clamp in three studies using an identical protocol (23). A priming dose of human insulin (Novolin; Novo Nordisk, Clayton, NC) was given followed by infusion for 120 min at a constant rate (60 mU ⋅ m−2 ⋅ min−1) to establish steady-state hyperinsulinemia. Blood was sampled every 5 min, and the rate of 20% dextrose coinfused was adjusted to maintain plasma glucose concentrations at 95–100 mg/dL. M over the last 30 min of steady-state insulin and glucose concentrations reflects glucose uptake by all tissues of the body (primarily insulin-mediated glucose uptake in muscle) and is therefore directly correlated with tissue insulin sensitivity (23). The insulin sensitivity index was calculated as M/I, where I is the steady-state insulin level. To distinguish between insulin sensitivity and clearance in this study, we relied on M as an approximation for insulin sensitivity because the calculations of M/I and insulin clearance both use steady-state insulin in the denominator. MCRI was calculated as the insulin infusion rate divided by the steady-state plasma insulin level of the euglycemic clamp (9,23). DI, a measure of β-cell compensation for insulin resistance, was calculated as M/I × Δ insulin, where Δ insulin was calculated as the difference between insulin at 30 min and insulin at baseline from a 2-h oral glucose tolerance test.

Glucose homeostasis traits were measured by FSIGT in four studies, with two modifications. An injection of insulin was used (one study, TRIPOD, injected tolbutamide) to ensure adequate plasma insulin levels for computation of insulin resistance across a broad range of glucose tolerance (24). Additionally, the reduced sampling protocol [which requires 12 rather than 30 plasma samples (25)] was used to facilitate study of large numbers of individuals. A 50% glucose solution (0.3 g/kg) and regular human insulin (0.03 units/kg) were injected through an intravenous line at 0 and 20 min, respectively. Blood was collected at −5, 2, 4, 8, 19, 22, 30, 40, 50, 70, 100, and 180 min for plasma glucose and insulin concentrations. SI and SG were calculated by mathematical modeling using the MINMOD program (version 3.0 [1994]) (22). AIRg was calculated as the increase in insulin concentrations at 2–8 min above the basal (fasting) insulin level after the bolus glucose injection at 0–1 min. DI was calculated as the product of SI and AIRg. MCRI was calculated as the ratio of the insulin dose over the incremental area under the curve of insulin from 20 min to infinity (26) (Eq. 1) as follows:where Dose is the amount of insulin injected at 20 min. Ins(t) is the plasma insulin concentration in standard units (μU/mL) at each FSIGT sampling point, and Ins(0) is the fasting plasma insulin concentration determined before the FSIGT glucose injection.

Genotyping

All samples were genotyped on the Illumina HumanOmniExpress BeadChip, and alleles were called using GenomeStudio software (Illumina, San Diego, CA) (27,28). Samples with call rates >0.98 and single nucleotide polymorphisms (SNPs) with call rates >0.99 and minor allele frequency (MAF) >0.001 passed laboratory quality control by usual best practices (e.g., sufficient signal and cluster separation with no replicate errors) (29). Additionally, ∼22,000 SNPs were manually reviewed for clustering accuracy.

Statistical Analysis

Quality Control

Samples were removed from analysis if the overall call rate was <0.98, self-reported ethnicity was inconsistent with genetic data (i.e., admixture proportions) relative to other members of the cohort (i.e., a genetic outlier), self-reported sex was inconsistent with genotype data, the sample exhibited excess or insufficient heterozygosity relative to cohort expectations, or the genotype data were inconsistent with the genotype data from existing SNP data (i.e., fingerprinting). The primary inferential SNPs did not exhibit differential missingness by trait, had a SNP call rate >95%, and were consistent with Hardy-Weinberg expectation proportions. For family-based studies, pedigree structures were confirmed using standard procedures (e.g., KING [Kinship-Based Inference for GWAS], http://people.virginia.edu/∼wc9c/KING). Each SNP was examined for Mendelian inconsistencies using PedCheck (Program for Detecting Marker Typing Incompatibilities in Pedigree Data, http://watson.hgen.pitt.edu/register/docs/pedcheck.html), and inconsistencies were converted to missing. A maximum of 693,128 SNPs were meta-analyzed among the Discovery cohorts.

Population Stratification

Population substructure was estimated using ADMIXTURE version 1.21 (http://www.genetics.ucla.edu/software/admixture) at each study site based on SNPs that passed quality control (n = 117,347 linkage disequilibrium [LD]–pruned SNPs). Data from the HapMap Project (CEU [northern and western European ancestry], CHB/JPT [Han Chinese in Beijing, China/Japanese in Tokyo], YRI [Yoruba in Ibadan, Nigeria], and MEX [Mexican ancestry in Los Angeles, CA]; n = 591) were used as reference populations. Depending on the cohort, up to k = 5 subpopulations were identified based on low cross-validation error. In all tests for association, admixture proportions were included as covariates in the linear or variance component models such that the covariates were not collinear and tests of association did not exhibit evidence of inflation.

Association

Variance component models as implemented in the GWAF (Genome-Wide Association analyses with Family) (30) or SOLAR (Sequential Oligogenic Linkage Analysis Routines) (31) programs were used to test for association in family cohorts and linear regression models as implemented in QSNPGWA (http://github.com/guyrt/WFUBMC) in nonfamily cohorts. All models included age, sex, study site (in multicenter recruitment studies), and admixture proportions. Conditional analyses were performed for significant loci with multiple uncorrelated variants by including the most significant variant as an additional covariate. If necessary, winsorization or transformation was applied to best approximate the distributional assumptions of conditional normality (conditional on the covariates) and homogeneity of variance. For traits warranting transformation, the same transformation was calculated across all cohorts and included the natural logarithm of the trait plus a constant (SI), natural logarithm (MCRI derived from FSIGT), and square root (M, AIRg, and DI); MCRI derived from clamp and SG were not transformed. The primary inference was derived from the additive genetic model. However, we also tested for a lack of fit to additivity using the orthogonal contrast. If the lack of fit was significant (P < 0.05), we reported the “best” P value as the minimum of the three genetic models. It can be shown that this approach has an inflation factor of 1.3. For robust estimation purposes, the additive and recessive genetic models were not calculated if there were not at least 10 and 20 individuals homozygous for the minor allele, respectively. In addition to single-variant association tests, a genetic risk score was calculated; that is, risk allele load was determined by the number of previously reported T2D risk alleles (Supplementary Table 2) carried by each individual and analyzed for association with the primary traits of interest (SI + M, MCRI, and AIRg). Subsequently, an enrichment analysis was performed among these variants to determine whether an excess of nominally significant values was observed.

A nonparametric meta-analysis was calculated to combine the evidence of association across cohorts as implemented in METAL (http://www.sph.umich.edu/csg/abecasis/metal). For each genetic model and T2D-related quantitative trait, a weighted, fixed-effects meta-analysis was calculated, weighting by cohort sample size and not by the SE of the parameter estimate because the traits were not identical and studies had different designs and ascertainment criteria. Power for the association analysis in the Discovery cohorts accounting for the familial correlations, with stimulation-based estimations resulting in an effective sample size of 92%, was estimated to be 80% to detect SNP-quantitative trait associations that explain 1% and 0.56% of the variance at α = 5 × 10−8 and α = 1 × 10−4, respectively.

Translation

Evaluation of T2D-related quantitative traits is a potentially powerful approach to identify genetic variants contributing to defects in specific underlying pathways leading to T2D; however, the true impact can be gauged only through direct validation in a population with clinically defined disease.

Translation Cohorts

Six cohorts were included in the translation effort: Los Angeles Latino Eye Study (LALES), Multi-Ethnic Study of Atherosclerosis (MESA) Family, MESA, Starr County Health Studies, Women’s Health Initiative (WHI), and SIGMA (Supplementary Data). All cohorts were of self-reported Mexican origin and provided either look-ups of the index SNPs or a preselected proxy.

SNP Selection for Translation

Results from the Discovery GWAS were reviewed to generate the list of SNPs to be examined for translation to T2D. For each of seven traits (SI + M, SI, M, MCRI, AIRg, DI, and SG), SNPs associated in the GWAS (P < 1.00 × 10−4) or with both primary traits (SI + M and MCRI, P < 0.05) were included. After removal of correlated SNPs (r2 > 0.90 within a 500-kb window in the largest set of unrelated samples; n = 553 from BetaGene and TRIPOD), this yielded a total of 594 SNPs for translation to T2D.

Meta-analysis of Discovery and Translation Results

Discovery and Translation cohort genotype data were aligned with the positive strand for compatibility. After alignment, the same nonparametric meta-analysis approach was used to combine the association statistics. Lower values for the quantitative traits were hypothesized to be associated with T2D risk (2,32). Supplementary Fig. 1 estimates the power of the Translation cohort to detect various odds ratios for T2D over a range of MAFs.

Functional Database Validation

Queries of the Encyclopedia of DNA Elements (ENCODE) data were carried out using both the University of California, Santa Cruz (UCSC), genome browser (http://genome.ucsc.edu) and RegulomeDB (http://regulome.stanford.edu). The positions of associated loci were overlaid with DNase I hypersensitivity hot spots from ENCODE that identified regions of chromatin accessibility and transcription factor motifs in 125 diverse cell lines and tissues. We used the browsers set up by the Genotype-Tissue Expression (GTEx) project to determine whether any of our association signals represented expression quantitative trait loci (eQTL) (i.e., SNPs associated with mRNA transcript levels) (33).

Results

The Discovery sample included 4,176 Mexican Americans without T2D (Table 1 and Supplementary Table 1). Characteristics of the sample have been previously reviewed (13). The Translation sample comprised 6,463 T2D case and 9,232 control subjects (Supplementary Data).

Table 1

Clinical characteristics of the Discovery cohorts

	FSIGT cohorts				Clamp cohorts
	BetaGene	TRIPOD	IRAS	IRAS-FS	HTN-IR	MACAD	NIDDM-Athero
Sample size	1,202	125	187	1,034	694	752	182
Age (years)	34.6 ± 7.9	34.8 ± 6.3	58.8 ± 8.3	40.6 ± 13.7	37.4 ± 14.2	34.5 ± 8.8	31.8 ± 9.69
Women (%)	72.1	100.0	58.3	59.0	59.4	56.7	58.2
BMI (kg/m2)	29.5 ± 6.1	30.6 ± 5.4	28.9 ± 5.1	28.3 ± 5.7	28.8 ± 5.5	28.9 ± 5.1	28.6 ± 6.3
AIRg (μU ⋅ mL−1 min)	569 ± 480	488 ± 450	673 ± 702	760 ± 649	NA	NA	NA
SG (min−1)	0.0178 ± 0.0067	0.0157 ± 0.0041	0.0208 ± 0.0088	0.0202 ± 0.0091	NA	NA	NA
MCRI (L/min)	10.1 ± 5.7	NA*	4.2 ± 2.0	5.5 ± 2.4	NA	NA	NA
MCRI (mL ⋅ m−2 ⋅ min−1)	NA	NA	NA	NA	458.2 ± 111.8	471.8 ± 116.3	416.2 ± 140.3
SI (×10−4 min−1 ⋅ μU−1 ⋅ mL−1)	3.03 ± 1.63	2.57 ± 1.79	1.33 ± 1.24	2.14 ± 1.86	NA	NA	NA
M (μmol ⋅ m−2 ⋅ min−1)	NA	NA	NA	NA	1,273 ± 547	1,364 ± 646	1,255 ± 533
DI	1,409 ± 946	1,004 ± 724	1,245 ± 1,184	1,202 ± 1,236	NA	NA	NA
DI (μmol ⋅ m−2 ⋅ min−1)	NA	NA	NA	NA	NA†	136.8 ± 100.2	93.23 ± 54.4

Data are mean ± SD unless otherwise indicated. NA, not available.

*MCRI is not available for TRIPOD because of the use of tolbutamide in the FSIGT.

†DI is not available for HTN-IR because of the lack of 30-min insulin values from oral glucose tolerance testing.

Figure 1 displays associations with T2D-related quantitative traits in the Discovery cohorts with signals that were significant at P < 2.00 × 10−6 listed in Table 2. (Supplementary Table 2 lists nominally significant hits.) Results were broadly similar with the inclusion of BMI as a covariate (Supplementary Table 3). The top signal (P = 5.23 × 10−12) was the association of rs10830963 in MTNR1B (melatonin receptor 1B gene) with AIRg; this SNP was also associated with DI but not with SI (Fig. 2). Associations with insulin sensitivity (SI, M, or SI + M), MCRI, and SG did not reach genome-wide significance levels. One signal for M (rs11683087) was located near IRS1, a locus previously identified for T2D and deemed to act through insulin resistance based on association with HOMA of insulin resistance (34). These variants were not highly correlated (r2 = 0.04), and the previously described variant (rs2943641) failed to show evidence of association with M (P = 0.63) or reduce the level of significance at rs11683087 upon conditional analysis (P = 1.29 × 10−6) (Supplementary Table 4).

Figure 1

Genome-wide Manhattan plots for the GUARDIAN Discovery meta-analysis. A: MCRI. B: Insulin sensitivity (SI + M). C: SI. D: M. E: AIRg. F: DI. G: SG.

Table 2

Top Discovery hits from the GUARDIAN Consortium, ordered by trait

SNP	Chr	Position*	Gene	Alleles†	RAF	Trait	β	P value
rs2302063	19	3150418	GNA15	A/C	0.336	MCRI	−0.29	7.31E-08
rs1602084	4	128843480	MFSD8	G/A	0.041	SI + M	8.97	5.20E-07
rs896232	2	2732877	MYT1L/TSSC1	T/C	0.291	SI + M	−5.28	1.26E-06
rs6719442	2	2722295	MYT1L/TSSC1	A/G	0.184	SI + M	−5.03	1.53E-06
rs1978648	2	43371542	HAAO/ZFP36L2/THADA§	T/C	0.324	SI	0.20	5.31E-07
rs896598	15	74036629	C15orf59	A/G	0.116	SI	0.37	5.83E-07
rs4887140	15	74046663	C15orf59/TBC1D21	G/T	0.139	SI	0.31	6.91E-08
rs196701	6	80147187	HMGN3/LCA5	C/T	0.132	SI	−0.35	1.37E-06
rs10492494	13	74920186	KLF12/LINC00347	A/C	0.240	M	−22.02	5.04E-07
rs11683087	2	227586606	LOC646736/IRS1‖	G/A	0.412	M	20.54	7.42E-07
rs10830963	11	92708710	MTNR1B¶	G/C	0.220	AIRg	−2.76	5.23E-12
rs1387153	11	92673828	FAT3/MTNR1B#	T/C	0.220	AIRg	−2.55	2.21E-09
rs2206734	6	20694884	CDKAL1**	T/C	0.198	AIRg	−2.05	1.02E-06
rs3847554	11	92668826	FAT3/MTNR1B††	A/G	0.341	AIRg	−1.64	1.08E-06
rs9368222	6	20686996	CDKAL1‡‡	A/C	0.264	AIRg	−1.46	1.28E-06
rs6803803	3	180116563	PEX5L/TTC14	C/T	0.003	AIRg	17.53	1.64E-06
rs10830963	11	92708710	MTNR1B¶	G/C	0.230	DI	−3.40	1.03E-11
rs1387153	11	92673828	FAT3/MTNR1B#	T/C	0.220	DI	−3.20	1.32E-09
rs2149423	13	36772381	CCDC169-SOHLH2; SOHLH2	G/A	0.315	DI	2.18	3.67E-07
rs3812570	9	139275204	SNAPC4	A/C	0.461	DI	−1.87	1.72E-06
rs523079	3	187615862	BCL6/LPP	T/C	0.069	SG	0.25	1.53E-07
rs780093	2	27742603	GCKR§§	T/C	0.341	SG	0.14	1.12E-06
rs788338	19	50778543	MYH14	C/T	0.287	SG	−0.17	1.66E-06

Independent signals (r2 < 0.80) with evidence of association (P < 2.00 × 10−6) with the nearest annotated RefSeq genes listed. Chr, chromosome; RAF, reference allele frequency.

*Build hg19.

†Reference allele/other allele.

§Previously identified T2D locus (THADA rs7578597, r2 = 0.0079).

‖Previously identified T2D locus (IRS1 rs2943641, r2 = 0.04).

¶Previously identified T2D locus (MTNR1B rs1387153, r2 = 0.69).

#Previously identified T2D locus (MTNR1B rs1387153).

**Previously identified T2D locus (CDKAL1 rs7754840, r2 = 0.42).

††Previously identified T2D locus (MTNR1B rs1387153, r2 = 0.54).

‡‡Previously identified T2D locus (CDKAL1 rs7754840, r2 = 0.72).

§§Previously identified T2D locus (GCKR rs780094, r2 = 0.98).

Figure 2

Regional plot of the MTNR1B locus in the GUARDIAN Discovery cohort meta-analysis. A: DI. B: AIRg. C: SI. Genotyped SNPs passing quality control measures across all Discovery cohorts are plotted with their Discovery meta-analysis P values (as −log10 values) as a function of genomic position (hg19). In each panel, the index variant is represented by a purple diamond. Color of additional variants indicates correlation with the index SNP (red, r2 ≥ 0.80; orange, 0.60 ≤ r2 < 0.80; green, 0.40 ≤ r2 < 0.60; light blue, 0.20 ≤ r2 < 0.40; dark blue, r2 < 0.20; gray, no r2 value available) based on pairwise r2 values from HapMap. Estimated recombination rates (taken from HapMap) are plotted to reflect the local LD structure. Gene annotations were taken from the UCSC genome browser.

Within the Discovery cohorts, we evaluated the association of previously reported T2D susceptibility variants (n = 90) (Supplementary Table 5) with seven T2D-related quantitative traits. Using the reported variant (n = 76) or a HapMap MEX proxy (n = 14; r2 > 0.80), the most profound effects were observed with decreased AIRg for 17 of the SNPs evaluated (P = 2.3 × 10−8–0.049). The most significant association was at the MTNR1B locus (rs1387153) (35). Comparatively, SI + M (n = 9; P = 0.0019–0.041) and MCRI (n = 5; P = 0.0053–0.050) had markedly fewer nominal associations. Similarly, the cumulative genetic risk score (P = 1.11 × 10−8) and enrichment analysis (P < 0.00001) were significantly associated with AIRg. Of note, we also observed an enrichment for previously reported T2D SNPs with insulin sensitivity (SI + M; P = 3.6 × 10−4), although the significance was attenuated in comparison.

Meta-analysis of the Discovery and Translation cohorts identified multiple SNPs that met or approached genome-wide significance (Table 3, Supplementary Fig. 2, and Supplementary Table 6) and included novel and established T2D loci. Results were broadly similar with the additional inclusion of BMI as a covariate (Supplementary Table 7). The most significant association observed was at rs2237897 (P = 1.24 × 10−21) in KCNQ1 (potassium voltage-gated channel, KQT-like subfamily, member 1 gene). This variant was associated with DI in the Discovery cohort (P = 7.04 × 10−6) and after conditional analysis for previously associated T2D variants (rs2237892, P = 1.57 × 10−4; rs231362, P = 1.06 × 10−5) (Supplementary Table 8). Three established T2D genes, motivated by their association with AIRg, remained associated in the meta-analysis: rs10830963 within MTNR1B (P = 5.86 × 10−9); rs2206734 within CDKAL1 (CDK5 regulatory subunit associated protein 1-like 1 gene) (P = 1.11 × 10−8); and rs7018745 near CDKN2A/B (cyclin-dependent kinase inhibitor 2A and 2B gene cluster) (P = 7.3 × 10−8), which is a strong genetic susceptibility locus for cardiovascular disease (36) and linked to T2D (37) (Fig. 3).

Table 3

Top regions from GWAS of T2D-related quantitative traits with translation to T2D in Mexican-origin cohorts, ordered by trait from the Discovery stage

				Discovery cohorts					Translation cohorts
Marker	Chr	Position*	Gene	Trait	RA	RAF	β	P value	OR (95% CI)	P value	Discovery and Translation meta-analysis P value
rs7219451	17	38957002	KRT28/KRT10	SI + M	C	0.380	−6.80	3.92E-06	1.08 (1.00–1.16)	1.06E-02	3.97E-07
rs6815953	4	183109012	TENM3	SI + M	T	0.431	4.10	3.14E-05	0.92 (0.86–0.99)	1.25E-02	2.47E-06
rs7581057	2	115958079	DPP10	SI + M	G	0.056	−5.47	1.74E-04	1.29 (1.10–1.50)	6.67E-03	4.82E-06
rs322394	5	172157768	NEURL1B/DUSP1	M	C	0.360	−15.83	3.38E-05	1.17 (1.09–1.26)	7.87E-04	1.12E-07
rs17060946	9	77808519	OSTF1/PCSK5	M	G	0.068	37.86	2.56E-06	0.87 (0.75–1.01)	6.89E-02	3.99E-06
rs13252932	8	25198091	DOCK5	M	T	0.049	47.71	5.17E-06	0.89 (0.75–1.04)	6.61E-02	6.11E-06
rs10830963	11	92708710	MTNR1B‖	AIRg	G	0.220	−2.76	5.23E-12	1.08 (1.00–1.18)	1.83E-01	5.86E-09
rs2206734	6	20694884	CDKAL1¶	AIRg	T	0.194	−2.05	1.02E-06	1.19 (1.09–1.30)	1.41E-03	1.11E-08
rs7018475	9	22137685	CDKN2B-AS1#/DMRTA1	AIRg	G	0.323	−1.64	4.90E-06	1.09 (1.02–1.17)	2.34E-03	7.31E-08
rs1387153	11	92673828	MTNR1B**	AIRg	T	0.220	−2.55	2.21E-09	1.07 (0.98–1.16)	1.97E-01	2.72E-07
rs2129969	11	45586169	PRDN11/CHST1	AIRg	G	0.477	1.14	2.10E-05	0.91 (0.85–0.98)	2.38E-02	4.10E-06
rs2053797	2	46370892	PRKCE	AIRg	A	0.162	1.89	5.02E-05	0.91 (0.83–1.00)	1.42E-02	4.21E-06
rs9553849	13	27082326	CDK8/WASF3	AIRg	A	0.199	1.57	6.70E-05	0.86 (0.79–0.94)	1.29E-02	4.70E-06
rs10870202	9	139257411	DNLZ	AIRg	T	0.430	−1.04	5.37E-05	1.12 (1.05–1.20)	2.08E-02	7.10E-06
rs10898909	11	72952496	P2RY2	AIRg	A	0.404	1.41	4.03E-06	0.94 (0.86–1.02)	9.27E-02	8.64E-06
rs2237897	11	2858546	KCNQ1††	DI	A	0.251	1.75	7.04E-06	0.73 (0.68–0.79)	1.89E-19	1.24E-21
rs4266763	9	139289825	SNAPC4	DI	A	0.489	−2.00	3.46E-06	1.10 (1.03–1.18)	1.23E-02	4.34E-07
rs2064197	6	8998811	SLC35B3/TFAP2A	SG	G	0.164	−0.19	7.23E-06	1.18 (1.08–1.30)	7.01E-04	2.56E-08
rs2291004	19	37997952	ZNF793	SG	A	0.153	0.16	3.82E-05	0.89 (0.82–0.97)	9.88E-03	2.18E-06
rs1260326	2	27730940	GCKR‡‡	SG	T	0.335	0.13	4.04E-06	0.96 (0.89–1.03)	5.58E-02	3.99E-06

SNPs with evidence of association with the nearest annotated RefSeq genes listed. Chr, chromosome; OR, odds ratio; RA, reference allele; RAF, reference allele frequency.

*Build hg19.

‖Previously identified T2D locus (MTNR1B rs1387153, r2 = 0.69).

¶Previously identified T2D locus (CDKAL1 rs7754840, r2 = 0.42).

#Previously identified T2D locus (CDKN2B rs7018475).

**Previously identified T2D locus (MTNR1B rs1387153).

††Previously identified T2D locus (KCNQ1 rs2237892, r2 = 0.90; KCNQ1 rs231362, r2 = 0.038).

‡‡Previously identified T2D locus (GCKR rs780094, r2 = 0.91).

Figure 3

Regional plots of loci attaining genome-wide significance (P < 5.00 × 10−8) in the combined Discovery and Translation meta-analysis. A: KCNQ1 rs2237897 with DI/T2D. B: MTNR1B rs10830963 with AIRg/T2D. C: CDKAL1 rs2206734 and AIRg/T2D. D: 6p24.3 rs2064197 and SG/T2D. Genotyped SNPs passing quality control measures across all Discovery cohorts are plotted with their Discovery meta-analysis P values (as −log10 values) as a function of genomic position (hg19). In each panel, the index variant from the Discovery cohort is represented by a purple circle, and the Discovery and Translation meta-analysis is represented by a purple square. Color of additional variants indicates correlation with the index SNP (red, r2 ≥ 0.80; orange, 0.60 ≤ r2 < 0.80; green, 0.40 ≤ r2 < 0.60; light blue, 0.20 ≤ r2 < 0.40; dark blue, r2 < 0.20; gray, no r2 value available) based on pairwise r2 values from HapMap. Estimated recombination rates (taken from HapMap) are plotted to reflect the local LD structure. Gene annotations were taken from the UCSC genome browser. chr, chromosome.

Four novel associations were observed that reached or approached genome-wide significance. At 6p24, rs2064197 was associated with SG, and the meta-analysis with T2D reached genome-wide significance (P = 2.56 × 10−8). Other novel associations included rs322394 (M/T2D, P = 1.12 × 10−7) at 5q35, rs7219451 at 17q21 (SI + M/T2D, P = 3.97 × 10−7), and rs4266763 (DI/T2D, P = 4.34 × 10−7) in SNAPC4 (small nuclear RNA activating complex, polypeptide 4 gene).

Discussion

GUARDIAN conducted a GWAS in seven Mexican American cohorts of insulin sensitivity, insulin secretion, insulin clearance, and glucose effectiveness directly quantified by the euglycemic clamp and FSIGT. We posited that the measurements of insulin sensitivity and clearance obtained by detailed physiologic phenotyping procedures are closer to the gene products and would yield increased statistical power to detect SNPs influencing trait variation. Establishing these loci in the Mexican American population will inform diabetes risk in an ethnicity that experiences a disproportionately high diabetes burden (16) and may explain risk in other ethnicities either directly or through a deeper understanding of the relevant biological pathways.

The most significant association observed (rs10830963, P = 5.23 × 10−12) (Table 2) that translated to T2D (P = 5.86 × 10−9) (Table 3) was at MTNR1B, which was initially identified as a locus for fasting glucose (35). Two modestly correlated variants in MTNR1B, rs10830963 and rs1387153 (r2 = 0.68), were associated with AIRg (P = 5.23 × 10−12 and 2.21 × 10−9, respectively). These variants were also, but less significantly, associated with fasting glucose in the Discovery cohorts (P = 3.92 × 10−8 and 2.09 × 10−5, respectively). As suggested by ENCODE, rs10830963 resides in an FOX2A transcription factor binding site and has a lower RegulomeDB score (3a vs. 5, respectively), which corroborates the stronger evidence of association observed at rs10830963. MTNR1B is expressed in both rodent and human islets and colocalizes with insulin. Gene expression increases with each copy of the rs10830963 risk allele in human islets from nondiabetic individuals. Consistent with this observation, MTNR1B gene expression levels are higher in human islets from patients with T2D than those from individuals without diabetes. MTNR1B is hypothesized to inhibit glucose-stimulated insulin secretion through binding of its ligand, melatonin, and decreasing cAMP levels (38), consistent with the direction of effect observed in the present study.

Among novel variants that translated to disease risk was rs2064197 (P = 2.56 × 10−8), which was also associated with SG in the Discovery cohorts (P = 7.23 × 10−6) and located intergenically on 6p24.3 between the SLC35B3 (solute carrier family 35, member B3) and TFAP2A (transcription factor AP-2 α). SG is the ability of glucose to enhance its own disappearance and suppress its production at fasting insulin levels (39,40). The role of SG in the regulation of glucose tolerance is often ignored but may be physiologically significant (40,41). SG varies by both physiologic and pathologic state (42) and has been shown to be predictive of conversion to T2D (43). Of note, this variant resides distally (1.7 Mb) to the recently implicated T2D susceptibility locus RREB1 (ras responsive element binding protein 1), which was not associated with SG in the present analysis (Supplementary Table 5). Other variants more nominally associated with this phenotype and translation to T2D included rs1260326, a missense variant located in GCKR (glucokinase regulator gene). This association is supported biologically because the ATP-dependent phosphorylation of glucose, which is catalyzed by glucokinase, is the first and rate-limiting step in liver glucose metabolism (44). Because this step of glucose metabolism is independent of dynamic insulin response, it is believed that a large portion of SG results from the ability of the liver to take up glucose through the glucokinase pathway, independent of insulin.

Replication of GWAS results in independent samples is widely accepted as critical. However, we are unaware of additional Mexican ancestry cohorts with highly detailed glucose homeostasis phenotypes available in which to directly replicate the present findings. Given that these phenotypes predict the subsequent occurrence of T2D, we have taken a unique approach by translating the findings to the directly relevant clinical phenotype T2D. This approach supports that these loci are involved in deterioration from impaired glucose homeostasis to T2D. Not surprisingly, we observed that only some of the quantitative trait loci identified through the Discovery sample—as loci associated with regulation of glucose homeostasis—were associated with T2D. Although likely not attributable to power (we had 80% power to detect modest effect sizes [odds ratio 1.10–1.15] among common variants [MAF >0.15] at stringent significance levels [P = 5.00 × 10−8]) (Supplementary Fig. 1), this observation likely reflects the pleiotropic nature of quantitative intermediate phenotypes of glucose homeostasis and the common observation that not all individuals with impaired glucose tolerance transition to overt T2D. Alternatively, a lack of association with T2D could reflect association of higher values for the quantitative traits with T2D, which conflicts with our underlying hypothesis. Although requiring further verification, these loci are still of substantial interest and could aid in understanding specific physiologic pathways that may ultimately lead to disease or phenotypic variation within the normal range.

GWAS for T2D have identified >70 susceptibility loci; association studies with quantitative traits have identified disturbed insulin secretion as the most frequent observation. The inability to identify insulin resistance loci may be partially explained by the high frequency of insulin resistance in nondiabetic control subjects (45). Furthermore, cohorts included in GWAS for T2D generally do not have detailed measures of insulin resistance. Fasting insulin and the closely related HOMA of insulin resistance have been most commonly used to represent insulin resistance in large-scale genetic studies (11). These traits only partially reflect insulin resistance (46) and therefore may be inadequate for gene discovery. Other than a small pilot study (47), the present GWAS is the first to include detailed measures of insulin resistance.

A few of the previously reported diabetes loci appear to act through altered insulin sensitivity (FTO, PPARG, IRS1, KLF14, ADAMTS9, GCKR, and RBMS1/ITGB6) (48), suggesting the likely presence of other, as yet undiscovered loci. The discovery of additional such traits was a major goal of GUARDIAN. However, consistent with prior GWAS, we did not identify any insulin sensitivity loci at genome-wide significance levels. Although none of the more modestly significant insulin sensitivity loci translated to T2D, rs7219451 and rs322394 nearly reached genome-wide significance. It is possible that environmental or lifestyle factors have a relatively greater effect on insulin sensitivity than genetic factors. We do not believe that differences in phenotyping of this trait (euglycemic clamp or FSIGT) hampered our ability to discover insulin sensitivity loci, given that these methods produce highly correlated measures (49).

Failure to meet genome-wide significance does not necessarily indicate that the detected variants are not of importance; such variants have been found to be enriched in enhancer elements in relevant tissues (50). Whether this is the case for the variants described herein will require further experimentation. To gain insight on the functional potential of our association signals, we queried ENCODE and GTEx databases. The linked SNAPC4 SNPs rs3812570 and rs4266763 (r2 = 0.85), associated with DI/T2D in the translational meta-analysis, had RegulomeDB scores of 1f and 1b, respectively, indicating a high likelihood of functionality based on eQTL evidence, residence in transcription factor binding sites, and DNase hypersensitivity sites. These SNPs are associated in multiple tissues with not only mRNA levels of SNAPC4 but also the nearby genes INPP5E and CARD9. INPP5E codes for an inositol polyphosphate-5-phosphatase that has been implicated in Golgi-vesicular trafficking (51), alterations in which might affect β-cell insulin granule formation. Of interest, another variant with a putatively functional RegulomeDB score, rs10870202 (score 1f, which is associated with AIRg/T2D in Table 3), is also an eQTL for INPP5E as well as DNLZ, a gene adjacent to CARD9 and SNAPC4, suggesting that this region on chromosome 9 may be key to insulin secretion. Additionally, SNP rs1978648 was associated with SI (RegulomeDB score 2b) and resides in a DNase hypersensitive region harboring multiple transcription factor binding sites in HepG2 cells.

In summary, GUARDIAN has performed the first GWAS to explore the genetic architecture of T2D-related quantitative phenotypes in a large Mexican American cohort. Because defects in the maintenance of glucose homeostasis are postulated to contribute to the development of T2D, a direct translation of the findings was performed to identify possible new disease pathways and test whether these variants explain T2D risk. Consistent with the literature, the present results suggest a strong contribution for variants that affect insulin secretion pathways as assessed by AIRg and DI (e.g., CDKAL1, MTNR1B, KCNQ1). Although novel signals of association with insulin sensitivity traits were observed, they did not translate with statistical significance to the clinically relevant phenotype of T2D. Of note, a novel association with glucose effectiveness was observed, adding further to the complex pathophysiology underlying T2D.