The mutational landscape of melanoma brain metastases presenting as the first visceral site of recurrence

Brain metastases are a major cause of melanoma-related mortality and morbidity. We undertook whole-exome sequencing of 50 tumours from patients undergoing surgical resection of brain metastases presenting as the first site of visceral disease spread and validated our findings in an independent dataset of 18 patients. Brain metastases had a similar driver mutational landscape to cutaneous melanomas in TCGA. However, KRAS was the most significantly enriched driver gene, with 4/50 (8%) of brain metastases harbouring non-synonymous mutations. Hotspot KRAS mutations were mutually exclusive from BRAFV600, NRAS and HRAS mutations and were associated with a reduced overall survival from the resection of brain metastases (HR 10.01, p = 0.001). Mutations in KRAS were clonal and concordant with extracranial disease, suggesting that these mutations are likely present within the primary. Our analyses suggest that KRAS mutations could help identify patients with primary melanoma at higher risk of brain metastases who may benefit from more intensive, protracted surveillance.


Extraction of clinical details
Patient demographics, primary tumour characteristics (date of primary diagnosis, Breslow thickness, ulceration, mitotic rate, N stage), date of diagnosis of brain metastasis, date of neurosurgical resection, sites of surgically resected brain metastases as well as date of last followup or death, were extracted primarily from prospectively maintained clinical research databases and, if needed, through further review of the clinical record. None of the patients received systemic therapies (included targeted and immunotherapies) prior to the resection of brain metastases.

Extraction and quality assessment of DNA and RNA
For each tumour, a hematoxylin and eosin-stained (H&E) section was reviewed by a consultant histopathologist (PF, RAS, CT and PE) for confirmation of the histopathological diagnosis and to identify appropriate regions for DNA extraction. All samples were obtained as either 1.0 mm diameter cores or 4um tissue sections micro-dissected from the original FFPE block. Germline DNA was extracted from micro-dissected adjacent normal tissue where available (n=36).
Genomic DNA was extracted using the QIAamp FFPE Tissue kit from Qiagen according to manufacturer's instructions.

Whole exome-sequencing of the discovery cohort
Exome capture was performed using the Agilent SureSelect Human All Exon V5 platform. Pairedend sequencing was performed using the Illumina HiSeq (Illumina, San Diego, CA, USA) platform at the Wellcome Sanger Institute to generate 75 bp paired-end reads. Sequencing reads were aligned using BWA-MEM (v0.7.17-r1188) 1 to the human reference genome hs37d5. PCR duplicates, secondary read alignments, and reads that failed Illumina chastity (purity) filtering were flagged and removed prior to running variant and copy number calling. For tumour samples (n=16) where no matching germline DNA was available, a panel of 39 FFPE-extracted normals was used to filter germline variants as well as artefacts. To ensure that tumour and matched normal samples had the best reciprocal genotype match, SAMtools mpileup, followed by BCFtools gtcheck were run to detect sample concordance, potential sample swaps and contamination. The resulting median sequencing coverage in the brain metastatic samples from the discovery cohort (excluding PCR duplicates) was 48x (range 11-134x) in the tumour samples and 47x (range 6-149x) in the germline samples.
Prior to running Strelka, Manta (1.5.0) 4 was run and candidate large indels were used as input to Strelka. The minimum base quality score for somatic and germline variant calling was set to Phred 30. The Ensembl Variant Effect Predictor 5 was used to predict the effect of variants on genes and proteins, relative to the gene build in Ensembl version 97. To remove artefacts, MuTect variant calls were filtered using a tiered approach, such that only the SNVs that met the following criteria were reported (1) Lower coverage samples (<40x); Variant Allele Frequency (VAF) >0.1. (2) All other depth samples; requirement of at least 2 alternative bases on each strand, total depth >=5, and VAF>=0.1 Variant calls found in the gnomAD database 6 were removed if the global population variant allele frequency was greater than or equal to 0.01. Mutational load was calculated as the number of non-synonymous mutations per Mb. The alignments for all reported variants in melanoma driver genes (Fig. 1a) were visually inspected using JBrowse and IGV.
Of note extracranial tumour tissue was available and whole-exome sequenced for two hotpot KRAS mutant patients. Patient MBM_Disc_23 sample PD42097c -regional lymph node also carried a KRAS Q61R mutation, and patient MBM_Dis_40 samples PD31211d/PD31211eextracranial skin metastases also carried KRAS G12D mutations, both concordant with the KRAS mutation in the matched brain metastases (raw data deposited in Supplementary Data).

Copy number profiling of the discovery cohort
Sequenza (v2.1.2) 7 was used to estimate tumour cellularity and ploidy from 30 tumour-normal pairs in the discovery WES cohort (36/50 discovery cohort samples were paired with germline; minus 5 samples with coverage < 20x and one sample with a noisy CNV profile PD36788e), as well as to calculate allele-specific copy number profiles. For each sample, the default best-fit solution was used for the cellularity and ploidy estimates. To call copy number gain or loss in each sample, the neutral copy number was set as the weighted mean copy number of the segments, rounded to the nearest whole number.

Targeted panel sequencing of the external validation cohort
Panel sequencing of the 18 samples in the external validation was performed using custom pulldown and sequencing of 549 key melanoma and related cancer driver genes (Supplementary Table 3

MSK-IMPACT datasets
In order to compare the mutational rates and frequencies to our dataset (represented by one sample per patient), the clinical and mutation data from The Cancer Genome Atlas (SKCM-TCGA) 8 , were downloaded from the cBioPortal 9,10 . Samples were filtered to a single sample per patient giving a total of 358 samples from 358 patients, of which all samples had appropriate SNV data. The demographic characteristics of these 358 SKCM-TCGA patients were largely comparable to our cohort, whereby (n=220, 62%) were male and the median age was 57 years (95% CI 47-71 years) (Supplementary Table 2). The majority of these samples (n=283, 79%) were classified as 'metastasis', as opposed to 'primary' (n=75, 21%). The tumour sites were mainly from skin, subcutaneous or nodal metastatic sites, including 144 (40%) classified as from 'extremities', 134 (37%) from 'truncal' locations, 26 (7%) from the 'head and neck', 23 (6%) were 'regional lymph nodes' and the remainder were from other (less frequent) metastatic sites. Only 5 tumour samples in SKCM-TCGA were classified as from the 'brain' and these were excluded when we subsetted to the 'extracranial' metastatic melanoma comparator (n=274).
The MSK-IMPACT dataset was extracted from the publication by Zehir et al 11 . Samples were also filtered to a single sample per patient giving a total of 186 samples from 186 patients, all of which had appropriate SNV data. These samples were all labelled as 'cutaneous melanomas', and the majority (n=161, 87%) were classified as 'metastasis' as opposed to 'primary' (n=25, 13%). All primary tumours in this dataset (n=25) were classified as from 'skin' whereas the majority of metastatic samples were from 'regional lymph node' (n=25, 16%), 'lung' (n=24, 15%) and 'intransit' (n=17, 11%) and the remainder were from other (less frequent) metastatic sites.
Both SKCM-TCGA and SKCM-MSK-IMPACT datasets included only cutaneous melanomas, in particular any melanomas within these datasets from acral, mucosal and other rarer sites were excluded. All SNVs reported from these datasets were non-synonymous mutations. We utilized the mutation calls provided by these resources and did not perform any additional filtering. The somatic mutational rate was calculated as the number of non-synonymous mutations divided by 30Mb, assuming that an average exome has 30Mb in protein-coding genes with sufficient coverage.

Copy number profiling of the SKCM-TCGA cohort
Copy number calls were generated from single nucleotide polymorphism (SNP data) within the SKCM-TCGA dataset using allele-specific copy number analysis (ASCAT version), as previously described. 12

Statistical methods
All statistical analysis and graphics were generated using R version 3.

Supplementary data
All the whole exome and targeted sequencing data (including raw sequencing files, variant calls and copy number calls) have been deposited at the European Genome-Phenome Archive Orthogonal validation for discovery cohort, raw sequencing files: EGAD00001005984 External validation cohort, raw sequencing files: EGAD00001005985

Supplemental Table Legends
Supplementary Table 1

KRAS mutations (including both patients from the discovery and validation cohorts). Of note all
hotspot KRAS mutations were identified in patients with either thin (T1/T2) or no prior history of primary melanoma, and from primary tumours in chronically sun-exposed locations. Mutations in KRAS also had a high variant allele frequency, indicating that these likely represent clonal driver mutations, which is further supported by their concordance in extracranial metastases.

Supplementary Table 6. Clinical and mutational characteristics of the patients with hotspot
KRAS mutations in SKCM-TCGA (n=6) and SKCM-MSK-IMPACT (n=2). KRAS mutations were mutually exclusive to BRAF/NRAS/HRAS hotspot mutations in both these datasets.

Supplementary Table 7. Clinical characteristics of hotspot KRAS-mutant (n=4) vs KRAS-wild
type (n=46) patient within the discovery cohort. Supplementary Fig. 1. Odds ratio plot showing the odds of observing a KRAS mutation in the three reference datasets relative to the early brain metastases discovery cohort (baseline). Dots correspond to logistic regression odd ratio estimates and range to the corresponding 95% confidence interval. The odds of observing a KRAS mutation in a given sample within the discovery cohort was~4-fold higher than in these three reference datasets, this difference is significant.

Centre
Auckland region Queensland region Supplementary Fig. 2. Tile plot of mutated positions within the RAS signalling genes in the external validation cohort (n=18). The KRAS G13C mutation was mutually exclusive from both BRAF and NRAS hotspot mutations, in keeping with the findings from the discovery cohort (Fig. 1B). Orange line: survival from resection of brain metastasis, defined as the time from the resection of brain metastasis to last follow-up (censored) or death from any cause as a function of time. Median 12.0 months, 95% CI 6.0-62.0 months. Data from resection of brain metastasis was available on 47 patients. Blue line: Time to progression from primary tumour to brain metastases. B) Timeline summary for (n=40) patients. Patients are ordered according to the age at primary tumour resection (y-axis). C) Impact of KRAS mutational status on overall survival from primary tumour as a function of time.