A Search for Variants Associated With Young-Onset Type 2 Diabetes in American Indians in a 100K Genotyping Array Free

The present data come from participants in a longitudinal study conducted in the Gila River Indian Community in central Arizona (19). In this study, community residents aged ≥5 years are invited to a research examination every 2 years. These examinations include a 75-g oral glucose tolerance test. Diabetes is diagnosed if the fasting plasma glucose concentration is ≥7.0 mmol/l, the 2-h plasma glucose concentration is ≥11.1 mmol/l (20), or the diagnosis is made during routine clinical care (19). DNA has been extracted from blood leukocytes. To detect variants associated with young-onset diabetes, individuals were selected from the extremes of the age-of-onset distribution. Thus, 300 case subjects were selected who had developed type 2 diabetes before the age of 25 years; for comparison, 334 control subjects were selected who were nondiabetic and aged ≥45 years when last examined. All individuals were full-heritage American Indian, and all individuals fulfilling criteria were selected regardless of affection status of other family members.

Although the case-control approach is potentially powerful for detection of associated variants (21,22), it is liable to spurious results due to population stratification. Therefore, to allow for within-family association tests that are robust to such confounding, siblings of these case or control subjects who defined discordant sibling pairs were selected. Thus, 121 siblings of case subjects were included who were nondiabetic when last examined and whose age was older than that at which the youngest case of diabetes onset in the family occurred. Likewise, 140 diabetic siblings of control subjects with younger age of onset than that of the oldest control subject in the family were included. These individuals constituted 340 discordant sibling pairs in 172 potentially informative sibships.

To examine the potential importance of associated markers on a population basis, a population-based sample of individuals from the longitudinal study was selected for genotyping with selected markers. This sample consisted of all participants from the longitudinal study with available DNA whose heritage was full Pima and/or Tohono O'odham (a closely related tribe); 1,561 of these individuals had diabetes, while 1,940 were nondiabetic. There were 2,834 individuals who were not included in the genome-wide association study, and analyses in these individuals were used to provide a replication of results from the genome-wide study. Analyses in all 3,501 individuals from the population study (who constituted 1,880 sibships) were used to determine population-based parameters. Because of differences in selection criteria, there were 228 individuals in the genome-wide study who were not in the population-based study. Characteristics of individuals in the various studies are shown in Table 1.

DNA was isolated using a proteinase K high-salt ethanol precipitation method. Before genotyping for the genome-wide association study, DNA was purified using Montage PCR plates (Millipore). Individuals were genotyped using the Affymetrix 100K Human Mapping array (Affymetrix, Santa Clara, CA), which contains 115,810 SNPs with known positions on the autosomal and X-chromosomes, according to the manufacturer's protocol. Genotypes were generated using the BRLMM (Robust Linear Model with Mahalanobis distance classifier and the addition of a Bayesian step) algorithm (23). A total of 50 individuals were genotyped in duplicate, and 5,122 SNPs were excluded from analysis that had discrepancy rates of >2.7% or that produced valid genotypes in <85% of individuals. The mean genotypic concordance rate among the 50 pairs of duplicates was 99.5%. Since spurious associations may occur more frequently with very rare alleles, a further 28,215 SNPs with minor allele frequency <1% were excluded. Hardy-Weinberg equilibrium was tested among all genotyped individuals using a continuity correction to produce a better approximate statistic for rarer alleles (24). Since systematic genotyping errors may produce severe violations of Hardy-Weinberg equilibrium, a further 2,429 SNPs were excluded that deviated from Hardy-Weinberg equilibrium at P < 0.001. Thus, the present genome-wide association analyses include results for 80,044 SNPs. Family relationships (parents and siblings) were confirmed by comparison of the observed proportion of alleles identical by state for these markers with that expected. SNPplex (Applied Biosystems, Foster City, CA) was used to genotype individuals in the follow-up population-based studies. Physical positions are given according to NCBI Build 35.

Analyses were performed using the SAS package (SAS Institute, Cary, NC). The association between genotype and case-control status was assessed with logistic regression. Genotype was analyzed as a numeric variable representing the number (0, 1, or 2) of copies of a given allele. For X-chromosome markers outside the pseudoautosomal regions, men were coded as homozygous. To account for the inclusion of multiple individuals in the same sibship, data were analyzed using a modified regressive model in which, for each individual, the prevalence of case status among all of his or her siblings was included as a covariate (25). This produces an approximation to the regressive model of Bonney, in which covariates are used to model the residual phenotypic correlation among siblings (26). To provide for a specific test of within-family association in sibships discordant for diabetes, data were also analyzed with conditional logistic regression (27). All models included sex as a covariate, and the likelihood ratio test was used to assess statistical significance. Odds ratios (ORs) were calculated from regression coefficients. In the event that one allele is absent from case or control subjects, parameter estimates will not reliably converge, but the likelihood ratio statistic is still approximately correct. In these cases, the estimated OR is infinite; therefore, we report an OR of infinity for these SNPs.

To summarize the results of case-control and within-family analyses, P values for the two tests were combined with Fisher's method (28). To maintain robustness to population stratification, priority was given to the P value for the within-family test (P_WI). This was accomplished by calculating a one-sided P value for the case-control test (P_CC1) for an association in the same direction as that observed in the within-family analysis. The contribution of the case-control result (P_CC*) to the combined test statistic was taken as the maximum of P_WI and 1 − (1 − P_CC1)², where P_CC1 employs a correction to partially account for the fact that the two tests have been performed in some of the same individuals. The combined test statistic was calculated as follows:

Simulation studies show that this method augments power of the within-family test by using information from the general test while maintaining robustness to population stratification (29). If P_WI and P_CC* were independent, the P value associated with χ²_CC-WI on 4 d.f. would have a uniform distribution under the null hypothesis (28). However, because P_CC* is truncated to be ≥P_WI and because of correlation among the tests, this “nominal” P value for χ²_CC-WI does not have a uniform distribution. To generate a P value corrected for these distributional deviations, 50 million pairs of standardized test statistics, representing within-family and case-control logarithms of the OR divided by their SEs, were simulated from a bivariate distribution with correlation of 0.32 (the observed correlation among these statistics for 80,044 SNPs). In these simulated data, χ²_CC-WI was calculated for each pair of tests, and the P value for the value of χ²_CC-WI observed for a given SNP was taken as the proportion of replicates at which the test statistic for these simulated values exceeded the observed value.

To assess experiment-wide statistical significance, affection status was permuted across individuals and analyses of the association of genotypes, with permuted affection status repeated genome wide. A total of 200 permutations were thus analyzed. To maintain the familial nature of the study, all data for a sibship were permuted in tandem (across sibships of the same size), and affection status of individuals, including covariates, was then further permuted within the sibship. The proportion of permutations in which the maximum of χ²_CC-WI exceeded the value observed for a given SNP was taken as the experiment-wide P value for that SNP. To compare the global distribution of test statistics with that expected under the null hypothesis of no association with any marker, the signed Kolmogorov-Smirnov deviation (d) was calculated for the observed distribution in comparison with the null distribution from each permutation (30). This statistic is the maximum deviation, over the full distribution of P values, of the observed cumulative distribution function (f) for any given observed P value [f(P)_obs] from the value expected under the null [f(P)_null]:

where i is an indicator variable that is 1 if f(P)_obs ≥ f(P)_null and −1 if f(P)_obs < f(P)_null at the point of maximum deviation. Under the global null hypothesis, the expected value of d is 0; thus, the test that the mean d = 0 across all permutations provides an empirical test of the global null hypothesis. Values of d > 0 indicate a shift toward lower P values in the observed compared with the null distribution. Since this analysis assumes exchangeability among men and women, it was restricted to the 78,568 autosomal SNPs.

SNPs that had the strongest associations in the genome-wide study were genotyped and analyzed in the population sample to examine the association on a population basis. To assess replication in a separate group of individuals, the association of genotype with prevalence of diabetes at the last available examination was analyzed among the 2,834 individuals not included in the genome-wide association study. These analyses were conducted using logistic regression models fit with generalized estimating equations to account for correlation among siblings (31). Within-family tests of association were also conducted using a modification of the method of Abecasis et al. (32), which partitions the association into between- and within-family components. In this method, the sibship mean of the numeric variable representing genotype is used to assess the between-family effects, and each individual's deviation from this mean is used to assess within-family effects. The P value for this within-family coefficient (P_WI) was further combined with the general association result using the modification of Fisher's method described above to produce a summary test of these two effects, with the difference that P_WI was calculated in a one-sided fashion to ensure that claims of replication would only be made if the direction of association was the same as that observed in the genome-wide association study. The distribution-corrected P value was calculated as described above by simulation of 50 million pairs of test statistics from a bivariate distribution in which the correlation was 0.52.

Analyses of the potential importance of these SNPs in the population were conducted among the 3,501 individuals from the population study. The hazard rate ratio (HRR) for diabetes associated with each copy of the marker allele was estimated. In this model, individuals who developed diabetes were considered to have developed the disease at the age of onset observed in the longitudinal study, while nondiabetic individuals were considered to be at risk for diabetes until the age at last examination. To account for familial resemblance among siblings, these analyses were conducted with a generalized estimating equations model in which diabetes incidence rates were represented as a Poisson function. The HRR was used to calculate population-attributable risk (PAR) for each marker (33). Where P_LL, P_LH, and P_HH represent, respectively, the proportion of individuals homozygous for the low-risk genotype and heterozygous and homozygous for the high-risk genotype, the PAR, under a multiplicative model, is as follows:

The PAR represents the proportion by which diabetes incidence would decrease if all individuals had the same risk as individuals with the low-risk genotype.

The distribution of P values for all 80,044 SNPs is shown in Fig. 1, along with the distribution expected under the global null hypothesis. Overall, the observed distribution of P values was similar to that expected under the null hypothesis; however, there was a slight excess of low P values beyond that expected, and this is more apparent if the portion of the distribution at P < 0.05 is examined (Fig. 1B). The average Kolmogorov-Smirnov deviation comparing the observed and expected distributions in 200 permutations was significantly greater than 0 (average d = 0.012, SE 0.00036, and P < 0.0001). This indicates that the observed distribution contains a statistically significantly greater proportion of low P values than would be expected under the global null hypothesis of no association with any SNP. Similar results were obtained with the within-family test alone (d = 0.006, SE 0.00035, and P < 0.0001). To explore genotyping artifact as a potential source of bias, analyses were repeated with more stringent thresholds for Hardy-Weinberg equilibrium, genotype success rate, and minor allele frequency with similar results. For the 25,915 SNPs with P > 0.05 for Hardy-Weinberg equilibrium, genotype success rate >0.99, and minor allele frequency >0.1, results were very similar (e.g., for the within-family test, d = 0.008, SE 0.00058, and P < 0.0001).

Results of tests for individual SNPs for association with young-onset type 2 diabetes are shown in Fig. 2. There were several regions where ≥1 SNP showed fairly strong evidence for association. The SNPs with the lowest P values were rs686989 (P = 2.7 × 10⁻⁶, which corresponds with an experiment-wide P value of 0.11) and rs672849 (P = 1.5 × 10⁻⁵, experiment-wide P = 0.55), both on chromosome 11 at 113.54 Mb; rs2164000 on chromosome 9 at 18.75 Mb (P = 2.3 × 10⁻⁵, experiment-wide P = 0.69); and rs10520926 on chromosome 5 at 25.36 Mb (P = 2.6 × 10⁻⁵, experiment-wide P = 0.73). The two chromosome 11 SNPs were highly concordant (r² = 0.99). There were 128 SNPs with P < 0.001 (∼80 expected under the null), and these SNPs were further genotyped in individuals from the larger population. The 128 SNPs are listed in supplementary Table 1 (available in an online appendix at http://dx.doi.org/10.2337/db07-0462).

There were 11 of these SNPs (out of 119 successfully genotyped) that also showed evidence for association (P < 0.05) with diabetes in individuals from the population who were not included in the genome-wide study and for whom the direction of association was the same as that in the genome-wide study. One would expect approximately six such replications by chance at this level of significance. Results of analyses for these SNPs are shown in Table 2. The chromosome 11 SNPs at 113.54 Mb (rs686989 and rs672849) and the chromosome 5 SNP at 25.36 Mb (rs10520926) were among those replicated. In general, the HRR in the population study was lower than the OR from the original study; this is expected given the selection procedure and the fact that the OR only approximates the HRR under certain assumptions (e.g., where incident cases are sampled or the disease is rare [34]). The HRR for most alleles is modest, and those with high PAR are largely those for which the risk alleles are at high frequency.

To assess the extent to which SNPs identified as associated with type 2 diabetes in the present study replicated in other populations, data from three other studies conducted with the Affymetrix 100K array were analyzed for all 646 SNPs with P < 0.007 in the present study. These data were provided by investigators from the Amish Family Diabetes Study (35), the Framingham Heart Study (36), and the Starr County Diabetes Study of Mexican-Americans (37). P values for an association in the same direction observed in the present study were combined across all three other studies by Fisher's method (28). Of the SNPs, 88.2% had valid data for all three other studies, 6.0% for two other studies, 5.6% for one other study, and 0.2% (one SNP) for none. There were 30 SNPs, shown in Table 3, that replicated the results of the present study at P < 0.05 (∼32 expected under the null). None of these overlap with those shown in Table 2. The overall analysis strategy and the number of SNPs selected at various stages of analysis are summarized in Fig. 3.

Type 2 diabetes is partially genetically determined (1–6), but susceptibility variants conclusively identified to date explain only a small portion of the familial risk. It is therefore likely that there are additional susceptibility variants that have not yet been identified. Genome-wide association studies are a potentially powerful way to detect these variants. The present analysis represents a genome-wide association study in Pima Indians, a population with a high prevalence of type 2 diabetes and obesity (3,19). Diabetes in this population is highly familial, particularly when it occurs at young ages (3,17–18,38); therefore, the present study was designed to detect variants associated with young-onset diabetes. Although case-control and related association tests can be powerful, they are liable to confounding by population stratification (39,40). In contrast, within-family tests, though less powerful, are robust to such confounding (40). In the present study, both case-control and family-based designs were employed and results combined so that the case-control results augment the power of the within-family test while maintaining robustness to population stratification (29). Thus, markers identified with this approach are likely to reflect both association and linkage with diabetes susceptibility variants.

Many of the associations observed in the present study have P values that are quite strong by conventional criteria. However, the interpretation of P values in genetic association studies is problematic. Because the prior probability that any given variant is in linkage disequilibrium with a disease susceptibility allele is low, many statisticians recommend very stringent criteria for declaration of statistical significance, e.g., P values from 10⁻⁵ to 10⁻⁸ (41–43). The problem is compounded in genome-wide studies, where because of multiple tests, one would expect some strong associations to occur by chance. Classical methods of adjustment for multiple comparisons, such as Bonferroni, are too stringent for genome-wide association studies because they assume that tests are independent and ignore linkage disequilibrium among markers. Furthermore, while corrections for multiple testing proceed on the assumption that the global null hypothesis of no association with any marker is of interest (44), the corrections are applied on a marker-wise basis and, thus, may not efficiently utilize information from the full distribution of P values relevant to this global null.

In the present analysis, a permutation procedure was used to assess experiment-wide statistical significance. This procedure maintains the linkage disequilibrium structure among markers. In addition, the Kolmogorov-Smirnov test was used to compare the overall distribution of P values observed in the actual data with that observed in data permuted under the null hypothesis. Although no single marker achieved experiment-wide statistical significance, analysis of the distribution of P values resulted in strong rejection of the global null hypothesis of no association with any marker. This result is expected if there are multiple susceptibility variants because, in such a situation, a test of the full distribution can be more powerful than a test of a single variant. A systematic bias is an alternative possibility that is difficult to completely exclude, but the present techniques are robust to bias from population stratification. The present results were also unmodified when restricted to SNPs with increased stringency of Hardy-Weinberg equilibrium, allele frequency, and genotype success rate, and this suggests that they are not explained by bias due to genotypic artifacts that can be captured by these measures. The number of true functional susceptibility variants is difficult to determine, given that modest linkage disequlibrium may extend over long distances and that many associations may be due to chance. Thus, while the present results strongly suggest that some of the low P values observed reflect linkage disequilibrium with diabetes susceptibility loci, the distinction between true- and false-positive results is difficult.

Replication of the association in a separate group of individuals provides additional evidence that a marker is associated with diabetes. In the present study, SNPs with the strongest evidence for association in the genome-wide study were further evaluated in a separate set of individuals from the population. (However, some of these individuals were related to those in the genome-wide study.) Several SNPs showed nominal evidence for association (P < 0.05) in this replication set. This analysis is limited in power because the vast majority of individuals from the extremes of the age-of-onset distribution, who provide much of the statistical power, were excluded due to participation in the genome-wide study. In addition, one would expect some of the markers to show association by chance. Thus, some false-positives are likely to remain among the SNPs that replicated and some true positives among those that did not. However, the probability that SNPs showing replication are in linkage disequilibrium with diabetes susceptibility variants is enhanced above that for the SNPS for which replication was not observed. The result for rs686989 is of particular interest because it is in a region identified as linked to diabetes and obesity in a previous genome-wide linkage study in this population (15); however, the association with rs686989 itself does not explain the linkage signal (data not shown).

The PAR, calculated in a group representative of the full population, provides a measure of the potential importance of each associated SNP, and this information can be used to prioritize regions for follow-up studies. The PAR is calculated on the assumption that the observed association is causal; therefore, it may be underestimated if the marker is not highly concordant with a functional allele or overestimated if confounded by population stratification. In the present analysis, the PAR was calculated from longitudinal data observed in the population study. Thus, in contrast to estimates derived from case-control studies based on prevalent cases, the present results do not depend on the questionable assumptions of a rare disease or of sampling only incident cases (34).

Studies of these SNPs in other populations may also be relevant in prioritization of regions for follow-up studies. The present results were further compared with those obtained from three other genome-wide association studies of type 2 diabetes on the Affymetrix 100K array (35–37). Comparison among studies is complicated by the fact that all had different study designs; the present study focused on young-onset diabetes and gave priority to within-family tests. (Characteristics of each study are presented in supplementary Tables 2 and 3.) Nonetheless, several of the SNPs identified in the present study had some evidence for association in the other studies (overall P < 0.05). These regions are also high priority for follow-up studies. It is noteworthy that rs516415, which is also in the diabetes-linked region on chromosome 11q, showed replication, as did rs1886004, which is in a region on chromosome 1q linked to diabetes in Pima Indians and other populations (14–15). However, none of the SNPs with evidence for replication in the Pima Indians at P < 0.05 also had a combined P < 0.05 in these other studies. This may reflect low power across the studies to detect alleles of modest effect.

The power of genetic association studies depends on the frequency of the functional polymorphisms and on the magnitude of their effects (22). Given that case subjects represent the lower 10% of the age-of-onset distribution and control subjects the upper 45%, we estimate that the present sample size has ∼75% power to detect a common (minor allele frequency >0.1) functional allele at P < 0.001 that explains 3% of the variance in age at onset of diabetes. Power also depends on the extent to which one of the typed markers is strongly concordant with an allele that influences disease susceptibility (21,22). Since the 100K array used in the present study does not exhaustively capture allelic variation, it is possible that additional regions that contain diabetes susceptibility variants were not identified. The pattern of linkage disequilibrium among SNPs identified in the HapMap project suggests that ∼30% of common variants have r² > 0.8 with a marker on this array in non-African populations (45). American Indians are not included in the HapMap, and they tend to share relatively few common haplotypes with HapMap populations (46). This may be reflected in the fact that a larger proportion of markers were nearly monomorphic in the present study than in the other studies that used the 100K array (35–37), and this could result in lower power for this array to detect associations in American Indians. On the other hand, surveys of various populations suggest that linkage disequilibrium is higher among American Indians than in many other populations such that fewer variants are required to capture common haplotypes. (46). Thus, fixed marker sets, such as the 100K array, may capture common variation in American Indians to an extent similar to that in other non-African populations.

Recently, several high-density genome-wide association studies have been conducted in northern European populations, and these have identified six gene regions (apart from TCF7L2, PPARG and KCNJ11) that contain putative diabetes susceptibility variants (47–50). None of the SNPs most strongly associated in the present study was in any of these regions. Furthermore, none of the SNPs consistently associated across these other genome-wide studies is well-captured by the 100K array. However, as shown in Table 4, some SNPs in these regions had modest evidence for association in the present study (P < 0.05). These putative diabetes susceptibility variants have been largely identified in northern Europeans, and it is not clear whether they play a significant role in diabetes in American Indians or other high-risk populations. A more exhaustive survey of variation in these regions is required to quantify their role in diabetes susceptibility in the Pima Indian population. The variants in TCF7L2 that are most strongly associated with type 2 diabetes in other populations (9–13) have been typed in the Pima Indians, in whom there is little evidence that they influence diabetes risk (51).

Genome-wide mapping studies are typically only an initial step in the elucidation of susceptibility variants. The present analyses have identified several regions that may harbor genetic variants that influence susceptibility to young-onset diabetes in American Indians. Several of the associated SNPs are in or near genes, including ZNF659 (chromosome 3, 21.47 Mb), FANCF (chromosome 11, 22.60 Mb), ZBTB15 (chromosome 11, 113.54 Mb), and SENP1 (chromosome 12, 46.71 Mb). Fine-mapping studies of these regions are needed to confidently localize the signals to specific genes. Confirmation of the role of genes in the regions identified in the present study will require replication studies in other populations and, ultimately, functional studies.

This research was supported by the Intramural Research Program of the National Institute of Diabetes and Digestive and Kidney Diseases (NIDDK).

We thank Jose Florez for review of the manuscript. We also thank the staff of the Diabetes Epidemiology and Clinical Research Section and the Diabetes Molecular Genetics Section (NIDDK, Phoenix, Arizona) for assistance. We further thank members of the Gila River Indian Community who participated in the study.