Multiplex melanoma families are enriched for polygenic risk.

Cancers, including cutaneous melanoma, can cluster in families. In addition to environmental etiological factors such as ultraviolet radiation, cutaneous melanoma has a strong genetic component. Genetic risks for cutaneous melanoma range from rare, high-penetrance mutations to common, low-penetrance variants. Known high-penetrance mutations account for only about half of all densely affected cutaneous melanoma families, and the causes of familial clustering in the remainder are unknown. We hypothesize that some clustering is due to the cumulative effect of a large number of variants of individually small effect. Common, low-penetrance genetic risk variants can be combined into polygenic risk scores. We used a polygenic risk score for cutaneous melanoma to compare families without known high-penetrance mutations with unrelated melanoma cases and melanoma-free controls. Family members had significantly higher mean polygenic load for cutaneous melanoma than unrelated cases or melanoma-free healthy controls (Bonferroni-corrected t-test P = 1.5 × 10-5 and 6.3 × 10-45, respectively). Whole genome sequencing of germline DNA from 51 members of 21 families with low polygenic risk for melanoma identified a CDKN2A p.G101W mutation in a single family but no other candidate high-penetrance melanoma susceptibility genes. This work provides further evidence that melanoma, like many other common complex disorders, can arise from the joint action of multiple predisposing factors, including rare high-penetrance mutations, as well as via a combination of large numbers of alleles of small effect.

Introduction

Cutaneous melanoma (hereafter melanoma) is an aggressive cancer that originates in melanocytes, the pigment-producing cells of the skin. Melanoma can cluster in families; one formal definition of this is when the number of familial cases is greater than expected based on total members, their ages and the melanoma incidence in their birth cohort (1). Using this definition, Aitken et al. found that 4.7% of 1116 families had a case count higher than expected given the number, age and cohort of family members. Other, related, approaches have estimated that up to 9% of melanoma cases have a high familial risk. Within melanoma-dense families, high-penetrance segregating mutations have been found in CDKN2A (2,3), CDK4 (4,5), BAP1 (6), MITF (7,8), TERT (9,10), POT1 (11,12), ACD, and TERF2IP (13). However, these rare high-penetrance mutations underpin melanoma development in only approximately 50% of densely affected families (14).

While melanoma rates globally increase with ultraviolet radiation exposure, indicating a substantial environmental component to risk, even within a uniformly high-risk environment, there will still be genetic variation in risk. Using twin and family studies, melanoma’s heritable component has been estimated to be approximately 55% (15,16). In addition to rare high-penetrance mutations, genome-wide association studies (GWASs) have identified many loci harboring common genetic variants with a low to moderate impact on melanoma risk (17,18). The polygenic disease burden arising from these moderate and common GWAS variants can be estimated from genome-wide single nucleotide polymorphism (SNP) heritability (19). The SNP heritability for melanoma has been estimated to be 30% (20). The difference between the SNP and twin heritability estimates is likely due to poorly tagged SNPs and is likely to diminish in time as SNP genotyping arrays and coverage (e.g. through whole genome-sequencing) improves. Although not all of the genes underpinning SNP heritability have been identified, Lu et al. estimated that approximately a third of the SNP heritability can be explained by the specific genes identified in GWASs. These variants can be aggregated in a weighted sum for predicting risk—a polygenic risk score (PRS) (21,22).

When compared to non-high-risk melanoma families, high-risk families from Aitken et al.’s study had significantly higher proportions of melanoma risk phenotypes that are themselves polygenic traits (higher nevus count, poor tanning and fair skin) suggesting that these families may be enriched for melanoma polygenic risk (1,23,24). Thus, in addition to the carriage of high-penetrance mutations, and shared environment/behaviors, familial aggregation of melanoma could arise from high polygenic risk within the family. For example, if the familial clustering is due solely to an unidentified high-penetrance variant, the PRS would be expected to be lower than in sporadic cases and more similar to the melanoma-free population. In contrast, if the clustering is due to the low penetrance variants in the PRS, the PRS would be expected to be higher than in the sporadic cases. Finally, it is possible that having a high polygenic load may increase the penetrance from a melanoma risk mutation. However, the data may not always be clear cut. For example, application of a PRS for bipolar disorder to a family densely enriched for this disease found that affected family members were closer to the mean of unrelated GWAS cases than to healthy controls (25). That is, family members had a higher polygenic risk for bipolar disorders than healthy controls, suggesting that it may, in part, explain their familial aggregation. A recent study reported significantly lower melanoma PRS scores in families carrying high-penetrance mutations compared to those families without such a mutation (26). Here we report the distribution of polygenic risk for cutaneous melanoma within densely affected melanoma families compared to unrelated melanoma cases and healthy controls.

Results

The summed PRS was not significantly different between the three unrelated melanoma groups from Melanoma Institute Australia (MIA) and Princess Alexandra Hospital (PAH) (Supplementary Material, Fig. S1). Hence, for subsequent analyses, we combined the unrelated melanoma cases into a single group; the distribution of PRS for each group is reported in raw values Supplementary Material, Table S1a and as deltaPRS relative to controls in Supplementary Material, Table S1b. Family members with melanoma from the Queensland Familial Melanoma Project (QFMP) had significantly higher scores compared to both unrelated melanoma cases and healthy controls (Fig. 1). All genotyped members of the families were included in these analyses; to confirm that inclusion of multiple family members was not impacting our interpretation, we also restricted the analysis to one case sample per family, and used the average of the PRS of genotyped family cases. In this sensitivity analysis, melanoma family cases were still enriched for polygenic risk for melanoma compared to the unrelated cases (deltaPRS = 0.505, SE = 0.036, N = 658; Holm corrected pairwise t-test P-value = 0.0015). Unaffected members of melanoma families had significantly higher PRS than healthy controls and were not significantly different from either the unrelated or family-based melanoma cases, possibly due to the relatively small numbers in this group (Fig. 1).

Figure 1

Distributions of PRSs by group PRS are plotted using dot-density binning, with each orange point representing a single person’s PRS. The left y-axis reports the mean 0-centered PRS values while the right axis displays the melanoma OR for a 1 SD change in the PRS (see Materials and Methods). The distribution of the PRS have been summarized using a blue layered notched whisker plot, with the midpoint the median and the notches an approximation of the 95% confidence interval of the median (see Martials and Methods). Each boxed area covers the first to third quartile of the group, with the whiskers extending out to 1.5 × the interquartile range. Untransformed values relating to this distribution are tabulated in Supplementary Material, Table S1a and transformed values in Supplementary Material, Table S1b. Bonferroni-corrected pairwise t-tests P-values (6 comparisons in total) are displayed where the adjusted two-sided P-value is <0.05. Total sample numbers are reported for each group.

Scores by familial risk class

The distribution of PRS for affected family members from medium-risk families (two family members with melanoma) did not differ significantly when compared to those from high-risk families (three or more members with melanoma; Bonferroni-corrected t-test P = 0.94; unadjusted P = 0.11, Supplementary Material, Table S3a and b, Fig. S2). In turn, the distribution of PRS in unaffected family members did not differ significantly to melanoma cases from either medium or high-risk families (unadjusted t-test P = 0.45 and P = 0.07 respectively, corrected P = 1.00 and P = 0.70).

Whole genome sequencing of melanoma families

Whole genome sequencing of 51 people from 21 families with melanoma and a relatively low PRS (Materials and Methods) was also performed to investigate whether a high- or medium-penetrance mutation might be identified. In one family for whom three individuals with melanoma were available for sequencing, we identified a pathogenic founder CDKN2A variant, p.G101W. This variant, first described in 1994 (27), is one of the most frequently reported deleterious founder missense variants reported in CDKN2A (28) and results in impaired binding of p16 to CDK4 (29). No further segregating mutations or near segregating (to take into account phenocopies in the high UVR exposed Australian population) were identified in the remaining 20 families.

Discussion

In a previous bipolar disorder family study, affected members had disease PRS similar to the PRS of cases recruited from clinics without regard to their family history of disease (25). In contrast, in our melanoma study, we found that melanoma cases from melanoma-prone families had higher mean PRS than unrelated melanoma cases from clinics and hospitals (Fig. 2). This difference may reflect differences in the genetic architecture of two diseases or may simply be due to the use of only a single pedigree in the bipolar disorder study, and a larger sample of familial bipolar cases may find polygenic risk is enriched in both diseases. In terms of genetic architecture, melanoma has a number of common variants with a relatively large odds ratio (OR) for a complex trait (Supplementary Material, Table S4). However, our results were not driven by GWAS variants of large effect, with the mean of affected family members’ PRS higher than that of unrelated controls in the absence of genetic variants from the MC1R or CDKN2A GWAS loci (data not shown).

Figure 2

Model evaluation for PRS constructed from SNPs binned by ɑ P-value threshold. The base model is a logistic regression of the first 10 principal components on unrelated melanoma cases versus controls. For each subsequent model, the PRS derived from the specified melanoma meta-analysis P-value range of SNPs was added to the regression (see Materials and Methods). For each PRS, we plot the model Nagelkerke R2 and report the number of included independent variants. The logistic regression P-values for all models including genome-wide variants were <1 × 10−4. We also repeat the analysis without the MC1R (no MC1R) or CDKN2A (no CDKN2A) regions (see Materials and Methods); logistic regression P-values for each model including a PRS with these regions dropped were all <0.001. This data are tabulated in Supplementary Material, Table S5.

Melanoma-free QFMP members had a mean score that was higher than healthy, unrelated controls (deltaPRS 0.347, SE 0.104, P = 0.0085; Fig. 2, Supplementary Material, Table S1b). The mean of melanoma-free family members was essentially the same (P > 0.05) as the unrelated melanoma cases (deltaPRS 0.354, SE 0.0222), and not significantly different to affected family members (deltaPRS 0.523, SE 0.0282). This suggests that even (currently) unaffected family members are at an elevated genetic risk for melanoma (as might be expected).

The vast majority of melanoma families in the QFMP were pre-screened for known high-penetrance melanoma mutations in CDKN2A and CDK4. This study only includes those without identified mutations; 633 of the included 660 families were screened. We hypothesized that those families with a low polygenic risk for melanoma (i.e. similar to the distribution of PRS in healthy controls) might be more likely to harbor a (novel) high-penetrance mutation, and indeed recent work has shown that melanoma families with a high-penetrance mutation may have lower polygenic risk for melanoma compared to families without (26). Sequencing of 21 families selected from those with low mean-family PRS (that is, the average PRS across genotyped, affected, family members was more similar to healthy controls) identified a single family carrying a previously undetected high-penetrance deleterious CDKN2A mutation, p.G101W (27–29).

Possible explanations for these densely affected families with a low melanoma PRS but no detected high-penetrance mutation include the following: (i) there are shared environmental factors within the family, (ii) there are rare moderate-penetrance mutations that we were underpowered to find by sequencing, (iii) although the families were selected to have a low PRS, the PRS values are not very low when benchmarked against the PRS values seen in healthy individuals (Supplementary Material, Fig. S3) and (iv) since our PRS only indexed a minority of the genetic variation in melanoma risk, these families may be carrying a large number of common low-penetrance risk variants, which are not in our current PRS. As GWAS sample sizes increase, we will identify more low penetrance variants, and in the future, we will be able to better assess point (iv).

As a potential limitation, the unrelated melanoma cases were used to select the best PRS (score 3, Supplementary Material, Tables S4 and S5), which may have led to a slight inflation in the significance of their PRS compared to healthy controls. However, the R2 difference between the chosen score and genome-wide significant SNPs alone was negligible, and the mean PRS of melanoma cases from melanoma families was still significantly higher than unrelated cases. To further explore this, we repeated the analyses using score 1, limited to SNPs in the PRS reaching P < 5 × 10−8. As shown in Supplementary Material, Fig. S4, cases from melanoma families still had significantly higher scores relative to unrelated melanoma cases (Bonferroni-corrected P = 0.014) indicating that our findings are not driven by selection of PRS in this manner. As a further limitation, pruning to genotyped SNPs across included arrays while reducing potential artefacts and biases may have reduced our power; however, we were still able to detect significant differences between analysis groups (e.g. Fig. 1).

Familial clustering in the context of a higher polygenic risk for melanoma than seen in unrelated melanoma cases from clinics and hospitals suggests that some familial clustering may be due to many common variants of low and moderate effect. Indeed, previous analysis of the QFMP found that high-risk families were more likely to have melanoma risk phenotypes that are themselves polygenic (e.g. high nevus count) (1,23,24). We explored if the number of affected family members was related to degree of polygenic risk for melanoma. The mean PRS for high-risk families (three or more affected members) while higher was not significantly different from those drawn from families with only two (Supplementary Material, Table S3a and b, Fig. S2). A larger dataset may be required to determine if there is a relationship between the number of affected family members and degree of enrichment for PRS. Regardless, high polygenic load in families does not rule out a role for segregating mutations as a high polygenic load may increase the penetrance of a rare variant (e.g. of a moderate-penetrance variant proposed above in example (ii)) such that a greater number of family members are affected.

This work provides further evidence that cutaneous melanoma, like other common complex disorders, can arise from multiple pathways including rare high-penetrance mutations as well as via a combination of large numbers of alleles of small effect.

Materials and Methods

Samples

Available to us were two large samples of unrelated melanoma cases, one large family study of melanoma and a sample of healthy controls (Table 1).

Table 1

Sample counts pre quality control (QC) and post-QC for QFMP melanoma family cases and unaffected members

Samples	Family-based		Unrelated
	Melanoma cases	Unaffected	PAH	MIA 1	MIA 2	Healthy controls
Pre-QC	1342	127	236	1439	301	1799
Post-QC	1292	111	232	1373	287	1774
Male	596	51	118	826	182	823
Female	696	60	114	547	105	951

Unrelated melanoma cases from the PAH and MIA are also reported. MIA samples were genotyped across two batches, with the counts shown separately.

Unrelated melanoma cases—MIA

Patients with histopathologically confirmed melanoma (N = 1740) with germline DNA available from blood samples in the MIA Biospecimen Bank (protocol HREC/10/RPAH/530, HREC/11/RPAH/444) were identified. Included samples were collected between April 2000 and February 2016. Further inclusion and exclusion criteria are reported in the GWAS cleaning section. Ethical approval was managed by the Human Research Ethics Committee of the Royal Prince Alfred Hospital (Sydney, New South Wales, Australia), and all biospecimens and patient data were collected with written informed consent.

Unrelated melanoma cases—PAH

Patients whose locally invasive primary melanoma was confirmed by histopathology between 1994 and 2007 were recruited with informed consent from the Multidisciplinary Melanoma Clinic at the PAH and related clinics (N = 236). Ethical oversight and management were through the Metro South Hospital and Health Service Human Research Ethics Committee (protocol HREC/12/QPAH/349, Woolloongabba, Queensland, Australia).

QFMP

Index cases for 1912 melanoma families recruited to the QFMP were ascertained by contacting all individuals with a histologically confirmed melanoma in the Queensland Cancer Register diagnosed from 1 January 1982 to 31 December 1990 (30). An iterative and sequential sampling scheme was used where if a first degree relative had melanoma, then the first-degree relatives of that new case were also included. Informed consent was a requirement for study/cohort inclusion, and all work was overseen and approved by the Human Research Ethic Committee of the Queensland Institute of Medical Research (now QIMR Berghofer Medical Research Institute), University of Sydney, University of Melbourne and cancer registries of New South Wales, Victoria and Queensland.

At the time of recruitment, the majority of QFMP families (1403 of 1912 families) had only one case member. Recruited families were ranked as low-, medium- or high-risk using a risk index (T) that factored in the number of confirmed cases of melanoma versus number of unaffected family members, ages and year of birth (31). All high and medium T-index families were approached for blood samples, as were 250 of the low T-index families; in total, 738 families consented to providing blood. A total of 1469 participants from 680 families provided blood samples from which DNA was extracted. Further inclusion and exclusion criteria are reported in the GWAS cleaning section (Table 1). Subsequent data linkages identified additional cases of melanoma within these families, resulting in families increasing their T index class. As a result, we adopted a more simplified risk class based on simple counts of affected family members.

Genotyped QFMP families were classified for PRS analyses into three classes based on their current number of affected family members: low-risk (one affected member), medium-risk (two affected members) and high-risk (three or more affected members). A single family was dropped as it had missing total affected family count data. Following data cleaning (see GWAS cleaning section below), genotyped samples were from a single low-risk family (one unaffected family member genotyped), 207 medium-risk families and 453 high-risk families. The single low-risk family was dropped from all analyses. The final count of included unaffected family members was 111 (Table 1). For affected family members, 261 were from the 207 medium-risk families and 1031 from the 453 high-risk families. For the 111 unaffected family members, one was from a medium-risk family, while 110 were from 110 different high-risk families (i.e. one per family). A more detailed breakdown of the number of individuals genotyped from families is presented in Supplementary Material, Table S2.

Healthy controls without melanoma

A total of 1799 unrelated individuals who self-reported as being free from melanoma and had no family history of melanoma were drawn from the Brisbane Adolescent Twin Study (32,33). The majority (80%) of samples were parents of twins, with the remainder being a single individual from twin or sibling pairs. Ethics approval was overseen by the Queensland Institute of Medical Research (now QIMR Berghofer Medical Research Institute) Human Research Ethics Committee (Herston, Queensland, Australia).

GWAS genotyping

Unrelated melanoma cases from MIA and PAH and QFMP samples were genotyped in two batches utilizing the Illumina Oncoarray (San Diego, CA, USA). The first batch included the majority of the MIA samples (N = 1439) and the PAH samples (N = 236), and the second batch included the 1469 QFMP samples and a further 301 samples from MIA. Healthy melanoma-free controls were genotyped using the Illumina HumanHap610 arrays (San Diego, CA, USA) (34).

GWAS cleaning

Raw genotypes were cleaned and exported using Illumina GenomeStudio/BeadStudio (San Diego, CA, USA). Following protocols published by the Oncoarray Consortium (35), oncoarray genotype data were cleaned of poorly performing assays, and variant positions were updated and aligned where required. Subsequently, oncoarray and Illumina HumanHap610 GWAS data were separately aligned to 1000 genomes and then cleaned using PLINK (v1.90b5.4) (36,37). Variants were filtered out for having a minor allele frequency <0.01 or for Hardy–Weinberg equilibrium P-value <5 × 10−4 in controls or <5 × 10−10 in cases. Individuals were retained if they had calls for at least 97% of variants, and their heterozygosity values were between −0.05 and 0.05. Samples were also discarded if their genetically determined sex did not match phenotype records. Principal component analysis in PLINK including reference populations from the 1000 Genomes Project was used to identify and exclude individuals with non-European ancestry (>6 SD from the centroid of principal components 1 and 2). For the unrelated melanoma cases and healthy controls, identity by descent was used to remove related individuals such that no pair existed with pihat >0.15 within and across groups. Within families, identity by descent was used to confirm reported relationships.

Merging genotyped data

Following GWAS cleaning, there were 1892 unrelated melanoma cases genotyped on the Oncoarray, 1660 from MIA and 232 from PAH. From the melanoma families, there were 1292 cases with melanoma and 111 unaffected members, also genotyped on the Oncoarray. A total of 1774 healthy controls without melanoma genotyped on the Illumina HumanHap610 array remained after cleaning.

The genotype data for the Oncorray (SNPs = 397 083) and Illumina HumanHap610 (SNPs = 494 381) were filtered to the autosomal SNPs in common across the two arrays. The combined GWAS data were cleaned again (missing genotypes/calls <3%, Hardy–Weinberg equilibrium P-value >1 × 10−6) using PLINK v1.90b6.10. One further MIA case from the second batch was removed for missing genotype calls >3% in the merged set. The final, merged, cleaned set included 135 325 SNPs, 3184 individuals with melanoma, 113 unaffected family members, and 1774 healthy control samples (Table 1).

Melanoma risk effect sizes

The published melanoma risk meta-analysis included two GWASs that incorporated probands from QFMP families (17). These overlapping GWAS were removed from the meta-analysis, and it was re-run prior to PRS calculations; contributing datasets were otherwise as previously reported (17). A fixed effects meta-analysis was performed using GWAMA v2.1 (38). A total of 133 717 SNPs overlapped between the genotype data and meta-analysis results.

A total of 5000 randomly selected European Ancestry (as measured by principal component analysis) UK Biobank (39) individuals were used for linkage-disequilibrium (LD) calculations. Haplotype reference consortium-imputed best guess variants (imputation quality rsq and best guess threshold both ≥0.3) were retained if their call rate was ≥97%, had a minor allele frequency ≥0.001 and a Hardy–Weinberg equilibrium P-value ≥1 × 10−6. Of the 133 717 variants overlapping between the genotype and meta-analysis data, 131 957 were in the LD reference panel. To explore the impact of MC1R (chromosome 16q) and CDKN2A (chromosome 9p) regions, which have common variants with relatively large effect sizes for a GWAS, on the PRS, we also generated SNP lists excluding all SNPs from 85 megabase (mb) to the end of chromosome 16 (328 variants) and 18–25 mb of chromosome 9 (505 variants).

LD clumping using the LD reference panel was performed in two stages using PLINK (36,37). The first used 250 kilobase windows at LD r2 0.5 (100 373 clumps formed from the 131 957 top variants), and the second round used 10 mb windows at LD r2 0.025 (27 895 clumps formed from 100 373 top variants). There were proportionally smaller counts of variants when MC1R or CDKN2A were excluded. The PRS for each sample was calculated using the score function in PLINK (36,37).

Selecting P-value ranges for PRSs

Clumped meta-analysis SNPs were binned by ɑ P-value threshold to calculate a series of 10 scores. The best performing PRS (score 3, including independent SNPs with a P-value ≤ 5 × 10−6) was identified using their Nagelkerke R2 value from a logistic regression on melanoma status for unrelated melanoma cases versus healthy controls with the first 10 principal components included as covariates (Fig. 2, Supplementary Material, Table S5). The most strongly associated independent (lead) SNPs with a P-value <5 × 10−6 are reported in Supplementary Material, Table S4. Due to the use of variants only in common across the genotyping arrays, meta-analysis and LD panel, not all loci reported in (17) are present, e.g. SLC45A2.

Selection of melanoma families for sequencing

To explore whether families with melanoma and a relatively low PRS for melanoma (i.e. similar to the distribution of PRS in healthy controls) harbored a high-penetrance mutation, we calculated the mean PRS for each family, dividing the sum of the PRS of genotyped, affected, family members by the number of affected genotyped members. Starting at the lowest end of the mean PRS distribution (lowest being the least polygenic load for melanoma), we selected 51 individuals from 21 families who (i) had multiple members genotyped in this study and (ii) had sufficient excess DNA available for further analysis/sequencing in the future. The PRS for sequenced family members is displayed in Supplementary Material, Fig. S3.

Whole genome sequencing was performed by Macrogen (Seoul, South Korea) on the Illumina HiSeq 2000 platform (San Diego, USA), with a mean coverage of 60–96×. Sequence output was mapped to the UCSC human genome reference build 19 by the Burrows-Wheeler Aligner (40). SNPs were detected using bcftools and SAMtools mpileup with disabled BAQ computation (41), in/del were detected with pindel (42) and both were annotated to dbSNP144, including ExAc population frequencies (43) by ANNOVAR (44). Variants were then filtered on quality score (>70) and alternate read counts (>2 and >20% of all reads at a given position) for stringency of data. Sanger sequencing was performed to confirm identified variants of interest.

Additional methods

Unless otherwise noted, Bonferroni-corrected pairwise t-tests in R were used to test for significant differences between groups, with the significance threshold for multiple testing corrected P-values set to <0.05.

PRS distributions were plotted using the ggplot2 and ggpubr R packages (45,46). For plotting purposes, PRS distributions were rescaled to have mean 0 and with a standard deviation of 1 unit over the whole population by subtracting the population mean PRS from each PRS value, and then dividing by the standard deviation of the PRS for the entire dataset. The mean of the healthy controls was then subtracted from all transformed PRS values such that the plots are centered with the mean of the control population = 0. To summarize distributions, we used a notched whisker plot where the midpoint is the PRS median with the notches corresponding to 1.58 × the interquartile range divided by the square root of the sample number, which is an approximation of the 95% confidence interval of the median. The boxed area is the first to third quartile of the PRS distribution. Whiskers extend as far as 1.5 × the interquartile range with points beyond this plotted individually. In addition to plotting PRS values scaled to mean 0 as above on the left y-axis, on the right axis, we also display the OR for a 1 SD change in PRS. Transformed PRS together with the first 10 principal components were used in a logistic regression with unrelated melanoma samples versus healthy controls to determine a melanoma OR of 1.46 per SD of the PRS.

Click here for additional data file.