BAP1 Loss Is Associated with DNA Methylomic Repatterning in Highly Aggressive Class 2 Uveal Melanomas

Uveal melanoma (UM) is the most common primary cancer of the eye, leading to fatal hematogenous metastasis in up to half of patients (1). In recent years, considerable progress has been achieved in elucidating the molecular landscape of UM (2). Primary UMs can be divided into two prognostically significant categories based on their gene expression profile (GEP), which remains the most accurate prognostic method for this cancer (3–7). UMs with the Class 1 GEP have low metastatic risk, whereas those with the Class 2 GEP have high metastatic risk. To provide a standardized instrument for classifying UMs based on GEP, a 12-gene prognostic test was developed and validated for routine clinical use in a Collaborative Ocular Oncology Group (COOG) prospective study of 459 patients (8, 9). The transcriptome of Class 1 UMs corresponds to that of differentiated melanocytes and Class 2 UMs to that of neural crest–like progenitors from which uveal melanocytes arise (10, 11), suggesting that aberrations may accrue during progression to Class 2 UM that disrupt the transcriptomic program maintaining cells in a differentiated melanocytic cell fate. We previously found that Class 2 UMs harbor loss-of-function mutations in BAP1, located on chromosome 3p21, and that the other copy of BAP1 is lost as a result of monosomy 3 or isodisomy 3 (12), thereby fulfilling Knudson "two hit" hypothesis for tumor suppressors (13). The discovery of BAP1 mutations finally explained the association between monosomy 3 and poor prognosis in UM that had been known for almost three decades with no mechanistic explanation (14). However, in contrast to other BAP1-associated cancers such as renal cell carcinoma and mesothelioma, where BAP1 is frequently eliminated by regional deletions around the gene locus (15, 16), Class 2 UMs demonstrate complete loss of an entire copy of chromosome 3 (14, 17). Indeed, partial deletion of chromosome 3 is associated with low-risk Class 1 UMs (17).

A critical unanswered question is why complete loss of chromosome 3 is required in UM for the Class 2 GEP and high metastatic risk. One possibility is that haploinsufficiency for multiple genes across the chromosome may be required for malignant progression, but this is unlikely because isodisomy 3 (loss of one copy of chromosome 3 and duplication of the remaining BAP1 mutation–bearing copy of chromosome 3) confers the same Class 2 GEP and metastatic risk as monosomy 3 (17, 18). A second possibility is that other tumor suppressor genes on chromosome 3 are also mutated in Class 2 tumors, and that these must also be reduced to homozygosity in order for the Class 2 phenotype to emerge (19). However, intensive investigations in using next-generation sequencing have failed to identify any commonly mutated genes on chromosome 3 other than BAP1 (12, 17, 20–22). A third possibility is that the retained copy of chromosome 3 in Class 2 tumors contains not only a mutant copy of BAP1, but also epigenetic alterations that are required for the Class 2 phenotype. If this model is correct, it raises the question of whether BAP1 loss might itself contribute to this epigenetic reprogramming. Here, we used an integrated transcriptomic, methylomic, and functional approach to address this question. Our findings provide evidence in support of the third model. BAP1-mutant Class 2 UMs exhibit extensive nonrandom alterations in DNA methylation, many of which are located on chromosome 3 and some of which may be specifically triggered by loss of BAP1. This study identifies novel prognostic biomarkers and provides functional evidence linking BAP1 loss to specific methylomic alterations, suggesting new therapeutic strategies in BAP1-mutant tumors.

This study was approved by the Institutional Review Board and adhered to the tenets of the Declaration of Helsinki. Written informed consent was obtained from each patient. Primary UM samples were obtained at enucleation, snap-frozen, and stored at −80°C. DNA was extracted using QIAamp DNA Mini Kit (Qiagen). DNA methylation analysis was performed on the Infinium Human Methylation 450K BeadChip Kit (Illumina, Inc.). In addition, raw data files from the 80 TCGA UM samples, which were also analyzed on the Infinium Human Methylation 450K BeadChip system (22), were obtained from the Cancer Genomics Hub (CGHub). Methylation data underwent quality control, normalization, batch identification with singular value decomposition and correction with Empirical Bayes frameworks, and analysis using the ChAMP methylation pipeline in R (23). Unsupervised principal component analysis (PCA) and 3D visualization were performed on the top 20% most variably methylated probes using the stats and rgl packages in R and ellipsoids were plotted at 95% confidence intervals (CI). RNA-Sequencing data from TCGA was analyzed as described previously (24) and then assessed for differential expression using EdgeR (25) and DESeq2 (26). Significantly differentially expressed genes found with both EdgeR and DESeq2 were used for downstream analyses. Plots of the number of significant probe sites in CpG-island feature, chromosome, or promoter regions were normalized to the overall number of probes present within each of those regions contained in the Infinium Human Methylation 450K BeadChip Kit and statistical significance was assessed using binomial testing. Genomic locations are expressed using the hg19/GRCh37 assembly.

Gene expression and associated methylation probe sites that were both significantly different (FDR < 0.05) between Class 1 and Class 2 tumors from the TCGA dataset consisting of 80 UM samples were compared using custom scripts in R. For genes with multiple significantly methylated probe sites, downstream analyses were conducted on the basis of the average and sum of all significantly methylated probes, as well as sequential and reverse sequential numeric ordering of the Delta Beta (Class 2 Beta – Class 1 Beta) methylation values with removal of additional probe sites for each gene. Primary findings of downstream analyses using each strategy showed similar results and reverse sequential numeric ordering of the Delta Beta values with removal of additional probe sites for each gene was selected for determination of gene quadrant location. Results were plotted in a quadrant graph comparing methylation and gene expression using ggplot2 in R. Quadrants that showed significant (i) hypermethylation and downregulated gene expression and (ii) hypomethylation and upregulated gene expression were selected for subsequent analyses; further filtering included an absolute value (log2(RNA fold change)) > 1 and absolute value (methylation Delta Beta value) > 0.05. RNA fold change was defined as the Class 2 RNA normalized count divided by the Class 1 RNA normalized count with normalization conducted in DeSeq2 as described previously (26). Circular plots of DNA methylation from probes within those quadrants were generated using the Perl-based Circos graphical program (27). Methylation location with respect to CpG sites and expression data of specific chromosomal regions were plotted using ggpubr, ggplot2, and GViz in R (28). Predictions of transcription factor–binding sites within DNA regions were conducted using PROMO (29).

Differentially hypermethylated/downregulated genes and hypomethylated/upregulated genes were input into gene set enrichment analysis (GSEA; ref. 30), which uses hypergeometric distribution to look for significant enrichment (FDR < 0.05) within annotated gene sets in the Molecular Signatures Database (MSigDB), including chromosomal location, transcription factor target binding, and KEGG and Reactome functional pathway gene sets. Protein–protein interaction networks of significantly enriched functional pathway gene sets were generated using STRING (31) and modified in Adobe Illustrator (Adobe Systems Inc.).

For validation of selected sites that were highly differentially methylated in Class 2 tumors, we filtered by hypermethylated CpG sites (FDR < 0.05) that correlated with significantly downregulated gene expression (FDR < 0.05) in Class 2 tumors compared with Class 1. In addition, to filter out differentially methylated sites arising from infiltrating immune cells, we obtained a publicly available dataset of 6 whole-blood samples analyzed on the Infinium Human Methylation 450K BeadChip Kit (Illumina, Inc.; ref. 32) and analyzed this dataset as described previously. Class 1 and Class 2 samples were compared with the whole-blood methylation data, and probes were discarded if there was hypermethylation in the whole-blood cell data compared with Class 1 or Class 2 samples. Next, to control for whether chromosomal aberrations may have biased the findings, probes were selected from genes in chromosomal regions frequently gained or lost in UM tumors (chr 1, 3, 6, and 8). Further probe selection was based on identifying probes from genes within those chromosomes that showed the greatest difference in methylation between Class 2 UMs when compared with both Class 1 UMs and WBCs and/or finding genes with multiple top hit probes that could be targeted with the same primer pair. The number of probes chosen for validation was limited because of the available DNA from fresh tumor samples. The genes that were selected, their chromosomal locations, the CpG probes targeted, and the forward and reverse primers targeting the sites were designed using Zymo Bisulfite Primer Seeker (Zymo Research) for each region of interest and were optimized to ensure specificity to bisulfite-converted DNA.

Genomic DNA was extracted using the Wizard Genomic DNA Purification Kit (Promega) according to the manufacturer instructions from 14 primary UM tumors (five Class 1 and nine Class 2). Two micrograms of DNA from each sample was bisulfite-converted using EZ DNA Methylation-Lightning Kit (Zymo Research) according to the manufacturer's instructions and subsequently PCR-amplified with the appropriate primers (Supplementary Table S1) using Epimark Hot Start Taq DNA Polymerase (New England BioLabs). PCR products were purified by agarose gel separation and extracted with the QIAquick Gel Extraction Kit (Qiagen). Sanger sequencing was conducted for each respective sample and primer site. CLC Sequence Viewer 7 (Qiagen) was used to align sequences, and SnapGene Viewer (GSL Biotech LLC) was used to visualize the sequencing traces and make methylation calls. CpG methylation probe site calls within the sequencing region of each primer pair were compiled independently by two investigators who were masked to all clinical information. In most methylation probe locations, either an unmethylated (T) or methylated (C) peak clearly predominated. At sites with approximately 50% methylation, the call was made as a Y (representing the cooccurrence of both methylation and nonmethylation). No discordance in calls between investigators occurred. After all methylation calls were made, the results were used to predict molecular prognostic class assignments (Class 1 or Class 2), which was independently determined using a prospectively validated 15-gene expression profile prognostic test performed in a CAP-accredited, CLIA-certified laboratory (33).

shRNA vectors targeting CGTCCGTGATTGATGATGATA, CCCTGTATATGGATTTATCTT, and CCACAACTACGATGAGTTCAT of human BAP1 cDNA (shBAP1) were individually cloned into a pLKO-TET-puro vector (Addgene #21915) and packaged into lentiviral particles by transient transfection into H293T cells, as described previously (34). UM cell lines Mel202 and 92.1 [kindly provided by Drs. B. Ksander and M. Jager (Harvard Medical School, Boston, Massachusetts), respectively] were then transduced with shBAP1-containing lentiviral particles and selected with puromycin (2 μg/mL) 48 hours after initial transduction. Cells were clonally selected for optimal doxycycline-inducible BAP1 knockdown (>85% knockdown) and expanded for use in this study. Cells were grown to <80% confluency and induced with 1 μg/mL of doxycycline for 5 days. After 5 days of doxycycline induction, BAP1 knockdown (BAP1KD) was confirmed by Western blot using an anti-BAP1 primary antibody (Santa Cruz Biotechnology, H-300), anti–β-actin loading control antibody (Cell Signaling Technology, 4967S), IRDye and VRDye secondary antibodies (Li-COR Biosciences), and a LI-COR Odyssey CLx imaging system. Genomic DNA was isolated using QIAamp DNA Mini Kit (Qiagen) and RNA was isolated using TRIzol and the RNeasy Mini RNA isolation kit with RNAase-Free DNase treatment (Qiagen). Uninduced cells were used as a control.

Genomic DNA from control and doxycycline-induced BAP1KD Mel202 and 92.1 UM cells were analyzed using the Infinium Methylation EPIC BeadChip Kit (Illumina, Inc.), with two biological replicates for each condition. Methylation data underwent quality control, normalization, batch correction, and analysis using the ChAMP methylation pipeline in R (23). Significantly hypermethylated or hypomethylated probes (FDR < 0.05) were identified with cell lines pooled together when comparing BAP1KD to control cells. RNA-seq libraries were prepared using the TruSeq Stranded Total RNA prep kit with Ribo-Zero Gold to remove cytoplasmic and mitochondrial rRNA according to the manufacturer's recommendation (Illumina, Inc.). tRNA-seq libraries were run on an Illumina NextSeq 500 sequencing instrument according to the protocols described by the manufacturer (Illumina, Inc.). Raw RNA-Seq FASTQ files were assessed for quality using FastQC and aligned to the genome with count files generated using STAR (35). Read counts were normalized and batch-corrected then assessed for differences in expression (P < 0.05) between groups using EdgeR (25). After library prep and quality control, two biological replicates were obtained for each condition from 92.1 cells and one each for Mel202 cells and the data were pooled together. Similar to the TCGA data analysis, hypermethylated/downregulated and hypomethylated/upregulated genes (P < 0.05) were selected for subsequent analyses; further filtering included an absolute value (log₂(RNA fold change)) > 1 and absolute value (methylation Delta Beta value) > 0.01. RNA fold change was defined as the doxycycline-induced BAP1KD RNA normalized count divided by the uninduced RNA normalized count. Genes with significant hypermethylation/downregulation or hypomethylation/upregulation from the BAP1KD cells (P < 0.05) were compared with hypermethylated/downregulated and hypomethylated/upregulated genes from the TCGA dataset (FDR < 0.05). Overlap of individual probes between datasets was also assessed. Chromosomal location and pathway analyses were conducted using GSEA and MSigDB as described previously.

All methylation and RNA-Seq data generated in this study have been deposited in Gene Expression Omnibus under accession code GSE130357.

Primary UMs can be divided on the basis of DNA methylomic profiling into two groups that correspond to BAP1 mutation status and clinical GEP classification (Class 1 vs. Class 2; refs. 17, 22). To elucidate methylation changes associated with BAP1 mutations, we performed unsupervised PCA of the top 20% most differentially methylated probes on two datasets, including the 80 TCGA samples and 12 independent samples from our center, both analyzed on the Infinium Human Methylation 450K BeadChip array (Supplementary Table S2). As expected, tumor samples from both datasets clustered into two distinct groups corresponding to Class 1/BAP1–wild-type and Class 2/BAP1–mutant subtypes (Fig. 1A and B). Of the 20% most variable probes, 82% were shared between datasets, and out of the top 5,000 probes making the greatest contribution to the first principal component (PC1), 204 probes were shared between datasets (Supplementary Table S2). Next, we sought to obtain functional insights into the methylation changes. Hypermethylated probes in Class 2 UMs relative to Class 1 UMs (FDR < 0.05) were most enriched within “shore” regions 1,000–1,500 bp upstream of transcription start sites (TSS-1500), followed by “shore” regions in the 5′UTR and "shelf" regions within the TSS-1500 (Fig. 1C; Supplementary Table S2). Notably, methylation in these “shore” regions correlated inversely with gene expression (36). Hypomethylated probes in Class 2 UMs relative to Class 1 UMs (FDR < 0.05) were most enriched in open sea regions. Hyper- and hypomethylated probes (FDR < 0.05) within promoter regions were significantly enriched on chromosomes 3 and 8 (Fig. 1D), which commonly display copy number alterations in UM (17). In particular, chromosome 3, which is reduced to homozygosity in most Class 2 tumors, demonstrated the most significant enrichment for promoter hypermethylation in Class 2 tumors (Fig. 1D). These findings suggest that the sole copy of chromosome 3 in Class 2/BAP1–mutant UM harbors extensive epigenetic alterations that may play a role in tumor progression.

To further explore the potential functional relevance of these findings, we then focused on integrating methylation with gene expression. As hypermethylation in promoter regions is usually associated with gene silencing (37), we identified genes with hypermethylation associated with decreased gene expression (hypermethylated/downregulated) or hypomethylation associated with increased gene expression (hypomethylated/upregulated) in Class 2 UMs relative to Class 1 UMs (Fig. 2A). This analysis revealed 1,621 hypermethylated probes associated with 508 downregulated genes (FDR < 0.05) and 3,876 hypomethylated probes associated with 923 upregulated genes (FDR < 0.05; Supplementary Table S3). Out of the 12 genes in the GEP test, six of the eight downregulated genes were hypermethylated (FXR1, ID2, ROBO1, LMCD1, SATB1, and MTUS1) and all four of the upregulated genes were hypomethylated (HTR2B, ECM1, RAB31, CDH1). Similar to our global analysis, hypermethylated/downregulated genes in Class 2 tumors were enriched within the promoter and 5′UTR “shore” and “shelf” regions (Fig. 2B), and hypomethylated/upregulated genes in Class 2 tumors occurred mostly in open sea regions (38). Chromosomal regions that were significantly enriched (FDR < 0.05) for hypermethylated/downregulated genes included 6p21, 6p24, 19q13, 10q24, 4p14, as well as multiple regions on chromosome 3 (3p21-23, 3p25-26, 3q12-21, and 3q27). Regions enriched (FDR < 0.05) for hypomethylated/upregulated genes included 6p21, 17q21, 15q21, and 12p13, 22q13, 1q31, and several regions on chromosome 8 (8p21, 8q13, 8q21-22; Fig. 2A–C; Supplementary Table S4). Interestingly, 6p21 contains regions of both hypomethylated/upregulated and hypermethylated/downregulated genes in Class 2 tumors, but only the former is enriched for HLA genes, which are known to be expressed as part of the altered immune microenvironment in Class 2 tumors.

Both hypermethylated/downregulated and hypomethylated/upregulated regions were enriched for genes with functions related to developmental biology and tissue development, focal adhesion, immune function and axon guidance (Supplementary Table S4), and they were enriched for regulatory binding sites for the transcription factors TCF3, TFAP4, FOXO4, and LEF1 (Supplementary Table S4), which are associated with cell fate, differentiation, stem cell maintenance, tumor growth, and metastasis (39, 40). Of note, LEF1 and TCF3 are important for melanocyte lineage commitment through their interaction with the MITF promoter, and they are involved in WNT signaling and axon guidance (41, 42). Because melanogenesis was one of the top pathways for hypermethylated/downregulated genes, we inspected MITF, which did not meet our strict RNA-Seq logFC cutoff, although it did exhibit significant hypermethylation and decreased gene expression (FDR < 0.05 for both) in Class 2 UMs relative Class 1 UMs. Furthermore, axon guidance was one of the most significant pathways for both hypermethylated/downregulated and hypomethylated/upregulated genes, many of which are involved in migration of neural crest from which melanocytes arise (Fig. 3; refs. 43–45). The proteins encoded by these axon guidance and melanogenesis genes demonstrated extensive interactome connectivity (Fig. 3) and share many of the same transcription factor–binding sites, suggesting that they may be deregulated in a functionally significant manner.

Plotting all hypermethylated/downregulated and hypomethylated/upregulated probes along the length of chromosome 3 demonstrates the extensive hypermethylation/downregulation occurring in Class 2 UM, with the most significant probes and highest methylation fold changes occurring within the previously identified chromosomal loci (3p26-25, 3p23-21, 3q12-21, and 3q27) that are enriched (FDR < 0.05) in hypermethylated/downregulated genes (Fig. 4). Of particular interest, these regions on chromosome 3 contained many of the axon guidance and melanogenesis genes described above, including DVL3, RAF1, MITF, SATB1, PLXNB1, CHL1, ROBO1, and SEMA3B (Fig. 4). Hypermethylated probes were present in the promoter region of all of these genes with the exception of ROBO1, which had several hypermethylated sites located within the gene body. Several of these genes (MITF, ROBO1, and SEMA3B) had both hypermethylated and hypomethylated sites distributed throughout gene body exons, suggesting that epigenetic regulation of these genes is complex and may involve splicing, alternative exon usage, or other regulatory mechanisms. Consistent with previous work (12, 46), no significant differences in methylation of the BAP1 promoter region were found between Class 2 and Class 1 UMs or between monosomy 3 and disomy 3 UMs. Interestingly, however, BAP1 is located within one of the hypermethylated/downregulated loci on chromosome 3 (3p21) and contained two significantly hypermethylated CpG probe sites (FDR < 0.05), one of which was located within the gene body (cg16871520) and the other in the 3′UTR (cg21746711; Fig. 5). The cg16871520 probe demonstrated the most significant inverse Spearman correlation coefficient between BAP1 gene expression and methylation (R = −0.79; P < 0.001). This single probe accurately distinguished Class 1 UMs from Class 2 UMs in 79 of the 80 samples (Fig. 5).

Several of the most significantly hypermethylated probes from the hypermethylated/downregulated genes in Class 2 UMs were selected for orthogonal validation using the bisulfite conversion method, including IL12RB2 (chr1, 6 probes), SATB1 (chr3, 1 probe), SESN1 (chr6, 1 probe), and ENPP2 (chr8, 2 probes; Supplementary Table S1). CpG methylation was assessed for sites within the sequencing region of each primer pair in an independent dataset of 14 primary UM samples (distinct from the 12 independent UM samples used for methylation array analysis) including five Class 1 and nine Class 2 tumors. Class 1 tumors were readily distinguished from the Class 2 tumors by the methylation status of these probes (Fig. 6). The probes that most accurately differentiated Class 1 UMs from Class 2 UMs were the IL12RB2 probes and the SESN1 probe.

To determine whether methylomic repatterning in Class 2 UM is directly attributable, at least in part, to BAP1 loss, we engineered UM cell lines Mel202 and 92.1 to allow efficient BAP1KD using a Tet-inducible shBAP1 system (Supplementary Fig. S1A). Following BAP1KD, the BAP1KD 92.1 and Mel202 cells clustered together and away from their respective control cells (Supplementary Fig. S1B and S1C) with 11,023 hypermethylated and 3,864 hypomethylated probes (FDR < 0.05) occurring within 4,481 and 2,153 genes, respectively (Supplementary Table S5). Consistent with findings in Class 2 UMs, development, axon guidance, and melanogenesis pathways were significantly enriched (FDR < 0.05) for genes with hypermethylated probes (FDR < 0.05) in the BAP1KD UM cell lines, which also showed the most significant hypermethylated loci enrichment for chromosomal band 3p21, where BAP1 is located, as well as another locus on chromosome 3 at 3q21 (FDR < 0.05, Supplementary Table S5). Similar to the Class 2 UMs, immune function, focal adhesion, and axon guidance pathways were significantly enriched for hypomethylated genes (FDR < 0.05).

After integrating methylation data with matched RNA-Seq data, there were 417 significantly hypermethylated/downregulated genes and 110 hypomethylated/upregulated genes (P < 0.05) in the BAP1KD cells (Supplementary Table S6). Again, development, axon guidance, and melanogenesis were among the pathways most significantly enriched for hypermethylated/downregulated genes in BAP1KD UM cell lines (FDR < 0.05), whereas loci on chromosome 3 were no longer significantly enriched (Supplementary Table S6). Hypomethylated/upregulated genes from the BAP1KD cells were significantly enriched for pathways involved in apoptosis and cell death (FDR < 0.05). Of the 12 genes from the GEP test, ROBO1 and MTUS1 became hypermethylated after BAP1KD, but none of the genes showed corresponding changes in methylation and gene expression.

Next, we compared BAP1KD cells to Class 2 UMs from the TCGA dataset, and we found 31 hypermethylated/downregulated and 9 hypomethylated/upregulated genes shared in common. Four of the 31 genes were in the development (ACVR2B, ADCY6, DPYSL2, ROBO3), two in the axon guidance (DPYSL2, ROBO3) and one in the melanogenesis (ADCY6) pathway. For individual probes, only two exact probe sites (one on PC and one on PALMD) were shared for the hypermethylated/downregulated genes and one (ZEB1) for the hypomethylated/upregulated genes. Interestingly, a recent paper has correlated ZEB1 upregulation with metastatic progression of UM, showing functionally that ZEB1 promotes UM dedifferentiation, depigmentation, and increased invasiveness (47). Overall, these findings confirm that loss of BAP1 leads to DNA methylomic repatterning, as recently described in liver cancer (48), with the most extensive changes in UM occurring in the axon guidance and melanogenesis pathways.

Here, we found that Class 2/BAP1–mutant UMs share a common DNA methylation profile reflecting extensive methylomic repatterning compared with Class 1/BAP1–wild-type UMs. The stereotyped nature of the methylomic profile and its strong association with BAP1 loss suggest that it may represent a latent epigenetic program regulated by BAP1 within the melanocyte/neural crest lineage. Consistent with this possibility, melanogenesis was one of the most significantly enriched functional pathways among hypermethylated/downregulated genes in Class 2 UMs, a finding concordant with the impairment of melanocyte differentiation in this subclass. Among the downregulated melanogenesis genes, EDNRB mediates commitment of neural crest cells to the melanocyte lineage, MITF is the master regulator of melanocyte differentiation, and DCT is an enzyme in the melanin synthesis pathway and is one of the earliest markers specific to melanocyte differentiation (49–51). In addition, the axon guidance pathway was a top hit for both hypermethylated/downregulated and hypomethylated/upregulated genes, demonstrating an extensive protein–protein interactome between all members (Fig. 3). This connectivity, even among inversely correlated genes, is a manifestation of the complex interaction of both attractive and repulsive cues coordinating migration and differentiation in the neural crest lineage. For example, the hypermethylated/downregulated axon guidance genes ROBO1, PLXNB1, SEMA3B, and CHL1 are involved in neuronal/neural crest migration but have been implicated as tumor suppressors (45, 52–54), whereas the hypomethylated/upregulated SEMA3E is an axonal path finding gene that increases tumor invasiveness and metastatic spread (55). These findings are consistent with recent work showing that BAP1 regulates neural crest migration and melanogenesis during vertebrate development (Dawn Owens, PhD, personal communication) and could explain why Class 1/BAP1–wild-type UMs demonstrate characteristics of differentiated melanocytes, whereas Class 2/BAP1–mutant UMs resemble neural crest–like progenitor cells (10, 11).

In addition, the most significant and densely clustered hypermethylated/downregulated gene loci in Class 2 UMs were located on chromosome 3, which contained many of the axon guidance cue, neural crest specification, and melanocyte differentiation genes (e.g., ROBO1, PLXNB1, SEMA3B, CHL1, SATB1, MITF, DVL3, and RAF1). Because all of these genes undergo repressive methylation changes on the sole remaining copy of chromosome 3, this could explain why the other copy of chromosome 3 must be lost to acquire the metastasizing Class 2 UM phenotype, as has been hypothesized previously (56, 57). The cg16871520 probe, which detects methylation in the BAP1 gene body, demonstrated a highly significant association between increased methylation and Class 2 GEP. As such, methylation sites in this region could potentially represent valuable new biomarkers. Because all but three Class 2 tumors in this analysis had detectable inactivating BAP1 mutations, it is not clear how hypermethylation of the gene would be functionally relevant to conversion to Class 2. BAP1 may negatively regulate its own methylation, as it does with the neural crest and melanogenesis genes described above, such that loss of BAP1 could lead to hypermethylation of the remaining allele as a “passenger event” in the absence of selective pressure, although this was not seen in the BAP1KD cells. Furthermore, experimental BAP1KD did show an enrichment for hypermethylation of axon guidance and melanogenesis genes, and recapitulated hypermethylation of specific loci on chromosome 3, including 3p21 where BAP1 is located, as observed in primary UM tumor samples. This is consistent with a recent large study that identified BAP1 as one of several genes in which mutations were associated with global methylomic alterations (48). When pairing methylation and RNA-Seq datasets, the hypermethylated/downregulated genes in the BAP1KD cells maintained enrichment for axon guidance and melanogenesis pathways; however, enrichment for loci on chromosome 3 was not maintained. In terms of the GEP, 10 of the 12 genes showed corresponding methylation and gene expression changes in the TCGA dataset, implicating methylation in regulating their gene expression, but this was not recapitulated in the BAP1KD cells.

Of note, we would not expect the depletion of BAP1 in 92.1 or Mel202 cells to precisely recapitulate the methylation or GEP findings in actual human Class 2 UMs. First, it is well known that cultured cells do not exactly match the tumors they were derived from. Second, the 92.1 and Mel202 cells are heterozygous for chromosome 3 and have different genetic profiles than Class 2 UMs. Accordingly, others have suggested that different populations of normal uveal melanocyte precursor cells exist in the eye with differing methylation profiles and that the specific methylation profile present may make those cells more susceptible to development of a specific UM subtype (56). Finally, knocking down BAP1 is not equivalent to knocking down expression of the entire chromosome 3, which occurs in Class 2 UMs and is required for the Class 2 GEP. Thus, we would anticipate that the differentially methylated probes shared between the BAP1KD and primary UM samples would be enriched for those specifically related to BAP1 loss but not necessarily for those related to chromosome 3 loss, tumor microenvironmental influences, and other factors. Consequently, and keeping in mind the well-known differences between tumors and derivative cell lines, we would not expect an exact match of probes between these sources. Nevertheless, it is remarkable that axon guidance and melanogenesis are enriched pathways associated with BAP1 loss in both tumor samples and cell lines, and that some of the same genes overlap between systems (i.e., ADCY6, ROBO3, DPYSL2). Therefore, taken in context with previous work showing the evolution of genomic aberrations in UM (17, 22), we hypothesize that the initial event presaging the divergence of Class 2 UM from Class 1 UM is loss of one copy of chromosome 3, which is followed by mutation of BAP1 on the other copy of chromosome 3, which then leads to the methylomic repatterning characteristic of Class 2 UMs. Future studies profiling the methylome of normal uveal melanocytes may provide insights into how preexisting epigenetic states may drive malignant transformation along the observed evolutionary trajectories.

Finally, we confirmed that methylation probes from the IL12RB2, SATB1, and SESN1 genes accurately differentiated Class 1 from Class 2 UMs as an orthogonal validation and as a proof-of-concept for developing a methyl DNA–based liquid biopsy assay in the future. For this purpose, we filtered the probe list for hypermethylated sites found in whole-blood samples to eliminate potential contamination from circulating immune cells. While our findings suggest that such an approach is possible and could reduce the need for invasive biopsies, several technical challenges need to be overcome that are beyond the scope of this study.

J.W. Harbour is an inventor of intellectual property and patents related to uveal melanoma. He is a paid consultant for Castle Biosciences, licensee of this intellectual property, and he receives royalties from its commercialization. He is also a consultant/advisory board member for Aura Biosciences and Immunocore. No potential conflicts of interest were disclosed by the other authors.

Conception and design: M.G. Field, L.Z. Cai, J.W. Harbour

Development of methodology: M.G. Field, J.N. Kuznetsov, L.Z. Cai, S. Kurtenbach, J.W. Harbour

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M.G. Field, P.L. Bussies, L.Z. Cai, C.L. Decatur, J.W. Harbour

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M.G. Field, P.L. Bussies, L.Z. Cai, K.A. Alawa, S. Kurtenbach, J.W. Harbour

Writing, review, and/or revision of the manuscript: M.G. Field, P.L. Bussies, K.A. Alawa, S. Kurtenbach, J.W. Harbour

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M.G. Field, J.N. Kuznetsov, C.L. Decatur, J.W. Harbour

Study supervision: J.W. Harbour

The authors are grateful to the patients who generously contributed samples for this research. The authors acknowledge the support of the Biostatistics & Bioinformatics and Oncogenomics Shared Resources at the Sylvester Comprehensive Cancer Center (Miami, FL), and the University of Miami Center for Computational Science (Coral Gables, FL). This work was supported by NCI grants R01 CA125970 (to J.W. Harbour) and F30 CA206430 (to M.G. Field), Alcon Research Institute (to J.W. Harbour), Research to Prevent Blindness, Inc. Senior Scientific Investigator Award (to J.W. Harbour), the University of Miami Sheila and David Fuente Graduate Program in Cancer Biology (to M.G. Field), the Center for Computational Science Fellowship (to M.G. Field), and a generous gift from Dr. Mark J. Daily (to J.W. Harbour). The Bascom Palmer Eye Institute also received funding from NIH Core Grant P30EY014801 and a Research to Prevent Blindness Unrestricted Grant.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.