Study participants

Participants were drawn from the Pancreatic Cancer Cohort Consortium (PanScan) and the Pancreatic Cancer Case-Control Consortium (PanC4) and individuals were included from 16 cohort and 13 case–control studies genotyped in four previous GWAS phases, namely PanScan I, PanScan II, PanScan III, and PanC412,13,14,16. Samples from the PANDoRA case–control consortium were used for replication23. The details on cases (individuals with PDAC) and controls have been previously described12,13,14,16.

All studies obtained informed consent from study participants and Institutional Review Board (IRB) approvals including IRB certifications permitting data sharing in accordance with the NIH Policy for Sharing of Data Obtained in NIH Supported or Conducted GWAS. The PanScan and PanC4 GWAS data are available through dbGAP (accession numbers phs000206.v5.p3 and phs000648.v1.p1, respectively).

Genotyping, imputation, and association analysis

Genotyping for PanScan was performed at the Cancer Genomics Research Laboratory (CGR) of the National Cancer Institute (NCI) of the National Institutes of Health (NIH) using the Illumina HumanHap series arrays (Illumina HumanHap550 Infinium II, Human 610-Quad) for PanScan I and II, respectively, and the Illumina Omni series arrays (OmniExpress, Omni1M, Omni2.5, and Omni5M) for PanScan III12,13,14. Genotyping for the PanC4 GWAS was performed at the Johns Hopkins Center for Inherited Disease Research (CIDR) using the Illumina HumanOmniExpressExome-8v1 array. Imputation was performed using the 1000 Genomes (1000 G) Phase 3, Release 1 reference data set22 and IMPUTE2 (http://mathgen.stats.ox.ac.uk/impute/impute_v2.html)60 as previously described14,18. Because of the large overlap of variants on genotyping arrays for PanScan I and II, these data sets were imputed and analyzed together. The PanScan III and PanC4 GWAS data sets were each imputed and analyzed separately. For quality control, variants were excluded based on (1) completion rate <90%; (2) MAF <0.01; and (3) low-quality imputation score (IMPUTE2 INFO score <0.3). After quality control, 11,381,182 SNPs genotyped or imputed in 5107 pancreatic cancer patients and 8845 controls of European ancestry were included in the analysis for PanScan I-III and 3933 cases and 3651 controls for PanC4 (Supplementary Table 1). The association analysis was performed using SNPTEST (http://mathgen.stats.ox.ac.uk/genetics_software/snptest/snptest.html)61 based on probabilistic genotype values from IMPUTE260, with parallel covariate adjustments: study, geographical region, age, sex, and principal components (PCA) of population substructure as were used in PanScan12,13,14 and study, age, sex, and PCA population substructure as were used in PanC416. The score statistic of the log additive genetic association magnitude was used. Summary statistics from PanScan I and II, PanScan III, and PanC4 were used for a meta-analysis using the fixed-effects inverse-variance method based on β estimates and SEs (http://csg.sph.umich.edu/abecasis/metal/). Heterogeneity was not observed for the SNPs identified as GWAS significant or suggestive in the combined study (P heterogeneity ≥ 0.32). IMPUTE2 information scores were above 90% for SNPs (P < 1 × 10−6), except for rs13303010 in the PanScan I+II data (INFO = 0.42; Table 1). The estimated inflation of the test statistic, λ, was 1.002 for PanScan I+II, 1.051 for PanScan III, and 1.024 for PanC462.

A Polygenic Risk Score (PRS) was constructed for each individual by summing the number of risk alleles carried for all established pancreatic cancer risk loci identified by GWAS, weighted by their estimated effect size. Individuals were grouped by percentiles, and the association of the PRS (as percentile groupings) with pancreatic cancer was estimated using logistic regression.

DEPICT analysis was used to prioritize causal genes at currently known pancreatic cancer risk loci12,13,14,15,16 identify gene sets enriched across risk loci, and tissues in which genes at risk loci are highly expressed57. No genes or SNP–gene pairs were significant at false discovery rate (FDR) < 0.05 (data not shown) but significant tissue enrichment is shown in Supplementary Table 7. Additional gene set enrichment analysis was performed using GeneCodis3 (http://genecodis.cnb.csic.es/). Genes (n = 65) were located in the currently known pancreatic cancer risk loci identified in subjects of European descent (for genes located +/−100 kb from the most significant SNP at the 22 risk loci) based on KEGG, Gene Ontology (Biological Process and Molecular Function) annotations using GeneCodis3 with reporting of FDR-corrected hypergeometric P values (Supplementary Table 8)56.

Replication

Ten promising signals (P < 10−6) were selected for replication in samples from the PANDoRA consortium23. Genotyping was performed by custom TaqMan genotyping assays (Applied Biosystems) at the German Cancer Research Center (DKFZ) in Heidelberg, Germany, for 2770 pancreatic cancer patients and 5178 controls, of which 2737 cases and 4753 controls had complete age and clinical data and did not overlap with other study individuals. Duplicate quality-control samples (n = 607 pairs) showed 99.48% genotype concordance. SNP quality metrics were performed for each SNP by plate; plates with <80% genotype completion rates were dropped from the analysis. Individuals were excluded if they were missing data on two or more SNPs after excluding SNPs on plates with low genotype completion rates. The association analysis for PANDoRA was adjusted for age and study in the same manner as previously described14,15,16. Heterogeneity between studies was assessed using the Cochran’s Q-test. Association analysis was also performed for the set of variants previously replicated in PANDoRA as part of the PanScan III14 and PanC416 GWAS studies (Supplementary Table 3).

Using SequenceLDhot, recombination hotspots for association plots were generated as previously described12,13,14. Recombination hotspot inference was performed using the 1000 G CEU samples (n = 99). The LD heatmap was prepared using the 1000 G Phase 3 CEU data, and the snp.plotter R software package63.

eQTL analysis

The publicly available GTEx data64 (http://www.gtexportal.org/; version 6) were used to assess eQTLs in pancreatic tissue samples (n = 149). RefSeq genes located within +/−500 kb of the marker SNP for each GWAS significant locus were assessed for cis-eQTL effects. Nominally significant eQTLs from this analysis (P < 0.05) were then taken forward to further analysis in two additional pancreatic tissue sample sets58: (1) the LTG and (2) The Cancer Genome Atlas (TCGA) pancreatic adenocarcinoma (PAAD) samples.

The LTG set included 95 histologically normal pancreatic tissue samples from participants of European ancestry collected at three participating sites: Mayo Clinic in Rochester, MN (45 samples, adjacent to tumor); Memorial Sloan Kettering Cancer Center in New York City, NY (34 samples, adjacent to tumor); and Penn State College of Medicine, Hershey, PA (16 samples, from tissue donors via the Gift of Life Donor Program, Philadelphia, PA) as previously described58. Samples were confirmed to be non-tumorous with ≥80% epithelial component by histological review and macro-dissected when needed. The project was approved by the Institutional Review Board of each participating institution as well as the NIH.

In short, RNA (RIN >7.5) isolated from fresh frozen histologically normal pancreatic tissue samples (LTG samples) with the Ambion mirVana kit was poly-A-enriched and subjected to massively parallel paired-end sequencing (Illumina’s HiSeq2000/TruSeq v3 sequencing). MapSplice was used to align reads and RSEM (v1.2.14) for gene expression quantification (TPM)65,66 using the hg19/GRCh37-based UCSC “RefSeq” track for gene annotation. DNA for genotyping was isolated from blood (Mayo Clinic samples), histologically normal fresh frozen pancreatic tissue samples (Penn State samples), or histologically normal fresh frozen spleen or duodenum tissue samples (MSKCC samples) using the Gentra Puregene Tissue Kit (Qiagen). DNA samples were genotyped on the Illumina OmniExpress or Omni1M arrays at the CGR of the Division of Cancer Epidemiology and Genetics, NCI, NIH. After quality control, genotypes were imputed using the 1000 G (Phase 1, v3) imputation reference data set and IMPUTE 267. Pre-imputation exclusion filters of Hardy Weinberg Equilibrium P < 1 × 10−6, minor allele frequency (MAF) <0.01, genotype missing rate >0.05, A/T and G/C pairs on ambiguous DNA strand (MAF > 0.45), and significantly different allele frequency between sample data and the 1000 G reference data (P < 7 × 10−8: Fisher’s exact test) were used. Post-imputation variants (single-nucleotide variants (SNP) and small insertion-deletion (indel) polymorphisms) with MAF < 0.05 or imputation quality scores (INFO score) <0.5 were removed from the final analysis58.

The second sample set included RNA-sequencing (RNA-seq) and genotype data from tumor-derived pancreatic samples obtained from TCGA PAAD data set by permission through the TCGA Data Access Committee. We excluded samples of non-European ancestry, with history of neo-adjuvant therapy prior to surgery, or with histological subtypes other than PDAC, leaving a total of 115 tumor samples for analysis58. TCGA mRNA-seq data (level 1 read data, generated using Illumina’s HiSeq2000) for pancreatic cancer samples (TCGA PAAD) were processed in the same manner as the histologically normal LTG samples described above. Blood-derived DNA samples for TCGA PAAD samples were genotyped on Affymetrix 6.0 arrays and processed in the same manner as for the LTG samples58.

The eQTL analysis was performed separately in histologically normal (LTG) and tumor-derived (TCGA PAAD) pancreatic samples using the Matrix eQTL (http://www.bios.unc.edu/research/genomic_software/Matrix_eQTL/) software package68. We tested associations between genotyped and imputed SNPs and the expression of genes evaluated by mRNA-sequencing after upper quantile normalization within samples and normal quantile transformation for each gene across samples by regressing the imputed dosage of the minor allele for each variant against normalized gene expression values68. Linear models were adjusted for age, sex, study, and the top five principal components (PCs) each for genotypes and gene expression to account for possible measured or hidden confounders58. The T-statistics from the linear regression is reported. For the tumor samples, we further adjusted for tumor stage and sequence-based tumor purity as per information provided by TCGA.

Bioinformatic analysis of functional potential

Variants at the new risk loci were assessed for potential functionality by examining their location in open (DNase Hypersensitivity Regions, DHS) and active chromatin (as per promoter and enhancer histone modification marks) in the ENCODE data. For this, we used HaploReg 4.1 (http://www.broadinstitute.org/mammals/haploreg/haploreg.php) and the UCSC Genome Browser (http://genome.ucsc.edu/). Variants correlated with the most significant variant at each locus at r2 > 0.7 in 1000 G EUR populations were included, except for 1p36.33/NOC2L where P value LR > 1:100 was used.

Candidate functional variants at 1p36.33 were selected by comparing the likelihood of each variant from the association analysis with the likelihood of the most significant variant. Ten variants had likelihood ratios, LR > 1:100 relative to the most significant SNP: rs13303010 (P value rank 1), rs13303327 (rank 2), rs113491766 (rank 3), rs3935066 (rank 4), rs111748052 (rank 5), rs10465242 (rank 6), rs13303160 (rank 7), rs7524174 (rank 8), rs13302957 (rank 9), and rs4970445 (rank 10). They were all highly correlated with rs13303010 (r2 = 0.52–1.00, 1000 G EUR data). These 10 variants were considered the set of variants most likely to contain the functional variant(s) at 1p36.33. Possible allelic effects of these top 10 variants on TF-binding motifs were determined using PrEdict Regulatory Functional Effect of SNPs by Approximate P value Estimation (PERFECTOS-APE; http://opera.autosome.ru/perfectosape/) analysis that determines the probability of a TF motif (using position weight matrices, from HOCOMOCO-10, JASPAR, HT-SELEX, SwissRegulon, and HOMER databases) in the DNA sequence overlapping each variant. The fold change in probability of there being a TFBS present for each allele of a variant is then calculated69. Two dbSNP variants at 1p36.33 with the same bp location, rs111748052 (−/ATTTT) and rs10465241 (C/T), may be two independent variants (as indicated in dbSNP), or a single tri-allelic marker (alleles are shown as C/CATTTT/T in 1000 G). As PERFECTOS-APE does not analyze indel variants, we analyzed the two indels among the 10 potential functional variants, rs111748052 and rs113491766, using a different program, sequence TF Affinity Prediction. This program calculates the total affinity of a sequence for a TF (as given by TRANSFAC and JASPAR databases) on the basis of a biophysical model of the binding energies between a TF and DNA70. The probability for a given TFBS for each variant of the indel was then compared as in PERFECTOS-APE to determine the fold change effect of the indel on the presence of the TFBS.

Gene expression analysis

Gene expression was assessed for genes that are closest to the reported variants at chromosomes 1p36.33: NOC2L, KLHL17, and PLEKHN1; 7p12: TNS3; 8q21.11: HNF4G; 17q12: HNF1B; 18q21.32: GRP, as well as two additional genes at 1p36.33 exhibiting nominally significant eQTLs in GTEx (1p36.33/SAMD11/DVL1). We assessed differential expression of these genes in pancreatic tumor samples (PDAC, n = 8), histologically normal (non-malignant) pancreatic tissue samples (n = 10), and pancreatic cell lines (n = 9) by RNA-seq as described previously44. We compared gene expression in tumors (T) and cell lines (C) to histologically normal pancreatic tissue samples (N) by EdgeR analysis. P values for differential expression in tumor vs. normal (TvN) and Cell lines vs. normal (CvN) represent an exact statistic using the normalized pseudo-counts and tagwise dispersion estimates per gene.

Data availability

The PanScan and PanC4 genome-wide association data that support the findings of this study are available through dbGAP (accession numbers phs000206.v5.p3 and phs000648.v1.p1, respectively).