Analysis of SLC16A11 Variants in 12,811 American Indians: Genotype-Obesity Interaction for Type 2 Diabetes and an Association With RNASEK Expression

The Slim Initiative in Genomic Medicine for the Americas (SIGMA) Type 2 Diabetes Consortium recently reported a haplotype consisting of five common exonic variants, including rs75493593, in the monocarboxylic acid transporter gene SLC16A11 that was strongly associated with type 2 diabetes in Mexican and other Latin American populations (1). The diabetes risk haplotype had a frequency of 50% in Native Americans from Mexico but was rare in Europeans and Africans (1), suggesting it may contribute to the higher burden of type 2 diabetes among American Indians. Our previous study of 63 established type 2 diabetes variants in 3,421 individuals from a longitudinal diabetes study in a Southwestern American Indian (SWAI) population showed little evidence of association between rs75493593 and type 2 diabetes (2). To further assess whether SLC16A11 variants are associated with type 2 diabetes in American Indians, we expanded our analysis by genotyping three tag single nucleotide polymorphisms (SNPs) (rs75493593, rs148775056, and rs9303212) that capture all variants with a frequency >0.01 across the SLC16A11 locus in our population-based SWAI sample (N = 7,295). As one of the tag SNPs, rs75493593, captures the previously reported Mexican diabetes risk haplotype, it was further analyzed in two additional samples of American Indians from the Phoenix extension of the Family Investigation of Nephropathy and Diabetes (FIND) (N = 3,095) and the Strong Heart Study (SHS) (N = 2,421) for a combined total of 12,811 North American Indians. As the SIGMA study found that the SLC16A11 diabetes risk alleles were associated with leaner BMI in individuals with diabetes, we also examined whether SLC16A11 genotypes interact with BMI in their association with diabetes. To date, none of the SLC16A11 missense variants have been shown to alter protein function; therefore, we merged rs75493593 genotypic data with gene expression data to search for potential expression quantitative loci (eQTLs) that could provide an alternative mechanism affecting the risk for type 2 diabetes.

Brief descriptions of the study cohorts are provided below and their clinical characteristics are shown in Supplementary Table 1. Each study was approved by an institutional review board and participants gave informed consent.

From 1965–2007, SWAIs from a community in Arizona participated in biennial health exams that included a 75-g oral glucose tolerance test (OGTT), measures of height and weight to calculate BMI, and a blood draw for DNA isolation (3). Diabetes status was determined according to the 1997 criteria of the American Diabetes Association (fasting plasma glucose ≥7.0 mmol/L or 2-h plasma glucose ≥11.1 mmol/L) (4) or from a previous clinical diagnosis. Fasting glucose, 2-h glucose, and insulin concentrations were analyzed using data from the participant’s last nondiabetic exam. Estimates of insulin resistance (HOMA-IR) and β-cell function (HOMA-B) were calculated using the homeostatic model (5). The current study includes all individuals (N = 7,710) for whom diabetes status is recorded and DNA is currently available.

Participants from the Phoenix extension of FIND, a multicenter study designed to identify genes involved in diabetic nephropathy (6), were required to be of ≥50% American Indian heritage (self-reported) and most were urban-dwelling American Indians living near Phoenix, AZ. Individuals were examined once and diabetes status was determined by 2010 American Diabetes Association criteria based on fasting blood glucose concentrations (≥7.0 mmol/L) or HbA_1c (≥6.5%) (7) or from a prior clinical diagnosis (6). The current study includes all participants (N = 3,106) who were not also studied as part of the SWAI sample.

Participants are American Indians who underwent health exams at study sites located in Arizona, North and South Dakota, and Oklahoma (8). Diabetes status was determined by 1998 World Health Organization criteria based on 2-h blood glucose ≥11.1 mmol/L, fasting blood glucose ≥7.0 mmol/L, current use of diabetes medication, and/or if an individual was previously diagnosed with diabetes (9). The current study includes all participants (N = 2,451) who were not also studied as part of the SWAI sample.

Some individuals (N = 561, age 18–45 years) from the SWAI sample also participated in inpatient metabolic exams in our Clinical Research Center. These individuals did not have diabetes, as determined by a 75-g OGTT; were healthy based on medical history, physical examination, and routine laboratory tests; and were not taking any medications. The OGTT was performed after an overnight fast and blood was drawn before and at 30, 60, 120, and 180 min after glucose ingestion for determining plasma glucose and insulin concentrations. Insulin action was assessed using the two-step hyperinsulinemic-euglycemic clamp technique to measure insulin-stimulated glucose disappearance (uptake) rates (10). To determine acute insulin response, a 25-g bolus of glucose was injected intravenously over 3 min and blood samples were collected before infusion and at 3, 4, and 5 min after the injection to determine plasma glucose and insulin concentrations. Acute insulin response was calculated as one-half the mean increment in plasma insulin concentrations from 3 to 5 min (10). Percent body fat was estimated by underwater weighing or by total body dual energy X-ray absorptiometry (DPX-1; Lunar Radiation, Madison, WI) (11). A subset of the inpatient subjects (N = 423) also participated in our metabolic chamber studies to determine energy expenditure and substrate (carbohydrate, lipid, and protein) oxidation in the metabolic respiratory chamber (12).

Genotyping was done by allelic discrimination using TaqMan genotyping assays (Applied Biosystems, Foster City, CA). All genotypic data included in the analyses met our quality-control criteria, which requires a successful call rate of ≥95%, a lack of deviation from Hardy-Weinberg equilibrium at P < 0.001, and a discrepancy rate of ≤1% in masked duplicates (430, 120, and 88 blind duplicates for the SWAI, FIND, and SHS samples, respectively).

The association between genotypes and type 2 diabetes was determined with logistic regression modeling (additive model), where homozygotes for the major allele (A1/A1), heterozygotes (A1/A2), and homozygotes for the minor allele (A2/A2) were coded to a numeric variable for genotype (0, 1, and 2). In individuals from the SWAI and SHS community-based studies, the model was fit with the generalized estimating procedure to account for sibship. For the FIND sample, the genomic control procedure (13) was used to account for unspecified relationships based on 42 random markers from throughout the genome. Individual estimates of European admixture, which were used as covariates in analyses, were derived from 45 markers with large differences in allele frequency between populations (14) with the method of Hanis et al. (15). A meta-analysis for the SWAI, FIND, and SHS samples was conducted by the inverse variance method (16). Heterogeneity among the samples was quantified by the I² measure and statistical significance was tested by Cochran Q statistic (17). For longitudinal analyses of diabetes, SWAI individuals who did not initially have diabetes were followed from their first examination after age 15 years until development of diabetes or until their last examination. Cox proportional hazards regression was used to estimate the risk (hazard ratio [HR]) of developing type 2 diabetes. A similar analysis was performed for the SHS using data from phases I through III. Linear regression was used to analyze the associations between genotype and continuous variables, and the models were fit with the generalized estimating equation procedure to account for sibship. Statistical analyses were performed using the statistical analysis system of the SAS Institute (Cary, NC).

Genome-wide mRNA levels had previously been measured using the Affymetrix 1.0 Human Exon microarray in adipose tissue (N = 187) and skeletal muscle (N = 196) biopsies from a subset of the SWAI sample (all were without diabetes) (18). Genome-wide mRNA levels in whole blood were also available from a prior study in a subset of the FIND sample (N = 1,412, 23% with diabetes) who had been genotyped for 42 established type 2 diabetes–susceptibility variants including rs75493593. Whole-blood gene expression levels were measured using the Illumina HumanHT-12 v4 Expression BeadChip. A χ² tail test was used to assess associations between genotype and 15,854 transcripts expressed above background (19) that mapped to unique locations and that were free of variants located in probe set binding sites according to ReMOAT annotation (20). For variants associated with both type 2 diabetes (P < 0.05) and whole-blood transcript levels (P < 0.005), a formal mediation analysis was used to quantify the extent to which the pattern of associations was consistent with an effect on transcript that mediates the variant-diabetes association. Significance of mediation was assessed by the Sobel test (21), and the percentage of the variant effect on diabetes that was potentially mediated by the transcript was calculated (22). Additional details are described in the Supplementary Data online. As associations between SLC16A11 variants and type 2 diabetes have previously been established in other populations, we considered a P < 0.05 as statistically significant for variant-diabetes–related phenotype associations. In the whole transcriptome analyses, which lack a specific a priori hypothesis, we calculated the false discovery rate (23) to account for multiple comparisons.

Previously analyzed whole-genome sequence data (40× coverage) from 296 SWAI individuals were used to identify 14 variants with a minor allele frequency ≥0.01 in the genomic region spanning the SLC16A11 locus (Supplementary Fig. 1A). Linkage disequilibrium (LD) information (r² ≥0.85) for the 14 variants could be captured by genotyping three tag SNPs, rs75493593, rs148775056, and rs9303212 (Supplementary Fig. 1B). Most of the common variants, including all five exonic variants that make up the previously reported type 2 diabetes risk haplotype in the SIGMA study, fell into one major LD group tagged by rs75493593. The three tag SNPs were genotyped in our population-based SWAI sample (N = 7,295), but none were associated with type 2 diabetes (Supplementary Table 2). As rs75493593 tagged the diabetes risk haplotype reported in the SIGMA study, this variant was further genotyped in two additional samples of American Indians from the FIND (N = 3,095) and SHS (N = 2,421) studies. Meta-analysis of the three samples (N = 12,811) detected a nominal association between the rs75493593 A allele and higher risk of type 2 diabetes (P = 0.001, odds ratio [OR] 1.11 [95% CI 1.04–1.19]) (Table 1), which is directionally consistent with the Mexican study.

A significant genotype × BMI interaction was observed for the association of rs75493593 with type 2 diabetes (combined American Indian sample, P_interaction = 5.1 × 10⁻⁷) (Fig. 1 and Supplementary Table 3). After stratifying the individuals into BMI categories, a modest association was observed between rs75493593 and diabetes for the nonobese groups (in the combined American Indian sample, P = 4.9 × 10⁻³ and 7.5 × 10⁻⁶ for the BMI categories <25 and 25–30 kg/m², respectively), where the A allele was associated with increased risk for diabetes (OR 1.41 and 1.36, respectively), but among the heavier individuals (BMI 35–40 and >40 kg/m²), the A allele was associated with decreased risk for diabetes (OR 0.88 and 0.84, respectively). Similar trends were observed in each American Indian sample when analyzed separately (Fig. 1 and Supplementary Table 3). This same phenomenon is reflected in a strong diabetes × genotype interaction for the association of rs75493593 with BMI (in the combined American Indian sample, P_interaction = 3.4 × 10⁻¹⁰) (Table 2). Although an association between rs75493593 and BMI in the combined sample was observed (6.4 × 10⁻⁵), where individuals with diabetes who carry the risk A allele have lower BMIs, this association was due entirely to individuals with type 2 diabetes (with diabetes, P = 3.4 × 10⁻¹⁵; without diabetes, P = 0.77) (Table 2). The distribution of genotype across BMI categories and diabetes status is given in Supplementary Table 4. The association between the risk haplotype and lower BMI in individuals with diabetes was also observed in the Mexican study (P = 5.2 × 10⁻⁴) (1); in combination with the present data, rs75493593 is a replicated variant for BMI at genome-wide statistical significance in individuals with diabetes (combined Mexican and American Indian sample, P = 5.2 × 10⁻¹⁶).

To assess whether individuals with the rs75493593 A allele develop diabetes at a leaner BMI or lose more weight after diabetes onset, 4,088 SWAI and 1,073 SHS individuals who did not have diabetes at their first exam and had longitudinal follow-up data were analyzed. Among the SWAI individuals (mean age 22.0 years), 1,498 (37%) developed diabetes during a mean follow-up time of 13.5 years; for the SHS participants (mean age 55.1 years), 222 (21%) developed diabetes during a mean follow-up of 6.6 years. A proportional hazards regression analysis (HR) showed that rs75493593 did not predict type 2 diabetes in either sample (SWAI, HR 1.04 [95% CI 0.97–1.12], P = 0.31; SHS, 1.00 [0.81–1.23], P = 0.97) or when the SWAI and SHS samples were combined (1.03 [0.97–1.11], P = 0.35). However, when grouped by BMI, a nominally significant HR was detected for SWAI individuals with a BMI <25 kg/m² (1.18 [1.01–1.38], P = 0.03) (Fig. 2), where the rs75493593 A allele was associated with increased risk of diabetes, a finding consistent with the cross-sectional analyses. The HR for rs75493593 was similarly elevated in the SHS individuals with BMI <25 kg/m², and the association was nominally significant in the combined sample (1.20 [1.03–1.39], P = 0.02) (Fig. 2). However, a significant rs75493593 genotype × baseline BMI interaction was not observed in the SWAI, SHS, or combined samples (P_interaction = 0.94, 0.08, and 0.57, respectively) (Fig. 2). The longitudinal relationship with change in BMI, as measured by the slope, across nondiabetic exams and across exams conducted after diabetes onset was also analyzed. When SWAI participants did not have diabetes, there was no difference in BMI change by genotype (N = 4,480, P = 0.32) (Fig. 3A). In contrast, after type 2 diabetes onset, there was an association between rs75493593 and BMI slope (N = 1,632, P = 0.02) (Fig. 3A), where the A allele was associated with increased weight loss. The decrease in BMI after diabetes onset was 1.2% versus 0.5% per year for individuals homozygous for the A and C alleles, respectively. There were no significant associations between rs75493593 and weight change in the SHS individuals, although it should be noted that the weight changes in this older population were smaller in magnitude (Fig. 3B).

To evaluate the physiology underlying the associations with type 2 diabetes, glucose and insulin concentrations along with HOMA measures of insulin secretion and resistance were analyzed in SWAI individuals who had available data from their last nondiabetic exam (N ∼5,200). Nominal evidence for a rs75493593 genotype × BMI interaction was detected for various traits including fasting insulin, HOMA-IR, and HOMA-B (P_interaction = 0.03–0.04, Supplementary Table 5); however, these associations are too weak to be conclusive of the underlying physiology. In a subset of SWAI individuals without diabetes who had undergone metabolic testing in our Clinical Research Center, there was a nominal but consistent association between rs75493593 and basal carbohydrate and lipid oxidation when measured during a hyperinsulinemic-euglycemic clamp (N = 521) and when measured in a human respiratory chamber (N = 423) (Supplementary Table 6). Individuals homozygous for the A allele had higher mean carbohydrate oxidation rates and lower mean lipid oxidation rates compared with the other two genotypes (P = 0.007 and 0.006, respectively) (Fig. 4A and B), resulting in a slightly higher mean 24-h respiratory quotient (P = 0.01). These associations, although weak, are consistent between two independent in vivo measures of carbohydrate and lipid oxidation.

On the basis of SIFT prediction, only one of the four SLC16A11 missense variants, rs13342692 (D127G), is predicted to be damaging with a SIFT score of 0.05 (Supplementary Table 7), just reaching the cutoff for deleterious amino acid substitutions (SIFT score thresholds: tolerated >0.05, damaging ≤0.05 [24]). To search for additional mechanisms whereby these variants could have a functional effect, genotypic data for rs75493593 were merged with previously measured global adipose tissue (N = 187) and skeletal muscle (N = 196) gene expression data to determine whether variants tagged by rs75493593 could function as an eQTL. The genome-wide eQTL analyses identified RNASEK, which is located approximately 27 kb upstream of rs75493593 (SLC16A11), as one of the top eQTL genes in both tissues that were biopsied from different SWAI individuals (adipose tissue P = 2.1 × 10⁻⁴ and skeletal muscle 1.7 × 10⁻⁸) (Table 3). The microarray mean RNASEK expression levels based on rs75493593 genotypes are shown in Fig. 5A and B. The association between rs75493593 and adipose tissue RNASEK expression was replicated by quantitative RT-PCR (N = 182, P = 8.0 × 10⁻⁴) (Fig. 5C). There were no microarray expression data for SLC16A11 in either adipose tissue or skeletal muscle; however, directly measuring SLC16A11 mRNA levels by quantitative RT-PCR showed no association between rs75493593 and expression in adipose tissue (N = 182, P = 0.30, data not shown). The association between rs75493593 and RNASEK expression was further analyzed in whole blood from a subset of the FIND sample (N = 1,412, 23% with diabetes) for whom both genome-wide expression data and genotyping data for 42 established type 2 diabetes variants, including rs75493593, were available. In these individuals, the rs75493593 A allele was also associated with lower RNASEK expression in whole blood (P = 1.6 × 10⁻⁵) (Fig. 5D) and increased risk for diabetes (P = 3.0 × 10⁻³, OR 1.36) (Fig. 5E); lower RNASEK expression was associated with an 18% increase per standard deviation in the odds of type 2 diabetes (P = 6.0 × 10⁻³, OR 0.82) (Fig. 5F). To examine the extent to which the rs75493593-RNASEK eQTL could account for the association between genotype and diabetes, a formal mediation analysis was conducted and significant mediation was observed (P = 0.01). To compare this mediation effect with the potential effect of other type 2 diabetes susceptibility variants in a transcriptomic context, the remaining 41 established type 2 diabetes variants genotyped in the 1,412 FIND participants were also analyzed. Six of the 41 variants were associated with type 2 diabetes (P < 0.05, data not shown), and these 6 variants along with rs75493593 gave rise to 85 potential eQTL genes (P < 0.005 for association between variant and transcript levels) (Supplementary Table 8). Only 3 additional eQTLs displayed statistically significant (P < 0.05) mediation, and none of the 3 displayed a stronger mediation effect than the rs75493593-RNASEK eQTL (Supplementary Table 9).

In the SIGMA study, it was noted that among all individuals with type 2 diabetes, those who carry the risk haplotype had lower BMIs than those who do not carry the risk haplotype, suggesting that diabetes onset occurred at a lower BMI among carriers (1). Our study of American Indians identified a significant rs75493593 genotype × BMI interaction where individuals with a BMI <35 kg/m² who carry the diabetes risk A allele were at increased risk for type 2 diabetes and individuals with higher BMIs who carry the A allele were protected from diabetes. Consistent with the cross-sectional results, for the SWAI and SHS individuals with longitudinal data on BMI and diabetes status, we observed that HRs for type 2 diabetes by genotype were higher in the leanest BMI category (<25 kg/m²) and lower in the heaviest BMI category (>35 kg/m²). As analyses of diabetes incidence, by necessity, exclude prevalent cases at baseline, the HRs from the longitudinal analyses may not be reflective of lifetime diabetes risk. However, the cross-sectional and longitudinal results together suggest that BMI may have a modulatory effect on the development of type 2 diabetes in individuals who carry the rs75493593 A allele. We also detected a strong rs75493593 genotype × diabetes interaction for BMI in the American Indians, where an association between rs75493593 and BMI only occurred in individuals with type 2 diabetes. However, the cross-sectional analyses cannot determine whether individuals were leaner before or after the onset of type 2 diabetes (i.e., are only lean individuals at increased risk for type 2 diabetes if they carry the A allele or do individuals with the risk allele who develop type 2 diabetes lose more weight as a consequence of the disease?). It has previously been shown that American Indians, like individuals in other ethnic groups, often lose weight after they develop diabetes (25). Analysis of the SWAI individuals with longitudinal data suggests that individuals with the risk A allele lose weight after developing type 2 diabetes. There was no difference in BMI by genotype in individuals who had not developed type 2 diabetes, but a nominal association (P = 0.02) between BMI and genotype was observed in individuals following onset of the disease, where individuals carrying the diabetes risk allele had greater weight loss. It is unclear whether the weight loss is due to a more severe pathophysiology of diabetes or medical treatment. Nevertheless, the association between rs75493593 and weight loss after diabetes diagnosis observed in the SWAI sample is very modest, suggesting that there may be other genetic or physiological reasons for the strong association between the diabetes risk allele and lower BMI only in the individual with diabetes. Although analysis of in vivo metabolic data available for a subset of SWAI subjects did not provide a conclusive explanation of the physiological link between rs75493593 and type 2 diabetes, a modest but reproducible effect on lipid oxidation rates was observed, which is consistent with the suggestion that SLC16A11 may be involved in lipid metabolism (1). Our genome-wide eQTL analysis suggests that rs75493593 or a variant in high LD may also affect expression of RNASEK, which is located approximately 27 kb upstream of SLC16A11, where individuals carrying the type 2 diabetes risk A allele have lower RNASEK expression levels compared with individuals homozygous for the C allele. Given these results observed in three different tissues that were biopsied from different individuals, we propose that RNASEK expression may be mediating, at least in part, the association between rs75493593 and type 2 diabetes. Our mediation analyses demonstrate that the pattern of associations is consistent with this hypothesis; however, these results could be confounded by complex pleiotropic relationships among the traits or LD between causal variants that independently influence RNASEK expression and diabetes risk. Studies in additional populations and molecular mechanistic analyses are required to evaluate this hypothesis more rigorously. Very little is known regarding RNASEK except that it is expressed in most tissues and encodes an endoribonuclease that cleaves ApU and ApG phosphodiester bonds (26,27). Additional studies are needed to determine how RNASEK may play a role in the development of type 2 diabetes.

In conclusion, our study identified a modest association between the SLC16A11 variant rs75493593 and type 2 diabetes in 12,811 American Indians, where the effect on diabetes was much more pronounced in nonobese individuals. Analysis of longitudinal measures of BMI suggests that at least some of the difference in BMI is due to greater weight loss in those that carry the A risk allele and develop type 2 diabetes. rs75493593 was also associated with RNASEK gene expression, suggesting that either altered SLC16A11 protein function, RNASEK gene expression, or both may be contributing to the phenotypic associations observed in this study.

Acknowledgments. The authors thank all of the American Indian volunteers included in this study.

Funding. This research was supported by the Intramural Research Program of the National Institute of Diabetes and Digestive and Kidney Diseases, National Institutes of Health. The SHS was supported by cooperative agreement grants U01-HL41642, U01-HL41652, U01-HL41654, U01-HL65520, and U01-HL65521 from the National Heart, Lung, and Blood Institute. This study utilized the computational resources of the Biowulf system at the National Institutes of Health, Bethesda, MD (for a description, see https://helix.nih.gov/).

The opinions expressed in this article are those of the author(s) and do not necessarily reflect the views of the Indian Health Service who helped support the SHS.

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

Author Contributions. M.T., R.L.H., L.J.B., and C.B. contributed to the study design, researched the data, and wrote the manuscript. A.M., S.K., P.P., S.C., J.E.C., J.B., H.G., S.K., R.G.N., and B.V.H. contributed data. W.C.K. contributed to the study design and contributed data. C.B. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.