# Heterogeneity revealed by integrated genomic analysis uncovers a molecular switch in malignant uveal melanoma

Oncotarget. 2015; 6:37824-37835. https://doi.org/10.18632/oncotarget.5637

Metrics: PDF 1472 views  |   HTML 1938 views  |   ?

## Abstract

Mark J. de Lange1, Sake I. van Pelt1, Mieke Versluis1, Ekaterina S. Jordanova2, Wilma G.M. Kroes3, Claudia Ruivenkamp3, Sjoerd H. van der Burg4, Grégorius P.M. Luyten1, Thorbald van Hall4, Martine J. Jager1, Pieter A. van der Velden1

1Department of Ophthalmology, LUMC, Leiden, The Netherlands

2Department of Pathology, LUMC, Leiden, The Netherlands

3Department of Clinical Genetics Laboratory for Diagnostic Genome Analysis (LDGA), LUMC, Leiden, The Netherlands

4Department of Clinical Oncology, LUMC, Leiden, The Netherlands

Correspondence to:

Pieter A. van der Velden, e-mail: velden@lumc.nl

Keywords: uveal melanoma, tumor heterogeneity, digital PCR, immune profile, molecular etiology

Received: April 15, 2015     Accepted: September 28, 2015     Published: October 08, 2015

ABSTRACT

Gene expression profiles as well as genomic imbalances are correlated with disease progression in uveal melanoma (UM). We integrated expression and genomic profiles to obtain insight into the oncogenic mechanisms in development and progression of UM. We used tumor tissue from 64 enucleated eyes of UM patients for profiling. Mutations and genomic imbalances were quantified with digital PCR to study tumor heterogeneity and molecular pathogenesis. Gene expression analysis divided the UM panel into three classes. Class I presented tumors with a good prognosis and a distinct genomic make up that is characterized by 6p gain. The UM with a bad prognosis were subdivided into class IIa and class IIb. These classes presented similar survival risks but could be distinguished by tumor heterogeneity. Class IIa presented homogeneous tumors while class IIb tumors, on average, contained 30% of non-mutant cells. Tumor heterogeneity coincided with expression of a set of immune genes revealing an extensive immune infiltrate in class IIb tumors. Molecularly, class IIa and IIb presented the same genomic configuration and could only be distinguished by 8q copy number. Moreover, UM establish in the void of the immune privileged eye indicating that in IIb tumors the infiltrate is attracted by the UM. Combined our data show that chromosome 8q contains the locus that causes the immune phentotype of UM. UM thereby provides an unique opportunity to study immune attraction by tumors.

## INTRODUCTION

Uveal melanoma (UM) is a rare ocular neoplasm characterized by GNAQ and GNA11 mutations [1, 2]. Despite changes in treatment, the overall survival rate remains low and patients still die due to metastases which are usually found after the initial diagnosis [3, 4].

For years, monosomy of chromosome 3 has been the most described biomarker that predicts survival in UM patients [5]. Markers on chromosome 3 that underlie this unique subdivision are largely unknown but the chromosome 3 status divides UM genetically into two groups which have a good and a bad prognosis [6]. This finding is supported by clinical and histopathological markers [7]: UM containing epithelioid cells are associated with monosomy 3 and a bad prognosis while UM purely made up of spindle cells are rarely lethal [8, 9]. Larger tumors are correlated with a worse prognosis and tend to have an epithelial/mixed phenotype rather than a spindle cell phenotype. With SNP [1013], CGH [1417] and karyotype analysis [18, 19], recurrent chromosomal aberrations were further investigated in UM and shown to be correlated with disease progression. Gain of chromosome 8q and 6p is frequently detected while loss of 1p and 16q loss is less common but is still found to be correlated with prognosis [15, 20, 21]. Although these tumor characteristics can be applied to precisely predict disease outcome, we do not yet understand the sequence of genetic events in the development of UM, nor the relationship between the genetic changes and tumor behavior.

To better understand UM progression and develop targeted therapy, we performed genome-wide gene expression analysis and chromosome analysis of 64 tumors. We integrated expression and structural analysis to obtain insight into the molecular etiology of prognostic groups. We investigated prognostic genes and correlated them to chromosomal aberrations to visualize UM molecularly. Additionally, we validated structural aberrations with digital PCR (dPCR) which revealed the heterogeneity of UM. Cellular heterogeneity provides an insight into the micro environment of UM while detailed molecular heterogeneity is the basis for a molecular progression model, showing that increased copy number of 8q precedes loss of chromosome 3. Moreover, our data show that distinct genetic events give rise to different classes of UM, impacting prognosis and therapy.

## RESULTS

### Expression profiles segregate UM into class I, IIa and IIb

Unsupervised cluster analysis of gene expression data of ~16,000 unique genes in 64 primary UM subdivided the tumors into two distinct classes [22]. Class I tumors presented a good prognosis while class II tumors were correlated with a bad prognosis similar to what has been shown previously (Supplementary Figure S1) [6]. Out of the 4,000 differentially expressed genes, 237 genes met the criteria of a Log fold change (LFC) of at least 1.0 and a p-value less than 0.05. A total of 132 genes presented a lower expression in class I than in class II, while the remaining 105 displayed a higher expression in class I. When using these 237 differentially expressed genes for cluster analysis, the UM cohort split into 3 clusters. Class I remained the same as with unsupervised analysis, but class II was subdivided into class IIa and IIb (Figure 1A).

Figure 1: Gene expression analysis. Unsupervised clustering of gene expression of 64 UM divides tumors in two classes. Supervised cluster analysis with the 237 most differentially expressed genes resulted in 3 classes because class II is subdivided in class IIa and class IIb A. Supervised clustering of the 53 class IIb classifier genes B.

Supervised gene expression analysis for IIa and IIb revealed a select panel of genes (n = 53) that was differentially expressed between class IIa and IIb (Supplementary Table S2). Among these genes, only two genes showed a higher expression in class IIa compared to IIb while 51 genes showed a higher expression in IIb (Figure 1B). Survival analysis revealed no significant difference between class IIa and IIb (not shown).

### Chromosomal aberrations are specific for UM classes

In order to investigate which genetic mechanisms underlie the subdivision into 3 classes, we investigated genomic aberrations in UM. With SNP analysis, five recurring chromosomal aberrations were detected which were validated with dPCR (Table 1). Loss of chromosome 1p, loss of chromosome 3, gain of chromosome 6p, gain of chromosome 8q and loss of chromosome 16q were common events. The distribution of the chromosomal aberrations (Figure 2) as well as their copy numbers (Supplementary Table S3) were plotted for the three classes showing that the chromosomal aberrations are non-randomly distributed over the three classes.

Table 1: Copy number analysis

Class

Tumor no.

FA

Chr 1

Chr 3

Chr 6p

Chr 8q

Chr 16q

GNAQ/11

dPCR

SNP

dPCR

SNP

dPCR

SNP

dPCR

SNP

dPCR

SNP

I

06-046*

35.0**

1.7

2.0

2.0

1.9

2.7

2.6

2.7

2.7

1.4

1.3

I

01-074

48.9

1.2

1.4

2.2

1.9

2.9

2.0

3.5

2.7

2.1

2.3

I

04-075

37.1

1.8

1.7

2.0

1.9

2.8

1.9

3.0

2.9

2.0

2.0

I

02-167

45.1

1.3

1.5

1.5

1.6

2.3

2.1

2.1

2.1

2.0

2.0

I

07-007

36.8

1.7

1.9

1.1

1.0

2.9

2.8

1.9

1.9

2.0

1.8

I

07-003

37.9

1.1

1.0

2.2

2.0

2.2

2.0

2.2

1.9

1.9

1.8

I

20-125

38.5

2.0

1.9

2.2

2.0

4.8

2.1

3.2

2.7

1.9

1.9

I

05-020

34.3

1.0

1.2

1.9

1.9

3.5

3.4

4.6

4.5

1.0

1.1

I

02-158

40.2

2.0

2.0

2.3

2.0

4.5

1.6

2.2

2.0

1.5

1.6

I

04-103

44.8

1.8

1.9

1.9

2.0

2.9

1.3

1.9

2.0

1.9

2.1

I

06-033

46.7

0.9

1.1

2.1

2.1

3.0

3.1

3.0

3.2

1.0

1.0

I

99-187

43.2

1.9

1.9

2.2

2.0

2.2

2.0

2.2

2.0

1.1

1.6

I

02-199

42.2

2.0

1.7

2.0

1.8

2.4

2.2

2.0

1.9

2.0

1.9

I

03-086

53.8

1.9

1.9

2.2

2.0

2.3

2.0

2.0

1.9

2.1

2.0

I

06-036

47.4

1.9

2.0

2.2

2.0

3.2

3.1

2.9

3.1

2.1

1.9

I

01-042

50.0**

1.8

1.8

2.1

2.0

2.1

2.1

1.9

1.9

2.0

2.0

I

05-058

49

1.8

2.0

1.9

2.0

2.8

2.8

1.9

1.9

1.8

1.9

I

06-010

41.2

1.7

1.9

2.0

2.0

2.4

2.2

2.0

2.0

1.9

1.8

I

06-011

43

2.2

2.0

2.1

2.0

3.5

2.9

3.5

3.2

1.2

1.2

I

03-087*

44.4

1.7

1.9

2.0

2.1

2.4

2.1

1.9

2.0

1.6

1.8

I

04-074

48.5

1.8

2.2

1.9

2.0

3.0

2.1

2.0

2.0

1.9

2.0

I

07-004

33.6

1.9

1.9

2.1

2.1

3.2

3.2

2.0

2.0

1.8

1.8

I

07-034

46.7

1.7

2.0

2.0

1.9

3.0

2.8

2.1

2.0

1.9

1.9

I

03-120

52

1.9

1.9

2.0

2.0

3.1

2.8

1.9

1.9

1.9

1.8

I

03-129

29.5

1.9

1.9

1.4

1.5

2.0

2.1

3.1

2.7

1.9

1.6

IIa

05-034

36.4

1.6

1.6

1.6

1.5

2.2

2.2

3.1

3.1

1.6

1.5

IIa

06-009

49

1.1

1.1

1.0

1.1

2.5

2.3

3.0

3.2

1.9

2.1

IIa

04-018

55.4

1.9

2.0

2.0

1.9

2.7

2.0

2.7

2.3

1.9

1.9

IIa

07-047

24.3

0.9

1.0

1.1

1.0

2.0

2.0

2.0

2.0

1.9

1.8

IIa

02-174

47.3

0.9

1.3

1.1

1.4

2.1

2.2

2.7

2.6

2.0

2.1

IIa

05-033

45.7

1.4

1.5

1.2

1.1

2.0

1.9

2.4

2.1

2.2

1.9

IIa

06-042

34

1.9

1.9

1.3

1.1

2.2

2.1

2.9

3.0

2.0

1.9

IIa

05-061

45.2

1.0

1.1

1.1

1.0

2.2

2.0

5.8

5.7

2.1

1.9

IIa

05-046

43.1

1.1

1.1

1.1

1.1

2.1

1.9

3.5

3.7

1.9

1.9

IIa

01-131

35.5

2.1

2.1

1.5

1.5

2.3

2.1

4.6

3.4

1.4

1.6

IIa

20-005

45.6

1.9

2.0

1.2

1.4

2.3

2.1

3.3

2.9

2.1

2.2

IIa

07-005

42.9

1.7

1.9

1.3

1.3

1.9

2.3

4.2

3.9

2.0

1.9

IIa

20-173

43.7

2.0

1.9

1.3

1.3

2.4

2.1

2.5

2.3

2.2

2.3

IIa

07-030

44.5

2.1

2.0

1.2

1.1

2.1

2.0

3.5

3.2

2.1

1.9

IIb

99-184

36.7

1.8

2.0

1.4

1.5

2.1

1.9

4.4

3.3

1.4

1.7

IIb

08-008

39.5**

2.0

1.9

1.2

1.1

2.1

2.0

2.9

2.7

2.2

1.8

IIb

07-050

37.1

1.1

1.2

1.2

1.1

2.1

2.0

4.9

4.9

2.1

2.0

IIb

04-112

44.3

1.7

2.0

1.1

1.3

1.9

2.1

3.5

3.5

1.9

2.0

IIb

06-023

31.5

1.9

2.0

1.4

1.3

2.1

2.1

2.5

2.6

1.9

1.9

IIb

01-129

12.8

1.9

2.0

1.8

1.8

2.1

2.1

2.2

2.1

2.0

2.0

IIb

06-014

19.8

2.0

2.0

1.8

1.6

2.4

2.2

4.7

4.5

1.8

1.6

IIb

07-012

15.8

2.0

1.9

1.7

1.5

2.3

2.0

3.8

3.7

2.2

2.0

IIb

06-041

42

1.8

1.9

1.3

1.2

2.2

2.0

4.6

4.7

1.9

1.8

IIb

99-239

43.4

1.6

1.8

0.9

1.4

1.8

1.9

3.5

3.3

1.3

1.9

IIb

08-029

40.4

1.1

1.3

1.3

1.2

2.1

2.1

3.5

3.2

2.2

1.9

IIb

05-005

21.9

1.7

1.9

1.4

1.3

2.0

2.0

3.1

3.2

1.8

1.8

IIb

20-042

35.9

1.8

1.6

1.6

1.5

2.9

1.9

3.9

3.1

2.3

2.0

IIb

06-004

43.5

1.9

2.1

1.1

1.1

2.1

2.0

3.8

4.2

1.8

1.9

IIb

20-178

27.3

0.9

1.2

1.1

1.3

2.1

2.0

6.0

4.2

1.9

2.1

IIb

06-015

44

1.2

1.3

1.2

1.0

2.3

2.1

3.8

3.9

2.1

1.9

IIb

20-128

39.3

2.3

2.0

1.5

1.5

2.7

1.6

6.2

3.9

1.4

1.7

IIb

08-005

27.4

1.3

1.4

1.2

1.3

2.4

2.0

7.0

5.7

1.3

1.3

IIb

06-008

36.3

1.2

1.1

1.4

1.2

2.2

2.0

3.0

2.9

2.3

1.9

IIb

08-004

33.5**

1.7

1.9

1.2

1.2

2.8

2.8

2.8

2.8

2.0

1.8

IIb

04-035

43

1.8

2.1

1.1

1.3

2.1

2.2

3.0

2.7

2.0

1.9

IIb

06-047

35.8

2.0

1.9

1.4

1.2

2.2

2.0

2.7

2.8

2.1

1.8

IIb

06-038

39.6

2.0

2.1

1.1

1.1

2.2

2.1

3.3

3.4

1.5

1.4

IIb

01-091

44

1.9

2.0

1.2

1.4

2.1

2.1

4.6

3.8

1.5

1.7

IIb

06-045

31.1

1.8

1.8

1.6

1.5

2.2

2.0

4.1

3.8

1.7

1.5

Copy number values of 64 UM measured with both SNP and dPCR analysis. Fractional abundance (FA) denotes the fraction of mutant alleles.

*Chromosome 5 imbalance, TTC5 copy number used as reference.

**FA values were calculated with copy number analysis.

Figure 2: Chromosomal aberrations in expression classes. Five recurrent abnormalities distributed over the gene expression classes. Differences between class I and II were mostly seen in chromosome 3 and 6. Gain of chromosome 8q appears to be the best classifier for IIa and IIb subdivision. Non-significant trends were seen in chromosome 1p and 16q.

### Chromosome 3 loss and gain of 6p define tumor class

Loss of chromosome 3 was predominantly presented in class II (a/b) tumors and was only detected in 3 out 25 class I tumors (Supplementary Table S3). In contrast, gain of chromosomal part 6p was most prominent in class I tumors. Moreover, co-occurrence of 6p gain and monosomy 3 in UM was rare and 6p copy number was low in these instances, marginally exceeding the thresholds for gain (Figure 3).

Figure 3: Chromosomal anomalies support UM subdivision. Monosomy 3 and 6p gain divided UM in class I and class II. In the mixed tumors 6p copy number was low.

### Increased 8q copy number in class IIb tumors

Monosomy 3 and 6p gain molecularly characterize class II and class I tumors, respectively, but do not differentiate class IIa from class IIb tumors. However with analysis of the aberrations of 8q and to a lesser extent 16q this turned out to be possible. Class IIb tumors presented a higher 8q copy number than class IIa tumors (Figure 2). Moreover, the mechanism underlying gain of 8q was different in these two UM classes. Amplification of 8q observed in class IIb was almost exclusively caused by isochromosome formation. The limited gain of 8q that was observed in class IIa tumors was more often due to gain of the entire chromosome [23]. Loss of 16q was not significantly differentially distributed over the three UM classes.

### Functional and genetic annotation

To investigate the link between chromosomal aberrations and differentially expressed genes we analyzed potential genetic and functional correlations of the 237 top differentially expressed genes (Figure 1). Functional annotation showed that 6 terms were significantly overrepresented in the 237 genes (Table 2). These terms concerned the immune system (n = 4) and the translation machinery (n = 2).

Table 2: Functional annotation of classifier genes

Term

Count

Bonferroni

Class

Antigen processing and presentation

16

5.8E-11

IIa–IIb

Immune response

34

6.7E-09

IIa–IIb

Translational elongation

15

1.8E-08

I–II

Antigen processing and presentation of peptide antigen

10

3E-08

IIa–IIb

Antigen processing and presentation of peptide antigen via MHC class I

7

0.000039

IIa–IIb

Translation

19

0.00011

I–II

Functional annotation of the most differentially expressed genes between class I and II. After correction for multiple testing, 6 terms were found to be significantly overrepresented.

Among the 237 genes we identified 151 genes that define the difference between class I and class II in general (Supplementary Table S2). These genes were differentially expressed (LFC > 1 or LFC < −1; p < 0.05) between class I and class II but did not significantly differ between class IIa and IIb. Genes on chromosome 3 and 6 were overrepresented with respectively 16 and 12 genes. All of the 16 genes from chromosome 3 showed a lower expression in class II tumors compared to class I tumors corresponding to loss of chromosome 3 in class II. In contrast, the 12 genes on 6p showed either an elevated or a decreased expression in class II. Functional annotation of the differentially-expressed genes showed that the terms all concerned ribosomal and other translational machinery proteins (Table 2).

### The class IIb classifier genes are involved in the immune response

Next the genes differentially-expressed between class IIa and IIb tumors were annotated. In total, 53 genes were significantly differentially expressed (LFC > 1 or < −1; p < 0.05) between class IIa and IIb tumors but did not differ between class I and IIa tumors (Supplementary Table S2). There was no clustering of genes to chromosome 8q or 16q. Remarkably, 26 genes on chromosome 6, many of which are involved in the immune response, revealed a significant overrepresentation. Of the 26 genes located on chromosome 6, a significant part is involved in HLA/T-cell reactions. Moreover, many of the genes that define class IIb, appear to be targets of interferon. These data suggest that part of the different gene expression between class IIb and class I and IIa tumors may be due to tumor-resident non-cancer cells.

### Tumor heterogeneity in UM

Based on the assumption that every tumor cell within a mutant UM l either carries a GNAQ or GNA11 mutation we were able to calculate the fractional abundance (FA), which represents the ratio between cancer cells and non-cancer cells (e.g. fibroblasts, immune cells) within the tumor. Four tumors did not present one of the hotspot mutations in GNAQ/11 when analyzed with dPCR and in these cases the monosomy 3 status was used to calculate the FA. This was warranted by the positive correlation between tumor fraction calculated with GNAQ/GNA11 mutation and monosomy 3 (Supplementary Figure S2).

Figure 4A shows the fraction of tumor cells and their distribution over the three classes. No significant difference was found between class I and IIa but a significantly decreased tumor cell percentage was detected in class IIb tumors, sustaining the previous notion that the different chromosome 6 associated gene expression in class IIb tumors may be due to non-tumor cells. As the copy number of the aberrations in this study is the result of admixture of normal cells and tumor cells, we also calculated corrected copy numbers. Fractional abundance of tumor cells was used for a more precise calculation of chromosome 8q copy numbers in the tumor cells (Figure 4B). The mean 8q copy number after adjustment in class I, class IIa and class IIb, was 2.6, 3.5 and 5.0, respectively. For the other aberrations, a dosage effect is not clinically relevant, and therefore we did not calculate adjusted copy numbers.

Figure 4: Quantification of tumor heterogeneity. Tumor fractions based on GNAQ/11 mutation analyzed with dPCR. Class IIb contains tumors with a significantly lower tumor fraction A. After adjusting chromosome 8q copy number for tumor fraction, significant differences were also found between class IIa and class IIb besides the significant differences between class I and class II B.

### Genetic heterogeneity in UM revealed clonal evolution

Tumor heterogeneity was also determined with fractional abundance of the chromosomal imbalances. Moreover, quantification of tumor heterogeneity with dPCR provided information on the sequence of events. In UM 20-042 for example, 71.8% of the cells contained the GNAQ Q209P mutation while a chromosome 3 copy number of 1,64 indicated that 36% of the cells contained monosomy 3 and thereby revealed that monosomy 3 occurred after GNAQ mutation (Table 1). Because 8q copy number in UM varied widely we used karyotyping to order monosomy 3 and 8q gain instead (Figure 5A). Among 39 UM with monosomy 3 and 8q gain, three tumors (20-173, 06-015 and 06-023) displayed 8q gain in all tumor cells, but some cells contained two copies of chromosome 3. This indicates that 8q gain occurred before monosomy 3 in these tumors.

Figure 5: Chronology of genomic imbalances in UM development. Karyograms of 06-023 revealed monosomy 3 heterogeneity while homogeneous isochromosome 8q indicated that it preceded monosomy 3 A. UM 20-042 displayed heterogeneity for monosomy 3 and loss of 1p, indicated by circles that represent the fractions of the tumor containing a specific chromosomal aberration B. Integration of tumor heterogeneity with gene-expression classes suggested distinct chronology in UM development in class I and class II C. Two UM in class I did not contain 6p gain but presented L1p and M3/8q+ genotypes respectively. The class II UM without gain of 8q did present monosomy 3.

With dPCR calculated tumor fractions we furthermore deduced that loss of 16q seems to be more common as a quaternary event following monosomy 3. Ultimately, 1p loss occurs as the last event in this series of five recurrent aberrations (Figure 5B). For monosomy 3 and chromosome 16q loss, 14 tumors showed heterogeneity (ΔX > 0.2) and all showed a larger monosomy 3 tumor fraction. For the combination of loss of 16q and 1p, only six UM were informative and all showed a higher 16q tumor fraction compared to the 1p fraction. Temporal distribution of chromosomal aberrations can also be determined with the number of imbalances in a tumor and this was applied for 6p gain [24]. In class I, three UM presented only one imbalance and in all of these cases it concerned 6p gain. This indicated that in class I tumors, 6p gain most likely occurred first, following the GNAQ/11 mutation (Figure 5C).

## DISCUSSION

For UM excellent molecular markers exist to predict disease outcome [25]. Chromosomal imbalances and gene expression profiles accurately predict disease progression in UM but so far have not resulted in improved treatment. We set out to integrate expression and genomic profiles of UM to study the underlying mechanisms as this may facilitate a knowledge based treatment. By applying a rigorous data reduction to the gene expression analysis, we revealed three expression classes. These classes overlapped with classes that were identified by unsupervised cluster analysis but divided the UM with a bad prognosis into two groups that were designated class IIa and class IIb. The immune response and the translation machinery were the only two functions significantly enriched among the 237 genes that defined the three classes.

The five recurrent aberrations that were analyzed, were non-randomly distributed over the expression classes. In our tumor set, monosomy 3 and 6p gain clearly divided the tumors into class I and class II (a/b) (Figure 3). Based on the fact that monosomy 3 and 6p gain are present in different UM it is commonly assumed that these aberrations are mutually exclusive [2628]. Quantification with digital PCR however showed some degree of 6p gain in combination with monosomy 3 and thereby indicates that the mutual exclusivity is not absolute. Though, our observations suggest that tumor clones presenting both monosomy 3 and 6p gain are unsuccessful and remain small. Class IIa and class IIb can be distinguished with 8q copy number. The number of chromosome 8q copies is highest in class IIb and this was related to a different underlying mechanism between classes. Class IIa more commonly yielded tumors with gain of a complete chromosome 8 while class IIb tumors almost exclusively displayed 8q copy gains due to isochromosome formation [23]. The formation of chromosome 8q isochromosomes is a phenomenon earlier described in UM [2931].

Remarkably, most of the genes highly expressed in class IIb tumors are involved in immune regulation and located on chromosome 6 which is rarely gained in class II tumors. We therefore postulated that these differentially-expressed genes were not the result of differences in tumor cells, but expressed by stromal cells (e.g. fibroblasts and immune cells). This notion was sustained by our calculation of the tumor cell fraction in the total tumor mass. The tumor cell percentage was similar between class I and class IIa but was significantly lower in class IIb tumors suggesting a higher percentage of non-cancer cells in class IIb tumors. The role of the immune system in UM development has not yet been elucidated but correlations with survival and other prognostic markers have been made [3234]. In UM, a bad prognosis is correlated with the presence of immune infiltrate [33] and an immune infiltrate may be what is causing the expression profile of class IIb (Table 2). Importantly, class IIa tumors that present a high tumor percentage, leaving few space for immune cell infiltrate, are also correlated with a bad prognosis. This implies that that the immune cells are not a prerequisite for the class IIa tumors to progress and metastasize.

The calculated tumor percentages of each chromosome loss furthermore allowed us to deduce the sequence of events in UM. Based on dPCR analysis we calculated which percentage of cells contained specific chromosomal losses (Figure 2). Because chromosomal gains could not be used to determine fractional abundance we used alternative approaches for 6p gain and 8q gain. In class I tumors 6p aberrations appear to occur first based on the fact that 22 out of 25 UM presented 6p gain in this class and because the tumors with only one chromosomal imbalance all presented 6p gain. Previous reports suggest that loss of chromosome 3 precedes the addition of 8q [35] but cytogenetic analysis of informative UM in our dataset supported the opposite. Karyotyping revealed monosomy 3 heterogeneity in tumors that are homogeneous for 8q gain and thereby suggested that chromosome 8q imbalance occurs prior to monosomy 3 in class II tumors. With regard to the sequence of chromosomal losses, we deduced that monosomy 3 occurs first, followed by loss of 16q and loss of 1p respectively

Surprisingly for such an early event, 8q copy number made molecular discrimination of the three gene-expression classes possible when we adjusted chromosome 8q copy number for tumor fraction. In a dosage-dependent manner, significant differences were revealed between the three expression classes. Previous studies showed that 8q copy number correlates with survival but no significant survival differences were revealed between class IIa and class IIb even though 8q copy number differed between these classes (Figure 4B). Both expression differences and genomic differences nevertheless suggest that class IIa and IIb represent different entities that may require different treatments.

In short, with integrated analysis of gene expression and genomic aberrations we were able to distinguish three molecular classes in UM. Based on an apparent immune infiltrate in class IIb we hypothesize that class II can progress either with or without immune involvement. Defining these classes gives us the opportunity to investigate the mechanism behind immune involvement in UM. Chromosome 8q is likely to be a candidate locus to search for targets since it is differentially changed in class IIa compared to class IIb. Dividing UM on a molecular level reveals tumor diversity and can lead to new possibilities for therapeutic intervention.

## MATERIALS AND METHODS

We used tumor material from 64 enucleated eyes of uveal melanoma patients that had been enucleated at the Leiden University Medical Center, Leiden, The Netherlands, between 1999 and 2008. Written informed consent was obtained for all patient samples. None of the tumors had prior treatment and we only used tumors with a follow-up time of at least 4 years. The maximum follow-up was 14 years. The average age at enucleation was 60.6 years (range 13 to 88); 33 patients were male and 31 female.

Tumor material was snap frozen using 2-methyl butane and RNA and DNA was isolated using the RNeasy mini kit and QIAmp DNA minikit, respectively, (both Qiagen, Valencia, USA) from 20 sections of 20 μm according to the manufacturer’s guidelines.

### Gene expression profiling

Gene expression was determined using the Illumina HumanHT-12 v4 chip containing 47,000 probes across the whole genome that will be referred to as genes in the remainder of the text. Genes in the supervised cluster analysis with a log-fold change (LFC) larger than 1 or smaller than −1 and a p-value smaller than 0.05 were labelled “significantly differentially expressed”. For differences between subgroups i.e. I versus II, we corrected for differences between IIa and IIb classified as LFC smaller than –0.5 or greater than 0.5 and a p-value smaller than 0.05. These most differentially-expressed genes were annotated and biological processes were analyzed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) [36, 37].

### SNP analysis

SNP microarray analysis was used to determine chromosomal aberrations. Two types of SNP microarray chips were used, the Affymetrix 250K_NSP-chip, with ~250,000 probes across the genome and the Affymetrix Cytoscan HD chip, with ~750,000. The first 28 samples were analyzed with the Affymetrix 250K_NSP chip, and the remaining 36 samples with the Affymetrix Cytoscan HD chip. Analysis of the Affymetrix 250K_NSP chips was performed with the ‘Genotyping Console’ to determine the copy number values and the ‘GCT Browser’ to visualize the data (both from Affymetrix, Santa Clara, USA). Affymetrix Cytoscan HD chips were analyzed with ‘ChAS’. Different loci per chromosome were evaluated to adjust for partial gains or deletions. ~200 probes per gene locus were averaged to determine eventual copy number.

### dPCR

#### Copy number analysis

The copy number of chromosome 1p, 3, 6p, 8q, and 16q was determined using probes for CDC42, PPARG, NEDD9, PTK2, and NFAT5 respectively. For calculation of normalized copy numbers, out of 3 control probes located on chromosome 5, 7, and 14 (TERT, VOPP1 and TTC5),TERT (situated at chromosome 5), was selected based on stability in this tumor set. In the two cases with chromosome 5 gain we used TTC5 as reference to calculate the copy numbers. Thresholds for copy number analysis were: loss, <1.9: normal, 1.9–2.1: gain, >2.1− <3.1: amplification, >3.1 (method described by Versluis et al. [23]). In short, 50–60 ng of DNA of each sample was used in a 20 ul reaction volume. Reaction mixture consisted of 2x droplet PCR supermix (Bio-Rad Laboratories, Inc., Hercules, USA), 20x target probe (FAM), 20x reference probe (HEX). Sequence context is provided in Supplementary Table S1. Droplets were generated using a QX100 droplet generator and after the PCR, the plate was loaded into the QX100 droplet reader (both Bio-Rad Laboratories, Inc.). The following end point PCR protocol using a T100 thermal cycler was used: 95°C, 10 min; (94°C, 30sec; 60°C, 1min) 40x; 98°C, 10 min; 4°C, till end. Digital PCR (dPCR) software (QuantaSoft) reads the positive and negative droplets in each sample and plots the fluorescence droplet by droplet. Fractional abundance was calculated as . The positive droplets represent the concentration of the target allele in the sample. Corrections made for the fraction of tumor cells containing chromosomal aberrations (Z) based on FA were performed as follows: X represents the dPCR value for copy number of chromosome n, and Y represents the tumor fraction of the associated tumor which is calculated by multiplying the FA by 2.

Heterogeneity of imbalances was determined by subtracting copy numbers and a difference of 0.2 copy number was used as threshold (ΔX > 0.2).

$Z=X−(2(1−(Y100)))Y100$

#### GNAQ/11 mutation detection

GNAQ and GNA11 mutations were detected using hydrolysis probes in a multiplex dPCR. 10 ng of sample DNA was used in a 20 ul reaction volume. The reaction mixture consisted of 2x droplet PCR supermix (Bio-Rad Laboratories, Inc.), 20x target probe (FAM), and 20x wildtype probe (HEX). Proprietary probes and primers (Bio-Rad Laboratories, Inc.) were used and the sequence context is provided in Supplementary Table S1. Droplet generation, thermal cycling and reads were similar to the method described in the previous section. PCR to end point protocol used: 95°C, 10 min; (94°C, 30 sec; 55°C, 1 min) 40x; 98°C, 10 min; 4°C, till end.

### Karyotyping

Following enucleation, a small part of each tumor was sent out for cell culture. Following mechanical dissection of the tumor biopsy, cells were washed and placed into one flask with RPMI 1640 (15% fetal bovine serum [Invitrogen, Breda, The Netherlands]) medium and another flask with Amniochrome II (Cambrix Bio Science, Verviers, Belgium). The flasks were cultured at 37°C with 5% CO2 for up to 4 weeks and harvested when at least 75% of the surface was covered with cells (after a mean of 18 days; SD, 9.4 days). When cell culturing was successful, conventional karyotyping was performed, to determine the presence of chromosomal changes. Two independent observers assessed all evaluations and scores, each without knowledge of the results obtained by the other investigator, to ensure accuracy of quantification of the slides. In case of a difference, consensus was reached during a simultaneous session. Cytogenetic analysis was performed on GTG-banded (G-banding with trypsin and Giemsa) metaphases. In the case of a normal karyotype, at least 20 metaphases were analyzed. When an abnormal clone was detected in the first ten karyotyped cells, no further analysis was performed; when three cells with loss of 1 copy of chromosome 3 were observed, monosomy 3 was identified.

### Statistical analysis

For gene expression analysis, the statistical programming language R was used (R: A Language and Environment for Statistical Computing, R Core Team, R foundation for Statistical Computing, Vienna, Austria, 2014). Since data have been obtained in two batches, a batch effect correction was applied. The R packages used were: ‘ber’ for batch correction and ‘lumi’ for unsupervised clustering. To compare survival between UM patients in different classes Kaplan-Meier functions were plotted. Survival analysis was performed using the log-rank test. The Pearson’s correlation test was used to assess the correlation between the fraction of tumor cells containing GNAQ/11 and the fraction of tumor cells containing chromosome 3 aberrations. The Chi square test was used to test frequency differences in chromosomal aberrations in different subgroups. Likelihood ratios were used for cases in which the Chi square assumption was violated. α = 0.05 was used as threshold for significance in all tests. For statistical analysis SPSS V.20.0.1 (IBM SPSS statistics, IBM corporation, Armonk, New York, USA) was used.

## ACKNOWLEDGMENTS AND FUNDING

The authors thank Ronald van Eijk from the Department of Pathology of the LUMC for helpful discussion. Eddy van Collenburg and Roy Wiggers (BioRad Laboratories) for providing probes and primers. Guido J.E.J. Hooiveld (WU Agrotechnology & Food Sciences) for use of the digital PCR facility and Frank de Gruijl for his expertise in mathematics.

## CONFLICTS OF INTEREST

All authors have declared that no conflict of interest exists.

## GRANT SUPPORT

This study was supported by grant UL 2011-4991 from the Dutch Cancer Society (KWF).

## REFERENCES

1. Van Raamsdonk CD, Bezrookove V, Green G, Bauer J, Gaugler L, O’Brien JM, Simpson EM, Barsh GS, Bastian BC. Frequent somatic mutations of GNAQ in uveal melanoma and blue naevi. Nature. 2009; 457:599–602.

2. Van Raamsdonk CD, Griewank KG, Crosby MB, Garrido MC, Vemula S, Wiesner T, Obenauf AC, Wackernagel W, Green G, Bouvier N, Sozen MM, Baimukanova G, Roy R, Heguy A, Dolgalev I, Khanin R, et al. Mutations in GNA11 in uveal melanoma. NEnglJMed. 2010; 363:2191–2199.

3. Eskelin S, Pyrhonen S, Summanen P, Hahka-Kemppinen M, Kivela T. Tumor doubling times in metastatic malignant melanoma of the uvea: tumor progression before and after treatment. Ophthalmology. 2000; 107:1443–1449.

4. Singh AD, Turell ME, Topham AK. Uveal melanoma: trends in incidence, treatment, and survival. Ophthalmology. 2011; 118:1881–1885.

5. Prescher G, Bornfeld N, Hirche H, Horsthemke B, Jockel KH, Becher R. Prognostic implications of monosomy 3 in uveal melanoma. Lancet. 1996; 347:1222–1225.

6. Onken MD, Worley LA, Ehlers JP, Harbour JW. Gene expression profiling in uveal melanoma reveals two molecular classes and predicts metastatic death. Cancer Res. 2004; 64:7205–7209.

7. Scholes AG, Damato BE, Nunn J, Hiscott P, Grierson I, Field JK. Monosomy 3 in uveal melanoma: correlation with clinical and histologic predictors of survival. Invest OphthalmolVisSci. 2003; 44:1008–1011.

8. Seddon JM, Polivogianis L, Hsieh CC, Albert DM, Gamel JW, Gragoudas ES. Death from uveal melanoma. Number of epithelioid cells and inverse SD of nucleolar area as prognostic factors. ArchOphthalmol. 1987; 105:801–806.

9. McLean IW, Foster WD, Zimmerman LE. Uveal melanoma: location, size, cell type, and enucleation as risk factors in metastasis. HumPathol. 1982; 13:123–132.

10. Onken MD, Worley LA, Person E, Char DH, Bowcock AM, Harbour JW. Loss of heterozygosity of chromosome 3 detected with single nucleotide polymorphisms is superior to monosomy 3 for predicting metastasis in uveal melanoma. ClinCancer Res. 2007; 13:2923–2927.

11. Lake SL, Coupland SE, Taktak AF, Damato BE. Whole-genome microarray detects deletions and loss of heterozygosity of chromosome 3 occurring exclusively in metastasizing uveal melanoma. Invest OphthalmolVisSci. 2010; 51:4884–4891.

12. McCannel TA, Burgess BL, Nelson SF, Eskin A, Straatsma BR. Genomic identification of significant targets in ciliochoroidal melanoma. Invest OphthalmolVisSci. 2011; 52:3018–3022.

13. Singh AD, Aronow ME, Sun Y, Bebek G, Saunthararajah Y, Schoenfield LR, Biscotti CV, Tubbs RR, Triozzi PL, Eng C. Chromosome 3 status in uveal melanoma: a comparison of fluorescence in situ hybridization and single-nucleotide polymorphism array. Invest OphthalmolVisSci. 2012; 53:3331–3339.

14. Gordon KB, Thompson CT, Char DH, O’Brien JM, Kroll S, Ghazvini S, Gray JW. Comparative genomic hybridization in the detection of DNA copy number abnormalities in uveal melanoma. Cancer Res. 1994; 54:4764–4768.

15. Trolet J, Hupe P, Huon I, Lebigot I, Decraene C, Delattre O, Sastre-Garau X, Saule S, Thiery JP, Plancher C, Asselain B, Desjardins L, Mariani P, Piperno-Neumann S, Barillot E, Couturier J. Genomic profiling and identification of high-risk uveal melanoma by array CGH analysis of primary tumors and liver metastases. Invest OphthalmolVisSci. 2009; 50:2572–2580.

16. Speicher MR, Prescher G, du MS, Jauch A, Horsthemke B, Bornfeld N, Becher R, Cremer T. Chromosomal gains and losses in uveal melanomas detected by comparative genomic hybridization. Cancer Res. 1994; 54:3817–3823.

17. Ghazvini S, Char DH, Kroll S, Waldman FM, Pinkel D. Comparative genomic hybridization analysis of archival formalin-fixed paraffin-embedded uveal melanomas. Cancer GenetCytogenet. 1996; 90:95–101.

18. Naus NC, van Drunen E, de Klein A, Luyten GP, Paridaens DA, Alers JC, Ksander BR, Beverloo HB, Slater RM. Characterization of complex chromosomal abnormalities in uveal melanoma by fluorescence in situ hybridization, spectral karyotyping, and comparative genomic hybridization. Genes ChromosomesCancer. 2001; 30:267–273.

19. White JS, Becker RL, McLean IW, Director-Myska AE, Nath J. Molecular cytogenetic evaluation of 10 uveal melanoma cell lines. Cancer GenetCytogenet. 2006; 168:11–21.

20. Sisley K, Parsons MA, Garnham J, Potter AM, Curtis D, Rees RC, Rennie IG. Association of specific chromosome alterations with tumour phenotype in posterior uveal melanoma. BrJCancer. 2000; 82:330–338.

21. White JS, McLean IW, Becker RL, Director-Myska AE, Nath J. Correlation of comparative genomic hybridization results of 100 archival uveal melanomas with patient survival. Cancer GenetCytogenet. 2006; 170:29–39.

22. Tschentscher F, Husing J, Holter T, Kruse E, Dresen IG, Jockel KH, Anastassiou G, Schilling H, Bornfeld N, Horsthemke B, Lohmann DR, Zeschnigk M. Tumor classification based on gene expression profiling shows that uveal melanomas with and without monosomy 3 represent two distinct entities. Cancer Res. 2003; 63:2578–2584.

23. Versluis M, de Lange MJ, van Pelt SI, Ruivenkamp CA, Kroes WG, Cao J, Jager MJ, Luyten GP PA. Digital PCR Validates 8q Dosage as Prognostic Tool in Uveal Melanoma. PLoSOne. 2015; 10:e0116371.

24. Hoglund M, Gisselsson D, Hansen GB, White VA, Sall T, Mitelman F, Horsman D. Dissecting karyotypic patterns in malignant melanomas: temporal clustering of losses and gains in melanoma karyotypic evolution. IntJCancer. 2004; 108:57–65.

25. Singh AD, Shields CL, Shields JA. Prognostic factors in uveal melanoma. Melanoma Res. 2001; 11:255–263.

26. Bonaldi L, Midena E, Filippi B, Tebaldi E, Marcato R, Parrozzani R, Amadori A. FISH analysis of chromosomes 3 and 6 on fine needle aspiration biopsy samples identifies distinct subgroups of uveal melanomas. JCancer ResClinOncol. 2008; 134:1123–1127.

27. Parrella P, Sidransky D, Merbs SL. Allelotype of posterior uveal melanoma: implications for a bifurcated tumor progression pathway. Cancer Res. 1999; 59:3032–3037.

28. Young TA, Burgess BL, Rao NP, Gorin MB, Straatsma BR. High-density genome array is superior to fluorescence in-situ hybridization analysis of monosomy 3 in choroidal melanoma fine needle aspiration biopsy. MolVis. 2007; 13:2328–2333.

29. Prescher G, Bornfeld N, Friedrichs W, Seeber S, Becher R. Cytogenetics of twelve cases of uveal melanoma and patterns of nonrandom anomalies and isochromosome formation. Cancer GenetCytogenet. 1995; 80:40–46.

30. Horsman DE, Sroka H, Rootman J, White VA. Monosomy 3 and isochromosome 8q in a uveal melanoma. Cancer GenetCytogenet. 1990; 45:249–253.

31. Couturier J, Saule S. Genetic determinants of uveal melanoma. DevOphthalmol. 2012; 49:150–165.

32. de Waard-Siebinga I, Hilders CG, Hansen BE, van Delft JL, Jager MJ. HLA expression and tumor-infiltrating immune cells in uveal melanoma. Graefes ArchClinExpOphthalmol. 1996; 234:34–42.

33. Bronkhorst IH, Ly LV, Jordanova ES, Vrolijk J, Versluis M, Luyten GP, Jager MJ. Detection of M2-macrophages in uveal melanoma and relation with survival. Invest OphthalmolVisSci. 2011; 52:643–650.

34. Makitie T, Summanen P, Tarkkanen A, Kivela T. Tumor-infiltrating macrophages (CD68(+) cells) and prognosis in malignant uveal melanoma. Invest OphthalmolVisSci. 2001; 42:1414–1421.

35. Landreville S, Agapova OA, Harbour JW. Emerging insights into the molecular pathogenesis of uveal melanoma. FutureOncol. 2008; 4:629–636.

36. Huang dW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009; 37:1–13.

37. Huang dW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. NatProtoc. 2009; 4:44–57.