Epistasis Between HLA-DRB1*16:02:01 and SLC16A11 T-C-G-T-T Reduces Odds for Type 2 Diabetes in Southwest American Indians

In 2011, we reported the protective effect of HLA-DRB1*02 on the susceptibility to type 2 diabetes (T2D) via an effect on insulin secretion in a population of southwestern indigenous Americans that suggested a role for the immune system in T2D as well as type 1 diabetes (T1D) (1). In a recent multiethnic study of T2D that included a stratum of American Indians, it was reported that 41.3% of patients had cellular islet autoimmunity and 13.5% had humoral islet autoimmunity (2). The natural correlate of these findings is to search for risk and protection HLA alleles for T2D in our southwestern indigenous sample with a very high prevalence of the disease (3,4). Since our 2011 article (1), we have HLA typed nearly the entire study group with targeted next-generation sequencing and computer algorithms that revealed high-resolution alleles. Our HLA-DRB1*02 resolves to DRB1*16:02:01 (5–7). In the present article, we present an analysis of the role of this protective molecular allele and its interaction with the monocarboxylate transporter SLC16A11, a known T2D risk locus specific to American ancestral groups, and determine their combined contributions to the genetic immunological risk for T2D. We used binding prediction algorithms to identify 15-mer SLC16A11 peptides that might work synergistically with DRB1*16:02:01 in the mechanics of protection, develop polygenic statistical models that reflect their joint action, and test for interaction and epistasis.

A haplotype of SLC16A11, defined by five single nucleotide polymorphisms (SNPs), four missense and one synonymous, including rs117767867(T), rs13342692(C), rs13342232(G), rs75418188(T), and rs75493593(T), was reported by the Slim Initiative Genomic Medicine for the Americas (SIGMA) Type 2 Diabetes Consortium to be strongly associated with T2D in Mexicans and Latin Americans, with age and weight effects (8). In addition, this ancient human haplotype dates from the origins of modern humans and the genus Homo as much as 500,000–750,000 years ago (9). We subsequently used the G340S missense mutation and reported risk for T2D in leaner people and protection from the disease in very heavy individuals (10). The association was replicated in the Hispanic Community Health Study/Study of Latinos (HCHS/SOL) in a large Mexican American study but not in other Latino samples with little indigenous American admixture (11). The search for the association of the ancient human haplotype with other expressions of T2D and with clinical markers of the disease has continued to be an active area of research (12–21). The ancient human SLC16A11 haplotype shares features with the concept of a private allele that was first defined by James Neel (22).

A sample of 6,669 Southwest American Indians who have some portion of American Indian heritage, who participated at some time in a longitudinal study of T2D conducted from 1965 to 2007, and who had exome sequencing data available were chosen for this study.

Diabetes was diagnosed using 1997 American Diabetes Association criteria: fasting plasma glucose ≥7.0 mmol/L, 2-h plasma glucose ≥11.1 mmol/L, or a diagnosis made in routine clinical care as previously described (10).

The five coding SNPs that define the ancient SLC16A11 haplotype were genotyped for association analyses using the TaqMan Allelic Discrimination (AD) Assay (Applied Biosystems) on an ABI 7900HT (Applied Biosystems). For genotyping using AD-PCR, we followed the manufacturer’s instructions (Applied Biosystems). AD-PCR primers and probes were custom designed from a published sequence using the Custom TaqMan Assay Design Tool (Thermo Fisher Scientific) (23–25).

HLA typing was performed from exome sequencing data, as previously described (5–7). HLA alleles were resolved at the four-field level of resolution from whole-genome sequences and at the three-field level from exomes. All analyses were performed at three-field resolution.

Allele frequencies were calculated either by gene counting or the estimator-maximum (EM) method (26). Haplotype frequencies were computed by the EM algorithm (27). Haplotypes were defined for each functional DRA-DRB1 heterodimer assigned to genotypes by a weighted probability function based on the estimated EM haplotype frequencies (5–7). Descriptive statistics and general linear models for the association of T2D with covariates were calculated using SAS 9.4 software (28). Covariates for the linear models included age, sex, and the first 5 principal components (PCs) derived from a genome-wide association study (PC1–PC5) (29). Each observation in the logistic regression was corrected for sibship using the REPEATED option in the SAS PROC GENMOD procedure that corrects for the familial correlation of the members of the sibship. The statistical power of the interaction term was calculated by a 10,000-iteration bootstrap of the fully parametrized logistic regression, with the power defined by the proportion of iterations <0.05. The sensitivity analysis was calculated from logistic regressions in strata defined by BMI categories <BMI and ≥BMI in BMI range 25–45 kg/m².

MHC binding prediction for HLA-DRB1 was performed using the algorithm NetMHCIIpan 4.0 EL from the Immune Epitope Database and Analysis Resource (30–32). For input, the program uses 15-mer peptides from a FASTA file for a given protein and a vector of class II alleles and calculates for each a probability of strong binding score and a percentile rank (rank%). The rank% compares the peptide’s score against the scores of 5 million random 15-mers selected from the Swiss-Prot database. A small rank% indicates high affinity of the allele with the 15-mer peptide. The program moves over one amino acid at a time and calculates the two parameters for each MHC allele, and then, by this method, continues to the end of the amino acid sequence.

Data and resource availability will be considered following joint guidelines developed with the indigenous community.

In Supplementary Table 1, the EM-estimated frequencies of DRB1 alleles are presented. DRB1*16:02:01(DRB1*02), has an EM-estimated allele frequency of 0.077. The haplotype frequencies for SLC16A11 reveal the segregation of two common specificities, the wild type (C-T-A-C-G, the 1 allele) and the ancient human risk (T-C-G-T-T, the 2 allele), and results in a simple diallelic, additive genetic model (SLC16A11*11, *12, *22). In the all-people stratum, the frequency of the wild-type allele is 0.603 and the risk allele, 0.397.

The allele DRB1*16:02:01, in a dominant model, and the locus SLC16A11 were first incorporated in fully specified logistic models (LMs) with T2D as the dependent variable in three strata: all people, BMI <35 kg/m², and BMI ≥35 kg/m² (Table 1). Furthermore, four logistic regressions were run within each stratum: each locus alone, together, and with interaction (Table 2). For all people, and when included with SLC16A11, DRB1*16:02:01 is significantly associated with T2D and protective, with an odds ratio (OR) of 0.80 (95% CI 0.66, 0.97); when tested alone or with the DRB1 allele, SLC16A11 has ORs not significantly different from 1.0. In the BMI <35 kg/m² stratum, DRB1*16:02:01 has similar ORs to the larger sample, with and without adjustment for SLC16A11 (0.86 and 0.83, respectively) but with marginal statistical significance. However, carriers of the SLC16A11 haplotype have a strong risk for T2D, with and without adjustment for DRB1*16:02:01. Genotype SLC16A11*22 has an OR in both models of ∼1.8 (1.3, 2.3; P < 0.0001), while among heterozygotes, it is ∼1.3 (1.1, 1.6; P ∼0.01). In the interaction model LM8, the ORs for SLC16A11 have a similar magnitude and significance, but the interaction term has an OR not significantly different from 1.0.

In the BMI ≥35 kg/m² stratum, the early human SLC16A11 T-C-G-T-T haplotype switches from risk to strong protection for T2D (Table 2 and Fig. 1), with an OR of 0.58 (95% CI 0.43, 0.78) and 0.78 (0.64, 0.96) for SLC16A11*22 and SLC16A11*12, respectively, in both LM10 and LM11. Allele DRB1*16:02:01 shows T2D protection, without adjustment for SLC16A11 (OR 0.73; 0.56, 0.95) and with adjustment for SLC16A11 (OR 0.73; 0.56, 0.95). The surprising result is model LM12, in which the interaction term between DRB1*16:02:01 and the SLC16A11*12 is statistically significant (OR 0.49; 0.28, 0.86) and the DRB1 allele and the early human SLC16A11 haplotype together reduce the risk of T2D by 51% (Fig. 1 and Supplementary Figs. 1 and 2). In contrast, there is no significant interaction between the DRB1 allele and SLC16A11*22 in LM12.

We designed a different, joint model for the combined protective alleles at DRB1*16:02:01 and SLC16A11 in the BMI ≥35 kg/m² stratum where we observed the protective effect of SLC16A11. People with three or four protective alleles also had an OR of 0.49, with a model significance P < 0.0001 (Supplementary Table 2).

To estimate the power of the interaction term in LM12 (Table 2) with the SLC16A11*12 and DRB1*16:02:01+ we performed a 10,000-iteration bootstrap of the model. The sample size for the BMI ≥35 kg/m² stratum is 2,513. For the interaction term, 4,801 iterations had P < 0.05 (Supplementary Fig. 3). Therefore, the power estimate is 0.4801, or 48%. The distribution of these significant P values has a mean very close to that of the term in the logistic model of ∼0.01. The bootstrap also allows us to estimate the interaction term’s bootstrap-mean or moment for the β and OR with 95% confidence interval (OR 0.43; 0.22, 0.80). These are empirical CIs that are free of any assumptions about the shape of the distribution of the data.

We first chose BMI ≥35 kg/m² as the stratum because the mean in the entire sample is 34.8 and we had used this cutoff in an earlier article (10). However, we decided to perform a sensitivity analysis to further define the role of BMI in the associations and interaction. Such an analysis usually tracks a single event, i.e., risk, whereas we have two events to monitor, i.e., risk and protection. Therefore, we performed three separate sets of logistic regressions with BMI values from 25 to 45 kg/m²: protection (≥BMI strata in Table 3 and Supplementary Table 3) and risk (<BMI strata in Supplementary Table 4). In Table 3, between BMIs of 31 and 39 kg/m², both the SLC16A11*22 and the interaction term have statistical significance with the joint maximum approximately at BMI ≥35 kg/m². We performed a second set of regressions with SLC16A11 genotypes alone in the ≥BMI strata, which maximizes protection in stratum BMI ≥35 kg/m² but exhibits a range of protections in BMI strata 31–41 kg/m² (Supplementary Table 3 and Fig. 2). For the risk analyses, we performed regressions in the <BMI strata where the risk OR for SLC16A11*22 maximizes in the BMI <26 kg/m² stratum (OR 2.95; P = 0.0019) but has statistical significance in the BMI 26–39 kg/m² strata (Supplementary Table 4 and Fig. 3).

To explore potential mechanisms that may explain the observed interactions, we analyzed the predicted binding of HLA class II alleles with SLC16A11. For HLA class II, the seven polymorphic DRB1 alleles in Table 1 were applied to the NetMHCIIpan EL 4.0 algorithm with the 15-mers of the ancient human SLC16A11 protein, which yielded 3,192 records (Supplementary Table 5). The program provides two numbers: the probability of a strong binding score and rank% for each 15-mer peptide. A total of 89.5% of the binding scores were in the interval 0.0–0.1, whereas only 1.4% of the 15-mers had a strong binding >0.4 (Supplementary Table 6). To assess the predicted binding in a larger sample of proteins, we directed the seven polymorphic DRB1 alleles against 113 proteins, which yielded 624,358 records (Supplementary Tables 5 and 6). Predicted strong binding scores >0.3 represented only 2.4% of the peptides. The rank% of binding is based on a library of >5 million random 15-mers in the Swiss-Prot database. The top 2% of binding values are considered as strong binding, 2.0–10% as weak binding, and >10% as what we classify as null binding.

The presence of four missense mutations (V113I, D127G, G340S, and P443T) in the ancient human haplotype of SLC16A11 allows the comparison of the wild-type amino acid and its substitution. Table 4 presents the 15 15-mers at D127G for predicted binding results with DRB1*16:02 (the prediction algorithm allows two levels of allele precision) and contrasts the wild-type SLC16A11 protein with aspartic acid (D) with the ancient human mutation with a substituted glycine (G). For the core peptide FSAFASDLL the DRB1*16:02 exhibits a maximum binding of 0.4717 for peptide 6 that includes the D amino acid, which is within the top 2% of values in Supplementary Table 6. Compared with the ancient human mutation, the effect is amplified. With the G amino acid in the core peptide FSAFASGLL, the maximum binding increases to 0.5650, with a range of ∼0.52–0.57, which is now within the top 1% of scores (Supplementary Tables 5 and 6). In addition, DRB1*16:02 recognizes a consistent core of amino acids in 15-mers 3–8, FSAFASDLL in the wild type and FSAFASGLL in the ancient human haplotype. A search with the UniProt peptide finder for Homo sapiens found that the wild-type peptide is unique in the human genome for SLC16A11. A search for the G mutation returned no match, meaning that the minor allele and its amino acid substitution are not represented in the protein database. Patterns of peptide binding for DRB1*16:02 and the three remaining missense mutations are found in Supplementary Tables 7–9. Supplementary Table 10 presents the 442 15-mer peptides for the ancient human protein when tested against DRB1*16:02 in amino acid order, excluding the leader sequence.

There is increasing recognition of the complexity of the human immune system (33,34). However, a fundamental reaction among the cascades is the presentation of peptide by the class I and class II HLA heterodimers. In the past, many studies of the association of HLA alleles alone have shown a statistical relation, risk or protective, with diseases. We expanded this model to identify and include a potential peptide and its protein source and to search for an epistatic relationship between them.

Since the 1977 Oxford HLA Workshop and Conference, when the serological HLA-DR locus was first elaborated and applied to disease association studies, the antigen HLA-DR2 has been found to be protective for T1D (35). In 1985, Tiwari and Terasaki (36) compiled data on HLA disease associations and reported results for HLA-DR2 from 14 studies from seven countries with a total of 1,268 patients with T1D and 2,768 control individuals. The mean relative risk was 0.18. Antigen HLA-DR2, in later serology, split into HLA-DR15, DR16. In 2008, as part of the Type 1 Diabetes Genetics Consortium, Erlich et al. (37) identified in a sample of individuals primarily of European, Caucasian heritage the protective molecular haplotype from HLA-DR15 as DRB1*15:01-DQA1*01:02-DQB1*06:02, with an OR of 0.03 (95% CI 0.01, 0.07; P = 2 × 10⁻²⁹). People with this haplotype did not develop the disease.

In 2011, we reported a protective effect in this population for the A allele at rs9268852 that tagged HLA-DRB1*02, which we had earlier typed as DRB1*16:02 with an OR of 0.723 (P = 0.002) (1,5). We have now extended this finding of protection for T2D and identify the associated haplotype as DRA*01:01:01-DRB1*16:02:01-DQA1*05:03:01-DQB1*03:01:01 (6). In an earlier report, we also showed the interaction between SLC16A11 and BMI such that in people with BMI ≥35 kg/m², the ancient human haplotype at SLC16A11 was protective for T2D in Southwest American Indians (10). The natural question is, what is the relationship between the two loci? This is a three-part question: Is there epistasis in which the loci interact; if yes, can we identify a 15-mer and core peptide in the protective allele with a strong binding pattern for DRB1*16:02? Furthermore, with the HLA-DR2 protective effect being a property common to T1D and T2D, can we identify a potential common immunological component and link among HLA-DRB1, the two diseases, and SLC16A11?

Our sensitivity analysis shows that there is strong epistasis between HLA-DRB1 and SLC16A11 in the fully parametrized logistic regression with a power of 0.4801 and a bootstrap estimate OR of 0.43. It is not, however, only a property of the BMI ≥35 kg/m² stratum, but the interaction term for DRB1*16:02:01+ and SLC16A11*12 has statistical significance between BMI 30 and 41 kg/m², with mean BMI values of 38.75–47.65 kg/m² (Table 3).

As for the second question, our hypothesis is that this epistatic effect comes from the binding of 15-mer peptides of SLC16A11 by the HLA-DRB1*16:02:01 heterodimer. Evidence comes from the monocarboxylate transporter mutation being both risk and protection for T2D in this population, while the DRB1 allele has consistent protection and strong predicted binding to both the wild-type and ancient human haplotype. In Table 4, the core peptide FSAFASDLL at mutation site D127G has a maximum predicted binding probability of 0.4717 for the basic aspartic amino acid, while for the nonpolar glycine mutation and core peptide FSAFASGLL, the predicted binding score increases to 0.5650. This is the seventh largest binding probability over the entire set of 15-mers from SLC16A11 (Supplementary Table 10). HLA-DRB1*16:02:01 demonstrates protection for the all-people stratum in Table 2, while the increase in predicted binding with the glycine mutation might contribute to its enhancement of protection for the BMI ≥35 kg/m² stratum.

What we have, then, is the observation of epistasis for the two molecules and the observation of strong predicted binding for the DRB1 protective allele; what we lack is the detailed mechanics of risk and protection observed here. They are almost certainly independent, with SLC16A11 having both risk and protection mechanisms deriving from its action as a monocarboxylate transporter. Rusu et al. (14) suggested two mechanisms by which this proton-coupled monocarboxylate transporter’s function could be changed by mutations: 1) The expression of SLC16A11 could be reduced in the liver, and 2) its interaction with basigin (CD147, BSG), a multifunctional protein with two immunoglobulin-like domains that plays a role in the orientation of monocarboxylic acid transporters, could disrupt the expression and/or orientation of SLC16A11 on the cell surface. Further study is needed to make more precise the effects of the mutations. It likely plays a passive role in epistasis with HLA-DRB1*16:02:01 in that it supplies a 15-mer peptide to the molecule when the protein is being turned over in the cell. What immunological cascade is initiated by the binding of peptide and heterodimer that leads to protection is also unknown. When the protective HLA allele and mutant transporter both are present, then the 15-mer will be a self-peptide, and its strong recognition might play a role in maintaining self-recognition and the prevention of autoimmune antibodies that might otherwise contribute to T2D.

With the historic role of protection for HLA-DRB1*15 and T1D, and now HLA-DRB1*16 for T2D, is there a common pattern in their binding for SLC16A11? Table 5 presents the probability of a strong binding score for HLA-DRB1*15:01, the most protective of the alleles in European-derived, Caucasian populations among whom the ancient human mutation haplotype at SLC16A11 is mostly absent. For 15-mers 4–6, the allele has binding scores of 0.2783–0.3526 and recognizes the same core peptide, FSAFASDLL, as does DRB1*16:02:01 for the wild-type protein. The ranks% of binding are all in the top 5% of binding scores compared with the Swiss-Prot database. In addition, DRB1*15:01 has very strong binding for peptides 13 and 14, with scores of 0.3035 and 0.4140, that recognize the core peptide LHLYLGLGL and for which the ranks are also in the top 5% of scores; the aspartic acid at position D127 lies outside the vector of amino acids of this core.

There is virtually no T1D in this Southwest American Indian population, which is consistent with the absence of the risk alleles for the disease in full indigenous American heritage people: HLA-DRB1*04:01, *04:02, *04:05, and *03:01 (37) (Supplementary Table 1). Furthermore, in contrast with HLA-DRB1*15:01 and *16:02, the highest probability of binding for these T1D risk alleles at the D127G missense site, >60 15-mers, is only 0.1139 with a rank% of 11. However, outside of the four missense sites and in the parts of the molecule common to the wild-type and mutant forms of SLC16A11, there is strong predicted binding for the T1D risk alleles (Supplementary Table 11).

With respect to locus SLC16A11, the private ancient human haplotype and risk allele has among the highest reported frequencies in this Southwest American Indian sample, with 0.397 in all people and 0.421 in full indigenous American heritage people. Of more importance than the allele frequency is the proportion of people with at least one early human mutant haplotype of 64%. This community has one of the highest frequencies of T2D and a distribution of SLC16A11 risk haplotypes that blankets two thirds of the population. The SIGMA Type 2 Diabetes Consortium article (8) on the association of SLC16A11 with T2D in Mexico reported a risk haplotype frequency in the entire SIGMA sample of 30%, while within people who were >95% Native American ancestry (n = 290), 48%, although it is not clear whether these are allele or combined genotype frequencies. The 1000 Genomes Project reported a frequency of 0% for the African sample (n = 185), <2% for the European sample (n = 379), 12% for the East Asian sample (n = 286), and 28% for a small sample of Mexican Americans from Los Angeles (n = 66). The frequency of the ancient human risk haplotype in Mexican Americans and admixed Mexicans and Latin Americans is likely derived by genetic admixture from indigenous American populations.

Our observation that enhanced protection from T2D associated with HLA DRB1*16:02:01 and SLC16A11 in the most obese stratum may seem counterintuitive given that obesity increases the risk of T2D. However, the current analyses are based on cross-sectional associations and do not reflect longitudinal risk. Our previous longitudinal analyses in this population suggested that the ancient SLC16A11 haplotype associates with increased risk in the leanest stratum and that the lower prevalence of diabetes in the heavier stratum was partially the result of increased weight loss that occurs after diabetes onset (10). However, in the present analysis, the heavier stratum (BMI ≥35 kg/m²) has a higher percentage of T2D than the leaner one (39.6 vs. 33.8%) (Table 1). The risk-to-protective switch is a key piece of evidence of the mechanics of SLC16A11’s overall function when further epidemiological studies and in vitro or animal experiments address it (38–40).

In summary, the demonstration of epistasis between HLA-DRB1 and SLC16A11 in the amplification of protection in the heavier stratum suggests that the two loci together have a complicated interplay of functions. The ancient human haplotype is a strong risk allele for T2D in the leaner stratum of people, where the combined risk genotypes’ frequency closely mimics the prevalence of the disease in the population. We suggest that the likely missense mutation in this risk is D127G because of the amplification of HLA-DRB1*16:02 binding in protection for the disease, while the independent mechanism of risk for SLC16A11 has yet to be completely elaborated. Absence of T1D in Southwest American Indians, and in full heritage indigenous people throughout the New World, is likely because of the absence of the HLA-DRB1 risk alleles. However, the strong protection of HLA-DRB1*15:01, and its shared binding pattern with DRB1*16:02 and a shared core peptide, suggests that the epistasis might not be exclusive to T2D in protection from diabetes.

Acknowledgments. The authors thank the Southwest American Indian community in this study for their cooperation and participation and thank the staff of the Diabetes Epidemiology and Clinical Research Section, National Institute of Diabetes and Digestive and Kidney Diseases, for conducting the examinations.

Funding. This research was supported in part by the Intramural Research Program of the National Institutes of Health National Institute of Diabetes and Digestive and Kidney Diseases.

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

Author Contributions. R.C.W. designed the study, performed the statistical analyses, and wrote the manuscript. R.L.H. and W.C.K. performed statistical analyses and reviewed and edited the manuscript. B.P. and K.K. contributed to the design and reviewed and edited the manuscript. C.B. and L.J.B. reviewed and edited the manuscript. R.C.W. 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.