1. Genetics and Genomics
Download icon

Integration of human pancreatic islet genomic data refines regulatory mechanisms at Type 2 Diabetes susceptibility loci

  1. Matthias Thurner
  2. Martijn van de Bunt
  3. Jason M Torres
  4. Anubha Mahajan
  5. Vibe Nylander
  6. Amanda J Bennett
  7. Kyle J Gaulton
  8. Amy Barrett
  9. Carla Burrows
  10. Christopher G Bell
  11. Robert Lowe
  12. Stephan Beck
  13. Vardhman K Rakyan
  14. Anna L Gloyn
  15. Mark I McCarthy  Is a corresponding author
  1. University of Oxford, United Kingdom
  2. University of California, San Diego, United States
  3. Kings College London, United Kingdom
  4. University of Southampton, United Kingdom
  5. Barts and The London School of Medicine and Dentistry, United Kingdom
  6. University College London, United Kingdom
  7. Churchill Hospital, United Kingdom
Research Article
  • Cited 7
  • Views 2,088
  • Annotations
Cite this article as: eLife 2018;7:e31977 doi: 10.7554/eLife.31977

Abstract

Human genetic studies have emphasised the dominant contribution of pancreatic islet dysfunction to development of Type 2 Diabetes (T2D). However, limited annotation of the islet epigenome has constrained efforts to define the molecular mechanisms mediating the, largely regulatory, signals revealed by Genome-Wide Association Studies (GWAS). We characterised patterns of chromatin accessibility (ATAC-seq, n = 17) and DNA methylation (whole-genome bisulphite sequencing, n = 10) in human islets, generating high-resolution chromatin state maps through integration with established ChIP-seq marks. We found enrichment of GWAS signals for T2D and fasting glucose was concentrated in subsets of islet enhancers characterised by open chromatin and hypomethylation, with the former annotation predominant. At several loci (including CDC123, ADCY5, KLHDC5) the combination of fine-mapping genetic data and chromatin state enrichment maps, supplemented by allelic imbalance in chromatin accessibility pinpointed likely causal variants. The combination of increasingly-precise genetic and islet epigenomic information accelerates definition of causal mechanisms implicated in T2D pathogenesis.

https://doi.org/10.7554/eLife.31977.001

Introduction

T2D is a complex disease characterised by insulin resistance and reduced beta cell function. Recent GWAS have identified a large number of T2D susceptibility loci (Scott et al., 2017; Mahajan et al., 2014; Wellcome Trust Case Control Consortium et al., 2012; Voight et al., 2010), the majority of which affect insulin secretion and beta cell function (Dimas et al., 2014; Wood et al., 2017). However, most GWAS signals map to the non-coding genome and identification of the molecular mechanisms through which non-coding variants exert their effect has proven challenging. Several studies have demonstrated that T2D-associated variants map disproportionately to regulatory elements, particularly those which influence RNA expression and cellular function of human pancreatic islets. (Parker et al., 2013; Pasquali et al., 2014; van de Bunt et al., 2015; Olsson et al., 2014; Dayeh et al., 2014; Volkov et al., 2017; Varshney et al., 2017; Gaulton et al., 2015, 2010).

Characterisation of the islet regulome has until now been limited in scope. The use of DNA methylation and open chromatin data to further annotate ChIP-seq derived chromatin states has successfully uncovered novel biology for other diseases (Wang et al., 2016). Existing methylation studies in islets, however, have either profiled a very small proportion of methylation sites using methylation arrays (Olsson et al., 2014; Dayeh et al., 2014) or focused on T2D-associated disease differentially methylated regions (dDMRs) rather than the integration of DNA methylation status with T2D-relevant GWAS data (Volkov et al., 2017). At the same time, assays of open chromatin in human islets have been restricted to small sample numbers (limiting the potential to capture allelic imbalance in chromatin accessibility for example): these have focussed predominantly on the impact of clustered or ‘stretch’ enhancers (Parker et al., 2013; Pasquali et al., 2014; Gaulton et al., 2010; Varshney et al., 2017).

Most importantly, in part due to historical challenges in accessing human islet material or authentic human cellular models, reference annotations of the islet epigenome and transcriptome (in the context of projects such as GTEx, ENCODE and Epigenome Roadmap) have been largely absent. It is worth noting that islets constitute only ~1% of the pancreas, and islet epigenomes and transcriptomes cannot therefore be reliably assayed in analyses involving the entire organ. Previous islet epigenome studies have, therefore, had only limited ability to directly relate genetic variation to regulatory performance or to broadly characterise the role of DNA methylation in these processes.

In this study, we set out to expand upon previous studies of the islet regulome in several ways. First, we explored the human islet methylome in unprecedented depth using Whole-Genome Bisulphite Sequencing (WGBS) applied to a set of 10 human islet preparations. Second, we explored both basal and genotype-dependent variation in chromatin accessibility through ATAC-seq in 17 human islet samples. Third, we integrated these genome-wide data with existing islet regulatory annotations to generate a high-resolution, epigenome map of this key tissue. Finally, we used this detailed map to interpret GWAS signals for T2D (and the related trait of fasting glucose) and deduce the molecular mechanisms through which some of these loci operate.

Results

Characterising the DNA methylation landscape of human pancreatic islets

To characterise the human islet methylome and characterise the role of DNA methylation with respect to T2D genetic risk, we performed WGBS (mean coverage 13X) in human pancreatic islet DNA samples isolated from 10 non-diabetic cadaveric donors of European descent. Methylation levels across the genome were highly correlated across individual donors (mean CpG methylation Spearman’s rho across 10 individual WGBS donors = 0.71, Figure 1—figure supplement 1A): we pooled the WGBS results to generate a single high-pass (mean coverage 85X) set of pooled human pancreatic islet methylation data covering 23.3 million CpG sites (minimum pooled coverage 10X).

Most previous studies of the relationship between GWAS data and tissue-specific methylation patterns (including those interrogating the relationship between islet methylation and T2D predisposition [Dayeh et al., 2014; Olsson et al., 2014]) had used data generated on the Illumina 450 k methylation array (Hannon et al., 2016; Mitchell et al., 2016; Kato et al., 2015; Ventham et al., 2016). For comparative purposes, we generated 450 k array methylation data from 32 islet samples ascertained from non-diabetic donors of European descent (five overlapping with those from whom WGBS data were generated). As with the WGBS data, methylation levels were highly correlated across individuals (mean CpG methylation Spearman’s rho across 32 individual 450 k donor = 0.98, Figure 1—figure supplement 1B). After pooling 450 k array data across samples, methylation profiles generated from the 450 k array and WGBS were highly correlated at the small subset of total CpG sites for which they overlap: this was observed across pooled samples (pooled WGBS vs. 450 k Spearman’s rho = 0.89, Figure 1A) and across the five donors analysed by both methods (mean Spearman’s rho = 0.80, not shown).

Figure 1 with 2 supplements see all
Comparison of human pancreatic islet WGBS and 450 k methylation data across the genome.

(A) Smooth Scatter plot shows Spearman’s rho correlation between the 450 k array (y-axis) and WGBS (x-axis) at overlapping sites. Darker colour indicates higher density of sites. (B) Comparison of the 450 k array (orange) and WGBS (yellow) methylation levels (x-axis) of all CpGs genome-wide assayed by either method (y-axis shows density). The P-value shown is derived using a Kolmogorov-Smirnov (KS) test. (C) For each chromatin state from Parker et al. (2013) the methylation levels of all CpG sites independent of overlap (diamond indicates the median) are shown as violin plots (left y-axis) and the CpG probe percentage per state for the 450 k array (orange) and WGBS (yellow) are shown as bar-plot (right y-axis). The 450 k probes represent the percentage of the total number of CpG sites which is determined by the number of WGBS CpG sites detected (WGBS = 100%). (D) Distribution of GWAS Posterior Probabilities (Type 2 Diabetes and Fasting Glucose) captured by CpG sites on the 450 k array (orange), 850 k array (green) and WGBS (yellow/black line). (E) Locuszoom plot showing CpG density and credible set SNPs. SNPs are shown with P-values (dots, y-axis left), recombination rate (line, y-axis right) and chromosome positions (x-axis) while CpG and gene annotations are shown below. These annotations include CpGs identified from WGBS (yellow stripes), 450 k CpG probes (orange stripes), 850 k CpG probes (green stripes) and gene overlap (DGKB label). The highlighted region in blue captures the 99% credible set region plus additional 1000 bp on either side. At the very bottom the position on chromosome seven is shown in Megabases (Mb).

https://doi.org/10.7554/eLife.31977.002

WGBS and 450 k array data differed substantially in terms of genome-wide coverage. The 450 k array was designed to interrogate with high precision and coverage ~480 k CpG sites (approximately 2% of all sites in the genome), selected primarily because they are located near gene promoters and CpG-island regions. The focus of the 450 k array on these regions, which tend to be less variable in terms of methylation, explains the high 450 k array correlation levels between donors. In addition, this selective design results in marked differences in the distributions of genome-wide methylation values between WGBS and the 450 k array. Whilst the WGBS data revealed the expected pattern of widespread high methylation levels with hypomethylation (<50%) restricted to 11.2% (2.6M/23.3M CpG sites) of the genome, the array disproportionately interrogated those hypomethylated sites (218 k [46%] of all 450 k CpG probes) (Kolmogorov–Smirnov (KS) test for difference, D = 0.40, p<2.2×10−16) (Figure 1B). These differences in methylation distribution were also evident within specific islet regulatory elements from previously defined standard chromatin state maps (Parker et al., 2013) (Figure 1C, Figure 1—figure supplement 1C–D). We found significant (FDR < 0.05) differences between the methylation levels of CpG sites accessed on the array, and those interrogated by WGBS, across most islet chromatin states: the largest differences were observed for weak promoters (median WGBS = 0.71 vs. median 450k = 0.11, KS test D = 0.51, p<2.2×10−16,) and weak enhancers (WGBS = 0.87 vs. 450 k median = 0.76, D = 0.39, p<2.2×10−16, Figure 1—figure supplement 1D).

In terms of coverage, most chromatin states, apart from promoters, were poorly represented by CpG sites directly interrogated by the array: for example the array assayed only ~2.9% of CpG sites in strong enhancer states (2.7–3.8% depending on strong enhancer subtype, Figure 1C). Although methylation levels were previously reported to be highly correlated across short (0.1–2 kb) genomic distances (Zhang et al., 2015; Bell et al., 2011; Eckhardt et al., 2006; Guo et al., 2017), the observed significant differences in the methylation distribution (Figure 1C, Figure 1—figure supplement 1D) across chromatin states including weak promoter (median size 600 bp) and enhancer subtypes (median size ranges from 200 to 1200 bp) indicate that these correlative effects are not strong enough to counterbalance the low coverage of the 450 k array. These findings are consistent with 450 k array content being focused towards CpG-dense hypomethylated and permissive promoter regions. This highlights the limited capacity of the array to comprehensively interrogate the global DNA methylome, in particular at distal regulatory states such as enhancers.

To understand the value of these data to reveal molecular mechanisms at GWAS loci, where we and others had shown enrichment for islet enhancer states (Pasquali et al., 2014; Gaulton et al., 2015; Parker et al., 2013), we were interested to see how the selective coverage of the array might impact on its ability to interrogate methylation in GWAS-identified regions. We used the largest currently available T2D DIAGRAM GWAS data set (involving 26.7 k cases and 132.5 k controls of predominantly European origin, see dataset section for details) to identify the ‘credible sets’ of variants at each locus which collectively account for 99% of the posterior probability of association (PPA) (Scott et al., 2017; Maller et al., 2012).

To estimate the respective proportions of these T2D-associated variants captured by CpG sites assayed by the 450 k array and WGBS, we determined, for each locus, the combined PPA of all 99% credible set variants mapping within 1000 bp of any CpG site captured. This is based on evidence that across short distances CpG methylation is highly correlated (Zhang et al., 2015; Bell et al., 2011; Eckhardt et al., 2006; Guo et al., 2017) and may be influenced by genetic variants associated with altered transcription factor binding (Do et al., 2016). We found that coverage of this space of putative T2D GWAS variants by the 450 k array is low: across GWAS loci, the combined PPA attributable to variants within regions assayed by the array ranged from 0–99% with a median PPA per locus of 16% (compared to a WGBS median PPA per locus = 99%, KS-test p<2.2×10−16, Figure 1D, top). We estimated that the equivalent figure for a recently developed upgrade of the 450 k array, which captures ~850 k CpG sites and aims to provide better coverage of enhancer regions, would be ~39% (range 0%–99% Figure 1D, top). For instance, at the DGKB T2D locus (centred on rs10276674), CpG sites covered by the 450 k array interrogated less than 1% of the PPA of associated variants (vs. 99% captured by WGBS); the figure for the 850 k array would be 23% (Figure 1E). We obtained similar results when we performed equivalent analyses using GWAS data for fasting glucose (FG, from the ENGAGE consortium [Horikoshi et al., 2015]), another phenotype dominated by islet dysfunction (Figure 1D, bottom).

These data indicate that available methylation arrays provide poor genome-wide coverage of methylation status and are notably deficient in capturing methylation status around the distal regulatory enhancer regions most relevant to T2D predisposition. For this reason, we focused subsequent analyses on the WGBS data.

Integration islet methylation and other epigenomic annotations

Studies in a variety of other tissues have shown that hypomethylation is a strong indicator of regulatory function (Stadler et al., 2011). More specifically, continuous stretches of CpG-poor Low-Methylated Regions (LMRs, with methylation ranging from 10–50% and containing fewer than 30 CpG sites) denote potential distal regulatory elements such as enhancers, while stretches of CpG-rich UnMethylated Regions (UMRs, containing more than 30 CpG sites) are more likely to represent proximal regulatory elements including promoters (Burger et al., 2013). We detected 37.1 k LMRs, 13.6 k UMRs (Figure 2A) and 10.7 k Partially Methylated Domains (PMDs) (Materials and methods and Figure 2—figure supplement 1A–B). PMDs represent large regions of unordered methylation states associated with DNA sequence features (Gaidatzis et al., 2014). As anticipated, we found significant enrichment of LMRs with weak and strong enhancer states as defined by islet chromatin state maps derived from existing ChIP-seq data (Parker et al., 2013) (69.2% of islet LMRs overlapped islet strong and weak enhancer states, log2FE = 2.2–2.9, Bonferroni p<0.05, Figure 2B, Figure 1—figure supplement 1C). Similarly, UMRs were enriched for islet active promoter chromatin states (90.8% of UMRs overlapped islet active promoters, log2FE = 3.9, Bonferroni p < 0.05, Figure 2B).

Figure 2 with 1 supplement see all
Overlap of WGBS hypomethylation and ATAC-seq open chromatin peaks with regulatory annotation.

(A) Methylation levels in percent (y-axis) and log2 CpG density (x-axis) of UMR and LMR regulatory regions with the dashed line indicating the CpG-number (30 CpGs) that distinguishes LMRs and UMRs. (B) Log2 Fold Enrichment (log2FE) of LMRs (green shape), UMRs (blue shape) in various islet annotations is shown. These annotations include islet chromatin states, islet relevant TFBS (FOXA2, MAFB, NKX2.2, NKX6.1, PDX1), islet eQTLs, WGBS derived T2D-associated islet disease DMRs (dDMRs) and ATAC-seq open chromatin peaks. The dDMRs were derived from 6 T2D and eight non-diabetic individuals by Volkov et al. (2017) and dDMRs (orange shape) were also tested for enrichment in the aforementioned islet regulatory annotations. For all annotations, the empirically determined Bonferroni adjusted P-value is ≤0.00032 unless otherwise indicated by the shape: a dot corresponds to an Bonferroni adjusted p-value<0.00032 while the three triangles indicates Bonferroni adjusted p-values>0.00032: UMR enrichment adjusted P-value for weak enhancers = 1; dDMR enrichment adjusted P-value for MAFB = 0.006 and dDMR enrichment adjusted P-value for islet eQTLs = 0.01.

https://doi.org/10.7554/eLife.31977.005

To further characterise these hypomethylation domains, we overlapped information from analyses of islet cis-expression QTLs (eQTLs) (van de Bunt et al., 2015) and islet ChIP-seq transcription factor binding sites (TFBS) (Pasquali et al., 2014). We observed marked enrichment for eQTLs (LMR log2FE = 1.1, UMR log2FE = 2.7, Bonferroni p<0.05) and TFBS (LMR log2FE = 4.1–4.6; UMR log2FE = 2.4–3.9, Bonferroni p<0.05, Figure 2B). These observations confirm that islet LMRs and UMRs correspond to important tissue-specific regulatory regions, overlapping cis-regulatory annotations known to be enriched for T2D GWAS signals (Pasquali et al., 2014; Gaulton et al., 2015).

We also considered the relationship between LMR and UMR regions defined in our non-diabetic islet WGBS, and a complementary set of methylation-based annotations previously derived from WGBS of islets from 6 T2D and 8 control individuals (Volkov et al., 2017). In that study, comparisons between diabetic and non-diabetic islets had been used to define a set of 25,820 ‘disease differentially methylated regions’ (dDMRs, minimum absolute methylation difference 5% and p<0.02). We found only limited overlap between these dDMRs and the UMRs and LMRs from our data: of the 25,820 dDMRs, 2.2% overlapped LMRs and 2.4% UMRs. This overlap was slightly greater than expected by chance (Bonferroni p<0.05, LMR log2FE = 1.0 and promoter-like UMRs log2FE = 1.1, Figure 2B) but more modest than seen for the other regulatory annotations. Similarly, we also observed that dDMRs showed more modest (log2FE = 0.4–1.0), but still significant (Bonferroni p<0.05) levels of enrichment with respect to all other islet regulatory annotations (Figure 2B). The modest enrichment of dDMRs indicates that only a fraction of these regions correspond to islet genomic regulatory sites. Given that T2D risk variants preferentially map in islet regulatory sites, the corollary is that most dDMRs are unlikely to directly contribute to the mediation of genetic T2D risk.

Refining islet enhancer function using methylation and open chromatin data

To further characterise the regulatory potential of hypomethylated regions, including LMRs and UMRs, we combined the islet WGBS methylation data with chromatin accessibility data generated from ATAC-seq assays of 17 human islet samples (from non-diabetic donors of European descent; mean read count after filtering = 130M, Figure 2—figure supplement 1C). We identified a total of 141 k open chromatin regions based on read depth, peak width and signal-to-noise ratio (see Materials and methods). These regions of islet open chromatin showed substantial overlap (78%) with equivalent regions described in a recent study of two human islets (Varshney et al., 2017) (log2FE = 2.8 compared to random sites, not shown). In addition, our islet ATAC-Seq sites demonstrated substantial overlap with LMRs: 53% of LMRs overlapped 16% of all ATAC-seq peaks (LMR log2FE = 3.8 compared to randomised sites, Figure 2B). Almost all UMRs (98%) were contained within regions overlapping (13% of) ATAC-seq peaks (UMR log2FE = 3.4 compared to randomised sites, Figure 2B).

To fully leverage information across multiple overlapping islet epigenome assays, we generated augmented chromatin state maps, using chromHMM (Ernst and Kellis, 2012). These maps combined the WGBS methylation and ATAC-Seq open chromatin data with previously generated ChIP-seq marks (Figure 3A, Figure 3—figure supplement 1A). For these analyses, we initially used a single definition for hypomethylated regions (methylation <60%) that captured both UMRs and LMRs (see Materials and methods).

Figure 3 with 2 supplements see all
Integration of islet epigenetic data to refine chromatin regulatory states and enrichment of these states in T2D GWAS data.

(A) 15 chromatin states (y-axis) were derived from ChIP histone marks, DNA methylation and ATAC-seq open chromatin annotations (x-axis) using chromHMM. For each state the relevant marks characterising the state are shown. The colour is based on the chromHMM emission parameters and a darker colour indicates a higher frequency of a mark at a given state. Weak enhancers (marked by H3K4me1 alone, red) and strong enhancers (marked by H3K27ac and H3K4me1, green) were subdivided by the chromHMM analysis according to methylation and ATAC-seq status (highlighted in red and green box). The black bar at the x-axis highlights the most important marks for characterising enhancer subtypes. (B-C) FGWAS Log2 Fold Enrichment including 95% CI (log2FE, x-axis) of all chromatin states (y-axis) in T2D GWAS regions is shown which demonstrate differential enrichment amongst enhancer subclasses in single-feature enrichment analysis. In addition, log2FE of Coding Sequence (CDS) and Conserved Sequence (CONS) annotations are shown to include the effect of protein-coding and conserved regions. Significantly enriched annotations are shown in black while non-siginificant annotations are shown in grey. (C) T2D FGWAS maximum likelihood model determined through cross-validation. Log2FE and 95% CI (x-axis) of annotations included in the maximum likelihood model (y-axis) also demonstrate differential enrichment amongst enhancer subclasses. *Analysis for Genic Enhancers (state 10) did not converge and hence, only a point log2FE estimate is provided. (D) Single feature log2FE including 95% CI (x-axis) results are shown highlighting the differences in T2D GWAS enrichment of various annotations. These include ATAC-seq open chromatin peaks (red), WGBS methylation regions (including enhancer-like LMRs, promoter-like UMRs and Partially Methylated Domains, blue), ChIP-seq chromatin states (orange) and CDS and CONS (green) annotations. (E) Chi-square distribution (curved black line) with the indicated results of a maximum likelihood ratio test based on the maximum likelihood difference between a model including LMRs or ATAC-seq peaks compared to the ChIP-only model. The dashed red line indicates significance (p-value<0.05). For all FGWAS enrichment plots the axis has been truncated at −6 to facilitate visualisation and accurate values are provided in the supplementary tables.

https://doi.org/10.7554/eLife.31977.008

This augmented and larger set of 15 islet chromatin states retained the broad classification of regulatory elements that included promoters (positive for H3K4me3), transcribed and genic regions (H3K36me3), strong enhancers (H3K4me1; H3K27ac), weak enhancers (H3K4me1), insulators (CTCF) and repressed elements (H3K27me) (Figure 3A). The addition of islet methylation and open chromatin data expanded existing chromatin state definitions to provide new subclasses, particularly amongst enhancer elements. Here, we observed two subclasses of strong enhancers and three of weak enhancers (Figure 3A). We denote the strong enhancer subtypes as ‘open’ (n = 32 k genome-wide), characterised by open chromatin and hypomethylation, and ‘closed’ (n = 110 k) with closed chromatin and hypermethylation (Figure 3A). The three weak enhancer states we denote as ‘open’ (n = 38k: open chromatin, hypomethylation), ‘lowly-methylated’ (n = 78 k; closed chromatin, hypomethylation) and ‘closed’ (n = 206k: closed chromatin, hypermethylation). No equivalent class of ‘lowly-methylated’ strong enhancers was observed in the 15-state model. When comparing these chromatin states to those identified using only ChIP-seq marks ([Parker et al., 2013], Figure 1—figure supplement 1C), the two strong enhancer subclasses we identified subdivided the ‘strong enhancer 1’ state as described by Parker (defined by H3K27ac and H3K4me1). Additional comparison to ‘stretch’ enhancer clusters (Parker et al., 2013), showed that there was considerable overlap between the ‘open’ strong and weak enhancer states we identify here and previously-described ‘stretch’ enhancer states (16.1 k out of 23 k stretch enhancer overlapped 32 k out of 70.1 k ‘open’ enhancers). Even so, most (55%) ‘open’ enhancer states, and in particular ‘open weak enhancers’ (70%), were not captured within ‘stretch’ enhancer intervals, and we regard these as distinct islet enhancer subclasses.

To understand the relationship of these various state definitions to genetic variants influencing T2D risk, we applied the hierarchical modelling approach FGWAS (Pickrell, 2014) to the same sets of large-scale GWAS data for T2D (from DIAGRAM [Scott et al., 2017]) and FG (ENGAGE [Horikoshi et al., 2015]) described in section 2.1. FGWAS allowed us to combine GWAS and genomic data to determine the genome-wide enrichment within islet regulatory features for variants associated with T2D risk. These enrichment priors were then used to generate credible variant sets that are informed by both GWAS and genomic data, as described in section 2.4.

In single-feature analyses, we found significant enrichment (lower limit of Confidence Interval (CI) >0) limited to four enhancer states (open weak enhancers, both types of strong enhancer and H3K36me3 marked genic enhancers) (Figure 3B, Table 1). To take into account protein-coding variant and conserved sequence effects, we also included CoDing exon Sequence (CDS) (Carlson and Maintainer, 2015) and CONServed sequence (CONS) (Lindblad-Toh et al., 2011) as additional annotations which were previously found to be strongly enriched for T2D GWAS signal (Finucane et al., 2015). We observed significant enrichment for CDS and CONS sequence in the single state results (Figure 3B, Table 1). FGWAS multi-feature analyses for T2D, incorporating all annotations positive in single-element analyses, retained both subclasses of strong enhancer, the subclass of open weak enhancers, genic enhancers and CDS in the joint model (Figure 3C and Materials and methods). Conserved sequence annotations were not retained in the joint model.

Table 1
Single FGWAS annotation enrichment in T2D and FG GWAS data.

For each chromatin state annotation the total number of sites and the single state FGWAS log2 Fold Enrichment (log2FE) in T2D and FG is shown. In addition, log2FE enrichment is also shown for CDS and CONS annotation. 95% Confidence Intervals (CI) for log2FE are shown in brackets and significantly enriched states are highlighted in bold (lower CI limit >0). Parker enhancer states refer to enhancer states defined by Parker et al., 2013.

https://doi.org/10.7554/eLife.31977.016
Chromatin
States
Total number of statesT2D log2FE
(CI)
FG log2FE
(CI)
1. Active Promoter20 k 1.6 (-0.8 to 2.7) 2.7 (0 to 4.1)
2. Weak Promoter33 k 1.7 (-4.8 to 2.9) 2.7 (-0.1 to 4.2)
3.Transcriptional Elongation71 k-0.4 (-20 to 1.1)-26.1 (-46.1 to 1.0)
4. Low Methylation73 k-1.5 (-3.1 to -0.6)-1.7 (-4.2 to -0.3)
5. Closed Weak Enhancer206 k 1.2 (-0.1 to 2) 1.7 (0 to 2.9)
6. Lowly-methylated Weak Enhancer78 k-0.5 (-20 to 1.6)-26.7 (-46.7 to 1.6)
7. Open Weak Enhancer38 k 3.4 (2.5 to 4.2) 3.1 (-0.6 to 4.6)
8. Closed Strong Enhancer110 k 2.7 (1.8 to 3.4) 3.3 (2 to 4.4)
9. Open Strong Enhancer32 k 3.8 (3.1 to 4.5) 4.3 (2.8 to 5.5)
10. Genic Enhancer39 k 2.5 (1.3 to 3.4) 2.9 (0.8 to 4.3)
11. Accessible chromatin14 k-25.2 (-45.2 to 2.5)-28.4 (-48.4 to 3.7)
12. Insulator31 k 0.9 (-20 to 2.6)-0.6 (-20 to 3.6)
13. Heterochromatin216 k 2.3 (-20 to 3.9) 1.8 (-1.5 to 4.0)
14. Polycomb Repressed71 k-25.5 (-45.5 to 0.9)-33.2 (-53.2 to 1.5)
15. Quiescent State1.7 k-1.0 (-2.2 to -0.1)-28.6 (-48.6 to -0.6)
CDSNA 2.6 (1.2 to 3.5) 2.7 (-0.2 to 4.3)
CONSNA 2.1 (1.1 to 2.9) 1.9 (0.2 to 3.2)
Parker Weak Enhancer119 k 0.9 (-2.5 to 2.0)-2.0 (-20.0 to 2.4)
Parker Strong Enhancer (all)123 k 2.7 (2.0 to 3.3) 3.1 (2.0 to 4.4)
Parker Strong Enhancer (open)64 k 3.1 (2.4 to 3.7) 3.6 (2.3 to 4.8)
Parker Strong Enhancer (closed)59 k 1.9 (0.8 to 2.7) 2.3 (0.5 to 3.5)

We observed markedly different levels of enrichment for T2D association between and within open and closed enhancer states (Figure 3B–3C, Table 1). Using these augmented chromatin state maps, we demonstrated clear enrichment for T2D association for the subset of ‘open’ weak enhancers (12% of all weak enhancer sites) with no evidence of enrichment in the remaining subclasses (‘closed’ and ‘lowly-methylated’) (Figure 3B and Table 1). This concentration of enrichment amongst a relatively small subset of the weak enhancers was consistent with the lack of enrichment across all weak enhancers defined solely on the basis of H3K4me1 signal ([Parker et al., 2013], single state log2FE = 0.9, CI = −2.5 to 2.0, Table 1, Figure 1—figure supplement 1C). We also saw differences in enrichment signal between open and closed strong enhancers, with the most marked enrichment amongst open strong enhancers (22% of the total, Figure 3B–C, Table 1). This effect was particularly obvious in the joint-analysis (open strong enhancer joint log2FE = 4.1, CI = 3.3 to 4.8 vs. closed strong enhancer joint log2FE = 2.4, CI = 0.5 to 3.3, Figure 3C).

Hypomethylation and open chromatin are highly correlated, but the observed difference in T2D enrichment between the weak enhancer states (particularly between ‘lowly-methylated’ and ‘open’ which differ markedly with respect to chromatin status) points to a primary role for open chromatin. To test this further, we regenerated chromatin state maps using different subsets of the data (ChIP-only, with optional addition of methylation and/or open chromatin information, see Materials and methods and Figure 3-figure supplement 1A-B). These analyses confirmed that the T2D GWAS enrichment signal was predominantly driven by the distribution of islet open chromatin (Figure 3—figure supplement 1C).

We further evaluated the role of subclasses of DNA methylation regulatory region with respect to T2D GWAS enrichment. We divided hypomethylated (<60% methylated) sequence into enhancer-like LMRs (6.5% of all hypomethylated sequence), promoter-like UMRs (7.5% of hypomethylated sequence), as well as PMDs (61% of hypomethylated sequence). The remaining 25% of hypomethylated sequence did not fit any category. LMRs were significantly (lower CI limit >0) enriched (log2FE = 3.2, CI = 2.3 to 3.9) for T2D association signals consistent with their co-localisation with distal regulatory elements, and displayed modestly increased enrichment compared to enhancer states derived from ChIP-seq alone (Figure 3D, Figure 3—source data 1). In contrast, no significant enrichment was found for human islet (promoter-like) UMRs (log2FE = 1.4, CI = −0.6 to 2.5) or PMDs (log2FE = −0.8, CI = −1.7 to −0.1). We also found no evidence that recently-described regions of T2D-associated differential methylation (dDMRs: derived from comparison of WGBS data from islets of diabetic and non-diabetic individuals) were enriched for genome-wide T2D association signals (log2FE = −24.6, CI = −44.6 to 3.7) (Figure 3D, Figure 3—source data 1).

Finally, since the hypomethylation signal for T2D enrichment was concentrated in LMRs (Figure 3D, Figure 3—source data 1), we reran a FGWAS joint-analysis combining open chromatin peaks, LMRs and ChIP-only states using a nested model (Figure 3E, Figure 3—figure supplement 1D–E, see Materials and methods). This confirmed that the improvement in enrichment was mainly driven by open chromatin but showed that LMRs also contributed significantly and independently to the enrichment (Figure 3E, Figure 3—source data 2).

FGWAS analysis for FG corroborated the observations from T2D analysis. Despite reduced power of the FG GWAS data due to a lower number of significantly associated FG GWAS loci, both single feature and joint-model analyses of human islet epigenome data found significant enrichment in strong enhancer states with the strongest enrichment in enhancers with open chromatin and hypomethylation (Figure 3—figure supplement 2A–B and Table 1). In addition, evaluation of the relative contributions of ATAC-seq open chromatin and DNA methylation to FG GWAS enrichment across both single-feature (Figure 3—figure supplement 2C–D) and joint-model analysis (Figure 3—figure supplement 2E–F and Figure 3—source data 3) indicated that open chromatin was primarily responsible for the enhanced enrichment.

Overall, these analyses demonstrate that the addition of open chromatin and DNA methylation data to ChIP-seq marks enhances the resolution of regulatory annotation for human islets. In particular, it defines subsets of weak and strong enhancers that differ markedly with respect to the impact of genetic variation on T2D risk. Although DNA accessibility and hypomethylation status are strongly correlated and provide broadly similar enrichments, the effects of the former predominate. In line with the dominance of open chromatin status for T2D GWAS enrichment, we observed that T2D risk in relation to methylation status is primarily invested in hypomethylated LMRs (i.e. enhancers) rather than UMRs, dDMRs or PMDs.

Augmented chromatin maps and open chromatin allelic imbalance refine likely causal variants at ADCY5, CDC123, and KLHDC5

We next deployed the insights from the global FGWAS enrichment analyses to define the molecular mechanisms at individual T2D susceptibility loci, refining T2D causal variant localisation using the combination of genetic data (from fine-mapping) and the genome-wide patterns of epigenomic enrichment.

Specifically, we applied FGWAS to the T2D DIAGRAM GWAS data (Scott et al., 2017) under the joint model (Figure 3C) derived from the augmented chromatin state maps. We divided the genome into 2327 segments (average size 5004 SNPs or 1.2 Mb) and identified 52 segments significantly associated with T2D genome-wide (segmental FGWAS PPA >= 0.9 or single variant GWAS p<5×10−8, see Materials and methods for details). These corresponded to 49 known T2D associated regions representing that subset of the ~120 known T2D GWAS loci which passed those significance/filtering criteria in this European-only dataset. We then calculated reweighted PPAs for each variant within each segment and generated reweighted 99% credible sets. (Of note, in line with traditional GWAS nomenclature, locus names were defined based on proximity between the lead variant and the closest gene and does not, of itself, indicate any causal role for the gene in T2D susceptibility).

Consistent with the increased T2D GWAS enrichment of states including open chromatin and DNA methylation information, we found that analyses using enrichments from the augmented chromatin state model (combining ChIP-seq, ATAC-seq and WGBS data) were associated with smaller 99% credible sets (median of 17 SNPs) than those derived from FGWAS enrichment derived from ChIP-seq data alone (median 23). In parallel, the PPA for the best variant per locus increased (median 0.39 vs 0.31). Individual T2D GWAS locus results are shown in Figure 4A–B. We also expanded the FGWAS PPA analysis to investigate open chromatin and DNA methylation effects on fine-mapping and found that the reduction in 99% credible set size and increase in maximum variant PPA was driven predominantly by open chromatin (Figure 4—figure supplement 1, Figure 4—source data 1). This demonstrates that the inclusion of open chromatin maps helps to improve prioritisation of causal variants at many T2D GWAS loci.

Figure 4 with 1 supplement see all
Evaluating Posterior Probabilities (PPA) derived from the FGWAS maximum likelihood model at significant T2D GWAS loci.

(A) Per locus the difference in the number of 99% credible set variants between ChIP +ATAC + Meth and ChIP-only model is shown (positive values indicate a reduction in the number of 99% credible set variants in the ChIP+ATAC+Meth model). (B) Per locus the difference in the maximum single variant PPA between the ChIP +ATAC + Meth and ChIP-only model is shown (positive values indicate an increase in the maximum single variant PPA in the ChIP +ATAC + Meth model). (C) T2D GWAS loci were classified into insulin secretion (ISR), insulin resistance (IR) or unclassified loci based on genetic association with physiological traits derived from Dimas et al. (2014) and Wood et al. (2017). In addition, loci with known role in islet genomic regulation or function are highlighted in bold. These include loci with islet eQTLs (ZMIZ1, CDC123) and mQTLs (WFS1, KCNJ11). (D) Identification of T2D GWAS loci and variants enriched for enhancer chromatin states using FGWAS PPA. Per locus the highest PPA variant is shown (y-axis) and the number of variants with PPA >0.01 (x-axis). Loci with high PPA variants (min PPA >0.1, dashed horizontal line) that overlap one of the enhancer states (green) are highlighted and the high PPA variants (PPA >0.1) were tested for allelic imbalance in open chromatin.

https://doi.org/10.7554/eLife.31977.017

A subset of T2D GWAS signals are known to influence T2D risk through a primary effect on insulin secretion, whilst others act primarily through insulin resistance. We used previous categorisations of T2D GWAS loci based on the patterns of association with quantitative measurements of metabolic function and anthropometry (Wood et al., 2017; Dimas et al., 2014), to define a set of 15/48 loci most clearly associated with deficient insulin secretion (and therefore most likely to involve islet dysfunction). At 11 of these 15 loci, we found that islet ‘open strong enhancer’ states, and to a lesser extent ‘open weak enhancer’ and ‘closed strong enhancer’, captured more than 60% of the PPA (median 92%, Figure 4C). Variants in these islet enhancer subclasses also captured at least 95% of the PPA at 4 T2D GWAS loci that could not be classified according to physiological association data but which have been previously implicated in human islet genome or functional regulation based on islet eQTL (van de Bunt et al., 2015) or mQTL (Olsson et al., 2014) data (Figure 4C, genes highlighted in bold). In contrast, at 3/6 of the insulin resistance and all but five unclassified loci, the PPA was mostly (>50%) attributable to other non-islet enhancer states (across all insulin resistance and unclassified loci, DNA not overlapping islet enhancers and defined as ‘Other’ capture a median PPA of 64%). Thus, islet regulatory annotations are particularly useful for fine-mapping T2D GWAS loci that affect insulin secretion and beta-cell function.

To obtain additional evidence to support the localisation of causal variants, we tested for allelic imbalance in ATAC-seq open chromatin data. We selected 54 variants within 33 T2D-associated GWAS segments for testing of allelic imbalance on the basis of (a) a reweighted variant PPA >= 10% and (b) overlap with an enriched regulatory state within the FGWAS T2D joint-model (Figure 4D, Figure 4—source data 2). Of these, 20 variants (at 16 loci) had sufficient numbers of heterozygous samples (>2) and ATAC-seq read depth (depth >9 and at least five reads for each allele). After correcting for mapping bias using WASP, we observed the strongest evidence for allelic imbalance (FDR < 0.05) at 3 out of the 20 variants (rs11257655 near CDC123 and CAMK1D, rs10842991 near KLHDC5 and rs11708067 at ADCY5) (Table 2). All three overlapped refined islet open strong or open weak enhancer regions characterised by open chromatin and hypomethylation.

Table 2
T2D-associated variants with allelic imbalance in open chromatin.
https://doi.org/10.7554/eLife.31977.021
VariantLocusDIAGRAM
P-value
Fgwas
 T2D PPA
Allelic imbalance
Allele Ratio (Allele #)
Allelic imbalance WASP P-valueDirection of effect (T2D)
 rs11708067ADCY58.8E-130.920.29
(38 A VS 94 G alleles)
1.2E-06risk allele A closed
 rs11257655CDC1234.0E-080.950.39
(278 C VS 435 T alleles)
4.5E-09risk allele T open
 rs10842991KLHDC57.3E-070.130.64
(75 C VS 43 T alleles)
4.1E-03risk C allele open

Variant rs11257655 accounts for 95% of the reweighted PPA (compared to a PPA of 20% from genetic data alone) at the CDC123/CAMK1D locus, overlaps an ‘open strong enhancer’ region (Figure 5A) and the risk allele correlates with increased chromatin accessibility. The same variant is in high LD (r2 = 0.82) with the lead variant for a cis-eQTL for CAMK1D in islets (van de Bunt et al., 2015). In experimental assays (Fogarty et al., 2014), the T2D-risk allele has been shown to be associated with increased CAMK1D gene expression and enhanced binding of the FOXA1 and FOXA2 transcription factors. These data all point to rs11257655 as the causal variant at this locus.

Figure 5 with 1 supplement see all
Epigenome Landscape of selected loci with allelic imbalance.

For each locus (A) CDC123, (B) KLHDC5 and C) ADCY5 the following information is shown: Variant level information (depending on the region GWAS lead SNP red, credible set black, eQTL blue and high LD SNPs with r2 >0.8 black), WGBS methylation data (black, middle), four human islet ATAC-seq tracks (green, middle), islet chromatin states (from this study as well as Parker et al., 2013) and Pasquali et al., 2014) and Encode chromatin states from 9 cell types (bottom). For ADCY5 3 ATAC-seq Endoß tracks (top green) and the Capture C results in the Endoß cell line are shown as well (middle blue). Abbreviation for cell types: B-lymphoblastoid cells (GM12878), embryonic stem cells (H1 ES), erythrocytic leukaemia cells (K562), hepatocellular carcinoma cells (HepG2), umbilical vein endothelial cells (HUVEC), mammary epithelial cells (HMEC), skeletal muscle myoblasts (HSMM),normal epidermal keratinocytes (NHEK) and normal lung fibroblasts (NHLF).

https://doi.org/10.7554/eLife.31977.022

At KLHDC5, no clear causal variant emerged based on genetic fine-mapping data alone as the credible set contained 23 variants in high mutual LD (r2 >0.8, top variant PPA <5%, Figure 5B). Of these, variants rs10771372 (genetic fine-mapping PPA = 5%), rs10842992 (genetic fine mapping PPA = 5%) and rs10842991 (genetic fine-mapping PPA = 3%) overlapped ‘open strong enhancer’ regions (Figure 5B), such that their reweighted PPAs rose to 21% (rs10771372), 21% (rs10842992) and 13% (rs10842991), respectively. We observed allelic imbalance only at rs10842991 with the T2D-risk C allele showing greater chromatin accessibility (binomial p=4.1×10−3, Table 2). This variant further overlapped a predicted TFBS motif for PAX6 as determined by the software tool FIMO (Grant et al., 2011): the T2D-risk allele was predicted to enhance PAX6 transcription factor binding consistent with the allelic effects on increasing chromatin accessibility (Figure 5—figure supplement 1A). This strong enhancer region is almost exclusively found in islets, with strong enhancer H3K27ac states overlapping rs10842991 in only two non-islet (heart and smooth muscle) Epigenome Roadmap tissues (out of 99 tissues with 18-state chromatin state information, Figure 5B)(Kundaje et al., 2015). Islet eQTL data (Varshney et al., 2017) also links rs10842991 and close proxy SNPs (including rs7960190) to islet transcription with the risk allele increasing KLHDC5 expression. These data prioritise rs10842991 as the likely causal variant at the KLHDC5 T2D GWAS locus, and indicate a likely molecular mechanism involving modified PAX6 transcription factor binding and an impact on KLHDC5 expression and islet function.

The third example of allelic imbalance mapped to the ADCY5 locus. Fine-mapping based solely on genetic data could not prioritise a distinct causal variant due to multiple variants in high LD (range for top five variants = 12–26%, Figure 5C). However, reweighting of variants based on epigenomic annotation clearly prioritised variant rs11708067: this SNP overlapped an ‘open weak enhancer’ and captured most of the reweighted PPA (PPA = 92%). Allelic imbalance analysis also showed that the T2D-risk A allele was associated with decreased chromatin accessibility (binomial p=1.2×10−6, Table 2). The same lead variant maps to an islet cis-eQTL and methylation QTL (Figure 5C, Figure 5—figure supplement 1B) at which the T2D-risk allele is associated with reduced ADCY5 expression and increased ADCY5 gene body DNA methylation.

To further understand the role of the rs11708067 variant, we performed ATAC-seq and Next Generation Capture-C, in the glucose-responsive human beta-cell line EndoC-βH1 (n = 3). We targeted the ADCY5 promoter to define distal regions interacting with the promoter, and confirmed physical contact with the hypomethylated open chromatin enhancer region harbouring rs11708067 (Figure 5C, Figure 5—figure supplement 1C). To resolve the significance of the interaction between the restriction fragment encompassing rs11708067 and the ADCY5 promoter, we used the programme peakC (de Wit and Geeven, 2017) (https://github.com/deWitLab/peakC) to evaluate the interactions of 12 fragments covering the lead SNP rs11708067 and 15 SNPs in high LD (r2 >0.8) across a region of 47 kb. After adjusting for multiple testing using FDR correction, only two fragments yielded a significant normalised read number over background. This included the open-chromatin overlapping fragment containing rs11708067 and another fragment harbouring rs2877716, rs6798189, rs56371916 (Figure 5—figure supplement 1D). These SNPs fall into a region that did not show evidence of open chromatin.

These findings support rs11708067 as the likely causal variant affecting islet accessible chromatin (in line with another recent study [Roman et al., 2017]), and link the open and hypomethylated enhancer element in which it sits to regulation of ADCY5 expression in islets.

Discussion

A key challenge in the quest to describe the molecular mechanisms through which GWAS signals influence traits of interest, involves the identification of the causal variants responsible and, given that most lie in non-coding sequence, the characterisation of the regulatory elements which they perturb. This underpins efforts to define the effector genes through which these variants operate and to reconstruct the biological networks that are central to disease pathogenesis.

Genetic and physiological studies have highlighted the singular importance of pancreatic islet dysfunction in type 2 diabetes, but epigenomic characterisation of this tissue has been limited in large-scale community projects such as ENCODE and GTEx. The present study seeks to address this deficit by describing, in unprecedented detail, genome-wide patterns of methylation and chromatin accessibility in human islet material. We have combined these data with existing islet epigenomic marks to generate a refined regulatory map which, based on the evidence of improved enrichment for T2D association signals, offers more granular annotation of functional impact.

Our data show that, for DNA methylation, the signal of T2D predisposition is primarily associated with enhancer-like LMRs rather than other categories of methylation elements including UMRs, dDMRs or PMDs. We highlight the strong correlation between islet methylation status and chromatin accessibility but demonstrate that open chromatin predominantly contributes to defining the regulatory impact associated with genetic T2D risk. Finally, we demonstrate how these enhanced epigenomic annotations, when analysed in concert with genetic fine-mapping data and information from allelic imbalance in chromatin accessibility allow us to home in on likely causal variants at T2D association signals such as those near ADCY5, CDC123 and KLHDC5.

While previous studies had explored the candidacy of selected variants at the CDC123 (Fogarty et al., 2014) and ADCY5 (Olsson et al., 2014; Hodson et al., 2014; van de Bunt et al., 2015) loci with respect to islet regulation and T2D predisposition, our integrative analysis of T2D GWAS and epigenetic data has enabled a detailed and comprehensive analysis that considers the regulatory impact of all variants at these loci across multiple islet samples. Our analysis implicates the rs11257655 and rs11708067 variants as the most likely causal variants at the CDC123 and ADCY5 loci respectively and highlights their relationship to islet enhancer activity. The findings at ADCY5 are supported by a recent paper that found allelic imbalance in H3K27 acetylation involving the rs11708067 variant in a single human islet sample, and which observed that deletion of the relevant enhancer element led to reduction in both ADCY5 gene expression and insulin secretion (Roman et al., 2017).

At the KLHDC5 locus, local LD frustrated efforts to define the causal variant using genetic data alone, but the integration of genetic and epigenetic data pinpointed rs10842991 as the likely culprit based on its impact on chromatin accessibility in an open enhancer region. Evidence that this variant co-localises with an islet cis-eQTL signal points to KLHDC5 as the likely downstream target (Varshney et al., 2017). Overall, our integrative approach provides useful insights into the functional mechanisms through which T2D GWAS signals operate. Our findings mirror those from other studies, which have, in various ways, and for other complex traits, combined diverse epigenomic annotations to explore the basis of genetic risk (Wang et al., 2016).

The whole genome methylation data generated in the present study also allowed us to evaluate the likely contribution of previously identified T2D-associated dDMRs (Volkov et al., 2017) with respect to T2D predisposition. These dDMRs, defined on the basis of observed differences in methylation between islets recovered from diabetic and non-diabetic donors, cover a substantial part of the genome, but we were able to show that only a small minority of these overlap functional islet regulatory regions. As a consequence, dDMR regions as a whole had no significant enrichment for T2D association signals. This suggests that most of the dDMR signal involves stochastic effects and/or the secondary consequences on methylation of the diabetic state. However, we cannot exclude that some of the dDMR signals are causal contributors to the diabetic phenotype either because they reflect environmental rather than genetic predisposition, or because they accelerate further perturbation of islet dysfunction as diabetes develops.

Although we provide highly detailed functional fine-mapping of T2D genetic variants to uncover causal variants, the FGWAS approach applied in this study is limited in its ability to determine the effect of multiple variants at individual loci. Specifically, FGWAS relies on the assumption of a single causal variant within each region, which may not necessarily be true for all loci. This assumption could be violated where there are multiple independent signals at a given locus, or where there are multiple (small effect size) variants on a single risk haplotype which jointly impact the phenotype. Analysis methods that combine functional fine-mapping with conditional analysis and consider LD and haplotype patterns are likely to provide a more complete overview of the causal interactions at T2D GWAS loci.

In addition, while the present study characterises islet epigenome status and variability in chromatin accessibility in substantially larger numbers of islet samples than those previously reported (Gaulton et al., 2015; Parker et al., 2013; Pasquali et al., 2014; Varshney et al., 2017), the number of islet preparations for which these data were available was still limited. As a result, our power to detect allelic imbalance in chromatin accessibility was restricted to sites with common variants and relatively large effects. We anticipate that expansion of these sample numbers will extend our capacity to detect such allelic imbalance, and offer more granular insights into the relationships between genetic variation and methylation status. A further limitation is that the genomic data we analysed was generated only from islet samples from non-diabetic donors. Whilst causal inference is possible through the integration of basal epigenomic annotations with genetic data, addition of epigenomic data from islets recovered from diabetic donors has the potential to add a further dimension to such analyses, and to unravel what are likely to be complex causal relationships between genetic variants, epigenomic phenotypes and disease states (Gutierrez-Arcelus et al., 2013). Finally, future work should also focus on experimental validation of likely causal variants and mechanisms e.g. differential binding of the TF PAX6 could be tested at the KLHDC5 rs10842991 variant through electrophoretic mobility shift assays. Our ongoing research efforts are now concentrated on improving the fine-mapping analysis and expanding these genomic enrichment analyses in larger numbers of human islet samples from healthy and diabetic islets. By coupling the integration of these data with empirical functional studies, we expect to provide an increasingly complete description of the causal interactions between DNA methylation, chromatin state, RNA expression and T2D susceptibility.

Materials and methods

Human pancreatic islet samples

WGBS and 450 k array human pancreatic islet sample collection

Human islets were retrieved from deceased Caucasian non-diabetic donors from the Oxford DRWF Human Islet Isolation Facility (n = 34) and at the Alberta Diabetes Institute in Edmonton in Canada (n = 10). For the analysis only samples with a purity >70% were used as determined by dithizone labeling. The Human Research Ethics Board at the University of Alberta (Pro00001754), the University of Oxford's Oxford Tropical Research Ethics Committee (OxTREC Reference: 2–15), or the Oxfordshire Regional Ethics Committee B (REC reference: 09/H0605/2) approved the studies. All organ donors provided informed consent for use of pancreatic tissue in research.

For all WGBS (n = 10) and a subset of 450 k array samples (n = 18) human pancreatic islet DNA was extracted from 100,000 to 150,000 islet cells using Trizol-(Ambion from Thermo Fisher Scientific, Waltham, MA was used for islets processed in Oxford, UK while Trizol from Sigma Aldrich, St. Louis, MO was used for islets processed in Edmonton, Canada) as described previously (van de Bunt et al., 2015). For the remaining 23 samples islet DNA was extracted using the ReliaPrep gDNA Tissue Miniprep system (Promega, Madison, WI). Extracted DNA was stored at −80°C before further use.

ATAC-seq human pancreatic islet sample collection

Human pancreatic islets preparations (n = 18) were retrieved from 17 deceased non-diabetic donors of European descent from the Oxford DRWF Human Islet Isolation Facility and stored for 1–3 days in CMRL or UW media. The latter were reactivated in CMRL for 1 hr before processing them further. Approximately 50,000 islet cells per sample were hand-picked and immediately processed for ATAC-seq as described previously (Buenrostro et al., 2013), however, an additional round of purification was performed using Agencourt AMPure XP magnetic beads (Beckman Coulter, Brea, CA).

WGBS data generation

Bisulphite conversion

400 ng of DNA per human islet samples (n = 10) were sent as part of a collaborative effort to the Blizard Institute, Queen Mary University, London, UK and bisulphite- converted using the Ovation Ultralow Methyl-Seq DR Multiplex System 1–8 (Nugen, Manchester, UK) and purified using Agencourt AMPure beads (Beckman Coulter) as described previously (Lowe et al., 2013).

Library generation and processing of reads

The libraries were sequenced by members of the High-Throughput Genomics group at the Wellcome Centre for Human Genetics, University of Oxford, Oxford, UK. Samples were sequenced as multiplex libraries across 3 HiSeq2000 lanes with 100 bp paired-end read length (including a PhIX spike-in of 5%) to obtain high-coverage read data. The obtained reads were trimmed using a customized python3 script (10 bp at the start and 15 bp at the end) and aligned to hg19 using the software Bismark (settings: L,0,–0.6, version 0.12.5, RRID:SCR_005604)(Krueger and Andrews, 2011). Specifically, paired-end alignment of trimmed reads was performed and unmapped reads from read one were realigned using Bismark and merged with the paired-end alignment using samtools (Li et al., 2009) (version 0.1.19, RRID:SCR_002105) in order to increase mapping efficiency. Coverage for the merged paired-end and realigned HiSeq read alignments was estimated for the human mappable genome (NCBI hg19 2.8 billion base pairs excluding gaps and unmappable and blacklisted regions according to UCSC and Encode [ENCODE Project Consortium, 2012]) using bedtools (version v2.21.0) (Quinlan, 2014).

WGBS DNA methylation quantification and prediction of hypomethylated regulatory regions

CpG methylation levels were determined for each sample by calculating the ratio of unmodified C (methylated) and bisulphite converted T (unmethylated) alleles using BiFAST (first described here [Lowe et al., 2013]). High-pass pooled WGBS data was generated by adding methylated and unmethylated read counts across individual low-pass samples to then estimate the average beta methylation levels.

Regulatory regions were identified using the R package methylseek (RRID:SCR_006513) (Burger et al., 2013). After removing PMDs, which represent highly heterogenous methylation states determined by DNA sequence features (Gaidatzis et al., 2014), LMRs (<30 CpGs) and UMRs (>30 CpGs) were predicted in hypomethylated regions (<50%) at an FDR of 0.05. The methylation level and FDR parameter was inferred from the data as suggested by the methylseek workflow (Burger et al., 2013).

450 k DNA methylation array data generation

In total, 41 samples were processed for the Illumina Infinium HumanMethylation450 BeadChip (Illumina, San Diego, CA). Of these 18 samples were bisulphite-converted and processed as part of a collaboration at the UCL Cancer Institute, University College London, London, UK while the remaining 23 samples were processed in OCDEM, University of Oxford, Oxford, UK. The DNA was bisulphite converted using the EZ DNA MethylationTM Kit ( Zymogen Research Corp, Irvine, CA) and hybridised to the Illumina 450 k array and scanned with iScan (Illumina) according to the manufacturer’s protocol.

The resulting data was analysed using the Package minfi (RRID:SCR_012830) (Aryee et al., 2014) and custom R scripts ([R Development Core Team, 2011], R version 3.0.2, RRID:SCR_001905). Specifically, CpG sites with a detection p-value>0.01 were removed from the analysis and samples with >5% of CpG sites failing this threshold (n = 9) were also removed from the analysis.

Following separate quantile normalisation of signal intensities derived from methylated and unmethylated Type I probes and Type II probes, methylation levels (ß) were estimated, based on the intensities of the methylated (M) and unmethylated (U) signal in the following way: β = M/(M + U + 100). To correct for batch effects the ComBat function implemented in the sva (Johnson et al., 2007; Leek et al., 2007) package was used (Figure 1—figure supplement 2).

ATAC-seq data generation

Sequencing of ATAC-seq reads

ATAC-seq libraries were sequenced at the High-Throughput Genomics group which is part of the Wellcome Centre for Human Genetics, University of Oxford, Oxford, UK. Samples were sequenced as 4-6plex libraries across 1–3 Hiseq2500 lanes with 50 bp paired-end read length.

Processing of ATAC-seq reads

Raw FASTQ reads were processed with an in-house pipeline (first described in (Hay et al., 2016) and on the website http://userweb.molbiol.ox.ac.uk/public/telenius/PipeSite.html). Specifically, library and sequencing quality was checked with FASTQC (RRID:SCR_014583) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) and reads were mapped to the human genome (hg19) via bowtie (Langmead et al., 2009) (version 1.1.0, RRID:SCR_005476) with default settings but -m 2, and maxins 2000 which allows mapping of reads with a maximum number of 2 alignments and a maximum insert size of 2000 bp. For reads that could not be aligned the first time, adapters were removed with Trim Galore at the three prime end (RRID:SCR_011847, settings -length 10, -qualFilter 20, http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) to enhance the chance of mapping. The resulting trimmed reads were then mapped again with bowtie. Any remaining unmapped and trimmed reads were processed with FLASH (Magoč and Salzberg, 2011) (version 1.2.8, RRID:SCR_005531, settings -m 9 -x 0.125) which combines overlapping read pairs and reconstructs read pairs without overlap. These are then realigned a third time using bowtie. PCR duplicates are then removed from the mapped bam files using samtools rmdup function (Li et al., 2009). Additionally, all reads overlapping any of the ‘unmappable’ UCSC Duke blacklisted hg19 regions (ENCODE Project Consortium, 2012) are also removed from the final bam file.

Open chromatin peaks were called through the aforementioned in-house pipeline by applying sample-specific read depth and width parameters, which were chosen based on the signal to noise ratio of a given sample.

ChIP-seq data and identification of chromatin states

Processing of available ChIP-seq data

Human islet ChIP-seq histone mark and TFBS data were obtained from various sources: H3K4me1, CTCF and H3K27ac (Pasquali et al., 2014), H3K36me3 and H3K4me3 (Morán et al., 2012) and H3K27me3 (Kundaje et al., 2015). Available raw fastq files were aligned to hg19 using bowtie1 (version 1.1.1) with modified default settings (-m 1 which removes reads with more than one valid alignment and -n 1 which defines the maximum number of mismatches in seed) and PCR duplicates were removed from the aligned bam files with Picard tools (RRID:SCR_006525, v1.119, http://broadinstitute.github.io/picard/). The resulting reads were converted into bed format using the bedtools bamToBed function (Quinlan, 2014) (RRID:SCR_006646, version v2.21.0) and extended by 200 bp towards the 3’ end to account for fragment size.

Identification of chromatin states using chromHMM

Binarised 200 bp density maps from the bed files of the 6 ChIP-seq marks were created using a Poisson distribution implemented in the BinaryBed function of the ChromHMM software as described in (Ernst and Kellis, 2012; Ernst et al., 2011). From these epigenomic density maps, 11 ChIP-only chromatin states were derived using a multivariate Hidden Markov Model implemented in the Learnmodel function (standard settings, h19 genome) of the software ChromHMM (Ernst and Kellis, 2012).

To generate additional sets of chromatin states based on ChIP-seq, ATAC-seq and DNA methylation data, ATAC-seq open chromatin and DNA methylation status were binarised. Specifically, ATAC-seq peaks (presence/absence) and whole-genome CpG methylation status (hypermethylation/hypomethylation based on a threshold of 60% methylation) were binarised across 200 bp windows of the genome.

These binarised 200 bp ChIP-seq, ATAC-seq and DNA methylation maps were combined and used to generate 3 sets of chromatin states derived from ChIP and DNA methylation data (ChIP +Meth), ChIP and ATAC-seq data (ChIP +ATAC) or ChIP, ATAC-seq and DNA methylation data (ChIP +ATAC + Meth) using the Learnmodel ChromHMM function (Figure 3A and Figure 3—figure supplement 1A–B). As suggested by (Ernst et al., 2011), after evaluating models with up to 20 chromatin states, a 15 state model was chosen based on the resolution provided by identified states

ADCY5 capture C analysis and ATAC-seq in EndoC-ßH1

Next-generation Capture-C was performed in order to map physical chromatin interactions with the ADCY5 promoter in EndoC-ßH1 (RRID:CVCL_L909) cell lines (n = 3) (see protocol in Materials and methods in [Davies et al., 2016]).

In brief, chromatin conformation capture (3C) libraries were generated by formaldehyde fixation prior to DpnII restriction enzyme digestion and subsequent DNA ligation. Following cross-link reversal, DNA extraction and sonication, sequencing adapters were added to sonicated fragments (~200 bp). Library fragments were subjected to a double capture through hybridisation with a biotinylated oligonucleotide probes encompassing the ADCY5 promoter and enriched using streptavidin bead pulldown. PCR amplified fragments were then sequenced through paired-end sequencing (Illumina Next-Seq). An in silico restriction enzyme digestion was performed on the set of reconstructed fragments (from paired-end sequenced reads) using the DpnII2E.pl script (Davies, 2015)(https://github.com/Hughes-Genome-Group/captureC). Uncaptured reads and PCR duplicates were removed prior to mapping to the human genome (hg19) with Bowtie (Langmead et al., 2009)(v 1.1.0). Removal of PCR duplicates and classification of fragments as ‘capture’ (i.e. including the ADCY5 promoter) or ‘reporter’ (outside the capture fragment on exclusion region) was performed with the CCanalyser2.pl wrapper (Davies, 2015)(https://github.com/Hughes-Genome-Group/captureC). Unique mapped interactions were normalized to the total number of cis interactions (i.e. same chromosome) per 100,000 interactions. Significant chromatin interactions were determined from a rank-sum test implemented in the program peakC (de Wit and Geeven, 2017)(https://github.com/deWitLab/peakC). Specifically, we evaluated interactions involving all SNPs in high LD (r2 >0.8) with the lead rs11708067. The lead variant (rs11708067) was in high LD with 15 SNPs (mapping to 12 DpnII fragments) that spanned a region of 47 kb. We applied the Benjamini-Hochberg correction to control the false discovery rate for the set of p-values corresponding to each restriction fragment within the 47 kb region at the ADCY5 locus.

In addition, ATAC-seq was performed in 50,000 cells of EndoC-ßH1 cell lines (n = 3) and the data was analysed in the same way as described above for human islet samples.

Endo-βH1 cells were obtained from Endocells and have been previously authenticated (Ravassard et al., 2011). In addition, the cell line was tested and found negative for mycoplasma contamination.

Overlaying generated epigenomic datasets generated here with other genomic regulatory regions

CpG sites and/or hypomethylated regulatory regions identified from the WGBS and/or 450 k array data were overlapped with existing islet chromatin state maps (Parker et al., 2013), islet transcription factor binding sites (FOXA2, MAFB, NKX2.2, NKX6.1, PDX1), T2D-associated islet dDMRs (Dayeh et al., 2014) and eQTLs (van de Bunt et al., 2015). Similarly, ATAC-seq open chromatin peaks generated here were overlapped with publicly available ATAC-seq peaks (Varshney et al., 2017).

In addition, we also obtained the 850 k array manifest file to determine overlap of 850 k array CpG sites with GWAS credible set regions (https://support.illumina.com/downloads/infinium-methylationepic-v1-0-product-files.html).

Genetic datasets used in this study

Credible sets from the DIAGRAM (Scott et al., 2017)(involving 26.7 k cases and 132.5 k controls of predominantly European origin, imputed to the 1000G March 2012 reference panel) and ENGAGE (Horikoshi et al., 2015)(including 46.7 k individuals, imputed to the 1000G March 2012 reference panel) consortium were used to compare the ability of the 450 k, 850 k and WGBS methylation array to interrogate T2D and FG GWAS regions.

The DIAGRAM and ENGAGE GWAS SNP summary level data was used for the FGWAS analysis to determine enrichment of regulatory annotations in T2D and FG GWAS signal.

Furthermore, data from (Wood et al., 2017) and (Dimas et al., 2014) were used to categories T2D GWAS loci into physiological groups of insulin secretion, insulin resistance or unclassified loci.

Statistical and computational analysis

Enrichment analysis of identified regulatory annotations in other genomic annotations

Enrichment of hypomethylated regulatory regions (LMRs and UMRs, result section 2.2.) and ATAC-seq open chromatin peaks (result section 2.3) in the aforementioned genomic annotations (method section 4.6) was determined through 100,000 random permutations. P-values and fold enrichment was determined by comparing the true overlap results to the permuted overlap results. The resulting P-values were multiple testing corrected using Bonferroni correction (an adjusted p-value<0.05 was considered significant).

FGWAS enrichment analysis

FGWAS (Pickrell, 2014) (version 0.3.6) applied a hierarchical model that determined shared properties of loci affecting a trait. The FGWAS model used SNP-based GWAS summary level data and divided the genome into windows (setting ‘k’=5000 which represents the number of SNPs per window), which are larger than the expected LD patterns in the population. The model assumed that each window either contained a single SNP that affected the trait or that there was no SNP in the window that influenced the trait. The model estimated the prior probability of a window to contain an association and the conditional prior probability that a SNP within the window was the causal variant. These prior probabilities were variable, dependent on regional annotations and estimated based on enrichment patterns of annotations across the genome using a Bayes approach.

FGWAS single state analysis

FGWAS was used with standard settings to determine enrichment of individual islet chromatin states, LMRs, UMRs, PMDS and ATAC-seq open chromatin peaks, CDS and CONS sequence in DIAGRAM (setting ‘cc’ was applied for use with T2D-case-control GWAS data) and ENGAGE GWAS SNP summary level data.

For each individual annotation, the model provided maximum likelihood enrichment parameters and annotations were considered as significantly enriched if the parameter estimate and 95% CI was above zero.

FGWAS joint model analysis

To determine the maximum likelihood model the following approach suggested by (Pickrell, 2014) was used for each set of chromatin states (ChIP-only, ChIP +ATAC, ChIP +Meth and ChIP +ATAC + Meth), separately. In addition, CDS and CONS sequenced were used as well for each set of chromatin states in the joint analysis. Firstly, a model was fitted for each annotation individually to identify all annotations that were significantly enriched with the trait. Secondly, the annotation with the highest increase (and enrichment) in the maximum log-likelihood was added to the model and the analysis was repeated with the additional annotation. Thirdly, annotations were added as long as they increase the maximum log-likelihood of the newly tested model. Fourthly, a 10-fold cross-validation approach was used after determining a penalty parameter based on the maximum likelihood of a penalised log-likelihood function to avoid overfitting. Fifthly, each annotation was dropped from the model and the maximum cross-validation likelihood was evaluated. If a reduced model has a higher cross-validation maximum likelihood, additional annotations are dropped until the model cannot be further improved. This model was described as the best fitted model and used for the remaining analysis. The maximum likelihood enrichment parameters and 95% CI for each annotation of the best model were reported (independent of significance).

Comparing FGWAS enrichment parameter across chromatin states

Initially, similar enhancer chromatin states derived from the four different ChromHMM analyses (ChIP-only, ChIP + ATAC, ChIP + Meth, ChIP + ATAC + Meth) were compared. Similarity was determined based on shared histone chromatin marks according to the chromHMM emission parameters. Further comparisons between the ChIP-only and ChIP + ATAC + Meth model were performed based on the reweighted FGWAS maximum variant PPA and the number of reweighted 99% credible set variants per T2D locus (for details regarding FGWAS PPA see next section).

However, considering that the chromatin states were derived from distinct sets of annotations across different analyses of ChromHMM, a direct comparison was not fully possible. Hence, a nested model approach was used to further dissect the contribution of open chromatin and DNA methylation to the enrichment. Specifically, an FGWAS analysis was performed that combined the ChIP-only chromHMM states with raw LMRs (representing DNA methylation) and ATAC-seq peaks (representing open chromatin). After determining the best maximum-likelihood cross-validation model (combining ChIP-only, ATAC-seq and LMR states) a nested model and log-likelihood ratio test were used to determine the contribution of each annotation to the model (Figure 3—figure supplement 1D).

Reweighting of variant PPA and testing of allelic imbalance

The enrichment priors derived from the FGWAS maximum likelihood model were used as a basis for evaluating both the significance and functional impact of associated variants in GWAS regions; allowing variants that map into annotations that show global enrichment to be afforded extra weight.

Specifically, variants at significant GWAS regions with a high FGWAS PPA (PPA >= 10%) and overlapping open enhancer states were prioritised for further follow-up. Genome-wide significance of loci was determined based on P-values (p<5×10−8) or a regional FGWAS PPA >= 90% (representing the sum of the PPAs of all SNPs in a given region). The latter threshold is based on a recommendation from (Pickrell, 2014) who observed that a regional PPA of 90% or above can be used to identify sub-threshold GWAS loci.

Of the prioritised variants, only variants with at least two heterozygous samples and ATAC-seq read depth of at least nine reads (minimum five reads for each allele) were tested for allelic imbalance.

To avoid read-mapping and reference allele bias the software WASP (van de Geijn et al., 2015)(Version 0.2261) was used to remove reads associated with mapping bias. In short, reads of the unfiltered bam file that overlapped the variant of interest were identified. For each read overlapping an SNP, the genotype of that SNP was changed to the alternative allele and the read was remapped using bwa (Li and Durbin, 2009) (version 0.5.8 c). Any read that failed to realign in the same position in the genome was discarded. Ultimately, PCR duplicates were filtered using the WASP ‘rmdup_pe.py’ script, which removed duplicated reads randomly (independent of the mapping score) to avoid any bias.

Allelic imbalance was determined using a binomial test as implemented in R.

Identification of TFBS at SNPs that display allelic imbalance

The tool ‘Fimo’(Grant et al., 2011) implemented in the ‘meme’ software package (RRID:SCR_001783) was applied to identify TF motifs that significantly (FDR < 0.05) matched the sequence overlapping a SNP variant showing allelic imbalance (20 bp up and downstream).

Overlap of regulatory regions

Overlap between genomic regulatory regions was performed using bedtools intersectBed function (Quinlan, 2014) (version 2.21.0). Summary statistics across 200 bp windows were determined using bedtools mapBed function. Random permutations of regulatory regions were performed by applying the bedtools shuffleBed function.

Statistical analysis

All statistical analysis (unless otherwise stated) was performed using R (version 3.0.2) including Spearman’s correlation analysis to compare the 450 k and WGBS array, the KS-test to compare 450 k and WGBS DNA methylation distributions, the binomial test to evaluate allelic imbalance and principal component analysis to identify batch effects in the 450 k data. Significance is defined as p<0.05 unless otherwise stated.

Visualisation and figure generation

All figures unless otherwise stated were generated using R (version 3.0.2) and/or ggplot2(Wickham, 2009). Figure 1E was generated using locuszoom (Pruim et al., 2010). Chromatin state CHiP-seq enrichment maps (Figure 3A, Figure 3—figure supplement 1A–B) were generated using chromHMM (Ernst and Kellis, 2012). The genome-browser views (Figure 5) were generated using the UCSC genome browser tool (Kent et al., 2002).

Sequencing data

ATAC-seq and WGBS sequencing data has been deposited at the EBI hosted European Genome-phenome Archive (EGA, https://ega-archive.org/) and is accessible via the EGA accession numbers: EGAS00001002592, EGAD00001003946 and EGAD00001003947.

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
    TxDb.Hsapiens.UCSC.hg19.knownGene: Annotation Package for TxDb Object(s)
    1. M Carlson
    2. B Maintainer
    (2015)
    TxDb.Hsapiens.UCSC.hg19.knownGene: Annotation Package for TxDb Object(s), R package version 3.2.2.
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
    Genetic fine mapping and genomic annotation defines causal mechanisms at type 2 diabetes susceptibility loci
    1. KJ Gaulton
    2. T Ferreira
    3. Y Lee
    4. A Raimondo
    5. R Mägi
    6. ME Reschen
    7. A Mahajan
    8. A Locke
    9. NW Rayner
    10. N Robertson
    11. RA Scott
    12. I Prokopenko
    13. LJ Scott
    14. T Green
    15. T Sparso
    16. D Thuillier
    17. L Yengo
    18. H Grallert
    19. S Wahl
    20. M Frånberg
    21. RJ Strawbridge
    22. H Kestler
    23. H Chheda
    24. L Eisele
    25. S Gustafsson
    26. V Steinthorsdottir
    27. G Thorleifsson
    28. L Qi
    29. LC Karssen
    30. EM van Leeuwen
    31. SM Willems
    32. M Li
    33. H Chen
    34. C Fuchsberger
    35. P Kwan
    36. C Ma
    37. M Linderman
    38. Y Lu
    39. SK Thomsen
    40. JK Rundle
    41. NL Beer
    42. M van de Bunt
    43. A Chalisey
    44. HM Kang
    45. BF Voight
    46. GR Abecasis
    47. P Almgren
    48. D Baldassarre
    49. B Balkau
    50. R Benediktsson
    51. M Blüher
    52. H Boeing
    53. LL Bonnycastle
    54. EP Bottinger
    55. NP Burtt
    56. J Carey
    57. G Charpentier
    58. PS Chines
    59. MC Cornelis
    60. DJ Couper
    61. AT Crenshaw
    62. RM van Dam
    63. AS Doney
    64. M Dorkhan
    65. S Edkins
    66. JG Eriksson
    67. T Esko
    68. E Eury
    69. J Fadista
    70. J Flannick
    71. P Fontanillas
    72. C Fox
    73. PW Franks
    74. K Gertow
    75. C Gieger
    76. B Gigante
    77. O Gottesman
    78. GB Grant
    79. N Grarup
    80. CJ Groves
    81. M Hassinen
    82. CT Have
    83. C Herder
    84. OL Holmen
    85. AB Hreidarsson
    86. SE Humphries
    87. DJ Hunter
    88. AU Jackson
    89. A Jonsson
    90. ME Jørgensen
    91. T Jørgensen
    92. WH Kao
    93. ND Kerrison
    94. L Kinnunen
    95. N Klopp
    96. A Kong
    97. P Kovacs
    98. P Kraft
    99. J Kravic
    100. C Langford
    101. K Leander
    102. L Liang
    103. P Lichtner
    104. CM Lindgren
    105. E Lindholm
    106. A Linneberg
    107. CT Liu
    108. S Lobbens
    109. J Luan
    110. V Lyssenko
    111. S Männistö
    112. O McLeod
    113. J Meyer
    114. E Mihailov
    115. G Mirza
    116. TW Mühleisen
    117. M Müller-Nurasyid
    118. C Navarro
    119. MM Nöthen
    120. NN Oskolkov
    121. KR Owen
    122. D Palli
    123. S Pechlivanis
    124. L Peltonen
    125. JR Perry
    126. CG Platou
    127. M Roden
    128. D Ruderfer
    129. D Rybin
    130. YT van der Schouw
    131. B Sennblad
    132. G Sigurðsson
    133. A Stančáková
    134. G Steinbach
    135. P Storm
    136. K Strauch
    137. HM Stringham
    138. Q Sun
    139. B Thorand
    140. E Tikkanen
    141. A Tonjes
    142. J Trakalo
    143. E Tremoli
    144. T Tuomi
    145. R Wennauer
    146. S Wiltshire
    147. AR Wood
    148. E Zeggini
    149. I Dunham
    150. E Birney
    151. L Pasquali
    152. J Ferrer
    153. RJ Loos
    154. J Dupuis
    155. JC Florez
    156. E Boerwinkle
    157. JS Pankow
    158. C van Duijn
    159. E Sijbrands
    160. JB Meigs
    161. FB Hu
    162. U Thorsteinsdottir
    163. K Stefansson
    164. TA Lakka
    165. R Rauramaa
    166. M Stumvoll
    167. NL Pedersen
    168. L Lind
    169. SM Keinanen-Kiukaanniemi
    170. E Korpi-Hyövälti
    171. TE Saaristo
    172. J Saltevo
    173. J Kuusisto
    174. M Laakso
    175. A Metspalu
    176. R Erbel
    177. KH Jöcke
    178. S Moebus
    179. S Ripatti
    180. V Salomaa
    181. E Ingelsson
    182. BO Boehm
    183. RN Bergman
    184. FS Collins
    185. KL Mohlke
    186. H Koistinen
    187. J Tuomilehto
    188. K Hveem
    189. I Njølstad
    190. P Deloukas
    191. PJ Donnelly
    192. TM Frayling
    193. AT Hattersley
    194. U de Faire
    195. A Hamsten
    196. T Illig
    197. A Peters
    198. S Cauchi
    199. R Sladek
    200. P Froguel
    201. T Hansen
    202. O Pedersen
    203. AD Morris
    204. CN Palmer
    205. S Kathiresan
    206. O Melander
    207. PM Nilsson
    208. LC Groop
    209. I Barroso
    210. C Langenberg
    211. NJ Wareham
    212. CA O'Callaghan
    213. AL Gloyn
    214. D Altshuler
    215. M Boehnke
    216. TM Teslovich
    217. MI McCarthy
    218. AP Morris
    219. DIAbetes Genetics Replication And Meta-analysis (DIAGRAM) Consortium
    (2015)
    Nature Genetics 47:1415–1425.
    https://doi.org/10.1038/ng.3437
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
    Trans-ancestry genome-wide association study identifies 12 genetic loci influencing blood pressure and implicates a role for DNA methylation
    1. N Kato
    2. M Loh
    3. F Takeuchi
    4. N Verweij
    5. X Wang
    6. W Zhang
    7. TN Kelly
    8. D Saleheen
    9. B Lehne
    10. IM Leach
    11. AW Drong
    12. J Abbott
    13. S Wahl
    14. ST Tan
    15. WR Scott
    16. G Campanella
    17. M Chadeau-Hyam
    18. U Afzal
    19. TS Ahluwalia
    20. MJ Bonder
    21. P Chen
    22. A Dehghan
    23. TL Edwards
    24. T Esko
    25. MJ Go
    26. SE Harris
    27. J Hartiala
    28. S Kasela
    29. A Kasturiratne
    30. CC Khor
    31. ME Kleber
    32. H Li
    33. Z Yu Mok
    34. M Nakatochi
    35. NS Sapari
    36. R Saxena
    37. AFR Stewart
    38. L Stolk
    39. Y Tabara
    40. AL Teh
    41. Y Wu
    42. JY Wu
    43. Y Zhang
    44. I Aits
    45. A Da Silva Couto Alves
    46. S Das
    47. R Dorajoo
    48. JC Hopewell
    49. YK Kim
    50. RW Koivula
    51. J Luan
    52. LP Lyytikäinen
    53. QN Nguyen
    54. MA Pereira
    55. I Postmus
    56. OT Raitakari
    57. MS Bryan
    58. RA Scott
    59. R Sorice
    60. V Tragante
    61. M Traglia
    62. J White
    63. K Yamamoto
    64. Y Zhang
    65. LS Adair
    66. A Ahmed
    67. K Akiyama
    68. R Asif
    69. T Aung
    70. I Barroso
    71. A Bjonnes
    72. TR Braun
    73. H Cai
    74. LC Chang
    75. CH Chen
    76. CY Cheng
    77. YS Chong
    78. R Collins
    79. R Courtney
    80. G Davies
    81. G Delgado
    82. LD Do
    83. PA Doevendans
    84. RT Gansevoort
    85. YT Gao
    86. TB Grammer
    87. N Grarup
    88. J Grewal
    89. D Gu
    90. GS Wander
    91. AL Hartikainen
    92. SL Hazen
    93. J He
    94. CK Heng
    95. JE Hixson
    96. A Hofman
    97. C Hsu
    98. W Huang
    99. LLN Husemoen
    100. JY Hwang
    101. S Ichihara
    102. M Igase
    103. M Isono
    104. JM Justesen
    105. T Katsuya
    106. MG Kibriya
    107. YJ Kim
    108. M Kishimoto
    109. WP Koh
    110. K Kohara
    111. M Kumari
    112. K Kwek
    113. NR Lee
    114. J Lee
    115. J Liao
    116. W Lieb
    117. DCM Liewald
    118. T Matsubara
    119. Y Matsushita
    120. T Meitinger
    121. E Mihailov
    122. L Milani
    123. R Mills
    124. N Mononen
    125. M Müller-Nurasyid
    126. T Nabika
    127. E Nakashima
    128. HK Ng
    129. K Nikus
    130. T Nutile
    131. T Ohkubo
    132. K Ohnaka
    133. S Parish
    134. L Paternoster
    135. H Peng
    136. A Peters
    137. ST Pham
    138. MJ Pinidiyapathirage
    139. M Rahman
    140. H Rakugi
    141. O Rolandsson
    142. M Ann Rozario
    143. D Ruggiero
    144. CF Sala
    145. R Sarju
    146. K Shimokawa
    147. H Snieder
    148. T Sparsø
    149. W Spiering
    150. JM Starr
    151. DJ Stott
    152. DO Stram
    153. T Sugiyama
    154. S Szymczak
    155. WHW Tang
    156. L Tong
    157. S Trompet
    158. V Turjanmaa
    159. H Ueshima
    160. AG Uitterlinden
    161. S Umemura
    162. M Vaarasmaki
    163. RM van Dam
    164. WH van Gilst
    165. DJ van Veldhuisen
    166. JS Viikari
    167. M Waldenberger
    168. Y Wang
    169. A Wang
    170. R Wilson
    171. TY Wong
    172. YB Xiang
    173. S Yamaguchi
    174. X Ye
    175. RD Young
    176. TL Young
    177. JM Yuan
    178. X Zhou
    179. FW Asselbergs
    180. M Ciullo
    181. R Clarke
    182. P Deloukas
    183. A Franke
    184. PW Franks
    185. S Franks
    186. Y Friedlander
    187. MD Gross
    188. Z Guo
    189. T Hansen
    190. MR Jarvelin
    191. T Jørgensen
    192. JW Jukema
    193. M Kähönen
    194. H Kajio
    195. M Kivimaki
    196. JY Lee
    197. T Lehtimäki
    198. A Linneberg
    199. T Miki
    200. O Pedersen
    201. NJ Samani
    202. TIA Sørensen
    203. R Takayanagi
    204. D Toniolo
    205. H Ahsan
    206. H Allayee
    207. YT Chen
    208. J Danesh
    209. IJ Deary
    210. OH Franco
    211. L Franke
    212. BT Heijman
    213. JD Holbrook
    214. A Isaacs
    215. BJ Kim
    216. X Lin
    217. J Liu
    218. W März
    219. A Metspalu
    220. KL Mohlke
    221. DK Sanghera
    222. XO Shu
    223. JBJ van Meurs
    224. E Vithana
    225. AR Wickremasinghe
    226. C Wijmenga
    227. BHW Wolffenbuttel
    228. M Yokota
    229. W Zheng
    230. D Zhu
    231. P Vineis
    232. SA Kyrtopoulos
    233. JCS Kleinjans
    234. MI McCarthy
    235. R Soong
    236. C Gieger
    237. J Scott
    238. YY Teo
    239. J He
    240. P Elliott
    241. ES Tai
    242. P van der Harst
    243. JS Kooner
    244. JC Chambers
    245. BIOS-consortium
    246. CARDIo GRAMplusCD
    247. LifeLines Cohort Study
    248. InterAct Consortium
    (2015)
    Nature Genetics 47:1282–1293.
    https://doi.org/10.1038/ng.3405
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
    Sva: Surrogate Variable Analysis, version R package version 3.8.0
    1. JT Leek
    2. WE Johnson
    3. HS Parker
    4. AE Jaffe
    5. JD Storey
    (2007)
    Sva: Surrogate Variable Analysis.
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
    Genome-wide trans-ancestry meta-analysis provides insight into the genetic architecture of type 2 diabetes susceptibility
    1. A Mahajan
    2. MJ Go
    3. W Zhang
    4. JE Below
    5. KJ Gaulton
    6. T Ferreira
    7. M Horikoshi
    8. AD Johnson
    9. MC Ng
    10. I Prokopenko
    11. D Saleheen
    12. X Wang
    13. E Zeggini
    14. GR Abecasis
    15. LS Adair
    16. P Almgren
    17. M Atalay
    18. T Aung
    19. D Baldassarre
    20. B Balkau
    21. Y Bao
    22. AH Barnett
    23. I Barroso
    24. A Basit
    25. LF Been
    26. J Beilby
    27. GI Bell
    28. R Benediktsson
    29. RN Bergman
    30. BO Boehm
    31. E Boerwinkle
    32. LL Bonnycastle
    33. N Burtt
    34. Q Cai
    35. H Campbell
    36. J Carey
    37. S Cauchi
    38. M Caulfield
    39. JC Chan
    40. LC Chang
    41. TJ Chang
    42. YC Chang
    43. G Charpentier
    44. CH Chen
    45. H Chen
    46. YT Chen
    47. KS Chia
    48. M Chidambaram
    49. PS Chines
    50. NH Cho
    51. YM Cho
    52. LM Chuang
    53. FS Collins
    54. MC Cornelis
    55. DJ Couper
    56. AT Crenshaw
    57. RM van Dam
    58. J Danesh
    59. D Das
    60. U de Faire
    61. G Dedoussis
    62. P Deloukas
    63. AS Dimas
    64. C Dina
    65. AS Doney
    66. PJ Donnelly
    67. M Dorkhan
    68. C van Duijn
    69. J Dupuis
    70. S Edkins
    71. P Elliott
    72. V Emilsson
    73. R Erbel
    74. JG Eriksson
    75. J Escobedo
    76. T Esko
    77. E Eury
    78. JC Florez
    79. P Fontanillas
    80. NG Forouhi
    81. T Forsen
    82. C Fox
    83. RM Fraser
    84. TM Frayling
    85. P Froguel
    86. P Frossard
    87. Y Gao
    88. K Gertow
    89. C Gieger
    90. B Gigante
    91. H Grallert
    92. GB Grant
    93. LC Grrop
    94. CJ Groves
    95. E Grundberg
    96. C Guiducci
    97. A Hamsten
    98. BG Han
    99. K Hara
    100. N Hassanali
    101. AT Hattersley
    102. C Hayward
    103. AK Hedman
    104. C Herder
    105. A Hofman
    106. OL Holmen
    107. K Hovingh
    108. AB Hreidarsson
    109. C Hu
    110. FB Hu
    111. J Hui
    112. SE Humphries
    113. SE Hunt
    114. DJ Hunter
    115. K Hveem
    116. ZI Hydrie
    117. H Ikegami
    118. T Illig
    119. E Ingelsson
    120. M Islam
    121. B Isomaa
    122. AU Jackson
    123. T Jafar
    124. A James
    125. W Jia
    126. KH Jöckel
    127. A Jonsson
    128. JB Jowett
    129. T Kadowaki
    130. HM Kang
    131. S Kanoni
    132. WH Kao
    133. S Kathiresan
    134. N Kato
    135. P Katulanda
    136. KM Keinanen-Kiukaanniemi
    137. AM Kelly
    138. H Khan
    139. KT Khaw
    140. CC Khor
    141. HL Kim
    142. S Kim
    143. YJ Kim
    144. L Kinnunen
    145. N Klopp
    146. A Kong
    147. E Korpi-Hyövälti
    148. S Kowlessur
    149. P Kraft
    150. J Kravic
    151. MM Kristensen
    152. S Krithika
    153. A Kumar
    154. J Kumate
    155. J Kuusisto
    156. SH Kwak
    157. M Laakso
    158. V Lagou
    159. TA Lakka
    160. C Langenberg
    161. C Langford
    162. R Lawrence
    163. K Leander
    164. JM Lee
    165. NR Lee
    166. M Li
    167. X Li
    168. Y Li
    169. J Liang
    170. S Liju
    171. WY Lim
    172. L Lind
    173. CM Lindgren
    174. E Lindholm
    175. CT Liu
    176. JJ Liu
    177. S Lobbens
    178. J Long
    179. RJ Loos
    180. W Lu
    181. J Luan
    182. V Lyssenko
    183. RC Ma
    184. S Maeda
    185. R Mägi
    186. S Männisto
    187. DR Matthews
    188. JB Meigs
    189. O Melander
    190. A Metspalu
    191. J Meyer
    192. G Mirza
    193. E Mihailov
    194. S Moebus
    195. V Mohan
    196. KL Mohlke
    197. AD Morris
    198. TW Mühleisen
    199. M Müller-Nurasyid
    200. B Musk
    201. J Nakamura
    202. E Nakashima
    203. P Navarro
    204. PK Ng
    205. AC Nica
    206. PM Nilsson
    207. I Njølstad
    208. MM Nöthen
    209. K Ohnaka
    210. TH Ong
    211. KR Owen
    212. CN Palmer
    213. JS Pankow
    214. KS Park
    215. M Parkin
    216. S Pechlivanis
    217. NL Pedersen
    218. L Peltonen
    219. JR Perry
    220. A Peters
    221. JM Pinidiyapathirage
    222. CG Platou
    223. S Potter
    224. JF Price
    225. L Qi
    226. V Radha
    227. L Rallidis
    228. A Rasheed
    229. W Rathman
    230. R Rauramaa
    231. S Raychaudhuri
    232. NW Rayner
    233. SD Rees
    234. E Rehnberg
    235. S Ripatti
    236. N Robertson
    237. M Roden
    238. EJ Rossin
    239. I Rudan
    240. D Rybin
    241. TE Saaristo
    242. V Salomaa
    243. J Saltevo
    244. M Samuel
    245. DK Sanghera
    246. J Saramies
    247. J Scott
    248. LJ Scott
    249. RA Scott
    250. AV Segrè
    251. J Sehmi
    252. B Sennblad
    253. N Shah
    254. S Shah
    255. AS Shera
    256. XO Shu
    257. AR Shuldiner
    258. G Sigurđsson
    259. E Sijbrands
    260. A Silveira
    261. X Sim
    262. S Sivapalaratnam
    263. KS Small
    264. WY So
    265. A Stančáková
    266. K Stefansson
    267. G Steinbach
    268. V Steinthorsdottir
    269. K Stirrups
    270. RJ Strawbridge
    271. HM Stringham
    272. Q Sun
    273. C Suo
    274. AC Syvänen
    275. R Takayanagi
    276. F Takeuchi
    277. WT Tay
    278. TM Teslovich
    279. B Thorand
    280. G Thorleifsson
    281. U Thorsteinsdottir
    282. E Tikkanen
    283. J Trakalo
    284. E Tremoli
    285. MD Trip
    286. FJ Tsai
    287. T Tuomi
    288. J Tuomilehto
    289. AG Uitterlinden
    290. A Valladares-Salgado
    291. S Vedantam
    292. F Veglia
    293. BF Voight
    294. C Wang
    295. NJ Wareham
    296. R Wennauer
    297. AR Wickremasinghe
    298. T Wilsgaard
    299. JF Wilson
    300. S Wiltshire
    301. W Winckler
    302. TY Wong
    303. AR Wood
    304. JY Wu
    305. Y Wu
    306. K Yamamoto
    307. T Yamauchi
    308. M Yang
    309. L Yengo
    310. M Yokota
    311. R Young
    312. D Zabaneh
    313. F Zhang
    314. R Zhang
    315. W Zheng
    316. PZ Zimmet
    317. D Altshuler
    318. DW Bowden
    319. YS Cho
    320. NJ Cox
    321. M Cruz
    322. CL Hanis
    323. J Kooner
    324. JY Lee
    325. M Seielstad
    326. YY Teo
    327. M Boehnke
    328. EJ Parra
    329. JC Chambers
    330. ES Tai
    331. MI McCarthy
    332. AP Morris
    333. DIAbetes Genetics Replication And Meta-analysis (DIAGRAM) Consortium
    334. Asian Genetic Epidemiology Network Type 2 Diabetes (AGEN-T2D) Consortium
    335. South Asian Type 2 Diabetes (SAT2D) Consortium
    336. Mexican American Type 2 Diabetes (MAT2D) Consortium
    337. Type 2 Diabetes Genetic Exploration by Nex-generation sequencing in muylti-Ethnic Samples (T2D-GENES) Consortium
    (2014)
    Nature Genetics 46:234–244.
    https://doi.org/10.1038/ng.2897
  41. 41
  42. 42
  43. 43
  44. 44
    Large-scale association analysis provides insights into the genetic architecture and pathophysiology of type 2 diabetes
    1. AP Morris
    2. BF Voight
    3. TM Teslovich
    4. T Ferreira
    5. AV Segrè
    6. V Steinthorsdottir
    7. RJ Strawbridge
    8. H Khan
    9. H Grallert
    10. A Mahajan
    11. I Prokopenko
    12. HM Kang
    13. C Dina
    14. T Esko
    15. RM Fraser
    16. S Kanoni
    17. A Kumar
    18. V Lagou
    19. C Langenberg
    20. J Luan
    21. CM Lindgren
    22. M Müller-Nurasyid
    23. S Pechlivanis
    24. NW Rayner
    25. LJ Scott
    26. S Wiltshire
    27. L Yengo
    28. L Kinnunen
    29. EJ Rossin
    30. S Raychaudhuri
    31. AD Johnson
    32. AS Dimas
    33. RJ Loos
    34. S Vedantam
    35. H Chen
    36. JC Florez
    37. C Fox
    38. CT Liu
    39. D Rybin
    40. DJ Couper
    41. WH Kao
    42. M Li
    43. MC Cornelis
    44. P Kraft
    45. Q Sun
    46. RM van Dam
    47. HM Stringham
    48. PS Chines
    49. K Fischer
    50. P Fontanillas
    51. OL Holmen
    52. SE Hunt
    53. AU Jackson
    54. A Kong
    55. R Lawrence
    56. J Meyer
    57. JR Perry
    58. CG Platou
    59. S Potter
    60. E Rehnberg
    61. N Robertson
    62. S Sivapalaratnam
    63. A Stančáková
    64. K Stirrups
    65. G Thorleifsson
    66. E Tikkanen
    67. AR Wood
    68. P Almgren
    69. M Atalay
    70. R Benediktsson
    71. LL Bonnycastle
    72. N Burtt
    73. J Carey
    74. G Charpentier
    75. AT Crenshaw
    76. AS Doney
    77. M Dorkhan
    78. S Edkins
    79. V Emilsson
    80. E Eury
    81. T Forsen
    82. K Gertow
    83. B Gigante
    84. GB Grant
    85. CJ Groves
    86. C Guiducci
    87. C Herder
    88. AB Hreidarsson
    89. J Hui
    90. A James
    91. A Jonsson
    92. W Rathmann
    93. N Klopp
    94. J Kravic
    95. K Krjutškov
    96. C Langford
    97. K Leander
    98. E Lindholm
    99. S Lobbens
    100. S Männistö
    101. G Mirza
    102. TW Mühleisen
    103. B Musk
    104. M Parkin
    105. L Rallidis
    106. J Saramies
    107. B Sennblad
    108. S Shah
    109. G Sigurðsson
    110. A Silveira
    111. G Steinbach
    112. B Thorand
    113. J Trakalo
    114. F Veglia
    115. R Wennauer
    116. W Winckler
    117. D Zabaneh
    118. H Campbell
    119. C van Duijn
    120. AG Uitterlinden
    121. A Hofman
    122. E Sijbrands
    123. GR Abecasis
    124. KR Owen
    125. E Zeggini
    126. MD Trip
    127. NG Forouhi
    128. AC Syvänen
    129. JG Eriksson
    130. L Peltonen
    131. MM Nöthen
    132. B Balkau
    133. CN Palmer
    134. V Lyssenko
    135. T Tuomi
    136. B Isomaa
    137. DJ Hunter
    138. L Qi
    139. Wellcome Trust Case Control Consortium
    140. Meta-Analyses of Glucose and Insulin-related traits Consortium (MAGIC) Investigators
    141. Genetic Investigation of ANthropometric Traits (GIANT) Consortium
    142. Asian Genetic Epidemiology Network–Type 2 Diabetes (AGEN-T2D) Consortium
    143. South Asian Type 2 Diabetes (SAT2D) Consortium
    144. AR Shuldiner
    145. M Roden
    146. I Barroso
    147. T Wilsgaard
    148. J Beilby
    149. K Hovingh
    150. JF Price
    151. JF Wilson
    152. R Rauramaa
    153. TA Lakka
    154. L Lind
    155. G Dedoussis
    156. I Njølstad
    157. NL Pedersen
    158. KT Khaw
    159. NJ Wareham
    160. SM Keinanen-Kiukaanniemi
    161. TE Saaristo
    162. E Korpi-Hyövälti
    163. J Saltevo
    164. M Laakso
    165. J Kuusisto
    166. A Metspalu
    167. FS Collins
    168. KL Mohlke
    169. RN Bergman
    170. J Tuomilehto
    171. BO Boehm
    172. C Gieger
    173. K Hveem
    174. S Cauchi
    175. P Froguel
    176. D Baldassarre
    177. E Tremoli
    178. SE Humphries
    179. D Saleheen
    180. J Danesh
    181. E Ingelsson
    182. S Ripatti
    183. V Salomaa
    184. R Erbel
    185. KH Jöckel
    186. S Moebus
    187. A Peters
    188. T Illig
    189. U de Faire
    190. A Hamsten
    191. AD Morris
    192. PJ Donnelly
    193. TM Frayling
    194. AT Hattersley
    195. E Boerwinkle
    196. O Melander
    197. S Kathiresan
    198. PM Nilsson
    199. P Deloukas
    200. U Thorsteinsdottir
    201. LC Groop
    202. K Stefansson
    203. F Hu
    204. JS Pankow
    205. J Dupuis
    206. JB Meigs
    207. D Altshuler
    208. M Boehnke
    209. MI McCarthy
    210. DIAbetes Genetics Replication And Meta-analysis (DIAGRAM) Consortium
    (2012)
    Nature Genetics 44:981–990.
    https://doi.org/10.1038/ng.2383
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
    A language and environment for statistical computing
    1. R Development Core Team
    (2011)
    R Foundation for Statistical Computing, Vienna, Austria.
  52. 52
  53. 53
  54. 54
    An expanded genome-wide association study of Type 2 diabetes in europeans
    1. RA Scott
    2. LJ Scott
    3. R Mägi
    4. L Marullo
    5. KJ Gaulton
    6. M Kaakinen
    7. N Pervjakova
    8. TH Pers
    9. AD Johnson
    10. JD Eicher
    11. AU Jackson
    12. T Ferreira
    13. Y Lee
    14. C Ma
    15. V Steinthorsdottir
    16. G Thorleifsson
    17. L Qi
    18. NR Van Zuydam
    19. A Mahajan
    20. H Chen
    21. P Almgren
    22. BF Voight
    23. H Grallert
    24. M Müller-Nurasyid
    25. JS Ried
    26. NW Rayner
    27. N Robertson
    28. LC Karssen
    29. EM van Leeuwen
    30. SM Willems
    31. C Fuchsberger
    32. P Kwan
    33. TM Teslovich
    34. P Chanda
    35. M Li
    36. Y Lu
    37. C Dina
    38. D Thuillier
    39. L Yengo
    40. L Jiang
    41. T Sparso
    42. HA Kestler
    43. H Chheda
    44. L Eisele
    45. S Gustafsson
    46. M Frånberg
    47. RJ Strawbridge
    48. R Benediktsson
    49. AB Hreidarsson
    50. A Kong
    51. G Sigurðsson
    52. ND Kerrison
    53. J Luan
    54. L Liang
    55. T Meitinger
    56. M Roden
    57. B Thorand
    58. T Esko
    59. E Mihailov
    60. C Fox
    61. CT Liu
    62. D Rybin
    63. B Isomaa
    64. V Lyssenko
    65. T Tuomi
    66. DJ Couper
    67. JS Pankow
    68. N Grarup
    69. CT Have
    70. ME Jørgensen
    71. T Jørgensen
    72. A Linneberg
    73. MC Cornelis
    74. RM van Dam
    75. DJ Hunter
    76. P Kraft
    77. Q Sun
    78. S Edkins
    79. KR Owen
    80. JRB Perry
    81. AR Wood
    82. E Zeggini
    83. J Tajes-Fernandes
    84. GR Abecasis
    85. LL Bonnycastle
    86. PS Chines
    87. HM Stringham
    88. HA Koistinen
    89. L Kinnunen
    90. B Sennblad
    91. TW Mühleisen
    92. MM Nöthen
    93. S Pechlivanis
    94. D Baldassarre
    95. K Gertow
    96. SE Humphries
    97. E Tremoli
    98. N Klopp
    99. J Meyer
    100. G Steinbach
    101. R Wennauer
    102. JG Eriksson
    103. S Mӓnnistö
    104. L Peltonen
    105. E Tikkanen
    106. G Charpentier
    107. E Eury
    108. S Lobbens
    109. B Gigante
    110. K Leander
    111. O McLeod
    112. EP Bottinger
    113. O Gottesman
    114. D Ruderfer
    115. M Blüher
    116. P Kovacs
    117. A Tonjes
    118. NM Maruthur
    119. C Scapoli
    120. R Erbel
    121. KH Jöckel
    122. S Moebus
    123. U de Faire
    124. A Hamsten
    125. M Stumvoll
    126. P Deloukas
    127. PJ Donnelly
    128. TM Frayling
    129. AT Hattersley
    130. S Ripatti
    131. V Salomaa
    132. NL Pedersen
    133. BO Boehm
    134. RN Bergman
    135. FS Collins
    136. KL Mohlke
    137. J Tuomilehto
    138. T Hansen
    139. O Pedersen
    140. I Barroso
    141. L Lannfelt
    142. E Ingelsson
    143. L Lind
    144. CM Lindgren
    145. S Cauchi
    146. P Froguel
    147. RJF Loos
    148. B Balkau
    149. H Boeing
    150. PW Franks
    151. A Barricarte Gurrea
    152. D Palli
    153. YT van der Schouw
    154. D Altshuler
    155. LC Groop
    156. C Langenberg
    157. NJ Wareham
    158. E Sijbrands
    159. CM van Duijn
    160. JC Florez
    161. JB Meigs
    162. E Boerwinkle
    163. C Gieger
    164. K Strauch
    165. A Metspalu
    166. AD Morris
    167. CNA Palmer
    168. FB Hu
    169. U Thorsteinsdottir
    170. K Stefansson
    171. J Dupuis
    172. AP Morris
    173. M Boehnke
    174. MI McCarthy
    175. I Prokopenko
    176. DIAbetes Genetics Replication And Meta-analysis (DIAGRAM) Consortium
    (2017)
    Diabetes 66:2888–2902.
    https://doi.org/10.2337/db16-1253
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
    Twelve type 2 diabetes susceptibility loci identified through large-scale association analysis
    1. BF Voight
    2. LJ Scott
    3. V Steinthorsdottir
    4. AP Morris
    5. C Dina
    6. RP Welch
    7. E Zeggini
    8. C Huth
    9. YS Aulchenko
    10. G Thorleifsson
    11. LJ McCulloch
    12. T Ferreira
    13. H Grallert
    14. N Amin
    15. G Wu
    16. CJ Willer
    17. S Raychaudhuri
    18. SA McCarroll
    19. C Langenberg
    20. OM Hofmann
    21. J Dupuis
    22. L Qi
    23. AV Segrè
    24. M van Hoek
    25. P Navarro
    26. K Ardlie
    27. B Balkau
    28. R Benediktsson
    29. AJ Bennett
    30. R Blagieva
    31. E Boerwinkle
    32. LL Bonnycastle
    33. K Bengtsson Boström
    34. B Bravenboer
    35. S Bumpstead
    36. NP Burtt
    37. G Charpentier
    38. PS Chines
    39. M Cornelis
    40. DJ Couper
    41. G Crawford
    42. AS Doney
    43. KS Elliott
    44. AL Elliott
    45. MR Erdos
    46. CS Fox
    47. CS Franklin
    48. M Ganser
    49. C Gieger
    50. N Grarup
    51. T Green
    52. S Griffin
    53. CJ Groves
    54. C Guiducci
    55. S Hadjadj
    56. N Hassanali
    57. C Herder
    58. B Isomaa
    59. AU Jackson
    60. PR Johnson
    61. T Jørgensen
    62. WH Kao
    63. N Klopp
    64. A Kong
    65. P Kraft
    66. J Kuusisto
    67. T Lauritzen
    68. M Li
    69. A Lieverse
    70. CM Lindgren
    71. V Lyssenko
    72. M Marre
    73. T Meitinger
    74. K Midthjell
    75. MA Morken
    76. N Narisu
    77. P Nilsson
    78. KR Owen
    79. F Payne
    80. JR Perry
    81. AK Petersen
    82. C Platou
    83. C Proença
    84. I Prokopenko
    85. W Rathmann
    86. NW Rayner
    87. NR Robertson
    88. G Rocheleau
    89. M Roden
    90. MJ Sampson
    91. R Saxena
    92. BM Shields
    93. P Shrader
    94. G Sigurdsson
    95. T Sparsø
    96. K Strassburger
    97. HM Stringham
    98. Q Sun
    99. AJ Swift
    100. B Thorand
    101. J Tichet
    102. T Tuomi
    103. RM van Dam
    104. TW van Haeften
    105. T van Herpt
    106. JV van Vliet-Ostaptchouk
    107. GB Walters
    108. MN Weedon
    109. C Wijmenga
    110. J Witteman
    111. RN Bergman
    112. S Cauchi
    113. FS Collins
    114. AL Gloyn
    115. U Gyllensten
    116. T Hansen
    117. WA Hide
    118. GA Hitman
    119. A Hofman
    120. DJ Hunter
    121. K Hveem
    122. M Laakso
    123. KL Mohlke
    124. AD Morris
    125. CN Palmer
    126. PP Pramstaller
    127. I Rudan
    128. E Sijbrands
    129. LD Stein
    130. J Tuomilehto
    131. A Uitterlinden
    132. M Walker
    133. NJ Wareham
    134. RM Watanabe
    135. GR Abecasis
    136. BO Boehm
    137. H Campbell
    138. MJ Daly
    139. AT Hattersley
    140. FB Hu
    141. JB Meigs
    142. JS Pankow
    143. O Pedersen
    144. HE Wichmann
    145. I Barroso
    146. JC Florez
    147. TM Frayling
    148. L Groop
    149. R Sladek
    150. U Thorsteinsdottir
    151. JF Wilson
    152. T Illig
    153. P Froguel
    154. CM van Duijn
    155. K Stefansson
    156. D Altshuler
    157. M Boehnke
    158. MI McCarthy
    159. MAGIC investigators
    160. GIANT Consortium
    (2010)
    Nature Genetics 42:579–589.
    https://doi.org/10.1038/ng.609
  61. 61
  62. 62
  63. 63
    Ggplot2: Elegant Graphics for Data Analysis
    1. H Wickham
    (2009)
    New York: Springer-Verlag.
  64. 64
  65. 65

Decision letter

  1. Marcelo Nobrega
    Reviewing Editor; University of Chicago, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Integration of human pancreatic islet genomic data refines regulatory mechanisms at Type 2 Diabetes susceptibility loci" for consideration by eLife. Your article has been reviewed by two peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Aviv Regev as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

The reviewers of your manuscript found that the publication of a comprehensive epigenomic map of human islets is of broad interest and a timely resource for the community. The study has significant strengths, including the important observation of the far superior performance of WGBS compared to even the more recent, expanded 800K methylation array in comprehensively annotating methylation sites across the genome. Another significant strength highlighted by all reviewers is the integration of this dataset with other functional annotations and genetic fine mapping for T2D-associted loci, readily identifying a subset of variants that are putatively causal of the associations, and laying out a general analytical strategy that can be used by others, beyond those with immediate interest in T2D genetics. Overall, the reviewers found that this study significantly contributes to the effort of generating comprehensive functional annotations of a tissue of great medical and biological importance and thus far not properly characterized yet.

Below is a summary of essential points to be addressed:

1) The data presented starting in Figure 2B and shown throughout the manuscript suggest that virtually all enrichments are significant using FDR (masked by reporting only "FDR < 0.05" in main text). A more appropriate representation would either be quantitative (size proportional to log10p) or using a more stringent (Bonferroni) correction. The enrichments shown could be depicted in a more quantitative manner.

2) Causal GWAS variants can be enriched in regions with other functional annotations, such as coding regions or regions with strong evolutionary conservation. Such regions have not been considered here. The authors chose to reweight their PPAs solely using chromatin data, but this would down-weight causal variants with other functional annotations in the same region.

3) How much have the segmentations from the HMM improved with the additional data generated (without using GWAS enrichment as the benchmark)? It appears as though FGWAS enrichments for the actual annotations themselves (shown in Figure 3—figure supplement 1 and Figure 3—figure supplement 2) are just as good, if not better, than the HMM analysis (could also look at simple non-HMM two-annotation enrichments). Another way to look at this in the authors' framework would be to compare the median # of credible set SNPs (and top PPA) using only these annotations for FGWAS.

4) Both the Bayes Factor and FGWAS approaches rely on the assumption that at most one SNP in each region is causal, which may not always be true. This limitation should be discussed.

5) The addition of the Capture-C data in Figure 5C does not appear to help resolve rs11708067 as the likely causal variant as a several LD variants are similarly (by eye) overlapping fragments of high Capture-C density. Could this be quantified?

6) Given the known advantages of Stratified LD Score Regression (Finucane et al., 2015) over FGWAS to uncover true enrichments, have the authors considered using this approach?

7) What proportion of the narrow sense heritability is captured by the PPA variants?

8) The nomenclature used to define individual loci emerging from GWAS has followed the historical approach of using the closest gene to the associated SNPs. It has recently been recognized that this is often an imprecise approximation, as pointed out by reviewer 1. While there is no expectation that the authors would come up with a solution to this nomenclature conundrum, it was felt that the authors might address this issue on the paper, highlighting that the gene names associated with each locus may not accurately reflect the gene target of the association.

9) It was also noted that, despite the depth and complexity of the analysis carried out in this work, key data are missing which would greatly enhance the value of this report. An example is provided by Figure 4D: why not present a table giving the identified variants at each of the loci (or at least the top half?) with their causal probability? Presumably TCF7L2 refers to rs7903146? This should be clarified (as for the other loci included). It would also be nice to provide eQTL data for this and a few other chosen loci.

10) Finally, we would suggest that the authors consider bringing up the possibility that the association signal of several of the GWAS loci analyzed, including the ones highlighted in this work (see comment 5, above), might result not from the phenotypic impact of a single causal variant, but rather by the combined effects of multiple variants within the associated haplotype. This is certainly a topic that will become more obvious as we perform more detailed and high throughput computational and experimental analysis of GWAS loci, and it would be important for this work to come out already pondering about this possibility.

https://doi.org/10.7554/eLife.31977.046

Author response

Below is a summary of essential points to be addressed:

1) The data presented starting in Figure 2B and shown throughout the manuscript suggest that virtually all enrichments are significant using FDR (masked by reporting only "FDR < 0.05" in main text). A more appropriate representation would either be quantitative (size proportional to log10p) or using a more stringent (Bonferroni) correction. The enrichments shown could be depicted in a more quantitative manner.

Since, in the original analysis, only 1000 permutations had been performed, the lowest empirical P-value (0.001) did not lend itself to Bonferroni correction for the number of tests performed. Given the reviewers’ comment, we have increased the number of permutations to 100k to be able to estimate a more accurate empirical P-value and then applied a Bonferroni correction. This had no relevant impact on the findings since the empirical P-values calculated for most annotations were still less than 1/100,000 (equating to a minimum Bonferroni corrected P-value of 0.00032). All annotations that were significant using an FDR <0.05 remained significant using the Bonferroni correction (adjusted P<0.05) and the newly calculated enrichment estimates were in the same range as the original ones.

As requested we also represented the P-values more appropriately by providing the P-values in Figure 2B and the associated figure legend. In addition, we added Figure 2—figure supplement 1D that has been changed to show both the enrichment (x-axis) and associated –log10 Bonferroni adjusted P-value (y-axis) for a more quantitative presentation of the enrichment.

We updated the results in section 2.2 and modified the Materials and methods section 4.9.1 as well as Figure 2B and Figure 2—figure supplement 1D.

2) Causal GWAS variants can be enriched in regions with other functional annotations, such as coding regions or regions with strong evolutionary conservation. Such regions have not been considered here. The authors chose to reweight their PPAs solely using chromatin data, but this would down-weight causal variants with other functional annotations in the same region.

We agree with the reviewers’ comments that it is important to include a wide range of functional annotations. In fact, and contrary to the reviewer comments, we had already included coding regions in the original FGWAS analysis.

In response to the reviewers’ suggestion, we have now repeated the analyses adding in conserved region annotation. We used conserved sequence annotations from Lindblad-Toh et al., 2011 and updated the results of the paper accordingly.

Consistent with previous observations (e.g. Finucane et al., 2015) we found that conserved sequence was significantly enriched for T2D association in the single-state analysis. However, in the joint model, conserved sequence annotation was not retained presumably because the signal arising from conserved sequence was adequately captured by the chromatin state and coding annotations. Therefore, the addition of conserved sequence annotation does not change the major findings and conclusions.

To include the updated analysis results the manuscript has been revised as follows: We included the updated results in Figure 3, Figure 4, Figure 3—figure supplement 1, Figure 3—figure supplement2 as well as Table 12 and Figure 3—source data 13. In addition, we updated the results described in the main text (in particular, section 2.3 Refining islet enhancer function using methylation and open chromatin data).

3) How much have the segmentations from the HMM improved with the additional data generated (without using GWAS enrichment as the benchmark)? It appears as though FGWAS enrichments for the actual annotations themselves (shown in Figure 3—figure supplement 1 and Figure 3—figure supplement 2) are just as good, if not better, than the HMM analysis (could also look at simple non-HMM two-annotation enrichments). Another way to look at this in the authors' framework would be to compare the median # of credible set SNPs (and top PPA) using only these annotations for FGWAS.

While high levels of enrichment do not necessarily reflect high levels of sensitivity across a large number of GWAS loci, we agree that this is an interesting question. To compare the effect of different annotations to the baseline model, we have generated the requested metrics (median number and distribution of 99% credible set SNPs and median top variant PPA and top variant PPA distributions) for the following datasets as suggested:

- ChIP-only

- ChIP+Meth

- ChIP+ATAC

- ChIP+ATAC+Meth

- ATAC-only

- LMR-only

As shown in the new Figure 4—source data 1 and Figure 4—figure supplement 1 (see below for additional information), we found that ATAC-seq on its own (ATAC-only) or in combination with ChIP-seq data (ChIP+ATAC), leads to a notable reduction in 99% credible set size and an increase in top variant PPA compared to the baseline ChIP-only model.

In contrast, DNA methylation did not improve fine-mapping resolution: DNA methylation on its own (LMR-only) or in combination with ChIP-seq data, actually led to a reduction in top variant PPA compared to the baseline ChIP-only model. Credible set size increased when only DNA methylation (LMR-only) was included in the FGWAS fine-mapping analysis.

Despite the opposing effects of ATAC-seq and DNA methylation, the combination of ChIP-seq, ATAC-seq and DNA methylation (ChIP+ATAC+Meth) also performed well. The ChIP+ATAC+Meth model showed the best results in terms of median credible set size, identified a higher number of significant segments/loci and was almost identical in performance to the ChIP+ATAC dataset in terms of median maximum variant PPA.

Overall these results indicate that the chromatin accessibility data provided by ATAC-Seq (either in its own or in combination with islet ChIP-seq annotations) has the most value for refining fine-mapping and homing in on likely causal variants at T2D GWAS loci.

These results are highly consistent with our original analysis described in section 2.3, which showed that the effect of open chromatin on genetic T2D risk predominates over the effects of DNA methylation. The additional results are shown in the new and updated Figure 4—figure supplement 1, Supplementary Figure 4—source data 1 and referred to in the main text:

“We also expanded the FGWAS PPA analysis to investigate open chromatin and DNA methylation effects on fine-mapping and found that the reduction in 99% credible set size and increase in maximum variant PPA was driven predominantly by open chromatin (Figure 4—figure supplement 1). This demonstrates that the inclusion of open chromatin maps helps to improve prioritisation of causal variants at many T2D GWAS loci.”

4) Both the Bayes Factor and FGWAS approaches rely on the assumption that at most one SNP in each region is causal, which may not always be true. This limitation should be discussed.

We agree with the reviewer. We have been working to implement additional analyses steps to circumvent this limitation. Specifically, we have been testing an adaptation to FGWAS that allows us to analyse regions with multiple signals, after conditional decomposition. This is something we hope to deploy in future analyses, but it is not yet ready for inclusion in this paper. Meantime, and in line with reviewers’ comment 10, we have included the following paragraph in the Discussion highlighting this limitation of FGWAS and how future studies could address this:

“Although we provide highly detailed functional fine-mapping of T2D genetic variants to uncover causal variants, the FGWAS approach applied in this study is limited in its ability to determine the effect of multiple variants at individual loci. […] Analysis methods that combine functional fine-mapping with conditional analysis and consider LD and haplotype patterns are likely to provide a more complete overview of the causal interactions at T2D GWAS loci.”

5) The addition of the Capture-C data in Figure 5C does not appear to help resolve rs11708067 as the likely causal variant as a several LD variants are similarly (by eye) overlapping fragments of high Capture-C density. Could this be quantified?

This is correct. In our experience using Capture-C, it is typical that promoter “baits” capture multiple (enhancer) sequences in the adjacent region. This is not unexpected given the strong correlation of functional activity within adjacent sequence (as reflected in TAD structure for example). We agree, therefore, that it is generally not possible to use capture-C in isolation to definitively resolve the mechanistic interactions at a GWAS signal. At the ADCY5 locus, the capture-C data we present demonstrate physical approximation between rs11708067 and the ADCY5 promoter in an appropriate cell type: this provides additional evidence (to be considered in tandem with other complementary data) to support their relationship. At ADCY5, we are able to push the interpretation a little further in terms of resolving the causal variant because rs11708067 is the sole variant in the enhancer region that shows this physical approximation to the promoter.

To support more rigorous assessments of Capture-C data, we have been implementing a method for quantitative detection of Capture-C peaks which takes into account the background distribution of interaction “noise”. We have quantified the evidence for interaction for a total of 12 fragments across the ADCY5 region which contain rs11708067 or other close T2D-associated proxies. In this analysis, only two fragments (one of which contains rs11708067) are called as “peaks” on the basis of significant read number over background.

These updated results are shown in the newly added Figure 5—figure supplement 1D and in the main text as shown below:

“To resolve the significance of the interaction between the restriction fragment encompassing rs11708067 and the ADCY5 promoter, we used the programme peakC (de Wit and Geeven et al., 2017) (https://github.com/deWitLab/peakC) to evaluate the interactions of 12 fragments covering the lead SNP rs11708067 and 15 SNPs in high LD (r2 > 0.8) across a region of 47kb. […] These SNPs fall into a region that did not show evidence of open chromatin.”

In addition, we also added the following paragraph to the Materials and methods section to include the additional analysis:

“Significant chromatin interactions were determined from a rank-sum test implemented in the program peakC (de Wit and Geeven et al., 2017)(https://github.com/deWitLab/peakC). […] We applied the Benjamini-Hochberg correction to control the false discovery rate for the set of p-values corresponding to each restriction fragment within the 47kb region at the ADCY5 locus.”

6) Given the known advantages of Stratified LD Score Regression (Finucane et al., 2015) over FGWAS to uncover true enrichments, have the authors considered using this approach?

We agree with the authors that Stratified LD Score Regression could be a valuable method to run on our dataset, in particular, since it can include multiple causal signals at a locus. However, after considering the use of multiple methods (Garfield, Stratified LD Score regression, FGWAS) we ultimately decided to use FGWAS for the following reasons:

- The FGWAS framework provides a stringent and straightforward way to determine a maximum likelihood model with multiple annotations using both forward and reverse mode selection.

- The framework provides a simple way of calculating updated posterior probabilities based on the enrichment priors of the best maximum likelihood annotation model. Importantly, the joint maximum likelihood model provides relative enrichment estimates for each category that allow us to up-weight the PPA of variants overlapping a given annotation in a proportional way dependent on the annotation.

- FGWAS is not dependent on reference datasets for LD-scores and can use all variants tested in the 1000G imputed DIAGRAM dataset including low frequency and rare variants.

- The observed enrichment patterns (highest enrichment for active enhancer states) are highly consistent with previous studies (e.g. Parker et al., 2013, Pasquali et al., 2014 and Gaulton et al., 2015).

We will continue to review these different approaches and will seek to implement other methods – such as LD score regression – in future work.

7) What proportion of the narrow sense heritability is captured by the PPA variants?

We have estimated the narrow sense heritability for T2D using UK Biobank data from ~350k individuals. To avoid inflation of the heritability estimates by variants in high LD, we restricted the analysis to the (single) top PPA variant for each of the 52 segments that were significantly associated with T2D genome-wide. The narrow sense heritability for T2D was estimated at 1.8%. To provide a reference point, we also estimated that all currently published T2D GWAS variants (Scott et al., 2017, Mahajan et al., 2014, Morris et al., 2012, Voight et al., 2010) explain 5.7% of variation in T2D risk. We have also used the top associated FG variants at 20 genome-wide significant segments to calculate narrow sense heritability estimates for FG, using Oxford BioBank data from 3.6k individuals. The narrow sense heritability for FG was estimated at 1.9%.

8) The nomenclature used to define individual loci emerging from GWAS has followed the historical approach of using the closest gene to the associated SNPs. It has recently been recognized that this is often an imprecise approximation, as pointed out by reviewer 1. While there is no expectation that the authors would come up with a solution to this nomenclature conundrum, it was felt that the authors might address this issue on the paper, highlighting that the gene names associated with each locus may not accurately reflect the gene target of the association.

We agree! We have added the following sentence:

“(Of note, in line with traditional GWAS nomenclature, locus names were defined based on proximity between the lead variant and the closest gene and does not, of itself, indicate any causal role for the gene in T2D susceptibility).”

9) It was also noted that, despite the depth and complexity of the analysis carried out in this work, key data are missing which would greatly enhance the value of this report. An example is provided by Figure 4D: why not present a table giving the identified variants at each of the loci (or at least the top half?) with their causal probability? Presumably TCF7L2 refers to rs7903146? This should be clarified (as for the other loci included). It would also be nice to provide eQTL data for this and a few other chosen loci.

- As requested we have created an additional table (Figure 4—source data 2) that provides information for all 190 variants from Figure 4D which overlap one of the annotations included in the FGWAS joint-model.

- Per variant, the following T2D/FGWAS information is provided: rsID; FGWAS PPA; T2D GWAS P-value; FGWAS segment number; T2D locus name; tested for allelic imbalance (Yes/No).

- If available, the following eQTL information (from Varshney et al., 2017) is also provided:

eQTL allele1 (effector), eQTL allele 2, eQTL q-value, eQTL effect and eQTL gene

10) Finally, we would suggest that the authors consider bringing up the possibility that the association signal of several of the GWAS loci analyzed, including the ones highlighted in this work (see comment 5, above), might result not from the phenotypic impact of a single causal variant, but rather by the combined effects of multiple variants within the associated haplotype. This is certainly a topic that will become more obvious as we perform more detailed and high throughput computational and experimental analysis of GWAS loci, and it would be important for this work to come out already pondering about this possibility.

In line with the reviewers’ comment 4 about the limitation of FGWAS to assume only single causal variants, we have added a paragraph in the Discussion describing the potential impact of multiple signals per locus and interacting causal variants on the phenotype and how this could be addressed. The part of the paragraph that highlights these effects is shown below:

“Specifically, FGWAS relies on the assumption of a single causal variant within each region, which may not necessarily be true for all loci. […] Analysis methods that combine functional fine-mapping with conditional analysis and consider LD and haplotype patterns are likely to provide a more complete overview of the causal interactions at T2D GWAS loci.”

https://doi.org/10.7554/eLife.31977.047

Article and author information

Author details

  1. Matthias Thurner

    1. The Wellcome Centre for Human Genetics, University of Oxford, Oxford, United Kingdom
    2. Oxford Centre for Diabetes, Endocrinology and Metabolism, University of Oxford, Oxford, United Kingdom
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing—original draft
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-7329-9769
  2. Martijn van de Bunt

    1. The Wellcome Centre for Human Genetics, University of Oxford, Oxford, United Kingdom
    2. Oxford Centre for Diabetes, Endocrinology and Metabolism, University of Oxford, Oxford, United Kingdom
    Contribution
    Conceptualization, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-6744-6125
  3. Jason M Torres

    The Wellcome Centre for Human Genetics, University of Oxford, Oxford, United Kingdom
    Contribution
    Formal analysis, Validation, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  4. Anubha Mahajan

    The Wellcome Centre for Human Genetics, University of Oxford, Oxford, United Kingdom
    Contribution
    Resources, Formal analysis, Supervision, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  5. Vibe Nylander

    Oxford Centre for Diabetes, Endocrinology and Metabolism, University of Oxford, Oxford, United Kingdom
    Contribution
    Data curation, Formal analysis, Validation, Investigation, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  6. Amanda J Bennett

    Oxford Centre for Diabetes, Endocrinology and Metabolism, University of Oxford, Oxford, United Kingdom
    Contribution
    Resources, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  7. Kyle J Gaulton

    Department of Pediatrics, University of California, San Diego, San Diego, United States
    Contribution
    Conceptualization, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  8. Amy Barrett

    Oxford Centre for Diabetes, Endocrinology and Metabolism, University of Oxford, Oxford, United Kingdom
    Contribution
    Resources, Methodology
    Competing interests
    No competing interests declared
  9. Carla Burrows

    Oxford Centre for Diabetes, Endocrinology and Metabolism, University of Oxford, Oxford, United Kingdom
    Contribution
    Resources, Methodology
    Competing interests
    No competing interests declared
  10. Christopher G Bell

    1. Department of Twin Research and Genetic Epidemiology, Kings College London, London, United Kingdom
    2. MRC Lifecourse Epidemiology Unit, University of Southampton, Southampton, United Kingdom
    Contribution
    Data curation, Writing—review and editing
    Competing interests
    No competing interests declared
  11. Robert Lowe

    Centre for Genomics and Child Health, Blizard Institute, Barts and The London School of Medicine and Dentistry, London, United Kingdom
    Contribution
    Data curation, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  12. Stephan Beck

    Department of Cancer Biology, UCL Cancer Institute, University College London, London, United Kingdom
    Contribution
    Conceptualization, Resources, Data curation, Writing—review and editing
    Competing interests
    No competing interests declared
  13. Vardhman K Rakyan

    Centre for Genomics and Child Health, Blizard Institute, Barts and The London School of Medicine and Dentistry, London, United Kingdom
    Contribution
    Conceptualization, Resources, Data curation, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
  14. Anna L Gloyn

    1. The Wellcome Centre for Human Genetics, University of Oxford, Oxford, United Kingdom
    2. Oxford Centre for Diabetes, Endocrinology and Metabolism, University of Oxford, Oxford, United Kingdom
    3. Oxford NIHR Biomedical Research Centre, Churchill Hospital, Oxford, United Kingdom
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Project administration, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1205-1844
  15. Mark I McCarthy

    1. The Wellcome Centre for Human Genetics, University of Oxford, Oxford, United Kingdom
    2. Oxford Centre for Diabetes, Endocrinology and Metabolism, University of Oxford, Oxford, United Kingdom
    3. Oxford NIHR Biomedical Research Centre, Churchill Hospital, Oxford, United Kingdom
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Project administration, Writing—review and editing
    For correspondence
    mark.mccarthy@drl.ox.ac.uk
    Competing interests
    Senior editor, eLife
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4393-0510

Funding

Novo Nordisk

  • Martijn van de Bunt

Horizon 2020 Framework Programme (HEALTH-F4-2007-201413)

  • Vibe Nylander
  • Anna L Gloyn

Royal Society

  • Stephan Beck

National Institute for Health Research

  • Anna L Gloyn
  • Mark I McCarthy

National Institutes of Health (U01-DK105535)

  • Anna L Gloyn
  • Mark I McCarthy

Wellcome Trust (099673/Z/12/Z)

  • Matthias Thurner

National Institutes of Health (R01-DK098032)

  • Anna L Gloyn
  • Mark I McCarthy

Wellcome Trust (106130)

  • Jason M Torres
  • Anna L Gloyn
  • Mark I McCarthy

National Institutes of Health (R01-MH090941)

  • Anna L Gloyn
  • Mark I McCarthy

Wellcome Trust (095101/Z/10/Z)

  • Anna L Gloyn

Wellcome Trust (200837/Z/16/Z)

  • Anna L Gloyn

Wellcome Trust (090532)

  • Anna L Gloyn
  • Mark I McCarthy

Wellcome Trust (098381)

  • Mark I McCarthy

National Institutes of Health (U01-DK085545)

  • Mark I McCarthy

Wellcome Trust (090367)

  • Mark I McCarthy

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank the High-Throughput Genomics Group at the Wellcome Centre for Human Genetics (funded by Wellcome grant reference 090532) for the generation of the Sequencing data. MT was supported by a Wellcome Doctoral Studentship. MvdB was supported by a Novo Nordisk postdoctoral fellowship run in partnership with the University of Oxford. JMT is supported by Wellcome Strategic Award. VN is funded by the European Union Horizon 2020 Programme (T2DSYSTEMS). SB was supported by a Royal Society Wolfson Research Merit Award (WM100023). ALG is a Wellcome Senior Fellow in Basic Biomedical Science (095101/Z/10/Z and 200837/Z/16/Z) and MIM is a Wellcome Senior Investigator. The research was supported by the National Institute for Health Research (NIHR) Oxford Biomedical Research Centre (BRC). This work was also supported by EU (HEALTH-F4-2007-201413), Wellcome (090367, 090532, 106130, 098381 and 099673/Z/12/Z) and NIH (U01-DK105535, U01-DK085545, R01-DK098032 and R01-MH090941) grants. The views expressed are those of the author(s) and not necessarily those of the NHS, the NIHR or the Department of Health.

Ethics

Human subjects: The Human Research Ethics Board at the University of Alberta (Pro00001754), the University of Oxford's Oxford Tropical Research Ethics Committee (OxTREC Reference: 2-15), or the Oxfordshire Regional Ethics Committee B (REC reference: 09/H0605/2) approved the studies. All organ donors provided informed consent for use of pancreatic tissue in research.

Reviewing Editor

  1. Marcelo Nobrega, University of Chicago, United States

Publication history

  1. Received: September 13, 2017
  2. Accepted: February 6, 2018
  3. Accepted Manuscript published: February 7, 2018 (version 1)
  4. Version of Record published: February 27, 2018 (version 2)

Copyright

© 2018, Thurner et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 2,088
    Page views
  • 366
    Downloads
  • 7
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Developmental Biology
    2. Genetics and Genomics
    Anne C Ferguson-Smith, Deborah Bourchis
    Feature Article
    1. Computational and Systems Biology
    2. Genetics and Genomics
    David G Hendrickson et al.
    Tools and Resources