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

Thirteen common susceptibility loci have been reproducibly associated with cutaneous malignant melanoma (CMM). We report the results of an international 2-stage meta-analysis of CMM genome-wide association studies (GWAS). This meta-analysis combines 11 GWAS (5 previously unpublished) and a further three stage 2 data sets, totaling 15,990 CMM cases and 26,409 controls. Five loci not previously associated with CMM risk reached genome-wide significance (P < 5 × 10(-8)), as did 2 previously reported but unreplicated loci and all 13 established loci. Newly associated SNPs fall within putative melanocyte regulatory elements, and bioinformatic and expression quantitative trait locus (eQTL) data highlight candidate genes in the associated regions, including one involved in telomere biology.

Cutaneous malignant melanoma (CMM) primarily occurs in fair-skinned individuals; the major host risk factors for CMM include pigmentation phenotypes[1-4], the number of melanocytic nevi[5,6] and a family history of melanoma[7].

Six population-based genome-wide association studies (GWAS) of CMM have been published[8-13] identifying 12 regions that reach genome-wide significance. Some of these regions were already established melanoma risk loci, for example through candidate gene studies[14] (for review see[15]). A 13th region in 1q42.12, tagged by rs3219090 in PARP1, that was borderline in the initial publication (P = 9.3 × 10−8)[12] was confirmed as genome-wide significant by a recent study (P = 1.03 × 10−8)[16]. As might be expected for common variants influencing CMM risk many of these loci contain genes that are implicated in one of the two well-established heritable risk phenotypes for melanoma, pigmentation (SLC45A2, TYR, MC1R and ASIP) and nevus count (CDKN2A/MTAP, PLA2G6 and TERT) (Supplementary Table 1)[17]. The presence of DNA repair genes such as PARP1 and ATM at two loci suggests a role for DNA maintenance pathways, leaving four loci where the functional mechanism is less clear (ARNT/SETDB1, CASP8, FTO and MX2).

Of particular interest is TERT, which is involved in telomere maintenance; SNPs in this region have been associated with a variety of cancers[11,18-22]. Further, ATM and PARP1’s DNA repair functions extend to telomere maintenance and response to telomere damage[23,24]. Longer telomeres have been associated with higher nevus counts and it has been proposed that longer telomeres delay the onset of cell senescence, allowing further time for mutations leading to malignancy to occur[20,25]. There is evidence that longer telomeres increase melanoma risk[20,26,27] and that other telomere-related genes are likely involved in the etiology of melanoma, but none of these loci has yet reached genome-wide significance (or even P < 10−6)[28].

In addition, two independent SNPs at 11q13.3, near CCND1, and 15q13.1, adjacent to the pigmentation gene OCA2, have been associated previously with melanoma, but did not meet the strict requirements for genome-wide significance, either not reaching P = 5 × 10−8 in the initial report, or not replicating in additional studies[10,11,29]. This meta-analysis has resolved the status of these two loci, as well as identified novel melanoma susceptibility loci.

Results and Discussion

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 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).

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 MDACC studies (see Online Methods). A fixed effects (P) or random effects (P) 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 < 1 × 10−6 and 2,543 reached P < 5 × 10−8. For reference we provide a list of SNPs that reached a P, or P if I2 > 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). Forest plots of the individual GWAS study results can be found in Supplementary Figure 1. 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 (Pfixed = 1.6 × 10−5) is 1.1 kb from TERC and rs75691080 in 20q13.33 (Pfixed = 1.0 × 10−6) is 19.4 kb from RTEL1. In total 18 SNPs were carried through to Stage two (Online methods).

Figure 1

Manhattan plot of the Stage one meta-analysis of GWAS of CMM from Europe, the USA and Australia

The P Stage one value for all SNPs present in at least two studies have been plotted using a log10(−log10) transform to truncate the strong signals at MC1R (P < 10−92) on chromosome 16 and CDKN2A (P < 10−31) on chromosome 9. The total Stage one meta-analysis included 11 CMM GWAS, totaling 12,874 cases and 23,203 controls. P < 5 × 10−8 (genome-wide significance) and P < 1 × 10−6 are indicated by a red and a blue line respectively. 18 loci reached genome-wide significance in Stage one. The 2 newly-confirmed loci 11q13.3 (CCND1) and 15q13.1 (HERC2/OCA2) are indicated by * and the 5 novel loci 2p22.2, 6p22.3, 7p21.1, 9q31.2 and 10q24.33 are highlighted by **. 2p22.2 (RMDN2/CYP1B1) and 10q24.33 (OBFC1) were genome-wide significant only in the Overall meta-analysis (Supplementary Table 3).

Table 1

Genome-wide significant results from a two-stage meta-analysis of GWAS of CMM from Europe, the USA and Australia

For each region we report the chromosomal location, nearest gene, and any other promising candidate gene in brackets for the top SNP. We also report the 1000 Genomes European population minor allele frequency (MAF) and minimum imputation quality across all studies (min INFO). The Stage one meta-analysis field reports the effect size estimate (beta) and P-value for the minor allele from the meta-analysis of 11 CMM GWAS, totaling 12,874 cases and 23,203 controls. Following their genotyping in three additional datasets (total 3,116 cases and 3,206 controls) we provide the Stage two meta-analysis results. Finally we provide the Overall meta-analysis of all available data. The results for the top SNP in each region that reached P < 1 × 10−6 in Stage one and so was genotyped in Stage two, per study results and evidence of heterogeneity of effect estimates across studies (I2) can be found in Supplementary Table 3. Where I2 values were below 31% fixed effects meta-analysis was used, otherwise random effects, and all genome-wide significant SNPs had low heterogeneity (I2 < 31%) in both Stage one and Overall. Regions previously confirmed as associated with melanoma (e.g. MC1R) are not shown.

				Stage onemeta-analysis	Stage twometa-analysis	Overallmeta-analysis
SNP	Region	Gene	MinorAllele:MAF(min INFO)	Beta (P)	Beta (P)	Beta (P)
rs6750047	2p22.2	RMDN2(CYP1B1)	A:0.43 (0.96)	0.088 (2.9 × 10−7)	0.113 (6.0 × 10−3)	0.092 (7.0 × 10−9)
rs6914598	6p22.3	CDKAL1	C:0.32 (0.88)	0.11 (2.6 × 10−8)	0.037 (0.63)	0.10 (3.5 × 10−8)*
rs1636744	7p21.1	AGR3	T:0.40 (0.96)	0.11 (1.8 × 10−9)	0.032 (0.38)	0.091 (7.1 × 10−9)
rs10739221	9q31.2	TMEM38B(RAD23B,TAL2)	T:0.24 (0.94)	0.12 (9.6 × 10−9)	0.145 (1.7 × 10−3)	0.12 (7.1 × 10−11)
rs2995264	10q24.33	OBFC1	G:0.088 (0.94)	0.14 (8.5 × 10−7)	0.206 (0.088)	0.16 (2.2 × 10−9)
rs498136	11q13.3	CCND1	A:0.32 (0.97)	0.12 (1.0 × 10−10)	0.124 (4.0 × 10−3)	0.12 (1.5 × 10−12)
rs4778138	15q13.1	OCA2	G:0.16 0.82	−0.18 (3.1 × 10−9)	−0.156 (1.7 × 10−3)	−0.17 (2.2 × 10−11)

Not genome-wide significant given after formal multiple testing correction e.g. P < 3.06 × 10−8 as in Li et al, (2012)[93].

Figure 2

Regional association plots for novel genome-wide significant loci 2p22.2, 6p22.3, 7p21.1, 9q31.2, 10q24.33 and the newly-confirmed region, 15q13.1 (OCA2)

The negative log10 of P values for SNPs from the Stage one meta-analysis of 12,874 cases and 23,203 controls have been plotted against their genomic position (Mb) using LocusZoom[89]. The rs ID is listed for the peak SNP in each region (purple diamond). The P-values and effect sizes for listed SNPs can be found in Supplementary Table 3. For the remaining SNPs the color indicates linkage disequilibrium r2 with the peak SNP. Note FAM82A1’s alternative gene ID is RMND2. Neither rs2995264 in 10q24.33 nor rs6750047 in 2p22.2 are genome-wide significant in Stage one, but reach this in the Overall meta-analysis. The plot for 11q13.3 (CCND1) can be found in Supplementary Figure 4.

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). 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.

2p22.2

While rs6750047 in 2p22.2 was not genome-wide significant in the Stage one meta-analysis it reached genome-wide significance (P = 7.0 × 10−9, OR = 1.10, I2 = 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 r2 =1.0, 98 bp apart; P = 3.91 ×10−6, I2 = 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 = 5.4 × 10−5, OR = 1.07, I2 = 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 = 1.83 × 10−5, OR = 0.90, I2 = 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 r2 = 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 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 (I2 = 0.00) is observed with or without the glaucoma samples.

6p22.3

rs6914598 (P = 3.5 × 10−8, OR = 1.11, I2 = 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 = 3.8 × 10−8, I2 = 0, in complete LD with rs6914598, r2 = 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].

7p21.1

rs1636744 (P = 7.1 × 10−9, OR = 1.10, I2 = 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 (GTEx 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. Evidence that the regions containing rs1636744 are 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 = 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.

9q31.2

The melanoma-associated variants at 9q31.2, peaking at rs10739221 (Overall P = 7.1 × 10−11; I2 = 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].

10q24.33

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 = 2.2 × 10−9; I2 = 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. 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 r2 <0.01) and in turn these SNPs are not associated with CMM in our meta-analysis (P > 0.2).

11q13.3

The CMM-associated variants at 11q13.3 peak at rs498136 (Overall P = 1.5 × 10−12, OR = 1.13, I2 = 0.00; Supplementary Figure 4) 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 = 4.6 × 10−4, I2 = 45.00). The previously-reported rs11263498 and the meta-analysis peak of rs498136/11q13.3 are in strong linkage disequilibrium (LD) (r2 = 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 = 1.8 × 10−9, OR = 1.12, I2 = 0.00) and rs868089 (Stage one P = 2.0 × 10−9 , OR = 1.12, I2 = 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 > 0.1, I = 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 > 0.0002, I = 0.00, LD r2=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.

15q13.1

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 r2 = 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 = 0.037 and 0.075 respectively, I2 > 77.00).

Amos et al., (2011) found that rs4778138 in OCA2, which is only in weak LD with rs12913832 (r2 = 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 = 2.2 × 10−11, OR = 0.84, I2 = 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.

Evidence of additional melanoma susceptibility loci

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 ATM, 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 = 0.013, OR = 0.93, I2 = 42.06). This is most likely due to the initial signal being a false positive, but may be due to a lack of power.

Conclusion

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.

Online Methods

Stage one array genotyping

The samples were genotyped on a variety of commercial arrays, detailed in the Supplementary Note.

Stage one genome-wide imputation

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 < 10−4 or missingness > 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 R2 > 0.95 were included in the final analysis.

Stage one genome-wide association 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 R2) 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].

Meta-analysis

Heterogeneity of per-SNP effect sizes in studies contributing to the_Stage one, Stage two and the Overall meta-analyses was assessed using the I2 metric[79]. I2 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 I2 was less than 31% a fixed effects model was used, with fixed effects P-values indicated by Pfixed; otherwise random effects were applied (P). The method of Dersimonian and Laird[80] was used to estimate the between-studies variance, . An overall random effects estimate was then calculated using the weights where ν is the variance of the estimated effect. 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. GIF-corrected P-values are provided in Supplementary Table 3.

Stage two genotyping

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 4) and is in linkage equilibrium with the genome-wide significant hit in this region (r2 = 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.

Statistical power 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).

Conditional Analysis

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 MAF <0.01 and > 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.

Proportion of Familial Relative Risk

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 (FRRsnp) as FRRsnp = (pr2 + 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 FRRmelanoma 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(product of FRRsnp)/ log(FRRmelanoma).

Association with nevus count or pigmentation

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] included pigmentation data 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, logged 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.

Bioinformatic annotation

As the SNP most associated with the phenotype is quite likely not the underlying functional variant[87] at each locus we considered SNPs with P if I2 < 31%, or P if I2 >= 31%, within a factor of 100 of the peak SNP for comprehensive bioinformatic assessment. To ensure we were not missing potentially interesting functional candidates, HaploReg was used to identify additional SNPs within 200kb and with LD r2 >0.8 using 1000 Genomes pilot data[37,77]. 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 r2 = 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 public databases ENCODE[39], the Roadmap epigenomics project[38] as well as a range of other functional tools. 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, H3K3me1 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 −log10P-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.

Additional methods

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].