Genome-wide meta-analysis identifies five new susceptibility loci for cutaneous malignant melanoma.

The samples were genotyped on a variety of commercial arrays, detailed in the Supplementary Note . All samples were collected with informed consent and with ethical approval (for details see Supplementary Note ). Imputation was conducted genome-wide, separately on each study, following a shared protocol. SNPs with MAF < 0.03 (MAF < 0.01 in AMFS, Q-MEGA_omni, Q-MEGA_610k, WAMHS, and HEIDELBERG), control Hardy-Weinberg equilibrium (HWE) P 0.03 were excluded, as were any individuals with call rates <0.97, identified as first degree relatives and/or European outliers by principal components analysis using EIGENSTRAT 71 . In addition, in each study where genotyping was conducted on more than one chip, any SNP not present on all chips was removed prior to imputation to avoid bias. IMPUTEv2.2 72 , 73 was used for imputation for all studies but Harvard, which used MaCH 74 , 75 and MDACC which used MaCH and minimac 76 . For GenoMEL, CIDRUK and MDACC samples the 1000 Genomes Feb 2012 data (build 37) was used as the reference panel, while for the AMFS, Q-MEGA_omni, Q-MEGA_610k, WAMHS, MELARISK and HEIDELBERG datasets the 1000 Genomes April 2012 data (build 37) was the reference for imputation 77 . In both cases any SNP with MAF < 0.001 in European (CEU) samples was dropped from the reference panel. The HARVARD data were imputed using MACH with the NCBI build 35 of phase II HapMap CEU data as the reference panel and only SNPs with imputation quality R 2 > 0.95 were included in the final analysis. Imputed genotypes were analyzed as expected genotype counts based on posterior probabilities (gene dosage) using SNPTEST2 78 assuming an additive model with geographic region as a covariate (SNPTEST v2.5 for chromosome X). MDACC imputed genotypes were analyzed using best guess genotypes from MACH and PLINK was used for logistic association test adjusting by the top two principle components. Only those with an imputation quality score (INFO/MaCH R 2 ) score >0.8 were analyzed. Potential stratification was dealt with in the GenoMEL samples by including geographic region as a covariate (inclusion of principal components as covariates was previously found to make little difference 9 ) and elsewhere by including principal components as covariates 71 . Heterogeneity of per-SNP effect sizes in studies contributing to the Stage one, Stage two and the Overall meta-analyses was assessed using the I 2 metric 79 . I 2 is commonly defined as the proportion of overall variance attributable to between-study variance, with values below 31% suggesting no more than mild heterogeneity. Where I 2 was less than 31% a fixed effects model was used, with fixed effects P -values indicated by P fixed ; otherwise random effects were applied ( P random ). The method of Dersimonian and Laird 80 was used to estimate the between-studies variance, τ ^ 2 . An overall random effects estimate was then calculated using the weights 1 ∕ ( v i + τ ^ 2 ) where v i is the variance of the estimated effect. τ ^ 2 = 0 for the fixed effects analyses. We report those loci reaching significance at > one marker incorporating information from > one study. Results for rs186133190 in 2p15 were only available from four studies; all other SNPs reported here utilize data from at least eight studies ( Supplementary Table 3 ). Per-study QQ plots of GWAS P -values are provided ( Supplementary Figure 3 ) and for the Stage one meta-analysis ( Supplementary Figure 3A ). We also provide the Stage one QQ plot with previously reported regions removed ( Supplementary Figure 2B ). While there was minimal inflation remaining following PC/region of origin correction, to ensure residual genomic inflation was not biasing our results the meta-analysis was repeated using the genome-wide association meta-analysis software, GWAMA v2.1 81 . Included studies were corrected by inflating SNP variance estimates by their genomic inflation (λ). As expected, given the low level of residual inflation, corrected and uncorrected results were very similar. Genome-wide inflation factor (λ)-corrected P- values are provided in Supplementary Table 3 . A single SNP for each of the 16 novel region reaching P < 10 –6 in Stage one was subsequently genotyped in 3 additional melanoma case-control sets in Stage two ( Supplementary Table 3 ). Any region that only showed evidence for association with CMM at a single imputed SNP and in only one study was not followed up. Included in the Stage two genotyping were rs75691080 in 20q13.33 which, while not quite reaching P < 10 –6 , lies 20 kb from RTEL1 ; and rs12696304 in 3q26.2 which lies 1 kb from TERC . Both these genes are known to be telomere-related and have been associated with leukocyte telomere length 35 . The 16 novel regions included rs2290419 at 11q13.3 which is 450 kb away from our primary hit in the region of CCND1 (rs498136, Supplementary Figure 5 and Supplementary Figure 6 ) and is in linkage equilibrium with the genome-wide significant hit in this region (r 2 = 0.002) so may represent an independent effect. The first Stage two dataset of 1,797 cases and 1,709 controls from two studies based in Cambridge, UK (see Supplementary Note for details of samples). These were genotyped using TaqMan® assays [Applied Biosystems]. 2 μl PCR reactions were performed in 384 well plates using 10 ng of DNA (dried), using 0.05 μl assay mix and 1 μl Universal Master Mix [Applied Biosystems] according to the manufacturer's instructions. End point reading of the genotypes was performed using an ABI 7900HT Real-time PCR system [Applied Biosystems]. The second was 711 cases and 890 controls from the Breakthrough Generations Study. These were genotyped in the same way as the Cambridge Stage two samples above. The third was 800 cases and 800 controls from Athens, Greece. Genomic DNA was isolated from 200 μl peripheral blood using the QIAamp DNA blood mini kit [Qiagen]. DNA concentration was quantified in samples prior to genotyping by using Quant-iT dsDNA HS Assay kit [Invitrogen]. The concentration of the DNA was adjusted to 5 ng/μl. Selected SNPs were genotyped using the Sequenom iPLEX assay [Sequenom]. Allele detection in this assay was performed using matrix-assisted laser desorption/ionization –time-of-flight mass spectrometry 82 . Since genotyping was performed by Sequenom, specific reaction details are not available. As it is described by Gabriel et al., 82 the assay consists of an initial locus-specific PCR reaction, followed by single base extension using mass-modified dideoxynucleotide terminators of an oligonucleotide primer which anneals immediately upstream of the polymorphic site of interest. Using MALDI-TOF mass spectrometry, the distinct mass of the extended primer identifies the SNP allele. Genotyping of 18 SNPs was attempted in Stage two; the rs186133190/2p15. rs6750047/2p22.2, rs498136/11q13.3 and rs4731742/7q32.3 assays failed in one or more Stage two datasets ( Supplementary Table 3 ). After QC (excluding individuals missing >1 genotype call, SNPs missing in > 3% of samples, SNPs with HWE P < 5 × 10 –4 ) there were 1,692 cases and 1,592 controls from Cambridge, 639 cases and 823 controls from Breakthrough Generations and 785 cases and 791 controls from Athens, Greece available for Stage two. We have estimated the power to reach P < 0.05 in the Stage two samples for all SNPs that reached genome-wide significance in the Stage one meta-analysis ( Supplementary Table 3 ). We converted ORs to genotype relative risks (as the SNPs are relatively frequent this is a reasonable assumption) and estimated power by simulating cases and controls (10,000 iterations) and conducting a Cochran-Armitage trend test (see Supplementary Table 5 ). Genome-wide Complex Trait Analysis (GCTA) 36 was used to perform conditional/joint GWAS analysis of newly identified or confirmed loci. GCTA allows conditional analysis of summary meta-analysis if provided with a sufficiently large reference population (2–5,000 samples) used in the meta-analysis to estimate LD. We used the QMEGA-610k set as a reference population to determine LD. QMEGA-610k imputation data for well imputed SNPs (INFO > 0.8) was converted to best guess genotypes using the GTOOL software (see URLs). Following best-guess conversion (genotype probability threshold 0.5), SNPs with a Hardy-Weinberg equilibrium P < 1 × 10 −6 , a MAF 3% missingness were removed. As per Yang et al., (2011) 36 we further cleaned the QMEGA-610K dataset to include only completely unrelated individuals (Identity by descent score ≤ 0.025 versus the standard 0.2 used in the meta-analysis), leaving a total of 4,437 people and 7.24 million autosomal SNPs in the reference panel. Within each new locus, Stage one fixed effects summary meta-analysis data for SNPs within 1 Mb either side of the top SNP were adjusted for the top SNP using the --cojo-cond option. As per Yang et al., (2011) we used the genomic control corrected GWAS-meta-analysis results. Where there remained an additional SNP with P < 5 × 10 –8 following adjustment for the top SNP we performed an additional round including both SNPs. If the remaining SNPs had P -values greater than 5 × 10 –8 no further analysis was performed. The results of this analysis are reported in Supplementary Table 6 . We have used the formula for calculating the proportion of familial relative risk (FRR) as outlined by the Cancer Oncological Gene-environment Study (see URLs). Given that CMM incidence is low, and the odds ratios reported small, we have assumed the odds ratios derived from the Stage one meta-analysis are equivalent to relative risks. With this assumption we have estimated the proportion of the FRR explained by each SNP (FRR snp ) as FRR snp = ( pr 2 + q )/ ( pr + q ) 2 Where risk allele and alternative allele frequency are p and q respectively, and r is the odds ratio for the risk allele Assuming a FRR melanoma for CMM of 2.19 83 and using the combined effect of all SNPs (assuming a multiplicative effect and a single SNP per loci), we computed the proportion of FRR is explained by a set of SNPs as log e (product of FRR snp )/ log e (FRR melanoma ). Pigmentation and nevus phenotype data were available for 980 melanoma cases and 499 control individuals from the Leeds case-control study 84 , 85 . Additional individuals from the Leeds melanoma cohort study 86 had pigmentation data available, giving a total of 1,458 subjects with melanoma and 499 control subjects. For the most significant SNP in each region reaching P < 1 × 10 –6 in the initial meta-analysis, log-transformed age- and sex-adjusted total nevus count was regressed on the number of risk alleles, adjusting for case-control status. A sun-sensitivity score was calculated for all subjects based on a factor analysis of six pigmentation variables (hair color, eye color, self-reported freckling as a child, propensity to burn, ability to tan and skin color on the inside upper arm) 19 . This score was similarly regressed on number of risk alleles and adjusted for case-control status. Full results can be found in Supplementary Table 1 . As the SNP most associated with the phenotype is quite likely not the underlying functional variant 87 , we performed a comprehensive bioinformatic assessment of SNPs with P fixed (if I 2 < 31%) or P random (if I 2 ≥ 31%) values within a factor of 100 of the P value for the peak SNP. To ensure we were not missing potentially interesting functional candidates, HaploReg was used to identify additional SNPs with r 2 >0.8 and no more than 200 kb away from those SNPs with P values within a factor of 100 of the peak SNP using 1000 Genomes Project pilot data. GCTA was used to confirm that SNPs carried forward for bioinformatic assessment derived from a common signal. Following adjustment for the locus’ top SNP, the SNPs selected for bioinformatic annotation at 6p22.3, 7p21.1, 10q24, 11q13.3 and 15q13.1 had CMM association P > 0.01. At 9q31.2 a single SNP rs1484384 retained a modest melanoma association ( P = 0.008) following adjustment for rs10739221; the rest were P > 0.01. At 2p22.2 those SNPs with P -values within 2 orders of magnitude of the peak SNP rs6750047 included rs1800440, a non-synonymous SNP with limited LD with rs6750047 (LD r 2 = 0.04). Following adjustment for rs6750047, the significance of rs1800440 was essentially unchanged ( P = 4.3 × 10 –4 ) and a second SNP rs163092 remained weakly associated with melanoma ( P = 0.008); all other SNPs had P > 0.01. Adjustment for both rs6750047 and rs1800440 removed the association between rs163092 and CMM ( P > 0.01). HaploReg 37 and RegulomeDB 40 were crosschecked to explore data reflecting transcription factor binding, open chromatin and the presence of putative enhancers. These tools summarize and collate data from the ENCODE 39 and Roadmap Epigenomics 38 public databases and from a range. ENCODE and Roadmap have assayed a large number of different cell types including keratinocyte and melanocyte primary cells, and, for a limited number of assays, melanoma cell lines; predicted functional activity in these cell types was given priority over cell types less likely to be involved in CMM risk. The summary results reported by HaploReg and RegulomeDB assign regions a putative function based on the combined results of multiple functional experiments and their position relative to known genes 38 , 39 . For example, ENCODE assigns the label of predicted enhancer to areas of open chromatin that overlaps a H3K4me1 signal, and binds transcription factors 39 . The Roadmap Epigenome uses as similar ranking system to ENCODE, and is summarized in the documentation for HaploReg 37 . For example, a weak enhancer will have only a weak H3K36me3 signal, while an active enhancer will have strong H3K36me3, H3K4me1 and H3K27ac signals. These labels are further divided into weak and strong depending on the quality of evidence. While these labels are predicted or putative, ENCODE reports that >50% of predicted enhancers are confirmed by follow up assays 39 , and these serve as a useful guide for interrogating CMM associated SNPs. Results from these tools were followed up in more detail using the UCSC genome browser 60 to explore the ENCODE 39 and the Roadmap epigenomics project 38 data. In addition, HaploReg uses genome-wide SNPs to estimate the background frequency of SNPs occurring in putative enhancer regions; this was used to test for enrichment in CMM-associated SNPs with an uncorrected binomial test threshold of P = 0.05 37 . The eQTL browser, the Genotype-Tissue Expression dataset (GTEx) 55 , and the Multiple Tissue Human Expression Resource (MuTHER 43 , 88 ) were further interrogated to attempt to resolve potential genes influenced by disease associated SNPs. For these databases we report only cis results; details of cell types and definition of cis boundaries can be found in Supplementary Tables 7 – 9 . The peak SNP for each locus, as well as other functionally interesting SNPs identified by HaploReg and RegulomeDB, were used to search listed eQTL databases. As the SNP coverage can differ for each database, where SNPs of interest were not present in the eQTL datasets we searched using high LD (>0.95) proxies. While priority was given to cell types more likely to be involved in CMM biology (e.g. sun-exposed skin from GTEx or skin from MuTHER) we reported eQTLs from other tissue types to highlight potential functional impact for identified SNPs. Regional plots of –log 10 P- values were generated using LocusZoom 89 . Where pairwise linkage disequilibrium measures are given, these were estimated from the 379 European ancestry 1000 Genomes Phase 1 April 2012 samples using PLINK 90 or the --hap-r2 command in VCFtools 91 unless otherwise indicated. To test for any overlap with published GWAS association, results reported in the NHGRI catalog for reported loci were extracted on 24/07/2014 and cross checked against the Stage one meta-analysis results. Manhattan plots were generated in R based on scripts written by Stephen Turner (see URLs). Forest plots were generated using the R rmeta package 92 .

We conducted a two-stage genome-wide meta-analysis. Stage one consisted of 11 GWAS totaling 12,874 cases and 23,203 controls from Europe, Australia and the USA; this includes all six published CMM GWAS and five unpublished ones ( Supplementary Table 2 ). In Stage two we genotyped 3,116 CMM cases and 3,206 controls from three additional datasets (consisting of 1,692 cases and 1,592 controls from Cambridge, UK, 639 cases and 823 controls from Breakthrough Generations, UK, and 785 cases and 791 controls from Athens, Greece; Online Methods) for the most significant SNP from each region reaching P < 10 –6 in Stage one and included these results in an Overall meta-analysis of both stages, totaling 15,990 melanoma cases and 26,409 controls. Details of these studies can be found in Supplementary Note . Given that the previous single-largest melanoma GWAS was of 2,804 cases and 7,618 controls 11 , this meta-analysis represents a fourfold increase in sample size compared to previous efforts to identify the genetic determinants of melanoma risk. Unless otherwise indicated we report the P -values from the Overall meta-analysis combining the two stages ( Supplementary Table 3 ). Forest plots of the individual GWAS study results can be found in Supplementary Figure 1 . All Stage one studies underwent similar quality control (QC) procedures, were imputed using the same reference panel and the results analyzed in the same way, with the exception of the Harvard and MD Anderson Cancer Center (MDACC) studies (see Online Methods). A fixed effects ( P fixed ) or random effects ( P random ) meta-analysis was conducted as appropriate depending on between-study heterogeneity. 9,470,333 imputed variants passed QC in at least two studies, of which 3,253 reached P fixed < 1 × 10 –6 and 2,543 reached P fixed 31%, value < 1 × 10 –7 ( Supplementary Table 4 ). The Stage one meta-analysis genome-wide inflation value (λ) was 1.032, and as λ increases with sample size we also adjusted the λ to a population of 1,000 cases and 1,000 controls 30 . The resulting λ 1000 of 1.002 suggested minimal inflation. Quantile-quantile (QQ) plots for the Stage one meta-analysis and individual GWAS studies can be found in Supplementary Figures 2 and 3 . To further confirm that our results were not influenced by inflation, the Stage one meta-analysis was repeated correcting for individual studies’ λ; P -values were essentially unchanged (Online Methods, Supplementary Table 3 ). All 13 previously-reported genome-wide significant loci (most first identified in one of the studies included here) reached P < 5 × 10 –8 in Stage one ( Figure 1 , Supplementary Table 4 ). In addition to confirming the two previously-reported sub-genome-wide significant loci at 11q13.3 (rs498136, 89 kb from CCND1 ) and 15q13.1 (rs4778138 in OCA2 ) we found three novel loci reaching genome-wide significance at 6p22.3, 7p21.1, and 9q31.2 ( Table 1 ; Figure 2 ). SNPs in another 11 regions reached P < 10 –6 ( Supplementary Table 3 ); notably, three were close to known telomere-related genes (rs2995264 is in OBFC1 31 in 10q24.33, rs11779437 is 1.1 Mb from TERF1 32 in 8q13.3, and rs4731207 is 66 kb from POT1 in 7q31.33, in which loss-of-function variants occur in some melanoma families 33 , 34 ). Given the importance of telomeres in melanoma we additionally genotyped two SNPs that did not quite reach our P < 10 –6 threshold but are close to telomere-related genes 35 : rs12696304 in 3q26.2 ( P fixed = 1.6 × 10 –5 ) is 1.1 kb from TERC and rs75691080 in 20q13.33 ( P fixed = 1.0 × 10 –6 ) is 19.4 kb from RTEL1 . In total 18 SNPs were carried through to Stage two (Online methods). Including the Stage two results in the Overall meta-analysis led to two new genome-wide significant regions, 2p22.2 and 10q24.33 ( Figure 2 ; Table 1 , Supplementary Table 3 ). The Stage two data also serve the purpose of independently confirming with genotype data the meta-analysis results from imputed SNPs. Five SNPs, rs4778138 ( OCA2 /15q13.1 ), rs498136 ( CCND1 /11q13.3 ), and the novel rs10739221 (9q31.2), rs6750047 (2p22.2) and rs2995264 (10q24.33) all reached P < 0.05 in the genotyped Stage two samples. We have estimated the power to reach P < 0.05 in the Stage two samples for all SNPs that reached genome-wide significance in the Stage one meta-analysis (Online Methods, Supplementary Table 5 ). rs6914598 (6p22.3) was only genotyped in the Athens sample and thus had a power of only 0.35. Of the remaining four SNPs that were genome-wide significant in Stage one, while the 7p21.1 SNP rs1636744 was well powered ( > 90%), the probability that all four of these well-powered SNPs would reach P < 0.05 in the analysis of Stage two data was only (0.916 × 0.736 × 0.787 × 0.955) = 0.51, so it is not surprising that one was not significant. The SNPs in 7p21.1 (rs1636744) and 6p22.3 (rs6914598) did not reach nominal significance in Stage two, but for both SNPs the confidence intervals for the effect estimates overlapped those from the Stage one meta-analysis. In terms of heritability the 13 loci that were genome-wide significant before this meta-analysis explained 16.9% of the familial relative risk (FRR) for CMM, with MC1R explaining 5.3% alone (Online Methods). Including the seven loci confirmed or reported here (2p22.2, 6p22.3, 7p21.1, 9q31.2, 10q24.33, 11q13.3, 15q13.1), an additional 2.3% of FRR is explained. In total, all 20 loci explain 19.2% of the FRR for CMM; this is a conservative estimate given the assumption of a single SNP per locus. We tested all new and known CMM risk loci for association with nevus count or pigmentation ( Supplementary Table 1 ). Aside from the known association between OCA2 and pigmentation, none of the newly-identified loci were associated ( P > 0.05). Following confirmation of the loci in the Stage two analysis, we performed conditional analysis on the Stage one meta-analysis results to determine whether there were additional association signals within 1 Mb either side of the top SNP at each locus using the Genome-wide Complex Trait Analysis (GCTA) software 36 (Online Methods; Supplementary Table 6 ). This indicated that while there are additional SNPs associated with CMM at each locus, for all but chromosome 7 and 11 the additional signals were not strongly associated with melanoma ( P < 1 × 10 –7 ; for more detail see Supplementary Note and Supplementary Figures 4 and 5 ). We then conducted a comprehensive bioinformatic assessment of the top SNP from each of the seven new genome-wide significant loci using a range of annotation tools, databases of functional and eQTL results and previously-published GWAS results (see Online Methods, Supplementary Tables 7 – 9 ). We applied the same analyses to each locus but, to limit repetition, where nothing was found for a given resource (e.g. NHGRI GWAS catalog) we do not explicitly report this. While rs6750047 in 2p22.2 was not genome-wide significant in the Stage one meta-analysis it reached genome-wide significance ( P fixed = 7.0 × 10 –9 , OR = 1.10, I 2 = 0.00; Table 1 , Supplementary Table 3 ) in the Stage two and Overall meta-analysis. The association signals for 2p22.2 ( Figure 2 ) span the 3' UTR of RMDN2 (also known as FAM82A1 ) and the entirety of the CYP1B1 gene, and as such there is a wealth of bioinformatic annotation for SNPs associated with CMM risk. Considering the 26 SNPs with P -values within two orders of magnitude of rs6750047 in 2p22.2 ( Supplementary Tables 7 – 9 ), HaploReg 37 reports a significant enrichment of strong enhancers in epidermal keratinocytes (4 observed, 0.6 expected, P = 0.003). The paired rs162329 and rs162330 (LD r 2 =1.0, 98 bp apart; P fixed = 3.91 ×10 –6 , I 2 = 11.23) lie approximately 10 kb upstream from the CYP1B1 transcription start site in a potential enhancer in keratinocytes and other cell types 37 - 40 . These two SNPs are eQTLs for CYP1B1 in three independent liver sample sets 41 , 42 . In addition several SNPs, including the peak SNP for 2p22.2, rs6750047, are strong CYP1B1 eQTLs in LCLs in the Multiple Tissue Human Expression Resource 43 (MuTHER; P < 5 × 10 –5 ). It is worth noting the overlap between the liver and lymphoblastoid cell line (LCL) eQTLs is incomplete; rs162330 and rs162331 are only weak eQTLs in MuTHER data ( P ~ 0.01). In terms of functional annotation the most promising SNP near rs6750047 is rs1374191 ( P fixed = 5.4 × 10 –5 , OR = 1.07, I 2 = 0.00); in addition to being a CYP1B1 eQTL in LCLs (MuTHER P = 6.9 × 10 –8 ); this SNP is positioned in a strong enhancer region in multiple cell types including melanocytes and keratinocytes 37 - 40 . In summary, SNPs in 2p22.2 associated with melanoma lie in putative melanocyte and keratinocyte enhancers and are also cross-tissue eQTLs for CYP1B1 . CYP1B1 metabolizes endogenous hormones, playing a role in hormone associated cancers including breast and prostate (reviewed in 44 ). CYP1B1 also metabolizes exogenous chemicals, resulting in pro-cancer (e.g. polycyclic aromatic hydrocarbons) and anti-cancer (e.g. tamoxifen) outcomes 44 . The former is of interest as CYP1B1 is regulated by ARNT , a gene at the melanoma-associated 1q21 locus 12 . The CYP1B1 promotor is methylated in melanoma cell lines and tumor samples 45 . CYP1B1 missense protein variants have been associated with cancers including squamous cell carcinoma and hormone associated cancers 44 , 46 . Of these only rs1800440 (p.Asn453Ser) is moderately associated with melanoma ( P fixed = 1.83 × 10 –5 , OR = 0.90, I 2 = 0.00), and it was included in the bioinformatic annotation ( Supplementary Tables 7 – 9 ). rs1800440 is not in LD with the CMM risk meta-analysis peak SNP rs6750047/2p22.2 (LD r 2 = 0.04) and adjusting for rs6750047 only slightly reduces its association with CMM ( P = 4.3 × 10 –4 , Online Methods). Truncating mutations in CYP1B1 are implicated in primary congenital glaucoma 47 and since glaucoma cases are used as controls in the contributing Western Australian Melanoma Health Study (WAMHS) melanoma GWAS, we considered the impact of excluding glaucoma cases; the SNP remains genome-wide significantly associated with CMM even after such exclusions ( Supplementary Note, Supplementary Table 10 ). While the association with melanoma in the WAMHS set is stronger without glaucoma cases (beta 0.05 vs. 0.19) both betas are within the range observed for other melanoma datasets and no heterogeneity (I 2 = 0.00) is observed with or without the glaucoma samples. rs6914598 ( P fixed = 3.5 × 10 –8 , OR = 1.11, I 2 = 0.00) lies in 6p22.3, in intron 12 of CDKAL1, a gene that modulates the expression of a range of genes including proinsulin via tRNA methylthiolation 48 , 49 . Bioinformatic assessment of the 35 SNPs with P -values within two orders of magnitude of the 6p22.3 peak rs6914598 by HaploReg 37 indicates that the most functionally interesting SNP is rs7776158 (Stage one P fixed = 3.8 × 10 –8 , I 2 = 0, in complete LD with rs6914598, r 2 = 1.0), which lies in a predicted melanocyte enhancer that binds IRF4 38 , 39 . IRF4 binding is of interest given the existence of a functional SNP rs12203592 in the IRF4 gene 50 , associated with nevus count, skin pigmentation and tanning response 51 - 54 . rs1636744 ( P fixed = 7.1 × 10 –9 , OR = 1.10, I 2 = 0.00; Figure 2 ) is in an intergenic region of 7p21.1 and lies 63 kb from AGR3 . rs1636744 is an eQTL for AGR3 in lung tissue (Genotype-Tissue Expression (GTEx) project P = 1.6 × 10 –6 ) 55 , 56 . AGR3 is a member of the protein disulphide isomerase family, which generate and modify disulphide bonds during protein folding 57 . AGR3 expression has been associated with breast cancer risk 58 and poor survival in ovarian cancer 59 . GTEx confirms that AGR3 is expressed in human skin samples. The region containing rs1636744 is not conserved in primates (UCSC genome browser 60 ), and RegulomeDB 40 indicates there is little functional activity at this SNP. More promising are rs847377 and rs847404 which, in addition to being both AGR3 eQTLs in lung tissue 55 and associated with CMM risk (Stage one P fixed = 3.89 × 10 –8 and 1.72 × 10 –7 ), are in putative weak enhancers in a range of cells including melanocytes and keratinocytes 37 - 40 . Adjusting for rs1636744 renders rs847377 and rs847404 non-significant ( P > 0.6) indicating that they are tagging a common signal. rs1636744, rs847377 and rs847404 are not eQTLs for AGR3 in sun-exposed skin. The melanoma-associated variants at 9q31.2, peaking at rs10739221 (Overall P fixed = 7.1 × 10 –11 ; I 2 = 0.00; Figure 2 ) are intergenic. The nearest genes are TMEM38B, ZNF462 and the nucleotide excision repair gene RAD23B 61 . While bioinformatic annotation did not reveal any putative functional SNPs, based on the importance of DNA repair in melanoma, RAD23B is of particular interest. rs10739221 is 635 Kb from the leukemia-associated TAL2 62 , and 1.2 Mb from KLF4 , which regulates both telomerase activity 63 and the melanoma-associated TERT 64 . While not genome-wide significant in Stage one, rs2995264 in 10q24.33 is strongly associated with telomere length 28 , 35 and was genotyped in Stage two. rs2995264 was significantly associated with CMM in the Cambridge study ( P = 0.046) and strong in the Breakthrough dataset ( P = 8.0 × 10 –4 ); in the Overall meta-analysis this SNP reached genome-wide significance ( P fixed = 2.2 × 10 –9 ; I 2 = 27.14). The melanoma association signal at 10q24.33 ( Figure 2 ) spans the OBFC1 gene and the promotor of SH3PXD2A. Given the strong telomere length association at this locus the most promising candidate is OBFC1, a component of the telomere maintenance complex 31 . HaploReg reports that SNPs within two orders of magnitude of rs2995264 in 10q24.33 are significantly more likely to fall in putative enhancers in keratinocytes than would be expected by chance ( P < 0.001). Promising candidate functional SNPs include the conserved rs11594668 and rs11191827 which lie in putative melanocyte and keratinocyte enhancers, and bind transcription factors 37 - 40 . The association observed at rs2995264/10q24.33 is independent of a recent report of a melanoma association at 10q25.1 65 . Our peak SNP for 10q24.33, rs2995264, and the 10q25.1 SNPs rs17119434, rs17119461, and rs17119490 reported in Teerlink et al., (2012) 65 are in linkage equilibrium (LD r 2 0.2). The CMM-associated variants at 11q13.3 peak at rs498136 (Overall P fixed = 1.5 × 10 –12 , OR = 1.13, I 2 = 0.00; Supplementary Figure 6 ) 5’ to the promotor of CCND1 . In the initial report of CCND1 11 rs11263498 was borderline in its association with melanoma ( P = 3.2 × 10 –7 ) and while supported ( P = 0.017) by the two replication studies exhibited significant heterogeneity and did not reach genome-wide significance (overall P random = 4.6 × 10 –4 , I 2 = 45.00). The previously-reported rs11263498 and the meta-analysis peak of rs498136/11q13.3 are in strong linkage disequilibrium (LD) (r 2 = 0.95). Bioinformatic assessment of the CCND1 region indicated the peak SNP rs498136/11q13.3 is in a putative enhancer in keratinocytes in both ENCODE and Roadmap data 37 - 40 . Considering other SNPs strongly associated with CMM, both the previously-reported 11 rs11263498 (Stage one P fixed = 1.8 × 10 –9 , OR = 1.12, I 2 = 0.00) and rs868089 (Stage one P fixed = 2.0 × 10 –9 , OR = 1.12, I 2 = 0.00) lie in putative melanocyte enhancers. Somatic CCND1 amplification in CMM tumors positively correlates with markers of reduced overall survival, including Breslow thickness and ulceration 66 , 67 . The CCND1 association with breast cancer has been extensively fine-mapped, revealing three independent association signals 68 . rs554219 and rs75915166 tag the two strongest functional associations with breast cancer 68 but are not themselves associated with CMM risk (Stage one P fixed > 0.1, I 2 = 0.00). While the third signal in breast cancer was not functionally characterized 68 , its tag SNP rs494406 is modestly associated with CMM (Stage one P fixed > 0.0002, I 2 = 0.00, LD r 2 =0.47 with rs498136/11q13.3). rs494406 is no longer significant after adjustment by rs498136 ( P = 0.53; Supplementary Table 6 ), suggesting that SNPs that are in LD in this region are associated with risk of both melanoma and breast cancer. Both OCA2 and nearby HERC2 at the 15q13.1 locus have long been associated with pigmentation traits 51 . rs12913832 in HERC2 , also known as rs11855019, is the major determinant of eye color in Europeans 69 , making this region a strong candidate for CMM risk. One of the studies contributing to this meta-analysis previously reported a genome-wide significant association between melanoma and rs1129038 and rs12913832 in HERC2 (in strong LD 10 reported as r 2 = 0.985), but this was not supported ( P > 0.05) by any of the three replication GWAS (final P = 2.5 × 10 –4 ) 10 . Stratification might be an issue for this locus as eye color frequencies vary markedly across European populations. Indeed, in our meta-analysis, which includes all four of these GWAS, both rs1129038 and rs12913832 showed highly heterogeneous effects in the CMM risk meta-analysis ( P random = 0.037 and 0.075 respectively, I 2 > 77.00). Amos et al., (2011) found that rs4778138 in OCA2 , which is only in weak LD with rs12913832 (r 2 = 0.12), exhibited a more consistent association across studies, albeit not genome-wide significantly. In our Overall meta-analysis we confirm rs4778138 in 15q13.1 is associated with CMM risk ( P fixed = 2.2 × 10 –11 , OR = 0.84, I 2 = 0.00; Figure 2 ). Following adjustment of the 15q13.1 signal by rs4778138 the effect size for the eye color SNP rs12913832 is reduced from beta = 0.12 to beta = 0.064. Conversely adjustment for rs12913832 reduces rs4778138's association with CMM (beta reduced from −0.178 to −0.114, corrected P = 1.6 × 10 –4 ). rs12913832 is poorly imputed across studies, reaching INFO > 0.8 in only 6 studies, and we are unable to conclusively exclude a role for rs12913832 at this locus. HaploReg indicates rs4778138 is within a putative melanocyte enhancer in Roadmap epigenetic data 37 - 40 . While it is not clear which gene(s) in 15q13.1 is/are influenced by melanoma-associated SNPs, the fact that rs4778138 is associated with eye colors intermediate to blue and brown 70 supports a role for OCA2 .

A further nine loci were associated with CMM risk at multiple SNPs with P < 10 –6 in Stage one but did not reach P < 5 ×10 –8 in the Overall meta-analysis ( Supplementary Table 3 ). Given that genome-wide significance is based on a Bonferroni correction assuming 1,000,000 independent tests, we would expect only one locus to reach P < 10 –6 and the probability that as many as nine loci reach this threshold is 1.1 × 10 –6 (exact binomial probability), so it is highly likely that several of these are genuine. Of the 16 regions that reached P < 10 –6 , three were near genes involved in telomere biology 7q31.33 (rs4731207 near POT1 ), 8q13.3 (rs11779437 near TERF1), and 10q24.33 (rs2995264 near OBFC1 ) ( Supplementary Table 3 ). Given the evidence for telomere biology in melanoma 20 , 25 - 28 and that previous genome-wide significant SNPs are near the telomere maintenance genes TERT , PARP1 and A TM, we included two further biological candidates: rs12696304, 1.1 kb from TERC in 3q26.2, and rs75691080 in 20q13.33 near RTEL1 . Of these five SNPs, rs2995264 (10q24.33/ OBFC1) attained genome-wide significance in the overall analysis while rs12696304 (3q26.2/ TERC ) was significant in Stage two ( P = 4.0 × 10 –3 ), and reached P = 2.8 × 10 –7 in the Overall meta-analysis ( Supplementary Table 3 ). While falling short of genome-wide significance this is nonetheless suggestive of an association at this locus. Neither rs4731207 (66 kb from POT1 in 7q31.33) nor rs75691080 (19.4 kb from RTEL1 in 20q13.33) were significantly associated with melanoma risk in Stage two, but in neither case did the estimated effect differ significantly from Stage one. In addition rs75691080 ( RTEL1 /20q13.33) is marginally associated with nevus count ( P = 0.058; Supplementary Table 1 ). Of the SNPs near telomere-related genes, rs11779437 in 8q13.3 was the most distant (1.1 Mb from TERF1 ) and was the only one to show a significantly different effect in Stage two (Overall P random = 0.013, OR = 0.93, I 2 = 42.06). This is most likely due to the initial signal being a false positive, but may be due to a lack of power.

This two-stage meta-analysis, representing a fourfold increase in sample size compared to the previous largest single melanoma GWAS, has confirmed all thirteen previously reported loci, as well as resolving two likely associations at CCND1 and HERC2 / OCA2 . The CCND1 association with melanoma only partially overlaps the signal observed for breast cancer 68 . The HERC2 / OCA2 association is with rs4778138/15q13.1, which may be a subtle modifier of eye color 70 , but we cannot rule out that the association at this locus is influenced by the canonical blue/brown eye color variant rs12913832. Our Stage one meta-analysis of over 12,000 melanoma cases identified three novel risk regions, with only rs10739221 formally significant ( P < 0.05) in Stage two ( Table 1 ). Two further loci (2p22.2 and 10q24.33) reached genome-wide significance with the addition of the Stage two data ( Figure 2 ; Table 1 , Supplementary Table 3 ). In total our Overall meta-analysis identified 20 genome-wide significant loci; 13 previously replicated, two previously-reported but first confirmed here and five that are novel to this report. The new loci identified in this meta-analysis explain an additional 2.3% of the familial relative risk for CMM. Overall, 19.2% of the FRR is explained by all 20 genome-wide significant loci combined. Except for the association at 9q31.2, the reported loci contain SNPs that are both strongly associated with melanoma and fall within putative regulatory elements in keratinocyte or melanocyte cells, with the nearby nucleotide excision repair gene RAD23B at 9q31.2 a promising candidate. eQTL datasets suggest that melanoma-associated SNPs at 7p21 regulate the expression of AGR3 albeit in lung tissue and not sun-exposed skin. AGR3 expression has been implicated in breast and ovarian cancer outcome. SNPs in 2p22.3 are associated with the expression of CYP1B1 . Although this gene is better known for its role in hormone-associated cancers it may influence melanoma risk through metabolism of exogenous compounds, a process regulated by ARNT at the 1q21 melanoma-associated locus. We have used the power of this large collection of CMM cases and controls to identify five novel loci, none of which are significantly associated with classical CMM risk factors and thus highlight novel disease pathways. Interestingly, we now have genome-wide significant evidence for association between CMM risk and a SNP in the telomere-related gene OBFC1 in 10q24.33, in addition to the established telomere-related associations at the TERT / CLPTM1L, PARP1 and ATM loci. We also have support, albeit not genome-wide significant, for TERC , the most significant predictor of leukocyte telomere length in a recent study 35 . Of the 20 loci that now reach genome-wide significance for CMM risk, five are in regions known to be related to pigmentation, three in nevus-related regions and four in regions related to telomere maintenance. This gives further evidence that the telomere pathway, with its effect on the growth and senescence of cells, may be important in understanding the development of melanoma.