Integrative Analysis Identifies Four Molecular and Clinical Subsets in Uveal Melanoma

Abstract

Comprehensive multiplatform analysis of 80 uveal melanomas (UM) identifies four molecularly distinct, clinically relevant subtypes: two associated with poor-prognosis monosomy 3 (M3) and two with better-prognosis disomy 3 (D3). We show that BAP1 loss follows M3 occurrence and correlates with a global DNA methylation state that is distinct from D3-UM. Poor-prognosis M3-UM divide into subsets with divergent genomic aberrations, transcriptional features, and clinical outcomes. We report change-of-function SRSF2 mutations. Within D3-UM, EIF1AX- and SRSF2/SF3B1-mutant tumors have distinct somatic copy number alterations and DNA methylation profiles, providing insight into the biology of these low- versus intermediate-risk clinical mutation subtypes.

Keywords: EIF1AX; GNA11; GNAQ; SF3B1; SRSF2; TCGA; molecular subtypes; monosomy 3; noncoding RNA; uveal melanoma.

Figure 1. Genomic Landscape of Primary UM

(A) Unsupervised clustering of somatic copy numberalterations(SCNAs) separated 80 primary UM into four clusters: 1 (n = 15),2(n = 23), 3(n = 22), and 4(n = 20), ordered by increasing chromosomal instability. The upper covariate tracks show SCNA clusters (1–4), chromosome 3 and 8q copy number, and ploidy level. The heatmap shows somatic copy number ratio (diploid = 0, white). Lower covariate tracks show (i) clinical outcome; (ii) BAP1 mRNA expression; (iii) unsupervised clusters for DNA methylation, mRNA, lncRNA, and miRNA; (iv) mutations in G-protein-signaling genes, splicing factors, and EIF1AX; (v) BAP1 alterations that include alternate splicing and rearrangements detected by assembly of DNA-seq and RNA-seq data. (B) BAP1 mRNA expression, grouped by SCNA clusters, with BAP1 alteration status determined by at least one method in (A). Dots show all data values. Box plots show median values, and the 25th to 75th percentile range in the data, i.e., the interquartile range (IQR). Whiskers extend 1.5 times the IQR. (C) Cancer cell fractions for chromosome 3 loss, BAP1 alterations, and other somatic mutations on chromosome 3, for tumors with BAP1 alterations detected either by standard SNP-indel algorithms or by local reassembly of WES data. Lines connect events that occurred in the same tumor. (D) Schematic depicting a probable sequence of somatic events resulting in those detected in the cluster 3 case V4-A9EO (M3, BAP1 mutation, WGD, then isochromosome 8q). See also Figure S1 and Table S1.

Figure 2. DNA Methylation Landscape in Primary UM

Unsupervised clustering of DNA methylation data, with the heatmap showing beta values ordered by DNA methylation clusters. CpG locus types (island, shore, and shelf) are indicated at the left border. Covariate tracks show unsupervised clusters for four other genomic data types, clinical outcomes, chromosome 3 and 8q copy number status, specific gene alterations, and gender. SF3B1 and EIF1AX mutations were statistically associated with the clusters (*p < 0.01, Fisher’s exact test). LOH, loss of heterozygosity.

Figure 3. Gene Expression Patterns in UM

The upper heatmap shows unsupervised consensus clustering for RNA-seq data of mRNA (left) or lncRNA (right) expression. Covariate annotation tracks show selected genomic and clinical features. The lower heatmap displays the expression profiles of 12 genes used in a prognostic test for the risk of developing metastasis (Harbour, 2014), with blue text highlighting genes that are on chromosome 3. High-risk primary tumors show low expression of eight of these genes and high expression of four genes (yellow versus green panels at the left). BAP1 structural alterations that include alternative splicing and rearrangements were detected by assembly of RNA-seq and DNA-seq data. Leukocyte fraction was estimated from DNA methylation data. LOH, loss of heterozygosity. See also Figures S2 and S3; Table S2. *, **, *** p value < 0.1, 0.01, 10-10, Fisher’s Exact or Chi-square test.

Figure 4. Immune Gene Expression in M3-versus D3-UM

Heatmap for 80 primary UM, highlighting mRNA expression levels of key immunological genes that represent the interferon-γ pathway, T cell cytolytic enzymes, chemokine factors, immunosuppressive factors, and macrophage markers, as well as individual immune checkpoint blockade genes (CD274, PDCD1LG2, PDCD1, CTLA4, IDO1, and TIGIT). Samples were separated by D3 versus M3 status, and sorted from lowest (left) to highest (right) CD8A expression level. Covariate tracks show mRNA, lncRNA, miRNA, PARADIGM, DNA methylation, and SCNA clusters. Leukocyte fraction was estimated from DNA methylation data. See also Figure S4.

Figure 5. Integrative Pathway Analysis of UM

(A) Heatmap of hierarchically clustered PARADIGM inferred pathway levels (IPLs) for 80 primary UMs. Samples are clustered into five groups (top horizontal track). Below this are cluster memberships for other platforms, and for chromosome 3 and 8q copy number, then IPL profiles for the MYC/MAX and MYC/MAX/MIZ1 complexes. The main heatmap shows PARADIGM features or nodes that have at least ten downstream regulatory targets and are differentially active in one-cluster-versus-othercomparisons; the annotation panel to the left indicates the cluster(s) in which a node satisfies these conditions. The vertical colored bars on the right highlight sets of pathway nodes that belong to common biological processes: MAPK/PI3K-AKT (purple), hypoxia (magenta), DNA damage repair/response (green), and immune response (blue). LOH, loss of heterozygosity. (B) Distributions of DDR pathway score and abundance for selected proteins, from RPPA data for M3/BAP1-aberrant versus D3/SF3B1-mutant UM(n = 11). PKC-α_pS657 denotes PKC-α phosphorylated at S657. Box plots show median values and the 25th to 75th percentile range in the data, i.e., the IQR. Whiskers extend 1.5 times the IQR. Dots show all data values. See also Figure S5 and Table S3.

Figure 6. Pathway and Regulators that were Differentially Active in Transcriptional Subtypes 3 and 4

Correlation network for transcriptional (lncRNA) subtype 3 (top) and subtype 4 (bottom), showing PARADIGM pathway features, (hierarchical) MARINa regulators, and lncRNAs. Red and blue lines indicate Spearman correlations (|rho| > 0.5) between the expression of a differentially expressed lncRNA and inferred activity of a differentially active PARADIGM or MARINa feature. The color of each node reflects differential expression for a lncRNA, and relative activity for a PARADIGM/MARINa feature (red for overexpressed/active, blue for underexpressed/inactive). See also Table S4.

Figure 7. Good-Prognosis D3-UM and Poor-Prognosis M3-UM Separate into Distinct Biological Subsets

(A) Kaplan-Meier plots and log-rank p values for the clinical event of UM metastasis for M3-versus D3-UM, then for unsupervised clusters for DNA methylation, SCNA, lncRNA, and mRNA. The number of cases and events in a cluster are shown on the plots. Median event times for clusters 3 and 4 were 10.8 versus 42.6 months for SCNA (p = 0.002, p = 0.01 with a Bonferroni correction [BC]); 13.0 versus >30 months for lncRNA (p = 0.19, p = 1.0 with BC); and 13.5 versus 30.0 months for mRNA (p = 0.43, p = 1.0 with BC). (B) Schematicof D3-UM and M3-UM molecularprognosissubtypes. D3-UM tumorswith EIF1AX versus SF3B1 mutations, which are known to be associated with low and intermediate risk of developing UM metastasis, respectively, correlated with distinct DNA methylation and SCNA profiles. D3-UM tumors also separated into two groups by transcription (mRNA, lncRNA, and miRNA) profile analysis. Loss of chromosome 3, followed by BAP1 alteration, results in bilallelic BAP1 loss. M3/BAP1 aberrancy is associated with a global DNA methylation profile that is not observed in D3-UM. Despite all M3/BAP1-aberrant UM sharing this common DNA methylation pattern, these tumors divide into two groups by SCNA and transcription profiles, with distinct pathway features indicative of hypoxia, DDR, MYC/MAX signaling, and proliferation. See also Figures S6 and S7; Tables S5 and S6.