DNA alteration‐based classification of uveal melanoma gives better prognostic stratification than immune infiltration, which has a neutral effect in high‐risk group

Abstract Background In uveal melanomas, immune infiltration is a marker of poor prognosis. This work intended to decipher the biological characteristics of intra‐tumor immune population, compare it to other established biomarkers and to patients' outcome. Methods Primary, untreated, and mainly large uveal melanomas with retinal detachment were analyzed using: transcriptomic profiling (n = 15), RT‐qPCR (n = 36), immunohistochemistry (n = 89), Multiplex Ligation‐dependent Probe Amplification (MLPA) for copy number alterations (CNA) analysis (n = 89), array‐CGH (n = 17), and survival statistics (n = 86). Results Gene expression analysis divided uveal melanomas into two groups, according to the IFNγ/STAT1‐IRF1 pathway activation. Tumors with IFNγ‐signature had poorer prognosis and showed increased infiltration of CD8+ T lymphocytes and macrophages. Cox multivariate analyses of immune cell infiltration with MLPA data delineated better prognostic value for three prognostic groups (three‐tier stratification) than two (two‐tier stratification). CNA‐based model comprising monosomy 3, 8q amplification, and LZTS1and NBL1 deletions emerged as the best predictor for disease‐free survival. It outperformed immune cell infiltration in receiver operating characteristic curves. The model that combined CNA and immune infiltration defined risk‐groups according to the number of DNA alterations. Immune cell infiltration was increased in the high‐risk group (73.7%), where it did not correlate with patient survival, while it was associated with poorer outcome in the intermediate risk‐group. Conclusions High degree of immune cell infiltration occurs in a subset of uveal melanomas, is interferon‐gamma‐related, and associated with poor survival. It allows for two‐tier stratification, which is prognostically less efficient than a three‐tier one. The best prognostic stratification is by CNA model with three risk‐groups where immune cell infiltration impacts only some subgroups.


Tumor series
For 89/91 tumors, clinico-radiological data and follow-up information were available and obtained from hospital and general practitioner records and histopathological details from haematoxylin and eosin (H&E) sections. Age at operation varied from 35 years to 90 years with a mean of 64.6 years and a median of 67 years. The male to female ratio was 1.1. The median basal diameter was 18mm and median thickness was 11.1mm with tumors mainly belonging to T3 and T4 (91%) categories. TNM staging was performed according to 7 th edition of American Joint Committee on Cancer guidelines 1 Table S1).

RNA Extraction
For the 15 snap-frozen samples, total RNA was extracted using guanidine isothiocyanate/cesium chloride procedure 3 . The first-strand cDNA was synthesized from 2µg of total RNA with M-MLV reverse transcriptase (catalog # 28025) (Invitrogen, Life Technologies) and random hexamer primer (Fermentas), as recommended by manufacturer.
For 21 samples stored in RNA later (-80°C), RNA was extracted using NucleoSpin RNA (catalog # 740955) (Macherey-Nagel, Düren, Germany). For each tumor, 20 cryosections of 20µm thickness were used for extraction. Isolation of total RNA was performed using NucleoSpin RNA (catalog # 740955) (Macherey-Nagel, Düren, Germany) as per the manufacturer"s guidelines. In brief, the cryosections were homogenized with RA1 buffer and beta-mercaptoethanol and filtered using NucleoSpin filters. The filtrate was mixed with 70% ethanol and bound to NucleoSpin RNA column. The column was washed and dried using buffers RAW2 and RA3, followed by elution using RNase-free water. Total RNA was treated with Turbo DNA-free TM kit (catalog # AM1907) (Ambion, Life Technologies, Carlsbad, CA) to remove the residual DNA. RNA was quantified using NanoDrop (ThermoScientific, Wilmington, DE). The A 260 /A 280 ratio for the samples ranged from 1.88-1.99. The melanin concentration in RNA was assessed by measuring the absorbance at 320nm and samples with absorbance >0.1 were re-purified on RNA column. The first strand cDNA was synthesized from 1µg of total RNA using RevertAid H minus first strand cDNA synthesis kit and random hexamer primer (catalog # K1632) (Fermentas, ThermoScientific), as per manufacturer's guidelines.

Array hybridization
RNA was quantified using NanoDrop (ThermoScientific, Wilmington, DE) and samples with A 260 /A 280 ratio between 1.70-1.83 were included. Also, only samples with intact RNA indicated by sharp, clear bands for 18S and 28S rRNA on agarose gels were utilized. Human Genome U133 Plus 2.0 arrays (Affymetrix, Santa Clara, CA) were used to generate gene expression profiles from 5µg of total RNA following manufacturer"s protocol. The intensity files were obtained by processing the images using Microarray suite 5.0 gene expression software. Probe sets present in <90% of the samples were excluded from analysis.

Array data analysis
Additional normalization was performed to reduce variation between individual hybridizations. Therefore, the mean intensity of each probe set was defined based on the intensity of all tumors. Then, individual intensities were divided by their respective mean intensity and converted to log2. The first 1500 probe sets with highest standard deviation (SD) were submitted for statistical analysis using TIGR MultiExperiment Viewer platform (MeV4.8.1, http://www.tm4.org/mev/) 4 . Subgroups of uveal melanoma were unraveled by performing hierarchical clustering with Pearson"s correlation. Estimation of sampling distribution was assessed by bootstrapping, using 100 iterations. Differentially expressed probe sets were identified by T-test applying high stringency (5000 permutations, p-value 0.001) and low stringency (1000 permutations, p-value 0.01) conditions, with adjusted Bonferroni method of p-value correction 5 .
The annotation of statistically significant probe sets was performed using InnateDB 6 by applying a false-discovery rate cut-off of 0.05 and p-value was correction by Benjamini-Hochberg method.

Immune genes
The immunologically relevant genes were compiled from three databases: 5998 genes from ImmPort (https://immport.niaid.nih.gov), 1773 genes related to "immune system process" from Gene Ontology database (http://amigo.geneontology.org) and 1326 immunologically important genes from Wiki for Immunology (http://wiki.geneontology.org/index.php/Immunology). There were a total of 6487 genes after accounting for common genes between databases and 4840 genes with available probe sets (n=11,103) were included for further analysis.

GBP1 immunohistochemistry
4µm sections were deparaffinized, endogenous peroxidase activity was blocked with 3% hydrogen peroxide and heat-induced epitope retrieval was performed in 0.1M, pH 6 in dark for 30 minutes. The slides were counterstained with haematoxylin and coverslipped using an aqueous mounting medium. All reactions were carried out at room temperature.

GBP1 quantification
Digital images of GBP1-stained slides were obtained at 20X magnification using a whole slide scanner, LeicaSCN400 (Leica Microsystems GmbH, Wetzlar, Germany). The images were stored and managed using a web-enabled image server, Digital Image Hub (DIH) (SlidePath, Leica). Quantification of GBP1 was performed using the image analysis software, Tissue IA 2.0 (SlidePath, Leica) using "measure stained area" algorithm. First, the total tumor area (reference area) to be analyzed was outlined manually using the annotation tool. Image segmentation was performed using two thresholds that were manually set. The first was a tissue segmentation threshold to separate the tumor tissue from background. The second was a color definition threshold to define the positively stained pixels. In order to ensure that the latter adequately reflected all positive staining, it was tested on multiple areas within the same slide and on multiple different slides. The output results included: total area of the annotated tumor (reference area), positively stained tumor area, average staining intensity and concentration of the stain within the tissue sample. The ratio of the positive tumor area with respect to the reference tumor area was computed and was referred to as the GBP1 labeling index 7 .

Quantitative real-time reverse transcription polymerase chain reaction (RT-qPCR)
The PCR mix was prepared as follows: PCR-grade water 3.4µl, forward and reverse primers 0.8µl each (5µM) and master mix 10µl. For one reaction, 15µl of PCR mix and 5µl of cDNA

Statistical methods
Continuous variables are reported as mean (standard deviation) and median (range), and are also categorized; binary and categorical variables are reported as frequencies and relative frequencies.
Cox proportional hazards regression model was used to investigate the association between each variable and disease-free survival (DFS).
DFS was defined as time from enucleation to radiological diagnosis of metastases for patients developing metastases or death. Patients still alive and metastases-free at the time of analysis were included only if at least 24 months post-enucleation follow-up was available.

i. Variable selection for Cox multivariate analysis
All categorical covariates were transformed into numeric codes before they were entered into Cox model. Potential prognostic biomarkers were first investigated in univariate analyses. ii. Assessment of model performance The performance of both univariate and multivariate Cox models were quantified using Harrell"s C-index, a measure of concordance. A C-index of 0.5 indicates that the model has no discriminative ability, whereas a C-index of 1 indicates that the model perfectly distinguishes between those with poor prognosis (early metastasis) from those with good prognosis (late metastasis). The models were also evaluated in a time-specific manner by constructing time-dependent ROC curves 10

iii. Prognostic indices
The prognostic score was based on the estimated regression coefficient of all variables in the final Cox model retained. The score of each variable was its estimated coefficient rounded to the first decimal point (using Microsoft Excel Round function) and multiplied by 10.
Individual patient score was obtained by adding up the scores of unfavourable factors present.
The samples were split into risk-groups based on tertiles (three groups) and quartiles (four groups). The KM plots of risk groups were compared by log-rank test (both global and pairwise comparisons with Benjamini-Hochberg adjustment of p-value) ( Table M3).

iv. Leave-one-out cross-validation
The misclassification error associated with each of the models was estimated by leave-oneout method. In each iteration, one observation was omitted and the prognostic score of the omitted observation was predicted based on the model built on the (n-1) other observations.
This was repeated for all observations, resulting in cross-validated scores. The samples were split into risk groups based on the cross-validated scores using the same cut-off as used for the original models. The prediction was considered as "true" if the predicted and original riskgroups were the same and as "false" if, they were different. The misclassification error rate was determined based on the number of false observations.

External validation of CNA and immune-CNA models on TCGA data
The copy number alterations, tumor-infiltrating-lymphocyte (TIL) density, tumor-associatedmacrophage (TAM) density, time to metastasis/last follow-up information were obtained from data supplement of Robertson, et al 12 . The combined TIL and TAM densities were taken into consideration to categorize tumors as immune-high and immune-low groups ( Table   M4). The LZTS1 deletion and NBL1 deletion status was obtained from cBioPortal 13,14 .