Whole genome landscapes of uveal melanoma show an ultraviolet radiation signature in iris tumours

Uveal melanoma (UM) is the most common intraocular tumour in adults and despite surgical or radiation treatment of primary tumours, ~50% of patients progress to metastatic disease. Therapeutic options for metastatic UM are limited, with clinical trials having little impact. Here we perform whole-genome sequencing (WGS) of 103 UM from all sites of the uveal tract (choroid, ciliary body, iris). While most UM have low tumour mutation burden (TMB), two subsets with high TMB are seen; one driven by germline MBD4 mutation, and another by ultraviolet radiation (UVR) exposure, which is restricted to iris UM. All but one tumour have a known UM driver gene mutation (GNAQ, GNA11, BAP1, PLCB4, CYSLTR2, SF3B1, EIF1AX). We identify three other significantly mutated genes (TP53, RPL5 and CENPE). Uveal melanoma has a propensity to metastasise. Here, the authors report the whole genome sequence of 103 uveal melanomas and find that the tumour mutational burden is variable and that two subsets of tumours are characterised by MBD4 mutations and a UV exposure signature.

U veal melanoma (UM) arises from melanocytes in the uveal tract, and though less common than cutaneous melanoma, a higher proportion of UM patients die from the disease [1][2][3] . Risk determination of metastatic spread can be obtained through assessment of specific chromosome copy number alterations (CNAs) 4 , gene expression profiles 5 and mutation status of known UM driver genes 6 .
Previous genomic studies have pointed to the existence of four UM categories, strongly linked to prognosis [7][8][9] . Similarly we segregate our tumours into four categories based on CNAs: category 1 are chromosome 3 disomy (D3) tumours lacking chromosome 8q copy-number gain and frequently possessing EIF1AX mutations; category 2 are D3 UM with chromosome 6p and 8q gain and a high proportion of SF3B1 mutations; category 3 are chromosome 3 monosomy (M3) tumours lacking chromosome 8q gain dominated by BAP1 mutations; category 4 UMs are M3 with chromosome 8q gain and BAP1 mutations. These genomic stratifications, while prognostic, are not indicative of treatment responses once progression has occurred.
To improve knowledge of UM genomics and to identify potential therapeutic targets, we conduct a whole-genome sequencing (WGS) study of 103 UM, comprised of 91 primary tumours and 12 metastases, with matched germline DNA. Eightyfour tumours originate from the choroid, eight from the iris, four from the ciliary body, and seven without known primary uveal site (Supplementary Data 1).
UV mutation signatures in iris tumours. Assessment of single base substitution (SBS) signatures revealed SBS5 to predominate in most cases, with strong representation of SBS3, SBS39 and SBS40 in some samples ( Fig. 1) 13 . Two samples were dominated by mutation signature SBS1 (associated with spontaneous deamination) and had correspondingly high TMB (>3 mutations per Mb). As previously observed 14,15 , these features corresponded to the presence of germline loss-of-function (LOF) MBD4 mutations; this takes the published tally of germline MBD4 mutant UM cases to six 14,15 , strengthening its role as a UM predisposition gene. The two UMs with germline LOF BAP1 mutations displayed no unique features. Strikingly, all iris melanomas displayed the genomic features associated with UVR damage (mutation signatures SBS7a, SBS7b and DBS1 13 combined with a high TMB). While exposure to UVR has been suggested as a cause of the elevated UM risk among arc-welders 16 , no molecular evidence of UVR as an aetiological factor has yet been observed in UM sequencing studies. The iris is located anteriorly within the uveal tract and is directly exposed to sunlight that breaches the cornea; we now show that UVR-associated DNA damage results from this exposure and is a unique genomic feature of iris UM.
Patterns of driver mutations and chromosomal aberrations. Assessment of known UM driver genes revealed an oncogenic driver mutation in 102 of 103 tumours: 51 in GNAQ (48 p. Q209P/L, two p.R183Q, one p.G48L), 46 in GNA11 (44 p.Q209L/ P, two p.R183C), five in PLCB4 (three p.D630Y, two p.D630N) and two in CYSLTR2 (p.L129Q). These mutations were generally mutually exclusive except for two PLCB4 p.D630 mutations that co-occurred with GNAQ/GNA11 p.R183H mutations. This cooccurrence between PLCB4 mutation and the minor hotspot p. R183, rather than the stronger oncogenic p.Q209 hotspot mutations, has previously been described in the UM TCGA data 8 . Though not previously highlighted, GNAQ p.G48L mutations have been reported in two UM samples from two separate studies 8,17 , as well as in two hepatic small vessel neoplasms, which are driven by activating GNAQ/GNA14 mutations 18 . This suggests that GNAQ p.G48L is another minor UM oncogenic hotspot mutation. Similar to previous observations 19,20 , BAP1 was the most altered gene in M3 samples (75%), including eight splice site mutations, two germline and 32 somatic LOF mutations, and three cases with disrupted BAP1 due to SV breakpoints. In addition, two D3 tumours carried BAP1 mutations, indicating that although BAP1 inactivation typically occurs after M3 8 , BAP1 aberration can also occur in D3 tumours, which may or may not later undergo loss of chromosome 3. Of note, one of these D3 tumours (MELA_0800) had a low BAP1 variant allele frequency (VAF = 9/80) suggesting it was only present in a subclone, and as copy number tools are not as sensitive as mutation callers, it is possible that the subclone had loss of heterozygosity (LOH) that was not detected by the algorithm. Five tumours had BAP1 mutations and copy-neutral LOH, suggesting that the mutations occurred before WGD in the two tetraploid UMs and before the LOH event in the three diploid UMs. SF3B1 mutations were present in 15 tumours, the majority occurring in category 2, in line with other studies [7][8][9] . EIF1AX hotspot mutations were observed in 19% of tumours. EIF1AX mutations were first discovered in D3 UMs 21 and in the TCGA cohort they were restricted to category 1 tumours (D3 and no 8q gain) 8 , while in the cohort presented by Royer-Bertrand and colleagues two of seven mutations were seen in tumours with M3 and/or 8q gain 9 . Similarly, here six of the 20 EIF1AX mutations (30%) were seen in UM with M3 (n = 5) or 8q gain (n = 1) (Fig. 3a).
Significantly mutated genes. In addition to these known UM genes, three other statistically significantly mutated genes (SMGs) were identified (CENPE, TP53, RPL5; Supplementary Table 1). Three missense mutations and two LOF SVs were identified in CENPE, along with one LOF germline mutation (late truncating) (collectively~5% of UM). An additional six samples were hemizygous for CENPE (Fig. 3a). CENPE is a plus end-directed kinetochore motor protein which plays a critical role in mitosis and chromosome segregation. Knockdown of CENPE has been shown to cause chromosome misalignment and lagging 22,23 . Two of the missense mutations (p.R14W and p.R251W) occurred in the kinesin motor domain (Fig. 3c) at highly evolutionarily conserved regions (ECRs) and the third (p.Q1098P) occurred in a reasonably conserved residue ( Supplementary Fig. 1). An additional two missense mutations were identified in the UM TCGA cohort at p.I1038T (weak ECR) and p.K1821N (reasonable ECR), both also in a coiled-coil domain 24 . CENPE has been shown to interact with CENPF, BUB1B and Aurora B, the latter two being critical in the activation of the spindle assembly checkpoint [25][26][27] . In the UM cohort described here, one sample had a BUB1B missense substitution (p.R691H) within the region reported to directly interact with CENPE 25 ; another had a p.D303E substitution in a highly ECR of the Aurora B catalytic domain. It is possible that disruption of this pathway is responsible for creating genomic instability allowing for chromosome aberrations to occur. Indeed, the twelve UM with CENPE alterations had significantly higher genome percentages with CNAs (Mann-Whitney, P = 0.028, median 23% vs 15%), though this association is confounded since tumours with high CNA generally have more genome-wide regions of LOH. Studies both in cell line and mouse models have demonstrated that CENPE/Cenpe functions in a haploinsufficient manner, with elevated levels of chromosome missegregation observed in heterozygous cells and animals 28,29 . Functional work on CENPE missense mutations is required to determine their impact in UM.
TP53 is commonly disrupted in uveal melanoma. High expression of p53 has been reported in some UM, often associated with histological and clinicopathological features correlated with poor prognosis, but the potential genetic basis of these observations was not assessed [30][31][32] . Somatic mutations in TP53 have been described in two UM. A hotspot mutation p.R175H (n = 1216 in IARC TP53 database) was observed in a hypermutated metastatic UM with deficient MBD4 33 and another hotspot mutation p.M237I (n = 196 in IARC TP53 database) was observed in a UM in a pan-cancer study of metastatic tumours 17 . Here we identified TP53 as an SMG and report six somatic TP53 mutations across four tumours in addition to eight cases of LOH (Figs. 3a, b). One LOH case overlapped with an LOF mutation (p.C277*) resulting in a double-hit in TP53. Another double-hit was seen in a sample with two missense mutations (p.H193R and p.T155I) confirmed as occurring on different alleles by assessing read pairs spanning both mutations. To evaluate the consequence of these mutations, we applied a computational prediction tool, FATHMM 34 , and assessed the results of two comprehensive characterisation studies of TP53 mutations (Table 1) 35,36 . p.H193R is a recurrent hotspot classified as pathogenic by PHANTM, RFS and FATHMM, while the consequence of p.T155I is more uncertain, as the variant is classified pathogenic by PHANTM but predicted to have neutral impact by RFS and FATHMM. Finally, one UM had a LOF mutation (p.R342*, COSM11073) and a missense p.R248Q mutation, both of which frequently occur in malignancies and are classified as pathogenic (Table 1). RNA-seq data revealed one read pair spanning both positions, which contained p.R248Q and was wildtype for p.R342; furthermore, there was a significantly lower VAF at p.R342 (4/78) than at p.R248 (25/68) (two-sided Fisher's exact test, P = 2 × 10 −6 ). These data suggest the two mutations occurred on different alleles, with the majority of the transcripts from the p.R342* allele undergoing nonsense mediated decay. TP53 is a tumour suppressor gene frequently deleted or mutated, resulting in either no production of p53 or the expression of a truncated and unstable protein. The spectrum of a few highly recurrent missense mutations, including p.R248Q, has, however, given rise to hypotheses that these hotspot mutations translate to mutant p53 with gained oncogenic functions 37 . For example, the p.R248Q mutation reported here has been shown to increase the migratory potential of cells in an in vitro model 38 .
RPL5 is significantly mutated in uveal melanoma. SMG analysis also identified RPL5, with truncating mutations in three cases and LOH in an additional 30 cases (Fig. 3a). No UM had a double hit in RPL5, in line with previous studies suggesting RPL5 is a haploinsufficient tumour suppressor, with heterozygous inactivation in glioblastoma (11%), breast cancer (34%) and cutaneous melanoma (28%) 39 . In UM the majority of RPL5 LOH occurred through loss of chromosome 1p, where the gene is located, but in four cases, focal loss occurred, suggesting that RPL5 may drive positive selection for 1p loss. Chromosome 1p loss in UM has been reported as a marker of poor prognosis, independent of M3 12 . RPL5 encodes ribosomal protein L5, which complexes with 5S rRNA and forms an important part of the impaired ribosome biogenesis checkpoint (IRBC). Together with RPL11 and ARF, RPL5 binds to and inhibits MDM2, resulting in p53 stabilisation in response to blocks in ribosome biogenesis and nucleolar stress [40][41][42] . Of note, one truncating mutation in RPL11 was observed in the UM TCGA cohort 8 further supporting the importance of this pathway in UM. Given the link between RPL5 and p53, we tested for an association between aberrations in RPL5 and TP53. Indeed, mutations in these genes were mutually exclusive, but when also considering copy-number loss there was no association. However, as p53 functions in multiple pathways, it may be inactivated in some tumours with defective RPL5 due to selective pressures outside the IRBC response. Interestingly, the IRBC is often triggered by oncogene-induced translational stress [40][41][42] , with oncogenic MYC being shown to significantly increase IRBC activation [41][42][43][44] . In addition to their IRBC role, RPL5 and RPL11 have also been shown to bind to MYC transcripts, mediating RNA-induced silencing and interfering with c-Myc driven transcription 45,46 . Overexpression of MYC is thought to be the driving force behind positive selection of chromosome 8q gains in UM. It is possible that c-Myc overexpression leads to ribosomal stress and IRBC activation, likely inhibiting tumour growth. To overcome this, it may be necessary to disrupt this pathway through mutation/loss of either RPL5, RPL11 or TP53. Supporting this notion, RPL5 and TP53 disruption (including mutations, SV breakpoints and chromosomal copy loss) was more common in cases with 8q gain (48%) than in cases without 8q gain (22%) (two-sided Fisher's exact test, P = 0.01). Tumours with alterations in this pathway were predominantly M3, making it difficult to disentangle the prognostic impact; however, these individuals had poorer prognosis ( Supplementary Fig. 2).
Genomic categories correlate with prognosis. To correlate UM genomic categories with prognosis, time from first presentation to metastasis was examined for all patients. Comparing the four TCGA categories, there was a trend that category 4 tumours had shorter relapse-free survival (RFS) (median: 2.5 years) than UMs in category 3 (median: 7.5 years) and similarly category 2 had shorter RFS (median: 7.0 years) than category 1 (median not reached), but the differences were not statistically significant (P 3vs4 = 0.11, P 1vs2 = 0.28) (Fig. 4a). However, M3 UM had significantly shorter RFS (median: 2.9 years) than D3 UM (median 7.0 years, log-rank test, P = 0.001, Fig. 4b), confirming the prognostic strength of M3/D3 status. Interestingly, while iris melanomas are associated with earlier detection and favourable prognosis, 4/8 iris cases had M3, with two already having progressed to metastatic disease. The high TMB in iris UM reported here suggests they are more likely to respond to immunotherapy, given the observations for MBD4 germline UM patients 14,15 , who have high TMB, a predictor of response to immunotherapy response in cutaneous melanoma and other cancers 47 . In combination with previous reports of metastatic iris UM 48 , these data suggest a subset of iris UM are at high risk for disease progression and will likely respond to immunotherapy in the event of progression and perhaps even in the adjuvant setting.   Somatic mutations. Somatic SNV and indels were detected using an established pipeline 53 in which SNVs were called with qSNP 54 and GATK HaplotypeCaller and indels were detected with GATK 55 . The contribution from different mutation signatures was inferred by approximating (minimising the squared error) the distribution of mutations as a linear combination of COSMIC signature v3 13 with the constraint that contributions were non-negative. SMGs were identified using MutSigCV 1.3.5 (via GenePattern) as well as Oncodrive-fm and OncodriveClust (via IntOGen) 56 . A Benjamini-Hochberg adjusted p-value (q-value) below 0.05 was considered significant. To avoid false negatives in hotspots of known UM genes, these regions (GNAQ p.48, p.183, and p.209; GNA11 p.183 and p.209; SF3B1 codons p.625, p.666 and p.700; EIF1AX codons p.1-20, PLCB4 p630; and CYSLTR2 p.129) were called with higher sensitivity. For each sample and genomic position the variant and reference read counts were compared with the variant and reference read counts in the pool of all 103 normal/germline samples at that specific position. Fisher's exact test was used to identify somatic mutations and a Bonferroni corrected p-value below 0.001 was considered statically significant.

Methods
Classification of TP53 mutations. To evaluate the TP53 mutations, they were compared with two comprehensive characterisation studies. The Phenotypic Annotation of TP53 Mutations (PHANTM) score v1.0 is a weighted sum of zscores for which common (i.e. benign) germline variants have values around 0 and recurrent somatic hotspot mutations have scores around 1 35 . The relative fitness score (RFS) is on average −2.50 for synonymous variants, while the average score for protein truncating variants is 0.42 36 . FATHMM is an in silico tool predicting the probability that variants are deleterious 34 . Scores above 0.5 are considered deleterious. The IARC TP53 mutation database (R20, July 2019) was used to evaluate how frequent mutations are 57 . Higher risk categories have shorter RFS. UM patients with monosomy of chromosome 3 has significantly shorter RFS than patients with disomy 3 (blue). Difference between curves were assessed using a two-sided log-rank (Mantel-Cox) test; p-values were not adjusted for multiple testing.
Copy number aberrations and SV. Copy number aberrations were identified using ascatNgs 58 . Copy number of at least 6 was considered high gain. The underlying model of ascatNgs assumes the data come from two clones of cells: the tumour and normal contamination. For a heterogeneous tumour it may therefore overestimate copy numbers; to distinguish heterogeneity from WGD, the distribution of VAF for somatic mutations within regions with allelic balance were used. For a tetraploid tumour, two peaks in VAF are expected corresponding to mutations occurring before and after the copy number gain, with the latter having half the VAF of the former. Statistical models for the data coming from a homogeneous tetraploid tumour and from a heterogeneous diploid tumour, respectively, were inferred using maximum-likelihood and if the heterogeneity was significantly more likely (P < 0.001), copy numbers were adjusted. Structural variants were identified using an in-house tool, qSV, as previously described 53 . Gene truncating breakpoints and consequence of the SVs were determined using in-house scripts and transcript annotation from Ensembl.
Survival analysis. Relapse-free curves were estimated using the Kaplan-Meier method. Difference between curves was assessed with the log-rank (Mantel-Cox) test.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The BAM files are deposited in the European Genome-phenome Archive ([https://www. ebi.ac.uk/ega/]) with accession number EGAS00001001552. The source data underlying Fig. 2 are provided as a data source file. All other data are available in the Article, Supplementary Information or available from the author upon reasonable request.