Germline variants are associated with increased primary melanoma tumor thickness at diagnosis.

Germline genetic variants have been identified, which predispose individuals and families to develop melanoma. Tumor thickness is the strongest predictor of outcome for clinically localized primary melanoma patients. We sought to determine whether there is a heritable genetic contribution to variation in tumor thickness. If confirmed, this will justify the search for specific genetic variants influencing tumor thickness. To address this, we estimated the proportion of variation in tumor thickness attributable to genome-wide genetic variation (variant-based heritability) using unrelated patients with measured primary cutaneous melanoma thickness. As a secondary analysis, we conducted a genome-wide association study (GWAS) of tumor thickness. The analyses utilized 10 604 individuals with primary cutaneous melanoma drawn from nine GWAS datasets from eight cohorts recruited from the general population, primary care and melanoma treatment centers. Following quality control and filtering to unrelated individuals with study phenotypes, 8125 patients were used in the primary analysis to test whether tumor thickness is heritable. An expanded set of 8505 individuals (47.6% female) were analyzed for the secondary GWAS meta-analysis. Analyses were adjusted for participant age, sex, cohort and ancestry. We found that 26.6% (SE 11.9%, P = 0.0128) of variation in tumor thickness is attributable to genome-wide genetic variation. While requiring replication, a chromosome 11 locus was associated (P < 5 × 10-8) with tumor thickness. Our work indicates that sufficiently large datasets will enable the discovery of genetic variants associated with greater tumor thickness, and this will lead to the identification of host biological processes influencing melanoma growth and invasion.

Introduction

Cutaneous melanoma (hereafter melanoma) is a potentially fatal skin cancer resulting from the uncontrolled growth of melanocytes, the pigment-producing cells of the skin. In 2019, in Australia, there were estimated to be over 38 000 new cases of melanoma (15 000 invasive) and nearly 2000 deaths, and in the United States an estimated 192 000 cases of melanoma (97 000 invasive) with over 7000 deaths (1,2).

Tumor thickness, measured from the top of the epidermal granular layer to the deepest point of tumor invasion, is the single strongest predictive factor for mortality in clinically localized primary cutaneous melanoma and widely used in disease staging and prognostication (3–5). It is not yet known whether primary melanoma thickness at the time of initial diagnosis is influenced by heritable factors. Beyond the time of presentation for biopsy/excision of the lesion, factors such as invasion potential, growth rate and immunosurveillance are likely to be relevant to tumor thickness at presentation, and these may be influenced by host germline characteristics (6–9).

Previous investigations have identified a limited number of potential loci associated with primary melanoma tumor thickness including at MMP1, FGFR4, VDR and loci previously associated with melanoma risk e.g. MC1R (10–16). No loci have been validated independently. Two genome-wide association studies (GWAS) of primary melanoma tumor thickness did not identify any genetic variants that were significant after multiple testing correction (17,18). These findings are summarized in Supplementary Material, Table S1.

Using over 8000 patients with primary cutaneous melanoma from eight cohorts (Table 1), we sought to determine whether heritable genetic factors influence melanoma tumor thickness and the extent of such influence. Due to the possible impact of shared environment on tumor thickness, and the increased surveillance from diagnosis of melanoma in a close relative, our analysis used distantly related individuals to remove this potential for confounding (19–24). Confirmation that germline genetic variation influences primary melanoma tumor thickness would imply that, with a sufficiently large dataset, specific genetic variants associated with greater tumor thickness can be identified. Identifying these genetic variants, and the genes they influence, will lead to a better understanding of the host biological processes that are important for the growth and invasion of melanoma.

Table 1

Population demographics

	AMFS	CAMBRIDGE	EPIGENE	LEEDS	MDACC	MIA	Q-MEGA 610 k	Q-MEGA omni	WAMHS
Genotyped	549	492	787	2192	1982	1745	912	656	1289
Age, sex	532	491	729	1898	1523	1240	543	457	1253
1 Tumor thickness	529	491	729	1895	1523	1177	534	451	1195
2Post-QC	528	488	728	1888	1522	1177	534	451	1189
3Thickness < 0.8 mm	350	223	539	200	523	335	351	315	711
3Thickness ≥ 0.8 mm and ≤ 1.0 mm	79	72	61	343	180	138	74	47	162
3Thickness > 1.0 mm and ≤ 2.0 mm	68	109	82	724	405	354	76	59	183
3Thickness > 2.0 mm and ≤ 4.0 mm	18	62	38	411	265	234	25	20	87
3Thickness > 4.0 mm	13	22	8	210	149	116	8	10	46
Mean tumor thickness mm (SD)	0.861 (1.088)	1.288 (1.252)	0.791 (1.184)	2.145 (2.013)	1.932 (2.548)	1.950 (2.399)	0.821 (0.875)	0.833 (1.005)	1.104 (1.770)
4Transformed tumor thickness mean (SD)	−0.479 (0.721)	−0.104 (0.838)	−0.593 (0.730)	0.476 (0.732)	0.190 (0.927)	0.255 (0.891)	−0.506 (0.739)	−0.518 (0.737)	−0.352 (0.855)
5Women (%)	336 (63.6%)	257 (52.7%)	231 (31.7%)	1082 (57.3%)	634 (41.7%)	474 (40.3%)	287 (53.8%)	223 (49.5%)	496 (41.8%)
Men	192	231	497	806	888	703	247	228	693
Genomic inflation λ	0.971	0.995	0.994	1.000	0.995	1.000	0.991	0.961	0.995
6IBD < 0.025	499	483	660	1864	1478	1079	504	420	1138

Number of individuals with tumor thickness as a continuous measurement in mm. 2Number of individuals following quality control and cleaning (Methods; Post-QC). 3 For reference, individuals for each cohort are grouped based on T categories as defined in the eighth edition American Joint Committee on Cancer Melanoma Staging System (4). 4 Continuous tumor thickness measurements in millimeter were natural log transformed. SDs are reported where indicated. 4 Sex ratio is reported for the post-QC set with all phenotypes. 5 Determining the genetic contribution to melanoma thickness requires distantly related participants so only those pairs with identify-by-descent (IBD) pi_hat < 0.025 are included; thickness distribution for the 8125 samples filtered to IBD 0.025 is reported in S3 Table. (Methods). Queensland study of Melanoma: Environment and Genetic Associations (Q-MEGA) samples GWAS were performed by the HumanHap610 (610 k) and Omni1-Quad (omni) genotyped arrays used (Methods). MIA. MDACC. WAMHS.

Results

Contribution of genome-wide germline genotype to variation in tumor thickness

Age- and sex- adjusted residuals of natural log transformed tumor thickness were generated, and GREML-LDMS-I (23) was used to determine the contribution of genome-wide germline variants with a minor allele frequency > 0.001 to tumor thickness variation (heritability, h2SNPs, Supplementary Material, Fig. S1, Supplementary Material, Fig. S2, Methods).

In the combined dataset (N = 8125) with the first six ancestry principal components and an individual cohort membership variable fitted as covariates, tumor thickness h2SNP was 0.266 (standard error (SE) = 0.119, P = 0.0128). Analyzing the age- and sex- corrected residuals of tumor thickness, and fitting a dataset membership covariate, may not have completely accounted for all differences across datasets. However, repeating analysis with rank normalized tumor thickness residuals gave a similar h2SNP = 0.248 (SE = 0.119, P = 0.0188, Supplementary Material, Fig. S3), indicating it is unlikely individual study differences are driving the observed heritability. The estimate from the random-effects meta-analysis of the h2SNP for each discrete dataset calculated separately, while losing power by not using genetic relationships across sample sets, was consistent with the overall result (h2SNP = 0.264) but with wider confidence intervals (SE = 0.158, P = 0.095; Methods, Supplementary Material, Fig. S4).

Acral lentiginous melanoma tends to be diagnosed later and have a greater thickness (25). Excluding 92 acral lentiginous melanoma cases (histology data were available for all but the MD Anderson Cancer Center (MDACC) cohort) did not meaningfully change the results (h2SNP = 0.281, SE = 0.121, P = 0.0104).

Genome-wide association study of tumor thickness

As the h2SNP was significantly different to zero, indicating a role for germline genetic variation in the thickness of primary tumors, we therefore performed a linear regression GWAS of primary melanoma tumor thickness in each dataset (Methods). Following the subsequent meta-analysis of the individual genome-wide results, there was no evidence of genomic inflation (N = 8505; genomic inflation λ = 1.00; Supplementary Material, Fig. S5). Two genetic variants are in linkage-disequilibrium (r2EUR = 0.69) at a single locus on chromosome 11 reached genome-wide significance (Fig. 1). The regional association plot is shown in Fig. 2. The lead variant was rs183471242 (P = 3.56 × 10−9; Table 2). rs183471242 is in an intron of the gene low-density lipoprotein receptor class A domain containing 3 (LDLRAD3; Fig. 2). The distribution of natural log transformed tumor thickness by genotype is displayed in Fig. 3; the distribution of effect sizes across GWAS is shown in Supplementary Material, Fig. S6. Both rs183471242 (G/A) and rs566382949 (C/A) are rare (HRC v1.1 minor A allele 0.0098 and 0.0082, respectively), with the minor allele associated with thicker tumors (Table 2,  Fig. 3). Each minor allele of rs183471242 translates to a 1.423-fold increase on the transformed residuals of tumor thickness.

Figure 1

Manhattan plot of P-values from the meta-analysis of genome-wide association studies of tumor thickness. Negative Log10 of observed fixed-effects meta-analysis P-values plotted by chromosome position. Red line indicates genome-wide significance (P = 5 × 10–8). Sample size is reported in Table 1.

Figure 2

Regional association plot for rs183471242. Negative Log10 of observed fixed-effects meta-analysis P-values plotted by chromosome position. Linkage disequilibrium r2 of plotted SNPs with the lead SNP rs183471242 is displayed. Plot generated using LocusZoom (26).

Table 2

Genetic variants associated with tumor thickness following multiple testing correction

CHR	BP	rsID	EA/NEA	EA FREQ	P	BETA	Q	I 2
rs183471242	11	36 019 025	A/G	0.0098	3.56 × 10−9	0.353	0.81	0
rs566382949	11	36 068 615	A/C	0.0082	1.58 × 10−8	0.382	0.97	0
rsID	AMFS	CAMBRIDGE	EPIGENE	LEEDS	MDACC	MIA	Q-MEGA 610 k	Q-MEGA omni	WAMHS
rs183471242	0.598	0.501	0.443	0.241	0.389	0.459	0.093	0.229	0.369
rs566382949	0.326	0.518	0.418	0.302	0.232	0.541	0.240	0.406	0.435

We report hg19 chromosome (CHR) and base pair (BP) positions for genetic variants (rsID). The effect allele (EA) and non-effect allele (NEA) are provided, as is the HRC frequency of the EA (EA FREQ). The fixed effects meta-analysis P (P), effect size (BETA) on the residuals of fitting age and sex on natural log transformed tumor thickness; the first six ancestry principal components were included as covariates in the regression (Methods). As there is no heterogeneity (Q, I2) the random effects meta-analysis values are identical. We also report effect size estimates for individual datasets. The distribution of these measures is displayed in Supplementary Material, Fig. S6. Q-MEGA samples were analyzed by HumanHap610 (610 k) and Omni1-Quad (omni) genotyped array used (Methods). Sample size is reported in Table 1.

Figure 3

Distribution of natural log transformed primary cutaneous melanoma tumor thickness by rs183471242 genotypes. Data are reported for the combined meta-analysis of all studies for this genetic variant (8505 individuals). For plotting purposes, we display the natural log transform of tumor thickness by genotype rather than the residuals of tumor thickness as used in the heritability estimation and GWAS. In total, there are 8310 homozygous GG, 194 AG and a single AA genotype (HRC v1.1 minor A allele 0.0098). AG and AA genotypes have been plotted together. Distribution of tumor thickness is described using a notched whisker plot (blue) where the midpoint of the notch is the median, and the 95% confidence interval of that median is represented by the notched region. The boundaries of the boxed area extend to the first and third quartiles. The whiskers represent the 1.5 × the interquartile range. The same data are displayed twice with differing secondary layers to display the distribution of tumor thickness residuals; the first is a violin plot, and the second displays the individual results. The individual with the AA genotype had a primary tumor thickness of 0.5 mm, and their position is indicated by an arrow in the second plot.

As a sensitivity analysis, we performed the regression of residual and rank normalized residual tumor thickness on rs183471242 and rs566382949 in the same combined dataset used for the GREML analysis (N = 8125), fitting the first six PCs and study covariates in the model. Both SNPs remained genome-wide significantly associated (P < 5 × 10−8) with tumor thickness.

While none of the previously reported genetic variants associated with primary tumor thickness reached genome-wide significance in this study, the IRF4 functional genetic variant rs12203592 was the most strongly associated (fixed P = 6.50 × 10−4; Supplementary Material, Table S1) (10–18).

Discussion

Given inconsistency in identifying specific germline variants associated with primary melanoma tumor thickness, we first determined whether tumor thickness is heritable (that is, a proportion of phenotypic variance can be explained by additive genetic variants) (10–18). Traditionally, twin- or family-based approaches are used to measure trait heritability, and do so by assessing whether more closely related individuals tend to have more similar phenotypes. However, shared environment and behaviors can confound heritability estimates. For example, diagnosis of melanoma in a relative can lead to increased surveillance and earlier detection, which may influence thickness (19,20,27–33). As a result, tumor thickness can become inversely correlated to degree of genetic relationship, biasing estimates of heritability. An effective alternative approach is to estimate the genetic contribution to thickness using distantly related individuals, removing the potential confounding between thickness and melanoma diagnosis in a relative and/or earlier diagnosis (22). This alternative approach relies on genome-wide single nucleotide polymorphism (SNP) arrays, yielding an estimate of the variation attributable to these genetic variants; this term is referred to as h2SNP. Since arrays do not genotype all genetic variants, h2SNP represents a lower bound of trait heritability.

Using a large collection of distantly related individuals with melanoma (N = 8125 following filtering to remove individuals such that there were no pairs where their relationship was closer than identity-by-descent pi-hat = 0.025, Methods, Table 1), we estimated the h2SNP for tumor thickness to be 0.266 (95% CI 0.033–0.500). This h2SNP result was robust to various sensitivity analyses (Methods, Results). As this result indicates that primary melanoma tumor thickness is heritable, it shows that with sufficiently well-powered cohorts, it will be possible to identify specific genetic variants associated with tumor thickness.

Following confirmation that tumor thickness is heritable, we used the complete dataset to perform a GWAS meta-analysis, identifying two genetic variants in a locus on chromosome 11 at genome-wide significance (P < 5 × 10−8). As tumor thickness values were natural log transformed prior to analysis, the effect size for each minor allele of rs183471242 is associated with a 1.42-fold higher tumor thickness. rs183471242 is rare (minor A allele frequency 0.0098) with 190 heterozygote GA and 1 homozygous AA samples (Fig. 3). This variant appears to be rarer in non-European populations with an MAF of 0.0018 in African populations and not observed in 780 East Asian samples (34). While these genetic variants are in an intron of LDLRAD3, the closest gene is not always the target of associated genetic variants. However, interrogation of public gene expression and annotation resources do not reveal an obvious functional target for these genetic variants (35–38). LDLRAD3 is a member of the LDL receptor family and may play a role in activating genes involved in protein ubiquitination (39). Neither of these two genetic variants were associated (P > 0.05) with ease of tanning, childhood sunburns or skin and hair color in the UK Biobank (data not shown), melanoma risk, nor nevus count (40,41). While none of the previously reported genetic variants associated with tumor thickness reached a P-value<5 × 10−8, the IRF4 functional genetic variant rs12203592 was the most strongly associated (P = 6.51 × 10−4, I2 = 36.4%, Supplementary Table 1). The direction of effect is the same as in the previous report from the Western Australian Melanoma Heath Study; this dataset is included in this analysis (10). The IRF4 SNP rs12203592 has been associated with risk of melanoma (42,43). It is not clear to what extent genetic variants associated with cancer risk influence outcome (e.g. there is no overlap between genetic variants associated with risk or survival for lung cancer (44) or breast cancer (45)), and as such, it is unclear if we would expect the other known melanoma risk variants (e.g. in MC1R) to associate with tumor thickness in a larger analysis.

The thickness of a tumor at diagnosis is the outcome of its duration and rate of growth. In practice, these two factors are hard to assess post detection and unknown before diagnosis, and as a result, the relative importance of either is unclear (6,27–29,33,46). Our work shows that germline genetic variation plays a role in thickness. While the specific mechanisms are yet to be determined, host germline genetic variants may influence tumor growth rate, immunosurveillance or characteristics relating to time to discovery such as appearance of the tumor or the assiduousness of personal skin checks. While different subtypes of melanoma grow at different rates, accurate heritability estimates require very large sample sizes; thus (and to avoid repeated multiple testing), we have not performed subtype specific analyses beyond excluding the acral subtype. Other factors may influence tumor thickness. Higher body mass index (BMI) has been associated with greater thickness; however, fitting a trait that is itself genetically controlled as a covariate can bias genetic association towards false positives and we hence did not include this in our analysis (47–50).

While both our heritability estimate and the genetic variants associated with tumor thickness following correction for multiple testing require replication in a sufficiently large dataset, these findings indicate germline variants impact tumor development, a strong predictor of melanoma outcome. Discovery of specific genetic variants will enable identification of host biological processes influencing melanoma growth and invasion.

Materials and Methods

Dataset descriptions

Overall demographics of contributing datasets are summarized in Table 1. Full descriptions of each contributing dataset can be found in the Supplementary Material, Note.

Quality control, cleaning and imputation of genome-wide genotype data

PLINK v1.9 and R 3.3.2 (51,52) were used for quality control and cleaning of genome-wide genotype data. Genotyped variants were filtered out if they had a minor allele frequency < 0.01, Hardy–Weinberg equilibrium P-value<5 × 10−4 in cancer-free individuals (where melanoma cases were genotyped/cleaned in concert with healthy individuals) or < 5 × 10−10 in those with melanoma. To remove samples with low-quality DNA, or other issues that may impact analysis (e.g. sample contamination or inbreeding) (53), individuals were excluded if they had missingness >0.03, heterozygosity more than three standard deviations (SDs) from the rest of the population, a mismatch between recorded sex and X chromosome determined sex or were considered non-European. European ancestry was determined by principal components analysis using 1000 Genomes European populations as a reference set (54). Individuals more than three SDs from the mean of principal component 1 or 2 were excluded. Relatedness across and within genotyped sets was measured by identity-by-descent pi-hat scores using PLINK (51). For pairs with pi-hat>0.15, the individual with the highest missing genotype rate was dropped.

The Michigan Imputation Server was used to impute individuals (Table 1) to the Haplotype Reference Consortium panel (HRC version 1), and genetic variants with an imputation quality score RSQ > 0.3, minor allele frequency > 0.001 and minor allele count >3 were retained for analysis (55,56). While the approach we have used to determine the genetic contribution to melanoma thickness (see below) is robust to imputation quality, RSQ > 0.3 excludes potentially poorly imputed genetic variants (23).

Cleaning and normalization of tumor thickness

Research participant demographics are presented in Table 1. Primary melanoma tumor thickness measurements were extracted from pathology reports. For participants with multiple melanomas, the first primary tumor was used. The distribution of tumor thickness is reported in Table 1. As tumor thickness is not normally distributed, measurements were natural log transformed. Age- and sex-adjusted residuals of transformed tumor thickness were used for analyses. The distribution of tumor thickness between cohorts is significantly different (P < 2 × 10−16 from ANOVA with age and sex fitted as covariates, see Supplementary Material, Methods). Pairwise comparisons are reported in Supplementary Material, Table S2. As a result, analyses were performed with a cohort variable fitted and as a sensitivity analysis repeated with rank transformation of tumor thickness.

Estimation of heritability for a complex trait

We sought to estimate the heritability (h2, defined as the proportion of phenotypic variance explained by additive genetic variants) for primary cutaneous melanoma tumor thickness. Heritability is traditionally estimated from family data, with high heritability inferred when individuals who are closely related have more similar phenotypes than those who are more distantly related. An alternative approach to family-based methods is to estimate the genetic contribution to thickness using only distantly related individuals to estimate heritability attributable to SNPs, h2SNP (22).

Early methods for estimating h2SNP such as genome-based restricted maximum likelihood as implemented in the Genome-wide Complex Trait Analysis software use only directly genotyped variants (22). However, imputation of genetic variants not present on genotyping arrays can improve discovery power and resolution for standard genome-wide association studies and heritability estimates (21,56). In aggregate, rare variants capture on average one third of the h2SNP, and their inclusion can yield more accurate h2SNP estimates (21,24).

The h2SNP was determined using an extension of the genome-based restricted maximum likelihood approach designed for imputed data, GREML-LDMS-I, as implemented in the Genome-wide Complex Trait Analysis software (21–23). GREML-LDMS was used recently to determine that ~50% of the heritability of height is accounted for by genetic variants with a minor allele frequency between 0.0001 and 0.1 (57). GREML-LDMS-I is made robust to the (unknown) underlying trait genetic architecture by dividing input genetic variants into bins based on their minor allele frequency and degree of linkage disequilibrium (23,24). This is important as incorrect modeling of the underlying genetic architecture can lead to under- or over-estimation of h2SNP (21,24).

Construction of a combined imputed dataset for array-based heritability estimates

Imputed dosage data from the Michigan Imputation Server in variant-call format were converted to best guess format (genotype dosage ≤0.5 as 0, 0.5–1.5 as 1 and >1.5 as 2) and merged into a single combined dataset using PLINK v1.91.4 beta3 (21,51). PLINK binary files were converted into a genetic relationship matrix by Genome-wide Complex Trait Analysis v1.91.4 (21,22). To ensure only distantly related individuals were included, the merged dataset was filtered such that no pair had identify-by-descent pi-hat>0.025 (21,22) (Table 1; final combined sample size of 8125).

We used GREML-LDMS-I to estimate h2SNP for genetic variants binned by minor allele frequency (0.4–0.5, 0.3–0.4, 0.2–0.3, 0.1–0.2, 0.01–0.1, 0.001–0.01) and linkage-disequilibrium quartile. Linkage-disequilibrium scores were estimated for individual genetic variants rather than regions of variants, as this approach produces unbiased estimates in the presence of all possible genetic architectures (21).

Spurious genetic similarities (e.g. ancestry, cohort or batch effects) can bias h2 estimates (21). To address this, we fitted the first six ancestry principal components and a dataset membership variable (a binary yes/no variable for membership in a given cohort) in the GREML-LDMS-I analyses. As a sensitivity analysis, we repeated analyses following rank normalization of tumor-thickness residuals within each cohort.

While the larger merged dataset increased power by leveraging distant genetic relationships across and between individual datasets, it may have introduced bias due to subtle differences in ancestry, tumor thickness or genotyping methods within each individual dataset. As an additional sensitivity analysis, we estimated h2SNP in the individual imputed datasets (filtered to identity-by-descent pi-hat<0.025 within the individual dataset rather than across all datasets) using standard genome-based restricted maximum likelihood (22). The binned GREML-LDMS-I approach was not applied to individual datasets as they were too small. Individual dataset’s h2SNP estimates were then combined using a random effects meta-analysis using the metafor package in R (58). metafor was also used to generate forest plots.

To identify specific genetic variants associated with tumor thickness, within each individual dataset, residual tumor thickness was regressed on imputed genome-wide genotype dosages with the first six ancestry principal components included as covariates. Individual dataset results were further filtered by removing variants with an extreme effect size estimate (>2 or < −2 on the rank-normalized residuals of natural log-transformed tumor thickness; as a reference, in the Australian Melanoma Family Study (AMFS) dataset, this removes genetic variants with an effect size estimate >8 SDs from an effect size of 0). Genome-wide association results for each individual dataset were combined by inverse variance-weighted fixed effects meta-analysis in PLINK v1.91.4 beta3 (51). The total number of genetic variants tested was 13 517 544. Genetic variants were deemed significant if their P-value was less than a genome-wide multiple testing corrected threshold of P < 5 × 10−8.

In a GWAS, the majority of genetic variants are not expected to be associated and their Χ2 distribution should match the null and have a median of ~0.456 (59). We report the genomic inflation λ, the median meta-analysis Χ2/0.456.

Additional analysis methods are reported in the supplementary note.

Ethics approval

Individual studies’ ethical approval details are reported in the Supplementary Material, Note. Overall, approval was managed by the Human Research Ethics Committee of the QIMR Berghofer Medical Research Institute.

Click here for additional data file.

Click here for additional data file.