Loss of BAP1 expression is associated with an immunosuppressive microenvironment in uveal melanoma, with implications for immunotherapy development.

Immunotherapy using immune checkpoint inhibitors (ICIs) induces durable responses in many metastatic cancers. Metastatic uveal melanoma (mUM), typically occurring in the liver, is one of the most refractory tumours to ICIs and has dismal outcomes. Monosomy 3 (M3), polysomy 8q, and BAP1 loss in primary uveal melanoma (pUM) are associated with poor prognoses. The presence of tumour-infiltrating lymphocytes (TILs) within pUM and surrounding mUM - and some evidence of clinical responses to adoptive TIL transfer - strongly suggests that UMs are indeed immunogenic despite their low mutational burden. The mechanisms that suppress TILs in pUM and mUM are unknown. We show that BAP1 loss is correlated with upregulation of several genes associated with suppressive immune responses, some of which build an immune suppressive axis, including HLA-DR, CD38, and CD74. Further, single-cell analysis of pUM by mass cytometry confirmed the expression of these and other markers revealing important functions of infiltrating immune cells in UM, most being regulatory CD8+ T lymphocytes and tumour-associated macrophages (TAMs). Transcriptomic analysis of hepatic mUM revealed similar immune profiles to pUM with BAP1 loss, including the expression of IDO1. At the protein level, we observed TAMs and TILs entrapped within peritumoural fibrotic areas surrounding mUM, with increased expression of IDO1, PD-L1, and β-catenin (CTNNB1), suggesting tumour-driven immune exclusion and hence the immunotherapy resistance. These findings aid the understanding of how the immune response is organised in BAP1 - mUM, which will further enable functional validation of detected biomarkers and the development of focused immunotherapeutic approaches.

Introduction

Uveal melanoma (UM) is the most common primary intraocular cancer in adults, accounting for 5% of all melanomas 1. Treatment options for primary UM (pUM) include radiotherapy and surgery 2, and usually achieve excellent local tumour control. Despite this, about 50% of UM patients develop metastatic disease, mainly in the liver 1. The average survival of patients with metastatic UM (mUM) is ∼12 months, as there are currently no proven effective treatments 3. Resection of isolated liver metastases may be attempted in selected cases, otherwise liver‐directed therapy (e.g. percutaneous perfusion with melphalan) or systemic chemotherapy. Recently, following the striking benefits in metastatic skin melanoma, immunotherapy using immune checkpoint inhibitors (ICIs) has been more widely used in cancer. However, in marked contrast to cutaneous melanoma, mUM is almost universally refractory to ICIs, mostly against CTLA‐4 and PD1/PDL‐1, with responses to single agents in the range of 3–8% 3.

While UMs have been partly ascribed to a low mutational burden 4, 5, evidence of specific TCR gene expression in tumour‐infiltrating lymphocytes (TILs) 6, promising responses to adoptive cell therapy using TILs 7, and encouraging results on targeting the melanocyte‐specific gp100 with the bispecific molecule tebentafusp (IMCgp100) 8 all suggest a specific immune response and that mutational burden is not the sole reason for the lack of response to ICI.

Monosomy 3 (M3) has long been known to be associated with increased risk of UM metastasis 9, 10, 11, and more recently, it has become apparent that this is primarily due to inactivating mutations of the BAP1 gene, which has been reported to be a stronger prognosticator than M3 12, 13. The Cancer Genome Atlas (TCGA) study of 80 pUMs demonstrated that patients with pUM at high metastatic risk [i.e. with UM characterised by M3 and loss of function of the tumour suppressor gene BAP1 (Chr 3p21.1)] could be further stratified, according to the presence of CD8+ T‐cell immune infiltrates and an altered transcriptional immune profile 4. The latter included elevated levels of HLA‐I molecules, which leads to natural killer (NK) cell suppression 14, TAM markers and expression of immune checkpoint regulators (ICRs), such as PD‐L1, indoleamine 2,3‐dioxygenase (IDO)‐1, and T‐cell Ig and ITIM domain (TIGIT) 4, 15.

Interestingly, previous work showed that loss of BAP1 in turn affects the expression of genes that impact the immune response 16. In this study, a comprehensive immune profiling of the 80 pUMs from the TCGA‐UM study revealed that several immune‐suppressive genes are significantly upregulated following BAP1 loss. We provide a novel and comprehensive understanding of UM immune evasion by profiling primary and metastatic UM at the transcriptomic and protein level using cutting‐edge approaches, including mass cytometry, NanoString, and digital spatial profiling of human patient tissues. Our findings suggest that UM cells, particularly those of BAP1‐negative (BAP1−) UM, shape the immune profile at both primary and metastatic sites, harnessing the expression of particular pathways and molecules to drive regulatory functions of myeloid cells and lymphocytes, and thus immunosuppression and immunotherapy resistance in advanced UM. These findings provide new insight for the functional validation of detected biomarkers for the further development of novel adjuvant immunotherapeutic approaches.

Materials and methods

Human subjects

This work was underpinned by the University of Liverpool (UoL) Ocular Oncology Biobank (OOB) and the Liverpool Bioinnovation Hub Biobank. Project specific approvals for work with pUM and mUM samples were obtained (REC‐18/LO/1027). Four fresh enucleated pUMs were included in this study for the CyTOF analyses.

TCGA analysis

mRNA expression and clinical data of The Cancer Genome Atlas (TCGA) GDC Ocular Melanomas dataset (UVM) were downloaded from the Xena Functional Genomics Explorer of University of California, Santa Cruz (https://xenabrowser.net/heatmap) 17. To provide understanding of the biological pathways involved in pUM pathogenesis via the expression of different immune genes, the nCounter PanCancer Immune Profiling gene set of 730 genes (NanoString Technologies, Seattle, WA, USA) was applied in the UCSC Cancer Genomics Browser to analyse the enrichment of immune genes sorted by BAP1 mRNA expression or chromosome 3 copy number variations. Generated data were extracted in comma‐separated values (CVS) format and analysed in GraphPad Prism 6 (GraphPad Software, Inc, San Diego, CA, USA) for correlation studies. Supervised clustering of immune genes of the TCGA RNA‐seq dataset was performed among those with significant Spearman's correlation to BAP1 expression or chromosome 3 copy number variation and sorted from the lowest rank (negative correlation) to the highest rank (positive correlation). The list of sorted genes was then uploaded in the Xena Browser for generation of heatmaps. Each of these genes was individually analysed as a prognosticator marker in Kaplan–Meier curves at the Xena Browser along the TCGA‐UM cohort. Those genes predicting significant survival differences (p < 0.05) were selected for further immune network analysis using the nCounter immune category list (NanoString Technologies) complemented by a custom‐built leukocyte functional immune response network collated by literature review (supplementary material, Table S1). Network plots were generated using the NodeXL Basic add‐in to Excel. In brief, immune genes were assigned to different immune categories in separated columns. In our analysis, we also considered the low‐ and high‐variance state of these genes along the TCGA‐UM cohort. This, in part, helped to define our hypothesis that the degree of variation in the expression of the genes associated with a particular network is indicative of the plasticity of that network 18. Therefore, high variance is associated with increased plasticity (higher thickness of network lines) and low variance with diminished plasticity (lower thickness of network lines) in response to BAP1 expression changes. We calculated the expression variance (σ2) of genes across the TCGA‐pUM cohort to predict how BAP1 loss impacts upon the expression of a particular gene by applying the following formula: σ2 = ∑(X − μ) 2/N, where X represents the RNA‐seq expression value of a particular gene, μ is the mean of the entire RNA expression for this particular gene in the cohort, and N is the distribution number (TCGA‐UM, N = 80). Therefore, the higher the effect of BAP1 on the gene expression, the higher the σ2 value of this particular gene in the cohort. In the network analysis, the highest variance value was limited to 5 units, assuming the CCL24 gene as reference for the highest variance (σ2 = 35.7). Sphere size represents the number of genes assigned to a given immune category (supplementary material, Tables S2 and S3). Box and whiskers analysis of specific genes according with different BAP1 expression levels was performed. RNA levels of BAP1 were defined as high (n = 27), mid (n = 26), or low (n = 27) according to a Kaplan–nMeier survival analysis of three groups in the TCGA‐UM cohort generated in the Xena Browser for BAP1 gene expression (p = 0.004843 and log‐rank test statistic = 10.66).

Immunohistochemistry

FFPE pUM and mUM samples were sectioned at 4 μm thickness and underwent antigen retrieval using the Dako pretreatment module (Agilent Technologies UK Ltd, Stockport, UK); slides were then incubated in a high‐pH bath containing Tris/EDTA buffer, pH 9.0 (Dako EnVision™ FLEX, Agilent) at 96°C for 20 min. IHC was performed using a Dako Autostainer PLUS machine, using the Dako Envision™ FLEX Kit (Agilent) according to the manufacturer's instructions. Slides were incubated with the following antibodies for 30 min: BAP1 (cat. No sc‐28 383/C‐4, dilution 1:200; Santa Cruz Biotechnology, Dallas, TX, USA), CD3 (cat. No IR503/polyclonal, ready to use; Dako Cytomation, CA, USA), CD4 (cat. No NCL‐L‐CD4/368, dilution 1:20; Leica Biosystems, Lincolnshire, IL, USA), CD8 (cat. No M7103/ C8/144B, dilution 1:200; Dako), CD163 (NCL‐L‐CD163/10D6, dilution 1:400; Leica Biosystems), and CD38 (NCL‐L‐CD38‐290/SPC32, dilution 1:100; Leica Biosystems).

The sections were counterstained with haematoxylin. Additional sections were treated with isotype controls at the same concentration as the primary antibodies.

Mass cytometry antibodies and reagents

All metal‐chelated optimised antibodies and reagents were purchased from Fluidigm (San Francisco, CA, USA). Full information for the antibodies and reagents used are provided in supplementary material, Table S4. The Maxpar Human Immune Monitoring Panel Kit was used as a reference antibody panel to immune profile primary uveal melanoma tumours, which includes the immune markers recommended by the Human ImmunoPhenotyping Consortium (HIPC) 19, with some modifications. The antibodies used cover the phenotype and functions of different subpanels of B cells, T cells, monocytes, dendritic cells, and NK cells. The MaxPar Panel Designer browser (Fluidigm) was used to predict and avoid metal spillover among tagged metals of the following additional markers included in the customised panel: CD74, LAG‐3, CD56, CD16, CTLA‐4, CD11b, and CD62L. The following markers were removed from the panel design in order to avoid spillover: CD194, TCRγδ, CD185, CD45RO, CD24, CD197, and CD20. Final spillover results were considered low between channels and this is shown in supplementary material, Figure S1A, top. The final wheel‐heatmap of the customised antibody panel shows the function that determines the best antibody–tag combinations to minimise background among channels, which ultimately contains targets with low tolerance of signal overlap (yellow‐green).

Mass cytometry of pUM

Four fresh histopathologically‐phenotyped BAP1− pUMs were manually minced prior to enzymatic digestion using collagenase A (cat. No C9722, 2 mg/ml; Sigma Aldrich, St Louis, MO, USA) and 40 units/ml DNase‐I (cat. No 79254; Qiagen, Germantown, MD, USA) in DMEM and incubated with agitation at 37°C for 60 min in a thermal mixer (Thermo Fisher, Waltham, MA, USA). Following incubation, digests were passed through a 70 μm filter to remove residual particulates. Cells were then pelleted (centrifugation at 1500 rpm for 5 min), washed in PBS, and viable cells were quantified using a Trypan Blue exclusion viability dye. Live cells were then washed twice with ice‐cold cell staining buffer (ic‐CSB; Fluidigm) and total cell concentration was determined using a Neubauer chamber. Up to three staining reactions of a maximum of 2.0 × 106 cells per sample were analysed. All samples were then incubated with 50 μl of 2% mouse serum in PBS with human TruStain FcX solution (Biolegend, San Diego, CA, USA) at 4 °C for 15 min. Samples were then processed for surface and intracellular staining with the panel described in supplementary material, Table S4 using the following protocol: 50 μl of a 2X surface antibody solution was made in ic‐CSB (final antibody dilution 1:100) and left on ice for 30 min. Cells were washed and fixed in 5 mm BS3 (Sigma) for 30 min, followed by fixation using 1x Fix‐I buffer according to the manufacturer's protocol (Fluidigm), and permeabilised in ice‐cold methanol for 10 min. Cells were washed and incubated with internal antibody cocktail (final dilution 1:100) for 30 min in ice. Then cells were washed and resuspended in intercalator‐Ir at 1:8000, and processed to be analysed using a Helios mass cytometer (Fluidigm).

Analysis of human tumour mass cytometry datasets

Data from mass cytometry were normalised to the EQ 4‐element bead signal using normalisation software version 2 (Fluidigm). Live Ir+CD45+ cells were manually gated as previously described 20 (supplementary material, Figure S1A, bottom), and FCS files were downloaded for concatenated analysis using Cytosplore V.2.2.1 for further downstream analysis by hierarchical stochastic neighbour embedding (HSNE) using a coefficient of 4 21, or individually processed for visualisation of t‐distributed stochastic neighbour embedding (viSNE) analysis in Cytobank. For accurate clustering and frequency calculations, a cut‐off of 1000 events was considered for the final gate. Eventually, Irhi CD45+ tumour infiltrated cells are detected among singlets and exhibit specific TAM markers, but not lymphocytic markers, excluding the possibility of doublets. These cells may often carry tumour‐derived DNA and melanin content given to phagocytosed tumour cells and are often observed in primary uveal melanoma tumours, classified as melanophages 22.

Mass cytometry and NanoString data (transcriptomic and DSP data) have been deposited at Flow CyTOF data, https://flowrepository.org/id/FR-FCM-Z2FD and GEO (gene datasets, GSE145782).

mRNA expression analysis using NanoString technology

For RNA immune gene expression analysis, four pUMs, six mUMs, and one normal liver (NL) formalin‐fixed, paraffin‐embedded (FFPE) samples were used. Only the tumour areas were selected for RNA extraction, or the entire normal liver tissue. The RNeasy FFPE Kit (Qiagen) was used for tissue dissociation, RNA extraction, and purification according to the manufacturer's instructions (Qiagen), as described in supplementary material, Supplementary materials and methods.

The NanoString nSolver 2.6 software was used for normalisation of expression counts using housekeeping genes following the manufacturer's recommendations 23. Data are displayed in expression count units of individual gene per patient compared with normal liver tissue, and internally normalised within each immune category (e.g. CTL suppression, M2 macrophage regulation, and immune checkpoint regulators). A full range of immune categories is displayed in supplementary material, Table S1.

Digital spatial profiling of mUM tissues

Digital spatial profiling analysis of one BAP1− mUM case was performed by NanoString's DSP technology platform to enable digital characterisation of protein distributed on the surface of FFPE tissue sections using the Human Immune Oncology panel (NanoString Technologies). In brief, 4‐ to 6‐μm‐thick FFPE mUM sections were stained for lymphocytes (CD3, red), macrophages (CD68, magenta), S100B (green), and DNA (blue) in order to detect the regions of interest (ROIs). The following workflow was used: deparaffinisation of FFPE unstained sections, antigen retrieval, antibody staining, ROI selection, DSP technology processing, nCounter analysis system. Data analysis and quality control were processed and normalised using positive and negative anti‐mouse and anti‐rabbit hybridisation control antibodies. S6 ribosomal protein and histone 3 were used as reference proteins. Area normalisation was applied between different ROI sizes varying from 100 to 650 μm in diameter. Results are displayed as absolute expression counts normalised with negative IgG controls.

Quantification and statistical analysis

All data were analysed using GraphPad Prism 6.0 and are presented as the means ± SD. Significant differences in the immune gene expression between the BAP1 lo, BAP1 mid, and BAP1 hi groups were estimated using one‐way analysis of variance followed by Bonferroni's multiple comparisons test. Survival analysis was performed in the Xena Browser using the Kaplan–Meier assay and was compared using the log‐rank test. The correlation between different mRNA expressions and overall survivals (OS) of TCGA‐UM patients was evaluated by non‐parametric Spearman's correlation, two‐tailed, where *0.01 < p < 0.05, **0.001 < p < 0.01; ***0.0001 < p < 0.001, and ****p < 0.0001 were considered to indicate significant differences.

Results

BAP1 loss significantly correlates with the modulation of immune genes and patient survival in UM

In the previous TCGA‐UM analysis, Robertson et al reported that the M3 phenotype in UM is associated with the upregulation of 30 immune genes 4, 24, and we therefore first sought to investigate whether this association could be the result of BAP1 loss. We found that an absence or reduced expression of BAP1 mRNA (cut‐off 19.54 for BAP1 expression and 2500 days as the default end‐point) also significantly correlated with decreased survival, similar to M3 status (cut‐off −0.3649 for Chr3 copy number and 2500 days as the default end‐point) (Figure 1A). Using Spearman's correlation analysis, we demonstrate that most of the investigated immune genes have a better expression correlation with BAP1 mRNA loss than with M3 (Figure 1B). In addition, only about 50% of these genes are significantly associated with patient survival, indicated by blue squares representative of Kaplan–Meier statistical test results (Figure 1B and supplementary material, Table S5).

Figure 1

BAP1 loss significantly correlates with the modulation of immune genes and patient survival in pUM. (A) Decreased mRNA expression of BAP1 is significantly correlated with a poor survival of primary UM patients similarly to the monosomy 3 (M3) status in the TCGA cohort. (B) Spearman's correlation analysis of specific immune genes compared to BAP1 expression and to chromosome 3 (Chr3) copy number variation. (C) Venn diagram depicting immune genes with significant correlation to Chr3 and BAP1 and with predictive survival significance. (D) Heatmap cluster analysis sorted by BAP1 expression showing upregulated and downregulated immune genes, including the P value profile of Kaplan–Meier survival scores.

We further expanded this analysis by interrogating the TCGA‐UM RNA‐seq data with a panel of 730 immune genes defined by the nCounter PanCancer immune panel (NanoString Technologies). One hundred and forty‐two immune genes exclusively correlated with BAP1 expression, but not with chromosome 3 copy number variation (supplementary material, Table S6); the expression of 117 genes negatively and another 25 genes positively correlated with BAP1 expression (r scores varying from −0.53 to 0.46). Among 181 immune genes that significantly correlated with both BAP1 expression and chromosome 3 copy number variation, 151 genes were negatively correlated and 30 genes were positively correlated, all with a higher correlation score to BAP1 expression than to chromosome 3 (r scores varying from −0.60 to 0.67) (supplementary material, Table S7). Among BAP1 correlated genes (n = 323), 168 genes that were negatively correlated with BAP1 expression were significantly associated with decreased survival, while 15 genes positively correlated with BAP1 expression were significantly associated with improved survival (Figure 1C and supplementary material, Table S8).

Independent of BAP1 expression, the M3‐UM genotype exclusively correlates with 43 immune genes and with 82 immune genes that also correlate with BAP1 to a lower degree (supplementary material, Tables S9 and S10), from which 48 are upregulated immune genes and 17 are downregulated immune genes, all significantly associated with patient survival (Figure 1C and supplementary material, Table S11). Supervised clustering analysis based on BAP1 mRNA expression is shown, including the P value profile of each gene related to survival outcome (Figure 1D). These findings suggest that loss of BAP1 expression is strongly associated with immune modulation of the microenvironment in pUM.

BAP1 loss correlates with immunosuppressive networks in pUM

We next analysed the group of immune genes upregulated following BAP1 loss, to predict the likely effects on the microenvironment. A scatter plot shows the gene expression and variance of 168 upregulated immune genes following BAP1 loss along the TCGA‐UM cohort, highlighting important immune genes involved in immune‐suppressive pathways, including LGALS3, CD74, CD38, PDCD1, IDO1, and HLA‐DR (Figure 2A). Differential expression analysis revealed that most of these immune genes are significantly upregulated in the TCGA‐UM cohort following BAP1 loss, as shown in the second quadrant of the volcano plot (Figure 2B and supplementary material, Table S8). PDCD1 and IDO1 have high significance with −log P values higher than 1.5 and log2FC lower than 0.2, and for that reason are not visible in the volcano plot.

Figure 2

BAP1 loss correlates with increased regulatory immune networks in primary uveal melanoma. (A) Gene expression profile of the TCGA‐UM cohort (n = 80) sorted from genes with the highest (left) to the lowest (right) gene variance expression. (B) Volcano plot depicting the most significantly upregulated immune genes with BAP1 loss (black arrows), which have potential immunosuppressive functions. (C) Immune network subcategory integrations with upregulated genes following BAP1 loss. The left panel shows general immune response networks and the right panel shows an expanded leukocyte effector immune response network. (D) Box and whiskers plots of selected upregulated and downregulated immune genes according to BAP1 expression levels (high, mid, and low). One‐way ANOVA was used for statistical analysis with Bonferroni's multiple comparisons test. ****p < 0.00001, ***p < 0.0001, **p < 0.001, *p < 0.05.

Importantly, all BAP1 significantly correlated immune genes simultaneously integrate different subcategories of the immune response, which were used to build an interactive transcriptomic network for visualisation of the predominant immune profile driven by pUM with BAP1 loss. Therefore, we performed a gene network analysis using two classification systems: a general immune response network based on major immune categories, and an amplified leukocyte functional immune network. In this analysis, the variance of gene expression (σ2) is represented by the thickness of the network lines, indicating the plasticity of the network toward BAP1 loss.

This analysis demonstrated that most of modulated immune genes correlate with leukocyte functions, and some of them with chemokine, interleukin, and cytokine expression; B‐cell functions; TLR and TNF superfamilies; and antigen presentation processes (Figure 2C). Within the leukocyte network, a dominance of immune‐suppressive pathways was observed, represented by red lines predominantly with higher expression variance (higher thickness), including Treg functions, Th1 suppression, Th2 activation, T‐cell tolerance responses, homing of Tregs, M2 macrophage functions, and ICRs (Figure 2C). Immune networks related to effective anti‐tumour immune responses are represented by green lines, predominantly with lower expression variance (lower thickness).

Importantly, some of the few immune genes that are downregulated following BAP1 loss (MICA, TNFSF13, and CD44) are important for the activation of anti‐tumour immune responses 25, 26, 27, as shown in box and whisker plots together with other immunosuppressive and exhaustion‐related genes (i.e. HLA‐DOB, CD74, CD38, LGALS3, IDO1, TIGIT, LAG3, and CD96) (Figure 2D). Although HLA‐DR has been classified as an immune response activation gene in disease given its importance in peptide presentation to CD4+ T cells 28, many regulatory functions have been attributed to HLA‐DR expression in the context of cancer 29. For that reason, in the network analysis, HLA‐DR was classified as immunosuppressive, although the generic immune activation classification was kept in Figure 2D. Immune genes that have a significant correlation with M3 status but do not correlate with BAP1 expression are probably regulated by different mechanisms that exclude BAP1 involvement. Among these genes, those related with immune response activation are downregulated, including IL12RB12, TLR1, and TLR5, and those involved with suppression of immune response are upregulated, including FN1, CD70, and CD73 (NTFE) (supplementary material, Figure S2A,B). All together, these findings show how different immune genes may integrate similar immune‐suppressive categories in high‐risk BAP1 − pUM, suggesting an importance in regulating the immune profile of mUM.

Transcriptomic analysis of mUM reveals a similar gene expression profile to BAP1 pUM

In order to compare the transcriptomic immune profiles of the primary and metastatic sites of UM, we performed a NanoString assay interrogating the expression profile of the nCounter PanCancer Immune Profiling panel using four pUMs and six mUMs, all lacking nuclear BAP1 expression. Unsupervised cluster analysis was performed revealing a high correlation between most pUM and mUM cases, with the exception of one mUM (mUM‐06). No significant correlation was observed between a normal liver control and tumour tissues (Figure 3A).

Figure 3

Transcriptomic analysis of BAP1‐negative metastatic uveal melanoma reveals similar immune profiles to BAP1‐negative primary tumours. (A) Heatmap of unsupervised clustering of all the samples [mUM (n = 6), pUM tissues (n = 4), and one normal liver] and all the transcripts (nCounter 730 immune genes panel). (B) Spearman's correlation analysis of gene expression between unmatched tumours from four pUMs and five mUMs (left), two matched tumours from one UM patient using Spearman's correlation rank (r = 0.92, p < 0.0001) (middle), and pUM/mUM NanoString data correlation analysis with TCGA‐UM Fpkm‐uq + 1 normalised RNA‐seq data (right). (C) Heatmap views of normalised RNA expression counts from six mUMs and one normal liver depicting the expression profile of CTL/NK suppression markers, immune checkpoint regulatory markers, and M2 macrophage regulation markers. Asterisks highlight selected highly expressed immune genes across the tissues.

Spearman's correlation test of the total gene expression counts revealed significant similarity between the gene expression of both primary and metastatic groups (r = 0.92, p < 0.0001) for BAP1 correlated immune genes (supplementary material, Table S8) (Figure 3B, left). A similar correlation score was also observed for one patient with matched primary and metastatic samples (Figure 3B, middle). Importantly, because the expression of samples of the same cancer generated using the same methodology and normalisation are often very similar, we also evaluated the correlation between our NanoString dataset (normalised counts) and the different samples of the TCGA‐UM cohort by applying the same genes. We observed that the normalised RNA‐seq data from the TCGA‐UM cohort are still highly correlative with the NanoString data, for both the primary and the metastatic tumours (Figure 3B, right).

Strikingly, when comparing the gene expression of six mUM patients with that of one human disease‐free liver normal biopsy (normal liver, NL), most mUM patients displayed upregulation of specific immune genes related to suppression of cytolytic T cells (CTLs) (Figure 3C), including HLA‐DRA, LGALS3, and CD38 partially 30, 31, 32, 33, 34, 35, 36, 37, 38; ICRs such as TIM‐3 (HAVCR2), HMGB1, IDO1, LAG3, and CD73 (NT5E) 39, 40, 41, 42, 43, 44; and TAM functional markers, such as ANXA1, CD74, CD9, INFAR2, MIF, PLA2G6, and CD163 36, 45, 46. In addition, transcript levels of NOS2, a predominant M1 macrophage marker 47, were similar to the levels found in normal liver without tumours. We also found high expression levels of CD74 and the macrophage migration inhibitory factor (MIF) across mUM tissues compared with normal liver. RNA levels of CD38 were increased compared with the levels found in normal liver. At the protein level, we found that CD38 was also positively expressed among regions of T‐cell infiltrates in one BAP1 mUM case, as evidenced by IHC staining (supplementary material, Figure S3H). B cells are nearly absent (low numbers of CD20+ cells) in mUM, and cytolytic effector cells seem to be at low activation states, given a paucity of TIA1+ cells among the TILs.

These findings suggest that the immune profile of BAP1 pUM is similar to mUM at the transcriptome level, suggesting an important role of CTL suppressive molecules, including HLA‐DRA and CD38, and TAM‐related pathways, where the CD74/MIF axis seems to play an important role driving the M2‐like phenotype and potential local tolerogenic responses. Therefore, some of these markers were further evaluated at the protein level of BAP1− pUM.

High‐resolution single‐cell analysis reveals regulatory T‐cell phenotype and a mixed macrophage phenotype within pUM with nuclear BAP loss

In order to phenotypically and functionally characterise pUM at the protein level, we performed a high‐resolution single‐cell analysis using mass cytometry in five pUMs, four of which showed nuclear BAP1 loss and one with normal nuclear BAP1 expression. Among infiltrating CD45+ cells, we observed a predominant cluster of macrophages, T lymphocytes (CD8+ and CD4+ T cells), B cells, and DCs, as evidenced in HSNE plots of all samples clustered together (Figure 4A,B). The frequency of each cell subtype was calculated for samples. The breakdown of the CD45+ infiltrating immune cells in BAP1− samples was as follows: CD4+ T cells (12%), CD8+ T cells (37%), DCs (10%), macrophages (30%), and B cells (11%) (Figure 4C, left). For one BAP1+ case, the breakdown of the CD45+ infiltrating immune cells was as follows: CD4+ T cells (14%), CD8+ T cells (59%), DCs (3%), macrophages (22%), and B cells (2%) (Figure 4C, right).

Figure 4

Mass cytometry analysis of infiltrated immune cells in primary UM. (A) Hierarchical stochastic neighbour embedding (HSNE) analysis showing the density of CD45+ infiltrated immune cells and selected phenotyping markers of concatenated pUM patients (BAP1−, n = 4 and BAP1+, n = 1). (B) Colour HSNE maps representing the phenotype of infiltrated immune cell subclusters. (C) Pie frequency charts of infiltrated immune cell subtypes detected by HSNE analysis for BAP1− and BAP1+ tumours. (D) Heatmap displaying normalised marker expression of each immune cell cluster for four concatenated BAP1− pUM samples. Analysis was generated in Cytosplore highlighting the most frequent clusters of CD45+ infiltrated immune cells, and an expanded analysis among tumour infiltrated monocytes, T lymphocytes, and B cells. (E) Pie chart showing the frequency of each cluster identified in HSNE analysis across the four merged BAP1− pUM tumours.

High‐dimensional HSNE single‐cell frequency clustering analysis and t‐SNE analysis (BAP1− cases) were performed among CD45+ cells for pUM cases in order to detect major clusters among the different cell subtypes (Figure 4D,E and supplementary material, Figure S1B–D). Across the T‐lymphocyte compartment, we observed high expression levels of CD28 receptor on both CD4+ (cluster A) and CD8+ T cells (clusters B and C), with low expression of CTLA‐4 and LAG‐3 checkpoint inhibitors (Figure 4D,E, and supplementary material, Figure S1D). Two out of four patients have low CD28 expression in CD8+ T‐cell clusters (supplementary material, Figure S1D), but no conclusions can be made, given the low number of tumours investigated. However, the frequency of CD28+CD8+ cells is still lower than CD28+CD4+ T cells in all four cases examined. CD4+ T cells mostly express CD25 and CD127 markers (cluster A, 5.1%), suggesting that they have a T‐regulatory phenotype 48, 49, 50, 51, 52. In addition, CD8+ T cells were positive for the proliferation marker Ki67 (clusters B–D), which had low expression among CD4+ T cells (cluster A) (Figure 4D). The high levels of Ki67 can be partially explained by the greater frequency of proliferating CD8+ T cells compared to CD4+ T cells. Cluster C (CD38+HLA‐DR+CD8+ T cells) is the most frequent CD8+ T‐cell cluster among BAP1− tumours (13.5%) and also showed increased expression of CD74 (Figure 4D,E).

In accordance with this suppressive phenotype, most CD8+CD28+ T cells (clusters B–D) were HLA‐DR+, a phenotype typical of regulatory CD8+ T cells 28. Importantly, CD8+ T‐cell clusters express high levels of CD38, recently reported to drive regulatory functions on CD8+ T cells 38, 53. Although B cells are not predominant in pUM, they could be divided into three subclusters: CD25+CD11b+ (cluster H), CD25lowCD11b+ (cluster I), and CD25−CD11b−CD74+ B cells (cluster J) (Figure 4D,E). Interestingly, cluster H showed increased expression of the LAG‐3 immune checkpoint regulator. In the macrophage and dendritic cell compartment (MOs and DCs), we found a mixed phenotype of M2‐like CD68+CD163+CD74+ macrophages (cluster E, 7.7%), M1‐like CD68+CD163−CD74−CD11c+CD11b+ macrophages (cluster F, 9.7%), and myeloid CD68−CD11b+CD11c+ dendritic cells (cluster G, 7.8%).

No significant differences were observed in the subclusters analysed in the pUM BAP1+ sample for the different immune cell subtypes compared with the BAP1− cases, apart from the CD8+ T‐cell compartment, which showed reduced levels of regulatory CD38+HLA‐DRhighCD8+ T‐cell cluster (cluster N,19.9%) compared with BAP1− tumours, and positive levels of the functional clusters CD38−HLA‐DRlowKi67+CD8+ T cells (cluster K 30.8%) and CD38−HLA‐DRlowKi67−CD8+ T cells (cluster M, 4.9%) and exhausted CTLA‐4+HLA‐DRlowCD8+ T cells (cluster L, 3.4%) (supplementary material, Figure S1B,C).

Taken together, we describe the regulatory nature of TILs in UM with nuclear BAP1 loss, particularly CD4+ and CD8+ T cells, and a mixed macrophage phenotype where M2‐like macrophages express higher levels of CD74.

Immune profile of mUM in regions of interaction between macrophages and lymphocytes

A digital spatial profiling assay (DSP, NanoString) revealed the protein expression profile of 31 immune markers in different regions of interest where macrophages (CD68) and lymphocytes (CD3) localised simultaneously in two mUM cases with nuclear BAP1 loss. Co‐localisation of macrophages and lymphocytes occurred both within and at the edge of the tumours (Figure 5A). Among the cancer‐related markers, we observed the expression of β2M, STAT3, STING, and β‐catenin (Figure 5B). The expression of β2M suggests that tumour antigens are presented via HLA‐A in these tumour areas, thus supporting the efficacy of ICI 54. Total levels of STAT3 were elevated, but not in its activated phosphorylated form (pY705), which regulates gene transcription for modulation of immunosuppressive factors 55. The expression of STING suggests a macrophage‐mediated hepatic inflammation and fibrogenic process 56.

Figure 5

Digital spatial profiling analysis of two mUM BAP1‐negative FFPE tissues using the NanoString immune oncology protein panel. (A) Regions of interest (ROI) to evaluate fibrotic areas with high infiltration of both macrophages (CD68) and lymphocytes (CD3). (B) Heatmap representation of different cancer‐related markers, immune phenotyping markers, and immune checkpoint and functional markers, all at the protein level among individual ROIs using normalised raw NanoString counts.

Importantly, high levels of β‐catenin were detected in both mUM patients, which is related to tumour‐induced immune exclusion mechanisms 57, 58, 59, suggesting an accessory mechanism by which tumours modulate infiltration and proliferation of lymphocytes in the metastatic site.

Among the immune phenotyping markers, we observed discrete, but positive, levels of CD163; high levels of CD68, HLA‐DR, and CD11c, and all macrophage and dendritic cell markers; and intermediate levels of the immune checkpoint PD‐L1 (Figure 5B). The neutrophil marker CD66b was not detected in the selected ROIs, suggesting that neutrophils are not involved at least in the crosstalk between macrophages and T cells in these particular BAP1− mUM cases. Lymphoid markers CD8A and CD4 were also found to be highly expressed. CD56 levels are relatively low compared with negative controls, suggesting absence of NK cells in the selected regions. However, T regulatory cells seem to be absent as intracellular levels of Foxp3 were not detected among these patients using this technique.

In addition, low positive levels of granzyme B were detected, together with high expression of B7‐H3, a checkpoint regulator of lymphocyte functions 60. We also observed the expression of IDO‐1, TIGIT, and VISTA (Figure 5B). IDO‐1 is known to induce adaptive resistance to anti‐PD1 and anti‐CTLA4 immunotherapies 61, 62. In addition, IDO and TIGIT were recently described to be expressed in pUM‐M3 with corresponding mUM tissues 63.

Considered altogether, these findings suggest alternative mechanisms of T‐cell exhaustion other than PD‐1 and CTLA‐4 engagement, as well as the involvement of mechanisms for immune exclusion that may undertake an important role to support tumour immune evasion and, consequently, immunotherapy failure.

Discussion

In this multi‐parametric immunophenotyping work in UM, we profiled the immune response of the 80 patients of the TCGA‐UM study, highlighting the involvement of BAP1 loss in the coordination of gene expression of several immune markers. Selected biomarkers were further investigated in a smaller number of primary and metastatic BAP1− UMs, at both transcriptome and protein levels, using cutting‐edge techniques. Recently, our group observed that different patterns of nuclear BAP1 expression in pUM provide insights into the prognostic significance of this tumour 64. Among UMs with an M3 status, the cumulative survival of patients with UM expressing nuclear BAP1 is significantly greater than that of UM patients whose tumours are M3 with nuclear BAP1 loss. These findings and previous reports associating BAP1 loss with a wide spectrum of cancers 65 underpin the molecular mechanisms behind the adverse prognostic effects of M3, supporting the importance of analysing immune gene expression from the aspect of BAP1loss in UM.

The genetic diversity of UM was recently described, including copy number variations (CNVs), somatic mutations, and BAP1 alterations 63. However, the diversity of immune gene expression is described by the tumour stroma (i.e. the features of the reactive cells in the tumour microenvironment), which shapes accordantly with the tumour phenotype (e.g. BAP1 loss). As a consequence of BAP1 loss in UM, tumour cells could therefore unleash metabolic mechanisms to secrete different factors that would induce the regulatory phenotype of T cells and macrophages in the tumour microenvironment (TME) to a more tolerogenic profile. The positive expression of MIF observed at the transcriptomic level in mUM is in accordance with previous reports showing that UM cells can secrete MIF as a mechanism of immune escape 66. The main receptor of MIF in different immune cells is CD74 36, 67. In addition, MIF was recently reported to induce M2 polarisation of TAMs, leading to immune suppression on several solid cancers 34, 36, 66, 68, 69, 70, 71, 72, and to downregulate the CTL responses 34, which may occur via its interaction with the CD74 receptor expressed on CD8+ T cells in UM, as observed for cluster C. Since CD74 expression was not observed in the CD8+ T‐cell clusters of one BAP1+ pUM, we hypothesise that tumours with reduced expression of MIF or its receptor CD74 may contribute to increase the frequency of more effective CD8+ T cells in the TME of UM.

Therefore, modulation of T lymphocytes and macrophages toward an immunosuppressive phenotype could be explained by the expression of CD74 on these cells, which can be affected by suppressive factors derived from tumour cells, including MIF. CD74, which was highly expressed at the transcriptomic level in pUM, was also found at the protein level across the regulatory CD8+ T‐cell cluster and CD163+ M2‐like macrophage cluster. CD74 is a chaperone involved in the trafficking of HLA‐DR molecules to the surface of immune cells, and while it remains expressed on the surface of the cells, it may bind to MIF secreted by tumour cells in the TME 46, 67, 73.

The pharmacological blockade of the MIF/CD74 interaction restores the TME immunogenic profile, as well as an effective anti‐tumour immune response against metastatic melanoma and gliomas 36, 46. The CD74 monoclonal blocking antibody milatuzumab is currently approved by the Food and Drug Administration (FDA) in the United States for the treatment of multiple myeloma, non‐Hodgkin lymphomas, and other CD74+ cancers 74, 75.

Changes in BAP1 expression have also been associated with immune transformation in mesothelioma, and became a predictive tool for immunotherapy of peritoneal mesothelioma 76, 77. The impaired ability of thymic development and the proliferative responses of T lymphocytes in the context of BAP1 inhibition are strong evidence that loss of BAP1 function is associated with immune suppression and systemic myeloid transformation 16, 78.

In this study, we also observed that increased transcriptome levels of CD38, HLA‐DRA, IDO1, and LAG‐3 are significantly correlated with BAP1 loss. These immune biomarkers are of extreme importance because they have been associated with different immune‐suppressive pathways, which suggests mechanistic insights for immune suppression and immunotherapy resistance using ICIs 30, 36, 37, 38, 41. In the protein single‐cell level, we show the functional state of UM infiltrating CD8+ T cells, which co‐express high levels of CD38, HLA‐DR, and CD28. The co‐expression of HLA‐DR/CD28 in CD8+ T cells suggests that these lymphocytes are distinct from cytolytic effector T cells 30, and can be classified as regulatory CD8+ T cells, with similar functions to classical CD4+Foxp3+ cells 28. In addition, higher levels of CD38 demarcates regulatory and memory status to CD8+ T cells in the context of IFN‐γ‐mediated immunosuppression, and was recently addressed to drive mechanisms of tumour‐mediated immune escape to immunotherapies using PD1/PD‐L1 blockade 37, 38. Indeed, IFN‐γ is upregulated in the context of BAP1 loss and is widely associated with several immune‐suppressive network categories, which is in accordance with recent reports showing the immune‐suppressive roles of IFN‐γ 79.Therefore, targeting CD38 in UM may be considered a suitable strategy to improve the efficacy of immunotherapy using ICI in metastatic UM. A recent study showed that targeting CD38 using isatuximab can preferentially block immunosuppressive T‐regulatory lymphocytes and therefore restore immune effector function against multiple myeloma 53.

The low expression of ICRs LAG‐3 and CTLA‐4 among the majority of T‐cell clusters suggests that these lymphocytes may not be exhausted but exist in a lower activation state in pUM 80, 81, 82, 83. Increased transcriptome levels of LAG‐3 in the TCGA‐UM study could be linked with LAG‐3 expression among CD25+ B‐cell clusters as evidenced by mass cytometry, suggesting a memory and natural regulatory phenotype for these cells 84, 85. Moreover, higher expression of IDO1 in both pUM and mUM suggests this molecule as an important adjuvant target for immunotherapy using ICIs, since IDO1 blockade has been shown to synergise the therapeutic effector of both CTLA‐4 and PD1/PD‐L1 inhibitors 61.

Our findings in this report also provide evidence that BAP1 − pUM could shape an immune response similar to that in mUM tissues, since BAP1‐loss‐correlated immune genes are similarly expressed in mUM, as observed using the NanoString approach. Furthermore, the DSP approach revealed that additional ICI‐resistant mechanisms not necessarily related to BAP1 changes may also be important in mUM‐induced exclusion of immune cells, such as the Wnt/β‐catenin axis 58, 59. A recent study showed that hepatocellular carcinoma patients displaying an altered Wnt/β‐catenin pathway were refractory to immune‐checkpoint blockade 86, which is aligned with evidence that melanoma‐intrinsic β‐catenin signalling prevents anti‐tumour immunity 87.

It is important to highlight that a weakness of this study is the low number of analysed BAP1 − UM samples for CyTOF studies. The reason behind that is the scarcity of the type of fresh tumour sample, not only because this type of tumour is very rare but also because the tissues must be sufficiently large to provide significant amounts of immune infiltrated cells for further downstream analysis, and thus reducing the sample size of this study. Despite this, we could not only reproduce and confirm previous data published regarding the higher frequency of infiltrated CD8+ T cells over CD4+ T cells at the transcriptomic and protein levels 4, 15, but also detect the expression of specific immune markers initially detected in the transcriptome analysis of the larger TCGA‐UM cohort (n = 80), and also observed in mUM tissues, expanding the impact of our pUM mass cytometry findings.

The present work shows an improved overview of the immune profile of pUM and mUM at both transcriptome and protein levels and suggests that immune modulation in UM may be driven by loss of BAP1 expression. Immunosuppressive networks found in BAP1 − tumours may not only influence the quality and quantity of local anti‐tumour immune responses but also affect immunotherapy outcomes using ICIs, leading to regulation or exclusion of T effector lymphocytes, as well as alternative polarisation of macrophages toward a tolerogenic phenotype in the TME. The relative importance of these findings will require further functional validation, and this study provides the solid ground to initiate these studies. Detected key immune biomarkers, such as CD38 and CD74, could be immediately investigated for functional validation in the adjuvant settings of ICI immunotherapies, since there are currently available FDA‐approved inhibitors against these targets 53, 74. Altogether, this work provides the most critical immune markers and pathways to consolidate the type of immune responses in the context of BAP1 loss in UM. This may help us to understand why this type of cancer is one of the most refractory to current immunotherapies using ICIs at present.

Author contributions statement

CRF performed the research and drafted the manuscript. HK assisted with methodology design, sample collection and preparation, conceptual advice, and manuscript drafting. JJS assisted with methodology design, conceptual advice, and manuscript drafting. RAA assisted with data analysis. AD provided support with mass cytometry data analysis and interpretation of results. SJR provided support to mass cytometry work. JMC assisted with conceptual advice and manuscript drafting. SEC was involved in methodology design and sample collection, supervised the research, and drafted the manuscript.

Figure S1. Effects of chromosome 3 loss

Click here for additional data file.

Figure S2. Immunohistochemistry analysis of one BAP1− mUM for selective immune markers: CD3, CD4, CD8, CD68, CD38, TIA1, and Tbet

Click here for additional data file.

Figure S3. Mass cytometry results

Click here for additional data file.

Table S1. Immune categories combined from the nCounter PanCancer immune profiling gene set with a literature review

Click here for additional data file.

Table S2. Immune categories for supervised network analysis

Click here for additional data file.

Table S3. Leukocyte immune response categories

Click here for additional data file.

Table S4. Mass cytometry antibodies and reagents

Click here for additional data file.

Table S5. Immune genes Spearman's correlation with BAP1 expression or Chr3 loss

Click here for additional data file.

Table S6. Spearman's rank correlation coefficient (r) and P values of immune genes exclusively correlated with BAP1 gene expression

Click here for additional data file.

Table S7. Spearman's rank correlation coefficient (r) and P values of immune genes with better correlation to BAP1 gene expression than Chr3 copy number variations

Click here for additional data file.

Table S8. List of Kaplan–Meier survival test scores and P values in the context of BAP1 loss

Click here for additional data file.

Table S9. Spearman's rank correlation coefficient (r) and P values of immune genes that exclusively correlate with Chr3 copy number variations (M3‐UM)

Click here for additional data file.

Table S10. Spearman's rank correlation coefficient (r) and P values of immune genes with better correlation to Chr3 copy number variations (M3‐UM) than BAP1 gene expression

Click here for additional data file.

Table S11. List of Kaplan–Meier survival test scores and P values in the context of Chr3 copy number variation

Click here for additional data file.