Multi-Modality Analysis Improves Survival Prediction in Enucleated Uveal Melanoma Patients

Abstract

Purpose: Uveal melanoma (UM) is characterized by multiple chromosomal rearrangements and recurrent mutated genes. The aim of this study was to investigate if copy number variations (CNV) alone and in combination with other genetic and clinico-histopathological variables can be used to stratify for disease-free survival (DFS) in enucleated patients with UM.

Methods: We analyzed single nucleotide polymorphisms (SNP) array data of primary tumors and other clinical variables of 214 UM patients from the Rotterdam Ocular Melanoma Study (ROMS) cohort. Nonweighted hierarchical clustering of SNP array data was used to identify molecular subclasses with distinct CNV patterns. The subclasses associate with mutational status of BAP1, SF3B1, or EIF1AX. Cox proportional hazard models were then used to study the predictive performance of SNP array cluster-, mutation-, and clinico-histopathological data, and their combination for study endpoint risk.

Results: Five clusters with distinct CNV patterns and concomitant mutations in BAP1, SF3B1, or EIF1AX were identified. The sample's cluster allocation contributed significantly to mutational status of samples in predicting the incidence of metastasis during a median of 45.6 (interquartile range [IQR]: 24.7–81.8) months of follow-up (P < 0.05) and vice versa. Furthermore, incorporating all data sources in one model yielded a 0.797 C-score during 100 months of follow-up.

Conclusions: UM has distinct CNV patterns that correspond to different mutated driver genes. Incorporating clinico-histopathological, cluster and mutation data in the analysis results in good performance for UM-related DFS prediction.

Uveal melanoma (UM) is the most common primary intraocular malignancy in adults, characterized by nonrandom chromosomal aberrations and a small set of recurrent mutations in the BAP1, SF3B1, and EIF1AX (BSE) genes. Prognostication has been performed initially only using clinico-histopathological parameters¹ and has gradually become more precise by using other genetic prognosticators, such as single chromosome copy number variation (CNV) status,^2,3 gene expression profiling,⁴ and mutation status,^5,6 and using them in different combinations.^7–11 In this report, we have used the 214 single nucleotide polymorphisms (SNP)-array set of our previous published study¹² and this dataset was expanded with in part unpublished clinico-histopathologic parameters and gene mutation data. In the 80 UM The Cancer Genome Atlas (TCGA) landmark study,¹³ four main molecular subtypes were identified using multiple genomic data sources including CNV. In a recent study,¹⁴ 658 UM fine-needle aspiration biopsies (FNABs) were categorized based on the TCGA classification and tumor clinical features and outcome was analyzed. Our findings correlate well with the aforementioned study results; however, due to our longer follow-up (up to 273 months in our cohort compared with approximately 60 months in the other studies) we are able to detect late-onset metastasis. In addition, using unsupervised CNV clustering, we were able to identify an extra fifth cluster of UM that shows a prolonged disease-free survival (DFS). Furthermore, the predictive value of BSE status compared with complete genomic CNV (interpreted as a cluster variable from 1 to 5) status has not yet been evaluated before and there is no literature that incorporates BSE mutation status with CNV and clinico-histopathological data in a Cox proportional hazard model. Therefore, we set out to incorporate clinico-histopathologic, CNV, and BSE mutation status variables that are known for their predictive value from literature and evaluate whether the combination of these data modalities yields improved accuracy for DFS prediction in our UM cohort with a long follow-up spanning up to 273 months. 

Materials and Methods

Study Population

Patients were selected from the Rotterdam Ocular Melanoma Study Group (ROMS) database, which contained 808 records based on the availability of whole-genome SNP array analysis. In total, 214 patients were included in SNP array analysis. However, patients with iris melanoma (n = 3) and one patient with iris melanocytoma (n = 1) were excluded from Cox proportional hazard analyses. We obtained tumor material from two (n = 2) UM patients from other hospitals (Groningen UMC and Nijmegen Radboud UMC) without clinical data, who were also excluded. All patients underwent primary enucleation (n = 212) or received a biopsy before radiation therapy (n = 2) between 1989 and 2015. An overview of patient and tumor characteristics is provided in Table 1 and mutation status of the analyzed UM is provided in Supplementary Table S1. Patients' follow-up data were updated until August 2017. Informed consent was obtained from all patients. This study was approved by the local ethics committee and performed according to the guidelines of the Declaration of Helsinki. 

Table 1

 

Clinical Data of the Patients and Histopathological Data of the UM

 

Histopathology

A histopathologic diagnosis of melanoma was confirmed by an experienced ophthalmic pathologist (R.M.V.) conforming to the Royal College of Pathologists guidelines (available at: https://www.rcpath.org/profession/guidelines/cancer-datasets-and-tissue-pathways.html). We used tumor thickness and largest basal diameter measurements during histological preparation. The T class of tumors (hereafter TNM) was determined using the American Joint Committee on Cancer eighth edition staging system.¹⁵ In short, hematoxylin and eosin staining was used to differentiate between epithelioid or spindle cells according to the modified Callender classification. Epithelioid cells were defined as present when more than 10% of tumor cells exhibited an epithelioid phenotype. Lymphocytic infiltration was defined as any obvious clusters of lymphoid inflammatory cells in the tumor. The mitotic rate was determined by counting the mitosis in 8 mm² equal to 50 high power fields. Microfoci of necrosis were accepted as positive. Extraocular extension was defined as tumor growth through the sclera and beyond the outer scleral surface. Extracellular matrix networks were visualized in tumor specimens stained with periodic acid-Schiff staining without hematoxylin staining and evaluated using a green filter, defined as at least three back-to-back closed loops. 

Immunohistochemistry (IHC) was performed with an automated IHC staining system (Ventana BenchMark ULTRA; Ventana Medical Systems, Tucson, AZ, USA) using the alkaline phosphatase method and a red chromogen. In brief, following deparaffinization and heat-induced antigen retrieval for 64 minutes, the tissue sections were incubated with a mouse monoclonal antibody raised against amino acids 430 to 729 of human BAP1 (clone sc-28383, 1:50 dilution; Santa Cruz Biotechnology, Dallas, TX, USA) for 1 hour at 36.1°C. A subsequent amplification step was followed by incubation with hematoxylin II counterstain for 8 minutes and then a blue-coloring reagent for 8 minutes according to the manufacturer's instructions (Ventana). Liver, tonsil, breast tissue, and the retinal pigment epithelium were used as positive controls for BAP1 expression. Loss of expression was defined as absent BAP1 expression in the tumor cell nucleus. Only nuclear expression was scored because only nuclear expression has proven to be of prognostic relevance in UM.^16–18 

Mutational Status

After enucleation or biopsy, the tumor was divided in three parts and snap-frozen in liquid nitrogen, formalin-fixed paraffin-embedded (FFPE), or directly used for DNA extraction. High-quality DNA was isolated from fresh tumor samples or from liquid nitrogen–stored tumor samples using the QIAmp DNA-mini kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. Mutation analysis on FFPE material was of sufficient quality to perform targeted sequencing using the Ion Torrent Personal Genome Machine system (Life Technologies, Carlsbad, CA, USA). 

Mutation analyses of GNAQ, GNA11, BAP1, SF3B1, and EIF1AX were performed using Sanger sequencing and targeted sequencing on the Ion Torrent Personal Genome Machine as described before.¹⁹ In 202 samples, BAP1 status was analyzed using BAP1 IHC, and in 118 samples BAP1 mutation analysis (44 extra samples and the 74 samples mentioned in the study of Koopmans et al.¹⁶). Sanger sequencing was used for validation. In 12 cases in which BAP1 IHC was positive and the BAP1 mutation was a missense or in-frame deletion, we used a cutoff of >20 for the combined annotation dependent depletion (CADD)²⁰ scores in order to reassign BAP1 IHC-positive samples (i.e., wild-type samples) to the BAP1^IHCneg/mut group in further analyses. In 202 samples, we sequenced exon 14 of SF3B1 to detect mutations in the R625 or K666 residues and in more than 150 samples, SF3B1 exons 12 to 16 (including the K700) were sequenced. EIF1AX exons 1 and 2 with flanking regions up to 25 base pairs into the 5′ untranslated region and intron 1 of EIF1AX (n = 200) were analyzed. We classified samples as having no recurrent mutation (NRM) if they had positive BAP1 immunostaining and had no mutations in BAP1, SF3B1, or EIF1AX. 

SNP Array Analyses

A total of 200 ng DNA input was used for whole-genome copy number analyses by SNP arrays. In all samples (n = 214), more than 90% of the SNPs were mapped, which were used for further analyses. For the entire set, four types of Illumina Human SNP array platforms were used: 18 UM were analyzed using the CytoSNP-12 v2.0 BeadChip; 58 UM were analyzed using the CytoSNP-12 v2.1; 46 UM were analyzed using OmniExpress-12 v1, and 92 UM were analyzed using CytoSNP-850K BeadChips (Bead Studio Illumina, San Diego, CA, USA) (Supplementary Fig. S1). Nexus Copy Number 8.0 (BioDiscovery, Inc., El Segundo, CA, USA) was used to calculate the total percentage aneuploidy of the entire genome, the total Copy Number (CN) events and for visualization of the whole-genome SNP array data. For these calculations, CN events smaller than 1 Mb were filtered out. Chromosomal aberrations of the sex chromosomes also were left out of the analysis to prevent gender-based bias. Cluster-wise comparisons were performed using the Mann-Whitney-Wilcoxon test with Holms correction for multiple testing. False discovery rate (FDR) correction was used for multiple testing, which was called significant when FDR ≤0.05 and P value ≤0.05 was called significant for single test comparisons.²¹ 

Nonweighted Hierarchical Clustering

A new dataset was obtained by combining samples across the four array comparative genomic hybridization (CGH) platforms. The dataset includes measurements of 215,138 SNPs that overlapped across the platforms on 214 samples (Supplementary Fig. S1). To remove germline-specific aberrations, tumor profiles in each platform underwent correction using the data from 10 genomic DNA samples isolated from blood of non-UM patients (R-package NoWaves version 0.6).²² Segmentation analysis was performed on the corrected data, which subsequently was used to calculate the called data using R-package CGHcall version 2.36.0.²³ The segmented and called data were combined to generate the regioned data, which has far fewer features than the former two and equal weights were applied for all regions. Finally, nonweighted hierarchical clustering analysis was performed using the soft calls (regioned call probability) data with symmetric Kullback-Leibler divergence as distance measure and Ward Linkage.²⁴ The analyses were performed using the dedicated R-package WECCA version 0.40 on the regioned data.²⁵ All analyses have been performed using the R statistical environment, version 3.3.1. 

Statistical Analysis

DFS was determined from the date of enucleation until the date of metastasis or until last follow-up. Patients were censored when they were lost to follow-up or when death from a cause other than UM occurred. The Kaplan-Meier method was applied to estimate the cumulative incidence of the study endpoint during follow-up, whereas differences between the patient groups were compared using the log-rank test. Up to 5% of data were missing for 24 variables. Taking this small percentage into account, we applied the simple method of mean imputation to complete the dataset for analysis. The relations between a broad range of clinico-histopathological, SNP array cluster, and mutation data, and the incidence of the study endpoint were evaluated by univariable Cox proportional hazard (PH) regression analyses. We report the corresponding unadjusted hazard ratios (HR) and 95% confidence intervals (CI). Inspection of the Schoenfield residuals showed no meaningful deviation of the PH assumption for all models. Subsequently, multivariable Cox PH models were fitted to study the predictive value of the combined clinico-histopathological data, and the added value of SNP array cluster and mutation data, which was evaluated by likelihood ratio tests. All variables included in the multivariable model reported in Table 2 were analyzed as categorical data. We report C-indices (indicating discriminating performance) for all models, and a calibration curve for the most extended model. 

Table 2

 

Multivariable Regression Models Combining Data From Three Different Domains for the Prediction of the Study Endpoint

 

The R statistical software package (version 3.5.0) was used for data analyses. Statistical significance was set at the 0.05 probability level. 

Results

Nonweighted Hierarchical Clustering Using SNP Data

The heat map obtained after hierarchical cluster (HC) analysis revealed five clusters, which are mainly driven by CNVs on chromosomes 1, 3, 6, and 8 (Fig. 1A). In almost all cases, chromosome 3 (disomy or monosomy) divides the three clusters (A, B, C) on the left major clustering branch from the right major clustering branch (D, E) (Supplementary Fig. S2A). 

Figure 1

 

Nonweighted hierarchical clustering of regional CNV and distinct distribution of UM mutational status. (A) Heat map of the CNVs constructed with unsupervised nonweighted hierarchical clustering. (B) Doughnut charts for every cluster showing the distribution of the BAP1, SF3B1, and EIF1AX mutation status, including several variables. (C) Chromosomal patterns for the five clusters constructed by nonweighted hierarchical clustering.

The clusters A, B, and C in the left major clustering branch of the HC tree are in general characterized by either a paucity of CNVs or specific CNVs in chromosome 6 (Figs. 1A, 1C). Samples in cluster A have no CNVs in chromosome 6q, whereas cluster B was characterized by entire chromosome 6 gain and cluster C is composed of UM with chromosome 6p gain and 6q loss in the tumor. Supplementary Fig. S2B and S2C show high discriminatory values for chromosome 6q when comparing clusters A versus B-C and clusters B versus C. The CNVs in chromosome 6q were also able to discriminate between cluster D and E (Fig. 1C, Supplementary Fig. S2D). 

Clusters A (n = 58) and B (n = 11) contained samples with only small CNVs or no CNVs at all. Typical for cluster A are chromosome 6p (58%) and 8q (28%) gain. Cluster B contains UM with gain of chromosome 6 (91%), 8p (36%), and 8q (45%) (Fig. 1A, 1B). In contrast to the clusters A and B, cluster C harbors UM tumors with multiple CNVs. Gains in chromosomes 6p (89%) and 8q (65%), and loss in chromosome 6q (91%) are observed in most cases. Moreover, these CNVs are in general larger than those in clusters A and B. Also, gain of the proximal part of chromosome 11q and loss of the terminal end of chromosome 11q are observed recurrently (26% and 43%, respectively) (Fig. 1C). 

The two clusters D (n = 88) and E (n = 34) of the right major clustering branch have losses in chromosomes 1p (31%; 50%), 3 (91%; 88%), and 8p (27%; 23%), and gain of chromosome 8q (78%; 88%) in common (Fig. 1C). Loss of 4q (20%), 6q (76%) and gain of chromosomes 1q (29%), 4p (23%), 6p (35%), and 8p (20%) are characteristic of cluster E and the entire right branch of cluster E consisted of UM cases with polyploidy. 

Next, aneuploidy and total CN events were calculated for the two major clustering branches and for each cluster separately. The median percentage of genomic aneuploidy in the left major clustering branch was 5.2% (interquartile range [IQR]: 1.0%–8.6%) compared with a median of 12.7% (IQR: 10.3%–18.9%) for the right major clustering branch. Also the median number of CN events in the left major clustering branch (seven events; IQR: 2.5–13.0) was significantly (P << 0.001) lower compared with the CN event number of right major clustering branch (12 events; IQR: 7.0–20.0) (Figs. 2A, 2B). 

Figure 2

 

Boxplots showing percentage aneuploidy and CN events. (A) Statistical comparison of the percentage of aneuploidy between samples allocated in the left major clustering branch versus the right major clustering branch, and (B) statistical comparison of the total number of CN events of samples in the left major clustering branch versus the right major clustering branch in the ROMS cohort. (C) The same analyses performed on samples allocated in the clusters for percentage of aneuploidy and (D) total number of CN events. Mann-Whitney-Wilcoxon test was used to compare groups with Holms correction for multiple testing. The continuous lines between the groups represent statistically significant differences (P < 0.05) between the clusters. The dashed line depicts no significant difference between the clusters.

Similar analyses were performed for the five clusters. The samples in the cluster A have the lowest median of 1.6% (IQR: 0.1%–5.2%) aneuploidy of the genome consistent with findings of Robertson et al.,¹³ whereas cluster E has the highest median aneuploidy of 23.7% (IQR: 16.1%–30.4%) (Fig. 2C). The absolute number of CN events was the lowest in cluster A with a median of four changes per sample (IQR: 1.0–8.0 events), whereas a median of 22.5 (IQR: 11.0–48.0 events) changes was observed in cluster E (Fig. 2D). 

Association of the Clusters With Mutation Status

In cluster A we observed UM cases with SF3B1 mutations (n = 13; SF3B1^mut), EIF1AX mutations (n = 18; EIF1AX^mut), SF3B1 and EIF1AX^mut UM samples (n = 2), or no BSE mutation and BAP1IHC positive (n = 18; NRM) (Supplementary Table S1). Six UM samples in cluster A were BAP1 IHC negative or contained a mutation in BAP1 (BAP1^IHCneg/mut) (Fig. 1B). Cluster B contained five EIF1AX^mut, one UM that is BAP1 IHC negative and harbors an EIF1AX^mut and four NRM UM, respectively. Cluster C predominantly harbors SF3B1^mut UM (n = 16; 70%) and in cluster D, BAP1^IHCneg/mut samples were enriched (n = 67; 76%). Enrichment of BAP1^IHCneg/mut (n = 24; 71%) was observed in cluster E. This distribution of mutational status across clusters is consistent with the distribution of UM samples found in the study of Robertson et al.¹³ except our analysis revealed an additional cluster. The NRM samples (n = 5; 15%) in cluster E were BAP1 IHC-positive polyploid UM cases. 

Association With Clinical and Tumor Data

Clinical and pathological data from 214 patients were revised and updated. The cohort consists of 107 women (50%) and 107 men (50%) with a mean age at diagnosis of 62.6 years (95% CI: 60.7–64.5). These cluster-wise and branch-wise clinical data are shown in Table 1. Mean age at diagnosis was higher in the right major clustering branch compared with the left (64.6 vs. 59.9 years, P = 0.016), respectively, whereas mean longest tumor diameters (LTDs) of UM in the left major clustering branch were smaller than LTD of UM in the right major clustering branch (12.1 vs. 13.8 mm, P = 0.001). 

Cluster-wise Kaplan-Meier survival curves are shown in Figure 3, and we used radiological-proven metastatic disease as study endpoint. Samples in the cluster B, mostly EIF1AX^mut and NRM, showed no events and have a median DFS of 72.2 (IQR: 55.9–175.0) months (Table 1). 

Figure 3

 

Survival curves of all patients included. Kaplan-Meier survival curves for each distinct CNV cluster of the ROMS with primary endpoint DFS (P < 0.001). Multiple black line types were used: “two dash” is used for cluster A, “long dash” for cluster B, “solid” for cluster C, “dotted” for cluster D, and “dot dash” for cluster E.

The clusters A and C correspond to intermediate survival (Cluster A: median DFS: 59.5 months, IQR: 43.7–111.9, Cluster C: median DFS: 49.1 months, IQR: 33.4–106.8), whereas the clusters D and E combined (right major clustering branch) correspond to worst survival with a median DFS of 35.6 (IQR: 20.4–60.7) months (Table 1). 

Prediction of DFS

Several factors in clinico-histopathological domain were associated with an increased incidence of the metastatic disease, including the presence of epithelioid cells, presence of closed extracellular matrix patterns, TNM (stages 3 and 4), longest tumor diameter, and lymphocytic infiltrate (Supplementary Table S2). With respect to SNP array cluster data: the incidence of the study endpoint was significantly higher in patients allocated in cluster D (HR: 10.05) or E (HR: 9.15) as compared with the reference cluster A. Finally, patients with a BAP1^IHCneg/mut tumor had a significantly increased risk (HR: 9.33) compared with those tumors that were BAP1 IHC-positive or had no mutation in BAP1, whereas SF3B1 mutated UM (HR: 0.48) patients and those UM with an EIF1AX mutation (HR: 0.07) or wild-type UM (HR: 0.19) had significantly reduced risk of the study endpoint. 

The regression model that included clinico-histopathological data yielded a C-index of only 0.728, indicating reasonable discriminative performance (Table 2). LTD and lymphocytic infiltrate lost significance in multivariate analysis and were not included in the regression models. The clinico-histopathological model was significantly improved by adding either SNP array cluster data (C-index 0.775, P < 0.001) or mutation data (C-index 0.794, P < 0.001). The SNP array cluster data model only was significantly improved by adding mutation data (P < 0.001), this was also the case when mutation data was added to the SNP array data model (P < 0.05). The model combining data from all three domains showed good discrimination (C-index 0.797), whereas the incidence of the study endpoint (at 100 months) was adequately predicted (Fig. 4). 

Figure 4

 

Calibration curve comparing the predicted probability of metastatic disease at 100 months follow-up with observed frequencies. Predicted probabilities, which are based on the most extended Cox model, are plotted against observed frequencies. The straight line describes the (optimal) scenario where predicted probabilities are equal to the observed frequencies.

Discussion

In this study, we observed that dividing our UM cohort in five clusters with distinct overall CNV patterns improved our proportional Cox regression analysis model in comparison with solely using clinico-pathological and/or mutation data. 

The added value of distinct CNVs of chromosomes 1, 3, 6, and 8 has been recognized and used before by other researchers.^8,10,26 However, using a cutoff of five clusters in the HC and our relatively long follow-up we defined in contrary to other studies, a sub cluster B in the disomy 3 UM with a typical gain of entire chromosome 6 and a prolonged DFS. After evaluation of the sample distribution in cluster B, we observed whole chromosome 6 gain almost exclusively in EIF1AX^mut tumors or UM without a recurrent mutation, except for two UM in cluster E that displayed a polyploid genome. Therefore, we suggest whole chromosome 6 gain without polyploidy is predictive of EIF1AX-mutated tumors or UM without a recurrent mutation. 

We validated our CNV clustering data results with a recently published study by Robertson et al.¹³ that performed CNV, hierarchical clustering analysis on UM data provided by TCGA.¹³ We observed major similarities in sample allocation (i.e., chromosome 3 loss samples were allocated in cluster 3 and 4 just as UM in our study were allocated in clusters D and E). Moreover, the clusters 3 and 4 contain nearly all BAP1 IHC negative/mutated samples. When observing the differences in survival between cluster 3 and 4 in the TCGA based on CNV cluster 4 UM develop metastasis earlier than cluster 3 UM. We also observed this phenomenon using a cutoff of 25 months (P < 0.05) of follow-up in our set. Therefore, cluster 4 from the TCGA is comparable to our cluster E, which is mainly characterized by chromosome 3 loss and chromosome 6 loss and 8q gain. The main difference between our cluster E and cluster 4 of the TCGA study is cluster 4 of the TCGA study consists of isochromosome 8q UM, which are distributed throughout our clusters D and E. A feature of our clustering algorithm on the other hand is the ability to detect polypoid samples that are located in our cluster E. However, if we take into account the follow-up of the TCGA cohort (60 months) and our maximum follow-up (273 months), clusters D and E from our set eventually show similar survival curves, suggesting no prognostic value in subdividing monosomy 3 UM based on CNV using follow-up longer than 273 months. Samples with partial chromosome 6p gain, partial chromosome 6q loss, and disomy, 3 CNV samples were combined in cluster 2 of the TCGA UM, which is comparable to cluster C from our study and contains most of the SF3B1-mutated UM. This cluster 2 UM appears to have a favorable prognosis in the TCGA study. However, UM in cluster C eventually shows late-onset metastasis with longer follow-up in our set. As an addition to the TCGA analysis, our extra fifth cluster over the TCGA set is cluster B, which is characterized by chromosome 6 gain, which is probably detected due to our relatively large SNP array sample size (n = 214) and these UMs show prolonged survival. 

The unsupervised nonweighted hierarchical SNP array clustering results corroborated our previous study¹² whereby we performed supervised clustering with the known mutational status of UM and combined both SNP array and cytogenetic data. Hence, whole chromosomal arm involvement in clusters D and E is explained by isochromosome formation in BAP1^IHCneg/mut samples, whereas cluster C was enriched for SF3B1-mutated samples with multiple distal and often smaller chromosomal structural variations (Fig. 1A). These findings were observed in both cytogenetic and in SNP array data.¹² Also, in other malignancies, the association of BAP1 mutations with entire chromosome arm CNVs are observed. For example, in hepatocellular carcinoma (HCC) analyses performed by TCGA Network, BAP1-mutated samples showed loss or gain of entire chromosomal arms suggestive of isochromosomes in 14 of 23 BAP1-mutated HCC samples, implying a similar underlying genetic mechanism.²⁷ 

We observed that the clusters are indicative of the samples' mutation status. BAP1^IHCneg/mut UM mainly clustered based on monosomy 3. UM with EIF1AX mutations or no detectable recurrent mutations tightly clustered together by either a very limited number of CNVs or the presence of chromosome 6 gain. SF3B1^mut UM cluster together based on recurring smaller CNVs in chromosome 6p, 6q, and 8q. These analyses indicate that the UM tumor cells, which are BAP1^IHCneg/mut, SF3B1^mut, or EIF1AX^mut, do have distinct CNV and chromosomal patterns. 

We used predominantly BAP1 IHC to define the BAP1 mutation status. This should not make a major difference, as we and others have shown a strong association between the loss of nuclear BAP1 protein expression and the presence of BAP1 mutations.^16,17 It might even be an advantage since BAP1 mutations are not always detected or evident from the sequencing data.²⁸ This could explain the fact that in the TCGA¹³ SNP array dataset all BAP1 mutated samples determined by sequencing were allocated in the chromosome 3 loss clusters, whereas this was not the case for nine samples in our set. Two of the nine BAP1^IHCneg/mut UM patients died of metastatic disease, both showing normal chromosome 3 status on SNP array (DFS of 34 and 33.4 months). However, one sample showed a small loss of heterozygosity of the BAP1 region and the other showed 21% FISH count signal of one copy of chromosome 3p, which was not consistent with SNP array results, possibly due to tumor heterogeneity or contamination of non-UM tissue in these particular cases. A third UM sample harbored the classical monosomy 3 CNV, which occurred in combination with gain of whole chromosome 6. In this particular UM patient, gain or overexpression of genes located on chromosome 6 such as the tumor suppressor genes PLAG1 or LATS1²⁹ may counteract the loss of BAP1 protein, resulting in a DFS of more than 266 months. None of the six other UM BAP1^IHCneg/mut patients with a DFS of 25.3, 43.9, 52, 67, 81, and 95 months, respectively, developed metastatic disease, suggesting loss of BAP1 protein expression or mutation of one allelic copy of BAP1 without concomitant loss of chromosome 3 appears not to correlate with increased metastatic risk. 

Our proportional Cox regression analysis model including clinico-pathological, cluster, and mutation data show significantly higher C-scores compared with independent analysis per data type. We also have shown there is additional prognostic value combining CNV data with mutation status data. This result implicates both SNP array and mutational status data should be acquired to provide a more accurate prognosis for UM patients when clinico-histopathological data are not available. Furthermore, our full model yields a C-score of 0.797, indicating a strong model and is comparable to the study of Eleuteri et al.³⁰ with a C-score of 0.79 and is higher than the validation study of this method on a cohort at the University of California San Francisco with a C-index of 0.67.³¹ However, by adding chromosome 8q CNV status in the model of Eleuteri et al.,⁸ the group successfully managed to improve the C-score to 0.861. 

One limitation of our and the TCGA study is that only tumors from enucleated eyes were used and therefore a bias toward larger tumors is expected. However, Vichitvejpaisal et al.¹⁴ published a study in which they applied the TCGA classification on FNABs and were able to corroborate the TCGA results in these, in general, smaller tumors. 

Nevertheless, we do not expect our CNV clustering results to be affected due to this size bias because chromosome 3 loss was a major discriminatory variable and approximately 90% of our BAP1^IHCneg/mut UM were allocated in the right major clustering branch. An overrepresentation of larger tumors most likely has a reinforcing significant effect on Kaplan-Meier curve due to high mortality in the right major clustering branch samples compared with samples in the left major clustering branch. In addition, our sample set of 214 is relatively small compared with the Eleuteri et al.⁸ cohort (n = 4161). The small sample size could be an explanation that the variable “extraocular extensions” is not significant in our Cox model. However, the 4161 samples of Eleuteri et al.⁸ included genetic data of 602 samples (i.e., chromosome 3 loss status and chromosome 8q gain), which is only 3-fold higher than our cohort, indicating we report a reasonable sample size. 

In conclusion, we show that patients with UM can be divided into molecular subclasses based on genetic CNV patterns, and above all, these groups correspond very well to the BAP1, SF3B1, and EIF1AX mutational status. Each group is characterized by recurring CNVs reflecting the different types of chromosomal aberrations (e.g., whole chromosome arm loss or gain as seen in the BAP1^IHCneg/mut tumors versus the smaller structural anomalies of chromosome 8 and 6 as seen in the SF3B1^mut subtype), indicating that different pathways are involved in the etiology of these UM subtypes and possibly also progression toward metastatic disease. However, mutation analysis or CNV detection using SNP arrays are not 100% sensitive, and combining these resulted in additional prognostic value in our proportional Cox regression analysis model. Last, BAP1, SF3B1, and EIF1AX have diverse cellular functions, and how exactly the changed proteins contribute toward tumorigenesis is not clear. Thus, by studying each subgroup separately and combining expression, protein, and epigenetic data, we will be able to pinpoint the essential downstream targets for tumor development in UM. Furthermore, because some of the driver genes are also found in more common malignancies (e.g., SF3B1 mutations in breast cancer or leukemia), targeted therapeutic options such as the use of spliceosome inhibitors might also be applicable for SF3B1^mut UM. 

Acknowledgments

The authors thank all patients for their contribution of sample material for this research. 

Supported by the Prof. dr. Henkes Foundation (AO, EK), the Combined Ophthalmic Research Rotterdam (CORR) (AdK, EK, DP), Stichting Coolsingel (WD, EK) and The Netherlands and the Dutch Cancer Society (EK). 

Disclosure: W. Drabarek, None; S. Yavuzyigitoglu, None; A. Obulkasim, None; J. van Riet, None; K.N. Smit, None; N.M. van Poppelen, None; J. Vaarwater, None; T. Brands, None; B. Eussen, None; R.M. Verdijk, None; N.C. Naus, None; H.W. Mensink, None; D. Paridaens, None; E. Boersma, None; H.J.G. van de Werken, None; E. Kilic, None; A. de Klein, None 

Appendix

Members of the Rotterdam Ocular Melanoma Study Group

Wojtek Drabarek,^1,2 Serdar Yavuzyigitoglu,^1,2 Askar Obulkasim,^1,3 Kyra N. Smit,^1,2 Natasha M. van Poppelen,^1,2 Jolanda Vaarwater,^1,2 Tom Brands,^1,2 Bert Eussen,² Robert M. Verdijk,^4–6 Nicole C. Naus,¹ Hanneke W. Mensink,⁵ Dion Paridaens,^1,5 Emine Kilic,¹ Annelies de Klein,² Ronald O. B. de Keizer,⁵ Jackelien van Beek,⁷ Caroline M. van Rij,⁸ Erwin Brosens,² Jolique A. van Ipenburg,⁴ and Daniel P. de Bruyn^1,2 

¹Department of Ophthalmology, Erasmus University Medical Center, Rotterdam, The Netherlands 

²Department of Clinical Genetics, Erasmus University Medical Center, Rotterdam, The Netherlands 

³Department of Pediatric Oncology/Hematology, Erasmus University Medical Center, Rotterdam, The Netherlands 

⁴Department of Pathology, section Ophthalmic Pathology, Erasmus MC Cancer Institute, Erasmus University Medical Center, Rotterdam, The Netherlands 

⁵The Rotterdam Eye Hospital, Rotterdam, The Netherlands 

⁶Department of Pathology, Leiden University Medical Center, Leiden, The Netherlands 

⁷Department of Ophthalmology, Albert Schweitzer Hospital, Dordrecht, The Netherlands 

⁸Department of Radiation Oncology, Erasmus Medical Center, Rotterdam, The Netherlands