Genetic variant effects on gene expression in human pancreatic islets and their implications for T2D

Abstract

Most signals detected by genome-wide association studies map to non-coding sequence and their tissue-specific effects influence transcriptional regulation. However, key tissues and cell-types required for functional inference are absent from large-scale resources. Here we explore the relationship between genetic variants influencing predisposition to type 2 diabetes (T2D) and related glycemic traits, and human pancreatic islet transcription using data from 420 donors. We find: (a) 7741 cis-eQTLs in islets with a replication rate across 44 GTEx tissues between 40% and 73%; (b) marked overlap between islet cis-eQTL signals and active regulatory sequences in islets, with reduced eQTL effect size observed in the stretch enhancers most strongly implicated in GWAS signal location; (c) enrichment of islet cis-eQTL signals with T2D risk variants identified in genome-wide association studies; and (d) colocalization between 47 islet cis-eQTLs and variants influencing T2D or glycemic traits, including DGKB and TCF7L2. Our findings illustrate the advantages of performing functional and regulatory studies in disease relevant tissues.

Fig. 1. Islet eQTLs and their activity in other tissues.

a Proportion of islet eQTLs active in GTEx tissues using p-value enrichment analysis (π₁ estimate for replication). b Comparison between eQTLs discovered in islets and their p-values in β-cells (top figure, N = 26) and whole pancreas tissue from GTEX (bottom figure, n = 149). The axes show the −log₁₀ p-value of the eQTL associations adjusted by the eQTL direction (positive or negative) of effect with respect to the reference allele. Source data are provided as a Source data file.

Fig. 2. Integration of islet eQTL with epigenomic information.

a Distribution of absolute effect sizes for islet eQTLs in each islet chromatin state. b Distribution of absolute effect sizes for islet eQTL in ATAC-seq peaks in three islet chromatin states. eQTL SNPs in ATAC-seq peaks in stretch enhancers have significantly lower effect sizes than SNPs in ATAC-seq peaks in active TSS and typical enhancer states. P-values were obtained from a Wilcoxon rank-sum test. c Fold enrichment for transcription factor footprint motifs to overlap low vs high effect size islet eQTL SNPs. d TF footprint motif directionality fraction vs fold enrichment for the TF footprint motif to overlap islet eQTLs. TF footprint motif directionality fraction is calculated as the fraction of eQTL SNPs overlapping a TF footprint motif, where the base preferred in the motif is associated with increased expression of the eQTL eGene. Significance of skew of this fraction from a null expectation of 0.5 was calculated using the binomial test. Source data are provided as a Source data file.

Fig. 3. GWAS SNPs in islet eQTLs.

a Enrichment of eQTL effect sizes in different GTEx tissues at T2D (all) and glycemic GWAS-associated variants. Numbers within square brackets denote the number of variants implicated for each the trait. Also shown a subset of T2D GWAS associated with reduced insulin secretion or islet β-cell dysfunction (T2D (BC)) and type 1 diabetes (T1D)-associated signals. b LocusCompare plot for the T2D GWAS p-values in the TCF7L2 locus. Plots on the right −log10 p-values for the GWAS (top) and for the the eQTL for TCF7L2, highlighting in both the GWAS lead SNP in the cis window tested for eQTLs. On the left it shows a comparison of the p-values in both analyses. Source data are provided as a Source data file.

Fig. 4. Functional assessment of DGKB eQTL locus.

a We show the two of the three independent islet eQTL signals that colocalize with identified independent GWAS variants near the DGKB gene locus (lead SNP rs17168486 referred at as the 5′ signal and lead SNP rs10231021 referred to as the 3′ signal). These signals colocalize with two independent T2D GWAS signals shown in b, where rs17168486 is referred to as the 5′ signal and lead SNP rs2191349 referred to as the 3′ signal. LD information was not available for SNPs denoted by (×). The third GWAS variant and the third eSNP are not shown as both are located outside this region and in opposite location with respect to DGKB, showing no evidence of colocalization. c Normalized DGKB gene expression levels relative to the T2D-risk-allele dosage at the 3′ islet eQTL for DGKB lead SNP rs10231021. eQTL p-value adjusted to the beta distribution is shown. d Genome browser view of the region highlighted in purple in a and b that contains the 3′ DGKB eQTL and T2D GWAS signals. Two regulatory elements (element 1 highlighted in green, element 2 highlighted in blue) overlapping ATAC-seq peaks in islet β-cells (islet single nuclei ATAC-seq) and bulk islets (islet track represents one islet sample from Varshney et al.) were cloned into a luciferase reporter assay construct for functional validation. All ATAC-seq tracks are normalized to 10 M reads and scaled from 0–15. e Log 2 luciferase assay activities (normalized to empty vector) are shown for in rat (832/13), mouse (MIN6), and human (EndoC-βH1) β-cell lines for the element 2 (cloned in the forward orientation), highlighted in blue in d. The risk haplotype shows significantly higher (p < 0.05) activity than the non-risk haplotype in 832/13 and MIN6, consistent with the eQTL direction shown in c. P-values were determined using unpaired two-sided t-tests. f EMSA for probes with risk and non-risk alleles at the four SNPs overlapping the regulatory element validated in e, using nuclear extract from MIN6 cells. Filled arrows, allele-specific binding; open arrows, non-allele-specific binding of proteins to probes. Source data are provided as a Source data file.