CDKN2A/B T2D Genome-Wide Association Study Risk SNPs Impact Locus Gene Expression and Proliferation in Human Islets

Type 2 diabetes (T2D) risk has a strong genetic component. Significant research investment has identified >100 genomic regions that influence T2D risk in human populations (1–3). Most T2D risk single nucleotide polymorphisms (SNPs) are noncoding, and the mechanism by which they impact local genome biology remains unclear for many loci (3). Risk alleles may act in multiple ways, interacting with other genes and polymorphism effects in a tissue-specific manner. Genome-wide expression quantitative trait loci studies seek to identify how polymorphisms impact biology at any given locus (1,4–7); however, depth of information at individual loci is limited in genome-wide studies. Most T2D SNPs influence risk by impacting islet biology (8), but the cost and inaccessibility of human islets, and poor utility of nonhuman models to study the human genome, have slowed progress in clarifying mechanisms.

SNPs at the CDKN2A/B genomic locus impact risk of T2D and related diseases, such as gestational diabetes mellitus, cystic fibrosis–related diabetes, and posttransplant diabetes, across ethnicities and cultures, suggesting a central diabetogenic mechanism (9). Multiple SNPs in different linkage blocks at the CDKN2A/B locus confer T2D risk (9); mechanisms impacting risk remain unknown. The CDKN2A/B locus encodes four genes (Fig. 1): MTAP, CDKN2A, CDKN2B, and a long noncoding RNA (lncRNA) named ANRIL. CDKN2A and CDKN2B are well studied, encoding cell cycle inhibitors (p14 and p16 are splice variants of CDKN2A, and p15 is encoded at CDKN2B) that impact aging, senescence, and tumorigenesis via regulation of retinoblastoma and p53 (10,11). Three T2D SNPs at this locus, rs10811661, rs2383208, and rs10757283, are noncoding, located downstream of known genes; rs2383208 and rs10811661 are in one linkage block, and rs10757283 is in a separate linkage block immediately downstream. A fourth SNP, rs564398, ∼100,000 base pairs upstream of these, falls within exon 2 of ANRIL. These SNPs were identified in large population studies seeking to identify genomic regions associated with T2D risk (12–14) (more details have previously been published [9]). The three downstream SNPs are mostly associated with T2D risk and not other diseases; the rs564398 SNP is also associated with coronary heart disease and glaucoma (15). The absolute magnitude of T2D risk is low with all identified SNPs (at this and other loci); for example, reported odds ratio for the linkage region containing rs10811661 and rs2383208 ranges from 1.18 to 1.46 (4,12,14,16,17). Weaker odds ratios were seen for rs564398 (1.12–1.26) (example previously published in 13); intriguingly, multiple studies show that rs564398 is associated with T2D risk in Caucasian but not Asian populations (18). A risk-risk combination of rs10811661/rs2383208 and rs10757283 conferred an odds ratio of 1.24, with stronger association than individual risk alleles (13). Although each T2D SNP is in linkage disequilibrium with multiple other SNPs, fine mapping has not identified linked SNPs with greater disease association than these genome-wide association study–identified SNPs (1,4). The causal SNP in any of these linkage blocks is not yet known.

In human populations, the rs10811661 risk allele is associated with reduced insulin secretory capacity after oral or intravenous glucose challenge (16,19–22). Insulin secretory capacity is a composite end point influenced by β-cell mass, insulin production, glucose sensing, and stimulus-secretion coupling (23), factors that cannot currently be effectively separated in living human subjects. Intriguingly, given the aging and senescence roles played by CDKN2A/B genes, the impact of rs10811661 on T2D risk was influenced by subject age (18). SNPs at this locus also influence insulin sensitivity and biology of other metabolic tissues, demonstrating the complexity of even a single genomic locus on T2D biology (9).

Since human studies suggest that CDKN2A/B locus SNPs impact T2D risk, at least in part, by reducing insulin secretory capacity, we hypothesized that locus SNPs influence pancreatic islet biology. Here, we present a detailed analysis of CDKN2A/B biology in nondiabetic human islets. We identified two overlapping coregulated gene sets: p14-p16-ANRIL and p15-MTAP. p14-p16-ANRIL expression, but not p15-MTAP expression, increased with donor age. Of the four T2D risk SNPs tested, rs2383208 and rs10811661 risk alleles were associated with inappropriate high expression of the ANRIL lncRNA in samples from younger donors. No other SNP-gene interaction was identified, but our data suggest certain SNP pairs that may impact locus gene expression in combinatorial fashion. Finally, risk alleles at rs564398 were associated with reduced β-cell proliferation index, suggesting a functional implication for this SNP, and perhaps the ANRIL lncRNA, in accrual or maintenance of human β-cell mass.

Human islets were obtained from the Integrated Islet Distribution Program (IIDP) at the City of Hope, supported by the National Institute of Diabetes and Digestive and Kidney Diseases (NIDDK), National Institutes of Health, or from a collaborative group headed at Vanderbilt University (24). Human islet studies were determined by the University of Massachusetts Institutional Review Board to not qualify for institutional review board review or exemption because they do not involve the use of human subjects. De-identified islet samples from 95 subjects without diabetes were live shipped in Prodo islet transport media. Donors (Supplementary Table 1) included 42 females, 48 males, and 5 without sex reported, with mean ± SD age 40 ± 16 years and ethnicity as follows: 1 Asian, 8 black or African American, 14 Hispanic/Latino, 66 white, and 6 unknown. Upon receipt, islets were plated in islet culture medium (RPMI, 10% FBS, 5 mmol/L glucose, and penicillin/streptomycin) and incubated at 37°, 5% CO₂, overnight for recovery from isolation and shipment. After recovery, 800 islet equivalents (IEQs) were handpicked, washed in PBS containing 100 nmol/L Na₃VO₄, and flash frozen at −80°C in 200-IEQ aliquots for future DNA and RNA analysis. Additional islets from a subset of donors were cultured as described below for glucose-stimulated proliferation.

DNA and RNA were extracted from flash-frozen 200-IEQ aliquots using the Norgen RNA/DNA/Protein Purification Kit (Norgen Biotek Corp. Thorold, Ontario, Canada) according to the manufacturer’s protocol. Genotyping for four CDKN2A/B SNPs (rs564398 [C/T], rs10811661 [C/T], rs2383208 [G/A], and rs10757283 [C/T]) was performed in duplicate with commercial (C_2618017_10, C_31288917_10, C_15789011_10, and C_31288916_10) TaqMan SNP genotyping assays (Thermo Fisher Scientific, Waltham, MA) on the C1000 Touch Thermal Cycler (Bio-Rad) or the RealPlex Cycler (Eppendorf) real-time PCR platforms, using 20 ng DNA in a 10-μL reaction volume under conditions recommended by the manufacturer. SNP determination was confirmed by both allelic discrimination and by manual cycle threshold (Ct) value assessment for all samples and all SNPs. Minor allele frequencies (MAFs) in our cohort were in agreement with expected MAF based on the 1000 Genomes Project (25) (Supplementary Table 2), and the observed frequency of SNP combinations predicted linkage disequilibrium similar to 1000 Genomes–reported values for these SNPs (Supplementary Table 3) (26).

Total RNA was reverse transcribed using a SuperScript IV VILO MasterMix kit (Thermo Fisher Scientific). The expression levels of target genes in human islets were quantitatively assessed in duplicate using Taqman-validated human gene expression assays (Thermo Fisher Scientific). Primers/probes used were as follows: ANRIL, Hs04259476_m1; p15, Hs00793225_m1; p14, Hs99999189_m1; p16, Hs02902543_mH; MTAP, Hs00559618_m1; KI67, Hs01032443_m1; PCNA, Hs00696862_m1; CCND2, Hs00153380_m1; ACTB, Hs01060665_g1; and GAPDH, Hs02758991_g1. ACTB and GAPDH were used as endogenous reference to normalize gene expression. Reproducibility of duplicate measurements was high, as assessed by the R² of the correlation between duplicates and by the absolute value of the relative percentage difference between the duplicates (Supplementary Fig. 1). Transcript expression levels were presented as log2-transformed expression (ΔCt).

Human islets cultured overnight in islet culture medium were dispersed to single cells using single-use–apportioned 0.05% trypsin and plated on uncoated glass coverslips (Fisherbrand) as previously described (27–29). Dispersed cells were cultured in islet culture medium containing either 5 mmol/L or 15 mmol/L glucose for 96 h, with 20 μg/mL BrdU included for the entire time. After culture, the islet cells were fixed for 10 min in 4% paraformaldehyde (Sigma-Aldrich). Immunofluorescence staining was performed after unmasking in 1 N HCl for 25 min at 37°C for insulin (ab7842, Abcam, or A056401-2; Dako), BrdU (ab6326; Abcam), and DAPI as previously described (27–29). β-Cell proliferation, defined as the percent of insulin-staining cells that were also BrdU labeled, was quantified on blinded images (30). Data were expressed as the proliferation index, calculated as the ratio of %BrdU+ β-cells in 15 mmol/L glucose divided by the %BrdU+ β-cells in 5 mmol/L glucose.

Univariate analyses were performed using GraphPad Prism, and data are expressed as mean ± SD. P values were determined by two-tailed Student t test for comparison of two conditions, with F test to compare variances, one-way ANOVA with Tukey posttest for correction for multiple comparisons for comparison of more than two conditions, or by linear regression for assessment of the relationship between two continuous variables. Multivariable linear models were performed to examine gene expressions (p14, p15, p16, ANRIL, and MTAP) simultaneously adjusted for donor sex, race/ethnicity, age (continuous), and BMI (continuous); additional models further adjusted for expression of the other gene products. Missing values were modeled with a missing indicator, replacing unknown values with sample means for linear variables. Infrequent or unknown race/ethnicity was grouped in a residual category. RNA expression associated with SNPs was estimated in linear multivariable models adjusting for demographics (as described above) in two fashions: first, by setting the population with no risk alleles as the common reference group, and second, by estimating the linear effect on a per-allele basis (treating the number of risk alleles as additive). Insulin secretion index was estimated as a function of demographic variables as described above and by including each SNP as a predictor of insulin secretion index. Interpretation of these models can be found in the Supplementary Data. P < 0.05 was considered significant, although this may be too generous for the exploratory analyses with multiple comparisons. For the SNP combination hypothesis-generating analyses, the false discovery rate (FDR), calculated by the original method of Benjamini and Hochberg, was set at 10%, based on our estimation that a hypothesis with 90% likelihood of being correct warranted experimental follow-up.

To understand the context of how T2D risk SNPs might impact biology at this locus in human islets, we first quantified expression of all locus genes (Fig. 1A). Validated Taqman probes were chosen that could independently quantify transcripts including MTAP, p14 (CDKN2A:ARF), p15 (CDKN2B:INK4B), p16 (CDKN2A:INK4A), and the ANRIL (CDKN2B-AS1) lncRNA. p14 and p16 are splice variants of CDKN2A, sharing exons 2 and 3 but with different first exons; exons 2 and 3 are in different reading frames, and thus p14 and p16 encode entirely different proteins with different functions (31). The ANRIL probe spans exons 5 and 6, thus detecting all known isoforms. In this cohort of islet samples from 95 unique donors without diabetes (Supplementary Table 1), RNA abundances of p14, p16, and ANRIL were highly correlated with each other (Fig. 1, normalized to ACTB, and Supplementary Fig. 2, normalized to GAPDH). In contrast, abundance of p15, despite being physically located within the first intron of ANRIL, was poorly (ACTB normalization) or not (GAPDH normalization) correlated with p14, p16, or ANRIL. On the other hand, p15, but not p14, p16, or ANRIL, was highly correlated with MTAP expression. When the data were examined in multivariable linear models, with integration of donor characteristics such as age into the model, again p14-ANRIL and p14-p16 were highly correlated, as were p15-MTAP (Supplementary Table 4). These results suggest two independent but overlapping coregulatory cassettes at the CDKN2A/B locus in human islets, with p14-p16-ANRIL in one and p15-MTAP in the other.

In many tissues, including islets, some CDKN2A/B locus genes increase with advancing age (9,10,32). In this cohort of human islets, expression of p14, p16, and ANRIL showed a modest positive correlation with donor age, whereas p15 and MTAP did not (Fig. 2A–E). Donor BMI could potentially confound the impact of age on gene expression; however, BMI was similar across donor ages (Fig. 2F). Furthermore, we observed no correlation between donor BMI, sex, or ethnicity and expression of any CDKN2A/B locus gene in univariate analysis (Supplementary Figs. 3–5). Multivariable linear models integrating age, sex, race, and BMI confirmed a positive correlation of p14, p16, and ANRIL with donor age and confirmed a lack of impact of sex, race, or BMI on locus gene expression (Supplementary Table 4). Scatterplots of gene expression versus age showed that some genes were expressed in very low abundance in islets from juvenile (age <10 years) donors, with points falling well below the linear regression curve (Fig. 2A, C, and D). Focused analysis of juvenile (<10 years) versus adolescent/adult (>10 years) islets (Fig. 2G–K) revealed that expression of p14, p16 , and ANRIL was markedly lower in juvenile islets, but expression of MTAP and p15 was not, again suggesting altered regulatory characteristics of these two genes relative to other locus genes. Interestingly, an F test showed that the variances were reduced in juvenile islets (see SD bars in Fig. 2G–K) for p14 (P < 0.0001), p16 (P < 0.001), and ANRIL (P < 0.0001) but not for p15 or MTAP (P = NS for both), despite the much smaller sample size, again suggesting fundamentally different biology of the juvenile samples. In sum, older age increased expression of p14, p16, and ANRIL, but not p15 or MTAP, and the youth-associated suppression was exaggerated in islets from very young donors.

We next tested whether T2D-related SNPs at CDKN2A/B impact locus gene expression. Validated Taqman genotyping procedures ascertained and confirmed the genotype of all n = 95 preps for four T2D SNPs: rs564398 (hg38 chr9:22029548), rs2383208 (hg38 chr9:22132077), rs10811661 (hg38 chr9:22134095), and rs10757283 (hg38 chr9:22134173). Measured MAFs (Supplementary Table 2) were similar to reported MAFs for ethnicity-matched populations, supporting genotyping accuracy. SNPs rs2383208 and rs10811661 were tightly linked, with only 2 of 95 samples differing in our cohort, consistent with the linkage disequilibrium reported in LDpair and HaploReg (26,33,34) (Supplementary Tables 3 and 5). In raw analysis across the entire cohort, no SNP genotype correlated with abundance of any locus transcript by univariate (Fig. 3) or multivariable (Supplementary Table 6) analysis. Since donor age impacted expression of p14, p16, and ANRIL, we assessed whether age interfered with the assessment of SNP effect on gene expression. Mean age was not significantly different between genotypes for any SNP (data not shown). However, expression of transcript abundance as a function of age revealed that for ANRIL, but not for p14 or p16, the age-dependent increase was genotype dependent, evident only in protective allele–carrying samples at rs2383208 (Fig. 4) and rs10811661 (Supplementary Fig. 6). Homozygous risk samples had high levels of ANRIL across all ages >10 years (Fig. 4B [flat slope of AA regression line even despite the influence of the juvenile samples]). Age-genotype interaction was not observed for any locus gene for rs10757283 or rs564398 (not shown). When the samples were reanalyzed using a different methodology, binning by quartiles, it was again evident that samples with a protective allele at rs2383208 or rs10811661 showed an age-dependent increase in ANRIL, but homozygous risk samples did not. In contrast, for p16 the slope of the age-dependent–gene expression regression lines (Fig. 4A) and binning analysis (Fig. 4C) were similar across genotypes. For testing of whether ANRIL abundance was inappropriately increased by homozygous risk at rs2383208 or rs10811661 in samples from young donors, samples of those between the ages of 10 years (to exclude juveniles, which were all suppressed independent of genotype) and 50 years (defined by the intersection of the linear regression curves in Fig. 4B), was stratified by genotype (Fig. 4E and F). ANRIL, but not p16, abundance was significantly increased in younger homozygous risk samples compared with protective-allele carriers. Taken together, T2D homozygous risk genotype at rs2383208 or rs10811661 prematurely increased ANRIL expression in islets of younger donors to older-donor levels.

This cohort was not powered for us to perform subgroup analyses to definitively detect gene expression impact of SNP combinations. For example, a sample size analysis using our ANRIL expression mean and SD for rs2383208 genotypes reveals that we would require n = 54 per subgroup to achieve a power of 80% and type 1 error of 0.05. Given the diminishing number of samples as we partition by SNP combination, we do not approach this sample size for subanalyses. Instead, we analyzed our data set using an FDR approach to prioritize hypotheses to test in future studies such as ex vivo promoter-enhancer experiments. We estimated that a risk tolerance of 90% likelihood that a hypothesis was correct would support future experimental investment. We then stratified our gene expression data by all SNP combinations and analyzed each comparison for likelihood of difference, defined by an FDR of 10% (Fig. 5 and Supplementary Fig. 7). By these criteria, we determined that genotype at rs564398 and rs10757283 may influence the impact of rs2383208 and/or rs10811661 on gene expression. For rs564398, in homozygous protective rs564398, but not risk allele–containing samples, protective alleles at rs2383208 and rs10811661 may decrease abundance of p16 compared with homozygous risk carriers. For rs10757283, in homozygous risk rs10757283 samples, but not protective allele samples, protective allele at rs2383208 or rs10811661 may decrease abundance of p15. These observations suggest that individual SNPs may contribute risk by impacting locus biology only in the presence of other locus SNP genotypes, support investment in future experiments to test specific combinations, and help narrow which combinations to target.

A subset of islet samples (n = 61) was tested for insulin secretion stimulation index at their respective islet isolation centers. Insulin secretion index was positively correlated with BMI but showed no relationship with donor age, sex, or islet isolation center by univariate (Supplementary Fig. 8) or multivariable (Supplementary Table 7) analysis. When insulin secretion index was stratified by SNP genotype, contrary to the hypotheses that CDKN2A/B T2D risk SNPs impair glucose sensing, insulin production, or stimulus secretion coupling, samples with T2D risk alleles did not show evidence for impairment in ex vivo insulin secretion in this cohort (Fig. 6 and Supplementary Table 8).

Since CDKN2A/B locus genes are best known for inhibiting the cell cycle, we assessed transcript markers of proliferation (KI67, PCNA, and CCND2) in this cohort, as well as the actual rate of S-phase entry in growth-stimulatory culture conditions, by BrdU labeling, in a subset of samples. Surprisingly, although PCNA and CCND2 showed a high degree of correlation with each other, KI67 did not correlate with either PCNA or CCND2 (Fig. 7A–C). No SNP genotype was correlated with abundance of PCNA, CCND2, or KI67 (not shown). Transcript level is only a surrogate for proliferation and lacks sensitivity in a tissue with a very low frequency of proliferation events and a mixture of cell types. To measure actual cell cycle entry in β-cells, we cultured a subset (n = 47) of islet preparations in 5 mmol/L or 15 mmol/L glucose and quantified nuclear BrdU incorporation in insulin-positive cells (Supplementary Fig. 9). BrdU incorporation rate in 5 mmol/L glucose was nominally correlated with basal PCNA abundance but not with KI67 or CCND2 (Supplementary Fig. 10). As previously observed (27,28,35), glucose increased human β-cell proliferation (P < 0.0001 [not shown]). To test whether any T2D SNP genotype impacted β-cell proliferation, the proliferation index (ratio of BrdU+ β-cells in 15 mmol/L compared with 5 mmol/L glucose) was stratified by SNP identity. Genotype at rs2383208, rs10811661, and rs10757283 did not influence the proliferation index (Fig. 7D–E). However, genotype at rs564398 was strongly associated with the human β-cell proliferation index, with homozygous protective alleles showing approximately doubled stimulation of proliferation by 15 mmol/L glucose compared with islet samples harboring risk alleles at this SNP (Fig. 7F).

We have performed a comprehensive analysis of how one T2D genome-wide association study locus associated with insulin secretory capacity in human populations influences human islet biology. In n = 95 islet samples, we quantified locus gene expression, SNP genotype, donor characteristics, β-cell function (insulin secretion), and β-cell proliferation. We have made several important observations. First, the locus contains two distinct gene cassettes that are physically overlapping but have different regulatory characteristics, with one age dependent and the other not. Juvenile islets have markedly suppressed expression of p14-p16-ANRIL but not p15-MTAP. Second, individual T2D SNPs at the locus do not substantially alter expression of locus genes, but subtle age-dependent influences are detected, which, contrary to expectations, impact the ANRIL lncRNA rather than the most prominent locus product, p16. SNPs may interact with each other to influence gene expression, increasing the complexity of genotype interpretation and raising intriguing mechanistic hypotheses. Third, genotype did not impact insulin secretion index. Finally, risk allele at rs564398, which is located within a transcribed exon of ANRIL, decreased the β-cell proliferation index. This work improves understanding of how CDKN2A/B T2D SNPs impact human islet biology and suggests that the influence of the locus on human insulin secretory capacity may be effected via β-cell mass rather than function.

How SNPs influence T2D risk is an important question in genetics today (3). Although the CDKN2A/B locus is associated with a number of different disease syndromes (36), the SNPs that we selected to study are most strongly related to T2D risk, with the exception of rs564398, which is also associated with coronary disease and glaucoma (15). CDKN2A/B SNPs are associated with a range of diabetes-related syndromes beyond T2D, such as gestational diabetes mellitus and transplant-related diabetes (9). Since this locus is active in many different cell types, and the coronary disease risk region is mostly nonoverlapping with the T2D risk region (36), it is likely that CDKN2A/B impact on diseases other than T2D is mediated by effects outside the islet. Given the association of rs10811661/rs2383208 with impaired insulin production capacity, to identify possible T2D risk mechanisms we performed our study in islets. Our data highlight the complexity of genetic inputs to human metabolism. Even starting with a genomic locus repeatedly associated with disease risk across ethnicities and T2D-related syndromes (9), with in vivo evidence that the islet is the risk-mediating tissue (16,19–21), abundant preclinical locus knowledge in model systems (9,32), and a fairly large sample size of the relevant tissue, we found the impact of risk SNPs to be subtle. As with many studies in human islets, our data illustrate the marked variability from one donor to the next, which reflects the variability of outbred human populations. We incorporated donor parameters such as age and BMI in our analyses but could not measure many potential other confounding premortem influences such as insulin sensitivity (liver, muscle, brain, and fat), coexisting diseases and medications, environmental effects (diet, stresses, and toxins), exercise history, prenatal events, and others. Islet stress related to donor demise, isolation, shipping, and culture may also introduce variability (37). In this context, subtle effects observed in this challenging system may reflect large effects in certain subpopulations or small effects present uniformly across variable conditions.

The presence of two gene expression cassettes at the locus, only one of which is age dependent, suggests interesting biology. p16 is well-known to increase with age in many tissues and organisms (10). Our observation that p15 abundance did not increase with age in human islets conflicts with published results in other tissues (38,39). The relationship between MTAP and CDKN2A/B locus genes has not been studied extensively, but these are contained in the same Hi-C–defined topologically associating domain in the human genome (Supplementary Fig. 11) (40,41). Whether coregulation of p15 and MTAP has functional importance in the islet remains unknown.

rs2383208/rs10811661 impacted ANRIL abundance in age-dependent fashion. Islets containing homozygous risk alleles showed a “premature aging” phenotype, with young risk allele islets having ANRIL levels similar to those of older protective-allele islets. The low MAF for these SNPs precluded comparison between homozygous protective and homozygous risk, which might have revealed a larger effect size. We do not currently know whether one of these, or another SNP in linkage with them, is causative. Fine mapping of this region has not revealed SNPs with greater impact on T2D risk (1,4), but islet ANRIL abundance was not the end point in those studies. These SNPs fall near a known regulatory region downstream of the 3′ end of ANRIL, which may regulate ANRIL transcription. Mechanistically, how higher ANRIL abundance in islets might increase T2D risk is unknown. In other cell types, ANRIL is proproliferative (42), an effect mediated by ANRIL-dependent suppression of locus cell cycle inhibitors p14, p15, and p16 (43). Our RNA analyses did not detect any hint of negative correlation between ANRIL and p14/p16 in islets; in fact, the strong positive correlation between these transcripts calls into question whether ANRIL negatively regulates other locus genes in human islets.

Beyond SNP regulation of ANRIL abundance, a second observation also points to a role for ANRIL in human islets: rs564398, which influences the β-cell proliferation index, is a transcribed polymorphism in this lncRNA. Although the causative SNP remains unknown, it is possible that rs564398 itself impacts lncRNA activity. ANRIL is a complicated gene, with 20 exons and at least 14 reported isoforms (42). Some ANRIL variants are circular (44). Exon 2 is not contained in all isoforms but is generally associated with linear variants (44). Whether rs564398 identity impacts ANRIL isoform production, splicing, stability, or interaction with other cellular DNA, RNA, or protein, to regulate human β-cell proliferation, remains unknown. Genotype at rs564398 did not correlate with expression of cell cycle genes, or BrdU incorporation, under basal conditions; rs564398 may be a marker for β-cell responsiveness to proliferation-inducing conditions rather than increased proliferation in unstimulated conditions. Importantly, there are many SNPs in linkage with rs564398, and any of these, or a combination of these, could be influencing biology instead of rs564398 itself. Also important, while rs564398 has been repeatedly confirmed to be associated with T2D in Caucasian populations, it has little to no relationship with T2D risk in Asian populations (45–47). Taken together with the “premature aging” influence of rs10811661/rs2383208 on ANRIL expression, CDKN2A/B SNPs may impact T2D risk by adversely impacting accrual of β-cell mass during early adulthood.

Combinatorial SNP regulation of gene expression increases the complexity of how genomic variation may impact cellular function. Whole-genome studies restricting locus analyses to single “lead” SNPs cannot detect biology related to two or more local SNPs interacting with each other. The mechanism by which SNPs interact to regulate CDKN2A/B locus gene expression in human islets is unknown. Interaction between rs10757283 and rs2383208 or rs10811661 linked SNPs to regulate p15 abundance could be via modulation of transcription factor occupancy or epigenetic regulation of the enhancer region in which they are located (48). Intriguingly, our observed interaction between these SNPs is supported by a complementary linkage disequilibrium block analysis, which revealed that a haplotype consisting of rs2383208/rs10811661 and rs10757283 was associated with T2D risk (13). There are multiple SNPs in linkage disequilibrium with rs10811661 and also with rs10757283; actual causal polymorphisms are unknown. A physical or functional interaction between the ANRIL lncRNA and this enhancer may mediate cooperation between rs564398 and rs2383208, rs10811661, or linked SNPs to regulate p16 expression. Our study cannot distinguish between in-cis and in-trans interaction. Focused, cell type–specific studies are needed to determine how SNP combinations influence locus gene expression.

This study adds to the body of knowledge debating the relative influence of β-cell mass versus function on T2D risk. Our study found that CDKN2A/B SNPs did not influence glucose-stimulated insulin secretion either in univariate analysis or in multivariable models incorporating donor age, sex, race, and BMI and islet isolation center. This is perhaps surprising given the in vivo data linking risk allele at rs10811661 with impaired insulin secretion (16,22). Since in vivo insulin secretion is a composite outcome of both mass and function, CDKN2A/B SNPs may impact β-cell mass but not β-cell function. This concept is in agreement with the widespread assumption that CDKN2A/B SNPs influence β-cell proliferation because of the known cell cycle–inhibitory effects of locus genes p14, p15, and p16 (49). Our rs564398 data are the first demonstration, to our knowledge, of a CDKN2A/B locus SNP impacting β-cell proliferation.

Our study has caveats. Multiple cell types are found in human islets. Other than the proliferation measurements, which were assessed in insulin-positive cells, all other studies were performed on whole islets. We did not assess the relative proportion of islet cells that were β-cells, and to the extent that gene expression may be cell type–specific, variable cellular makeup may have influenced results. In addition, there is considerable heterogeneity even among β-cells (50,51). Repetition of our current analyses on sorted β-cells, or on single cells, exceeds our current resources. The insulin secretion data have caveats; use of the IIDP-reported insulin secretion index introduces technical variability but benefits from freshly isolated, preshipment tissue. Our confirmation that insulin secretion correlated with BMI, but not with isolation center, is reassuring in this regard.

In sum, this work provides new information about how CDKN2A/B T2D SNPs impact islet biology, suggests that the ANRIL lncRNA may play a role in human islets, and uncovers a link between a T2D SNP and β-cell proliferation. Further studies into the CDKN2A/B locus to develop a mechanistic understanding of how these SNPs impact islet biology to influence T2D risk could one day open the door for using personalized genomic information to inform T2D subtype definitions and therapeutic choice.

Acknowledgments. Human pancreatic islets were provided by the NIDDK-funded IIDP at City of Hope and by Sambra Redick and David Harlan from the University of Massachusetts Medical School and Alvin C. Powers from Vanderbilt University. The authors are grateful to Ahmet Rasim Barutcu, from the Broad Institute and Harvard University, for helpful guidance with the topologically associating domain analysis. The authors thank the Beta Cell Biology Group at the University of Massachusetts Medical School for many helpful discussions.

Funding. This work was supported by grants from NIDDK, National Institutes of Health (R01-DK-095140 to L.C.A.; DK-104211 and DK-106755 to Vanderbilt University group, headed by Al Powers; and 2UC4-DK-098085 to the IIDP). L.C.A. received support from the American Diabetes Association grant 1-15-BS-003 in collaboration with the Order of the Amaranth.

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

Author Contributions. Y.K. performed the majority of the experiments. R.B.S., S.L., and R.E.S. also performed experiments. Y.K. and L.C.A. analyzed data. W.M.J. performed the multivariable linear modeling. L.C.A. devised and planned the experiments and wrote the manuscript. All authors viewed and had the opportunity to edit and approve the manuscript. L.C.A. 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.

Prior Presentation. Parts of this study were presented in abstract form at the 77th Scientific Sessions of the American Diabetes Association, San Diego, CA, 9–13 June 2017.