Research | Open | Published:
BAP1 haploinsufficiency predicts a distinct immunogenic class of malignant peritoneal mesothelioma
Genome Medicinevolume 11, Article number: 8 (2019)
The Research Highlight to this article has been published in Genome Medicine 2019 11:18
Malignant peritoneal mesothelioma (PeM) is a rare and fatal cancer that originates from the peritoneal lining of the abdomen. Standard treatment of PeM is limited to cytoreductive surgery and/or chemotherapy, and no effective targeted therapies for PeM exist. Some immune checkpoint inhibitor studies of mesothelioma have found positivity to be associated with a worse prognosis.
To search for novel therapeutic targets for PeM, we performed a comprehensive integrative multi-omics analysis of the genome, transcriptome, and proteome of 19 treatment-naïve PeM, and in particular, we examined BAP1 mutation and copy number status and its relationship to immune checkpoint inhibitor activation.
We found that PeM could be divided into tumors with an inflammatory tumor microenvironment and those without and that this distinction correlated with haploinsufficiency of BAP1. To further investigate the role of BAP1, we used our recently developed cancer driver gene prioritization algorithm, HIT’nDRIVE, and observed that PeM with BAP1 haploinsufficiency form a distinct molecular subtype characterized by distinct gene expression patterns of chromatin remodeling, DNA repair pathways, and immune checkpoint receptor activation. We demonstrate that this subtype is correlated with an inflammatory tumor microenvironment and thus is a candidate for immune checkpoint blockade therapies.
Our findings reveal BAP1 to be a potential, easily trackable prognostic and predictive biomarker for PeM immunotherapy that refines PeM disease classification. BAP1 stratification may improve drug response rates in ongoing phases I and II clinical trials exploring the use of immune checkpoint blockade therapies in PeM in which BAP1 status is not considered. This integrated molecular characterization provides a comprehensive foundation for improved management of a subset of PeM patients.
Malignant mesothelioma is a rare but aggressive cancer that arises from the internal membrane lining of the pleura and the peritoneum. While the majority of mesotheliomas are pleural in origin, the incidence of peritoneal mesothelioma (PeM) accounts for approximately 20–30% of all mesothelioma cases in the USA and possibly higher in other countries . Occupational asbestos exposure is a significant risk factor in the development of pleural mesothelioma (PM). However, epidemiological studies suggest that unlike PM, asbestos exposure plays a far smaller role in the etiology of PeM tumors . More importantly, the incidence of PeM is skewed towards young women of childbearing ages rather than in old patients  making PeM a malignancy often associated with many years of life lost.
Previous studies in mesotheliomas have revealed that over 60% of mesotheliomas harbor BRCA1 associated protein 1 (BAP1) inactivating mutation or copy number loss, making BAP1 the most commonly altered gene in this malignancy [3,4,5,6,7]. BAP1 is a tumor suppressor and deubiquitinase, localized to the nucleus, known to regulate chromatin remodeling and maintain genome integrity [8, 9]. Furthermore, BAP1 localized in endoplasmic reticulum regulate calcium (Ca2+) flux to promote apoptosis . Thus, the combined reduced BAP1 nuclear and cytoplasmic activity results in the accumulation of DNA-damaged cells and greater susceptibility to the development of malignancy. In addition, inactivating mutations of neurofibromin 2 (NF2) and cyclin-dependent kinase inhibitor 2A (CDKN2A) are also relatively common, while other mutations are rare. Previous studies in PeM [11,12,13,14,15,16,17,18] have only focused on genomic information; therefore, the downstream consequences of these genomic alterations are not well understood. Genome information coupled with transcriptome and proteome information is more likely to be successful in revealing potential therapeutic modalities.
Mesothelioma is typically diagnosed in the advanced stages of the disease. A combination of cytoreductive surgery (CRS) and hyperthermic intraperitoneal chemotherapy (HIPEC), sometimes followed by normothermic intraperitoneal or systemic chemotherapy (NIPEC), has recently emerged as the first-line treatment for PeM . However, even with this regime, complete cytoreduction is hard to achieve and death ensues for many patients. Actionable molecular targets for PeM critical for precision oncology remain to be defined. Immune checkpoint blockade therapy in PM has recently gained traction [7, 20] given that 20–40% of PM cases are reported to show an inflammatory phenotype . However, the role of immunostaining for PD-L1, the usual approach to predicting a response to immunotherapy for other tumor types, is controversial in PM, since positive stating has generally been associated with a worse prognosis, and it is unclear what marker should be used to predict tumors that may respond to immunotherapy.
Although, clinical trials typically lump PeM and PM together for immune checkpoint blockade [22,23,24,25,26], even less is known about PeM and immunotherapy. Thus, there has been no attempt to stratify PeM patients. In this study, we performed an integrated multi-omics analysis of the genome, transcriptome, and proteome of 19 PeM, predominantly of epithelioid subtype, and correlated these with tumor inflammation.
We assembled a cohort of 19 PeM from 18 patients (Table 1 and Additional file 2: Table S1) undergoing CRS at Vancouver General Hospital (Vancouver, Canada), Mount Sinai Hospital (Toronto, Canada), and Moores Cancer Centre (San Diego, CA, USA). We obtained 19 fresh-frozen primary treatment-naïve PeM tumor tissue and adjacent benign tissues or whole blood from the 18 patients. For 1 patient, MESO-18, 2 tumors from distinct sites were available. Immunohistochemical analyses using different biomarkers were evaluated by 2 independent pathologists (Additional file 1: Figure S1-S4). Both pathologists categorized all 19 tumors as epithelioid PeM with a content of higher than 75% tumor cellularity. To the best of our knowledge, this is the largest cohort of PeM subjected to an integrative multi-omics analysis.
Immunohistochemistry and histopathology
Freshly cut tissue microarray (TMA) sections were analyzed for immunoexpression using Ventana Discovery Ultra autostainer (Ventana Medical Systems, Tucson, AZ). In brief, tissue sections were incubated in Tris-EDTA buffer (CC1) at 37 °C to retrieve antigenicity, followed by incubation with respective primary antibodies at room temperature or 37 °C for 60–120 min. For primary antibodies, mouse monoclonal antibodies against CD8 (Leica, NCL-L-CD8-4B11, 1:100), CK5/cytokeratin 5 (Abcam, ab17130, 1:100), BAP1 (SantaCruz, clone C4, sc-28383, 1:50), rabbit monoclonal antibody against CD3 (Abcam, ab16669, 1:100), and rabbit polyclonal antibodies against CALB2/calretinin (LifeSpan BioSciences, LS-B4220, 1:20 dilution) were used. Bound primary antibodies were incubated with Ventana Ultra HRP kit or Ventana universal secondary antibody and visualized using Ventana ChromoMap or DAB Map detection kit, respectively. All stained slides were digitalized with the SL801 autoloader and Leica SCN400 scanning system (Leica Microsystems; Concord, Ontario, Canada) at magnification equivalent to × 20. The images were subsequently stored in the SlidePath digital imaging hub (DIH; Leica Microsystems) of the Vancouver Prostate Centre. Representative tissue cores were manually identified by two pathologists.
Whole exome sequencing
DNA was isolated from snap-frozen tumors with 0.2 mg/ml Proteinase K (Roche) in a cell lysis solution using Wizard Genomic DNA Purification Kit (Promega Corporation, USA). Digestion was carried out overnight at 55 °C before incubation with RNase solution at 37 °C for 30 min and treatment with protein precipitation solution followed by isopropanol precipitation of the DNA. The amount of DNA was quantified on the NanoDrop 1000 Spectrophotometer and an additional quality check done by reviewing the 260/280 ratios. Quality check was done on the extracted DNA by running the samples on a 0.8% agarose/TBE gel with ethidium bromide.
For Ion AmpliSeq™ Exome Sequencing, 100 ng of DNA based on Qubit® dsDNA HS Assay (Thermo Fisher Scientific) quantitation was used as input for Ion AmpliSeq™ Exome RDY Library preparation. This is a polymerase chain reaction (PCR)-based sequencing approach using 294,000 primer pairs (amplicon size range 225–275 bp) and covers > 97% of Consensus CDS (CCDS; release 12), > 19,000 coding genes, and > 198,000 coding exons. Libraries were prepared, quantified by quantitative PCR (qPCR), and sequenced according to the manufacturer’s instructions (Thermo Fisher Scientific). Samples were sequenced on the Ion Proton System using the Ion PI™ Hi-Q™ Sequencing 200 Kit and Ion PI™ v3 chip. Two libraries were run per chip for a projected coverage of 40 M reads per sample.
Somatic variant calling
Torrent Server (Thermo Fisher Scientific) was used for signal processing, base calling, read alignment, and generation of results files. Specifically, following sequencing, reads were mapped against the human reference genome hg19 using the Torrent Mapping Alignment Program. Variants were identified by using Torrent Variant Caller plugin with the optimized parameters for AmpliSeq exome-sequencing recommended by Thermo Fisher. The variant call format (VCF) files from all samples were combined using GATK (3.2-2) , and all variants were annotated using ANNOVAR . Only non-silent exonic variants including non-synonymous single nucleotide variations (SNVs), stop-codon gain SNVs, stop-codon loss SNVs, splice site SNVs, and In-Dels in coding regions were kept if they were supported by more than ten reads and had allele frequency higher than 10%. To obtain somatic variants, we filtered against dbSNP build 138 (non-flagged only) and the matched adjacent benign or blood samples sequenced in this study. Putative variants were manually scrutinized on the Binary Alignment Map (BAM) files through Integrative Genomics Viewer version 2.3.25 .
Copy number aberration (CNA) analysis
Copy number changes were assessed using Nexus Copy Number Discovery Edition version 9.0 (BioDiscovery, Inc., El Segundo, CA). Nexus NGS functionality (BAM ngCGH) with the FASST2 segmentation algorithm was used to make copy number calls (a circular binary segmentation/hidden Markov model approach). The significance threshold for segmentation was set at 5 × 10−6, also requiring a minimum of 3 probes per segment and a maximum probe spacing of 1000 between adjacent probes before breaking a segment. The log ratio thresholds for single copy gain and single copy loss were set at + 0.2 and − 0.2, respectively. The log ratio thresholds for the gain of 2 or more copies and for a homozygous loss were set at + 0.6 and − 1.0, respectively. Tumor sample BAM files were processed with corresponding normal tissue BAM files. Reference reads per CNA point (window size) was set at 8000. Probes were normalized to the median. Relative copy number profiles from exome sequencing data were determined by normalizing tumor exome coverage to values from whole blood controls.
Whole transcriptome sequencing (RNA-seq)
Total RNA from 100 μm sections of snap-frozen tissue was isolated using the mirVana Isolation Kit from Ambion (AM-1560). Strand-specific RNA sequencing was performed on quality controlled high RIN value (> 7) RNA samples (Bioanalyzer Agilent Technologies) before processing at the high throughput sequencing facility core at BGI Genomics Co., Ltd. (The Children’s Hospital of Philadelphia, PA, USA). In brief, 200 ng of total DNAse-treated RNA was first treated to remove the ribosomal RNA (rRNA) and then purified using the Agencourt RNA Clean XP Kit (Beckman Coulter) prior to analysis with the Agilent RNA 6000 Pico Chip to confirm rRNA removal. Next, the rRNA-depleted RNA was fragmented and converted to cDNA. Subsequent steps include end repair, addition of an “A” overhang at the 3′ end, and ligation of the indexing-specific adaptor, followed by purification with Agencourt Ampure XP beads. The strand-specific RNA library prepared using TruSeq (Illumina catalog no. RS-122-2201) was amplified and purified with Ampure XP beads. Size and yield of the barcoded libraries were assessed on the LabChip GX (Caliper), with an expected distribution around 260 bp. The concentration of each library was measured with real-time PCR. Pools of the indexed library were then prepared for cluster generation and PE100 sequencing on Illumina HiSeq 4000. The RNA-seq reads were aligned using STAR (2.3.1z)  onto the human genome reference (GRCh38), and the transcripts were annotated based on Ensembl release 80 gene models. Only the reads unique to one gene and which corresponded exactly to one gene structure were assigned to the corresponding genes by using HTSeq . Normalization of the read counts was conducted by DESeq . For a detailed description, see Additional file 1: Supplementary Methods.
Proteomics analysis using mass spectrometry
Fresh-frozen samples dissected from tumor and adjacent normal were individually lysed in 50 mM of HEPES pH 8.5, 1% SDS, and the chromatin content was degraded with Benzonase. The tumor lysates were sonicated (Bioruptor Pico, Diagenode, NJ, USA), and disulfide bonds were reduced with DTT and capped with iodoacetamide. Proteins were cleaned up using the SP3 method [33, 34] (Single Pot, Solid Phase, Sample Prep) then digested overnight with trypsin in HEPES pH 8, peptide concentration determined by Nanodrop (Thermo) and adjusted to an equal level. A pooled internal standard control was generated comprising of equal volumes of every sample (10 μl of each of the 100 μl total digests) and split into 3 equal aliquots. The labeling reactions were run as 3 TMT 10-plex panels (9+IS) then desalted, and each panel is divided into 48 fractions by reverse-phase HPLC at pH 10 with an Agilent 1100 LC system. The 48 fractions were concatenated into 12 superfractions per panel by pooling every fourth fraction eluted resulting in a total of 36 overall samples. These samples were analyzed with an Orbitrap Fusion Tribrid Mass Spectrometer (Thermo Fisher Scientific) coupled with EasyNanoLC 1000 using a data-dependent method with synchronous precursor selection (SPS) MS3 scanning for TMT tags. Based on ProteomeDiscoverer 18.104.22.168 (Thermo Fisher Scientific), we selected peptide-spectrum match (PSM) results with q value ≤ 0.05 and extract proteins from both high and medium confidence level after false discovery rate filtering for protein identification and quantification results. For a detailed description, see Additional file 1: Supplementary Methods.
Prioritization of driver genes using HIT’nDRIVE
Non-silent somatic mutation calls, CNA gain or loss, and gene-fusion calls were collapsed in gene-patient alteration matrix with binary labels. Gene expression values were used to derive an expression-outlier gene-patient outlier matrix using the generalized extreme studentized deviate (GESD) test. STRING ver10  protein interaction network was used to compute pairwise influence value between the nodes in the interaction network. We integrated these genome and transcriptome data using the HIT’nDRIVE algorithm . The following parameters were used: α = 0.9, β = 0.6, and γ = 0.8. We used IBM-CPLEX as the integer linear programming (ILP) solver.
Stromal and immune score
We used 2 sets of 141 genes (1 each for stromal and immune gene signatures) as described in . We used the “inverse normal transformation” method to transform the distribution of expression data into the standard normal distribution. The stromal and immune scores were calculated, for each sample, using the summation of standard normal deviates of each gene in the given set.
Enumeration of tissue-resident immune cell types using mRNA expression profiles
CIBERSORT software  was applied to the RNA-seq gene expression data to estimate the proportions of 22 immune cell types (B cells naive, B cells memory, plasma cells, T cells CD8, T cells CD4 naive, T cells CD4 memory resting, T cells CD4 memory activated, T cells follicular helper, T cells gamma delta, T cells regulatory (Tregs), NK cells resting, NK cells activated, monocytes, macrophages M0, macrophages M1, macrophages M2, dendritic cells resting, dendritic cells activated, mast cells resting, mast cells activated, eosinophils, and neutrophils) using LM22 dataset provided by CIBERSORT platform. Genes not expressed in any of the PeM tumor samples were removed from the LM22 dataset. The analysis was performed using 1000 permutations. The 22 immune cell types were later aggregated into 9 distinct groups.
Landscape of somatic mutations in PeM
To investigate the landscape of somatic gene mutations in PeM, we performed high-coverage whole exome sequencing of 19 tumors and 16 matched normal samples (Additional file 2: Table S1). We achieved a mean coverage of 180× for cancerous samples and 96× for non-cancerous samples (Additional file 2: Table S2). We identified 346 unique non-silent mutations affecting 202 unique genes (Additional file 1: Figure S5 and Additional file 2: Table S3). We observed an average of 0.015 protein-coding non-silent mutations per Mb per tumor sample.
We first identified driver genes of PeM using our recently developed computational algorithm HIT’nDRIVE . Briefly, HIT’nDRIVE measures the potential impact of genomic aberrations on changes in the global expression of other genes/proteins which are in close proximity in a gene/protein interaction network. It then prioritizes those aberrations with the highest impact as cancer driver genes. HIT’nDRIVE prioritized 25 unique driver genes in 15 PeM samples for which matched genome and transcriptome data were available (Fig. 1 and Additional file 2: Table S4). Six genes (BAP1, BZW2, ABCA7, TP53, ARID2, and FMN2) were prioritized as drivers, harboring single nucleotide changes.
BAP1 was the most frequently mutated gene (5 out of 19 tumors) in PeM. Among the 5 BAP1-mutated cases, 2 cases (MESO-06 and MESO-09) were predicted to have inactivated BAP1, whereas despite BAP1 mutation in 3 cases (MESO-18A/E and MESO-19), their mRNA transcripts were expressed in high levels (Fig. 2c and Additional file 1: Figure S6-S7). We identified that all variants of BAP1 (except a 42-bp deletion in MESO-09) were expressed at the RNA level (Additional file 2: Table S16). In addition, we identified mutations in genes such as TP53, SETD2, SETDB1, and LATS1 each present in just a single case (Fig. 1).
Copy number landscape in PeM
The aggregate somatic copy number aberration (CNA) profile of PeM is shown in Fig. 2a, b. We observed a total of 1281 CNA events across all samples (Additional file 2: Table S5). On average, 10% of the protein-coding genome was altered per PeM. Interestingly, the CNA burden in PeM was strongly correlated (R = 0.74) with its mutation burden (Additional file 1: Figure S9).
Using HIT’nDRIVE, we identified genes in chromosomes 3p21, BAP1, PBRM1, and SETD2, as key driver genes of PeM (Fig. 1 and Additional file 2: Table S4). This region was also identified as significantly recurrent focal CNAs using the GISTIC  algorithm (Additional file 1: Figure S9). Chromosome 3p21 was deleted (homozygous or heterozygous) in almost half of the tumors (8 of 19) in the cohort. Here, we call tumors with 3p21 (or BAP1) loss as BAP1del and the rest of the tumors with 3p21 (or BAP1) copy number intact as BAP1intact. Interestingly, BAP1 mRNA transcripts in BAP1del tumors were expressed at lower levels as compared to those in BAP1intact tumors (p value = 3 × 10−4) (Fig. 2c). We validated this using immunohistochemical (IHC) staining demonstrating a lack of BAP1 nuclear staining in the tumors with BAP1 homozygous deletion (Fig. 2d). Tumors with BAP1 heterozygous loss still displayed BAP1 nuclear staining (Additional file 1: Figure S10). We observed 3 BAP1-mutated cases (MESO-18A/E and MESO-19) among BAP1intact tumors. BAP1 mRNA transcripts in these 3 tumors were expressed at high levels (Fig. 2c). Furthermore, we found DNA copy loss of 3p21 locus to include 4 concomitantly deleted cancer genes—BAP1, SETD2, SMARCC1, and PBRM1—consistent with . Copy number status of these 4 genes was significantly correlated with their corresponding mRNA expression (Additional file 1: Figure S11), suggesting that the allelic loss of these genes is associated with decreased transcript levels. These 4 genes are chromatin modifiers, and PBRM1 and SMARCC1 are part of the SWI/SNF complex that regulates transcription of a number of genes.
The global transcriptome and proteome profile of PeM
To segregate transcriptional subtypes of PeM, we performed total RNA-seq (Illumina HiSeq 4000) and its quantification of 15 PeM tumor samples for which RNA were available. Using principal component analyses, we found that tumor samples in BAP1intact and BAP1del subtypes have distinct transcriptomic patterns; however, a few samples showed an overlapping pattern (Additional file 1: Figure S16A).
We performed mass spectrometry (Fusion Orbitrap LC/MS/MS) with isobaric tagging for expressed peptide identification and its corresponding protein quantification using Proteome Discoverer for processing pipeline for 16 PeM tumors and 7 matched normal tissues. We identified 8242 unique proteins in 23 samples analyzed. We were surprised BAP1 protein was however not detected in our MS experiment, likely due to inherent technical limitations with these samples and/or processing. Quality control analysis of in-solution Hela digests also has very low BAP1 with only a single peptide observed in occasional runs. Unlike in transcriptome profiles, the proteome profiles of tumor samples in BAP1intact and BAP1del subtypes did not group into distinct clusters (Additional file 1: Figure S16B).
Next, we identified differentially expressed genes and proteins between BAP1intact and BAP1del subtypes (see Additional file 1: Supplemental Methods). As expected, BAP1, PBRM1, SMARCA4, and SMARCD3 were among the top 500 differentially expressed genes. Many other important cancer-related genes were differentially expressed such as CDK20, HIST1H4F, ERCC1, APOBEC3A, CDK11A, CSPG4, TGFB1, IL6, LAG3, and ATM.
To identify the pathways dysregulated by the differentially expressed genes between the PeM subtypes, we performed gene set enrichment analysis (see Additional file 1: Supplementary Methods). Intriguingly, we observed high concordance between pathways dysregulated by the 2 sets (mRNA and protein expression data) of top 500 differentially expressed genes and proteins (Fig. 3a, b). The unsupervised clustering of pathways revealed 2 distinct clusters for BAP1del and BAP1intact tumors. This indicates that the enriched pathways, between the patient groups, are also differentially expressed. BAP1del tumors demonstrated elevated levels of RNA and protein metabolism as compared to BAP1intact tumors. Many genes involved in chromatin remodeling and DNA damage repair were differently expressed between the groups (Additional file 1: Figure S20-S21). Genes in DNA damage repair pathway—PARP1, ERCC1, and APC—were downregulated, and CHEK1/2, BRAC2, RAD50, and ATM were upregulated in BAP1del tumors. Genes involved in cell cycle and apoptotic pathways were observed to be highly expressed in BAP1del patients. Furthermore, glucose and fatty acid metabolism pathways were repressed in BAP1del as compared to BAP1intact. More interestingly, we observed a striking difference in immune system-associated pathways between the PeM subtypes, whereas BAP1del tumors demonstrated strong activity of cytokine signaling and the innate immune system; the MHC-I/II antigen presentation system and adaptive immune system were active in BAP1intact tumors.
BAP1 del subtype is correlated with tumor inflammation characterized by immune checkpoint receptor activation
Prompted by this finding, we next analyzed whether PeM were infiltrated with leukocytes. To assess the extent of leukocyte infiltration, we computed an expression-based (RNA-seq and protein) score (see the “Methods” section) using the immune cell and stromal markers proposed by . We discovered that the immune marker gene score was strongly correlated with the stromal marker gene score (Fig. 4a) suggesting possible leukocyte infiltration in PeM from the tumor microenvironment. Furthermore, using CIBERSORT  software, we computationally estimated the leukocyte representation in the bulk tumor transcriptome. We observed massive infiltration of T cells in majority of the PeM (Fig. 4b). A subset of PeM had massive infiltration of B cells in addition to T cells. Interestingly, when we group the PeM by their BAP1 aberration status, there was a marked difference in the proportion of infiltrated plasma cells, natural killer (NK) cells, mast cells, and B cells between the groups. Whereas the proportions of plasma cells, NK cells, and B cells were less in the BAP1del tumors, there was more infiltration of mast cells and T cells in BAP1del tumors as compared to BAP1intact tumors. We performed TMA IHC staining of CD3 and CD8 antibodies on PeM tumors. We observed that BAP1del PeM were positively stained for both CD3 and CD8 confirming infiltration of T cells in BAP1del PeM (Fig. 4c and Additional file 1: Figure S22-S23). Combined, this strongly indicates that PeM could be divided into tumors with an inflammatory tumor microenvironment and those without and that this distinction correlated with BAP1 haploinsufficiency.
Finally, we surveyed PeM for the expression of genes involved in immune checkpoint pathways. A number of immune checkpoint receptors were highly expressed in BAP1del tumors relative to BAP1intact tumors. These included PDCD1 (PD1), CD274 (PD-L1), CD80, CTLA4, LAG3, and ICOS (Fig. 4d and Additional file 1: Figure S30) for which inhibitors are either clinically approved or are at varying stages of clinical trials. Notably, the differential gene expression pattern of LAG3, ICOS, and CTLA4 between the PeM subtypes suggests potential opportunities for immune checkpoint blockade beyond conventional PD1/PD-L1. Moreover, a number of MHC genes, immuno-inhibitor genes, and immuno-stimulator genes were differentially expressed between BAP1del and BAP1intact tumors (Additional file 1: Figure S24). Furthermore, we analyzed whether the immune checkpoint receptors were differentially expressed in tumors with and without 3p21 loss in the pleural mesotheliomas (PM) from The Cancer Gene Atlas (TCGA) project . Unlike in PeM, we did not observe a significant difference in immune checkpoint receptor expression between the PM groups (i.e., BAP1del and BAP1intact) (Additional file 1: Figure S25). These findings suggest that BAP1del PeM tumors could potentially be targeted with immune checkpoint inhibitors while PM tumors may less likely to respond.
In this study, we present a comprehensive integrative multi-omics analysis of malignant peritoneal mesotheliomas. Even though this is a rare disease, we managed to amass a cohort of 19 tumors. Prior studies of mesotheliomas, performed using a single omic platform, have established BAP1 inactivation as a key driver event in mesotheliomas. Our novel contribution to PeM is that we provide evidence from integrative multi-omics analyses that BAP1 haploinsufficiency (BAP1del) forms a distinct molecular subtype of PeM. This subtype is characterized by distinct expression patterns of genes involved in chromatin remodeling, DNA repair pathway, and immune checkpoint activation. Moreover, BAP1del subtype is correlated with inflammatory tumor microenvironment. Our results suggest that BAP1del tumors might be prioritized for immune checkpoint blockade therapies. Thus, BAP1 is likely both prognostic and predictive biomarker for PeM enabling better disease stratification and patient treatment. Further corroborating our findings, BAP1 status has been recently shown to be correlated with perturbed immune signaling in PM .
Loss of BAP1 is known to alter chromatin architecture exposing the DNA to damage and also impairing the DNA repair machinery [9, 40]. The DNA repair defects thus drive genomic instability and dysregulate tumor microenvironment . DNA repair deficiency leads to the increased secretion of cytokines, including interferons that promote tumor-antigen presentation and trigger recruitment of T lymphocytes to destroy tumor cells. As a response, tumor cells evade this immune surveillance by increased expression of immune checkpoint receptors. The results presented here also indicate that PeM are infiltrated with immune cells from the tumor microenvironment. Moreover, the BAP1del subtype displays elevated levels of immune checkpoint receptor expression which strongly suggests the use of immune checkpoint inhibitors to treat this subtype of PeM. However, in a small subset of PM tumors in TCGA dataset, the loss of BAP1 did not elevate the expression of immune checkpoint marker genes. This warrants further investigation on the characteristics of these groups of PM.
The main challenge in mesothelioma treatment is that all current efforts made towards testing new therapy options are limited to using therapies that have been proven successful in other cancer types, without a good knowledge of underlying molecular mechanisms of the disease. As a result of sheer desperation, some patients have been treated even though no targeted therapy for mesothelioma has been proven effective as yet. For example, a number of clinical trials exploring the use of immune checkpoint blockade (anti-PD1/PD-L1 or anti-CTLA4) in PM and/or PeM patients are currently under progress. The results of the first few clinical trials report either a very low response rate or no benefit to the patients [22,23,24, 26, 42]. Notably, BAP1 copy number or mutation status was not assessed in these studies. Our study warrants further investigation of immune checkpoint molecules targeting beyond conventional PD1/PD-L1. We hypothesize based on this evidence presented that response rates for immune checkpoint blockade therapies in clinical trials for PeM will improve when patients are segregated by their BAP1 copy number status.
Our first-in-field multi-omics analysis of PeM tumors identified BAP1 haploinsufficiency as a distinct molecular subtype and a candidate for immune checkpoint blockade therapies. This is significant because almost half of PeM cases are now candidates for these therapies. BAP1 status is not currently taken into account in the ongoing phases I and II clinical trials exploring the use of immune checkpoint blockade therapies in PeM. Moreover, this is the first study to demonstrate evidence of inflammatory tumor microenvironment in PeM. Our findings identify BAP1 as a tractable prognostic and predictive biomarker for immunotherapy that refines PeM disease stratification and may improve drug response rates.
Binary alignment map
BRCA1-associated protein 1
Cyclin-dependent kinase inhibitor 2A
Copy number aberration
Generalized extreme studentized deviate
Hyperthermic intraperitoneal chemotherapy
Integer linear programming
Normothermic intraperitoneal chemotherapy
Polymerase chain reaction
SET domain containing 2
Single nucleotide variation
Synchronous precursor selection
Variant call format
Mao W, Zhang X, Guo Z, Gao Z, Pass HI, Yang H, Carbone M. Association of asbestos exposure with malignant mesothelioma incidence in Eastern China. JAMA Oncol. 2017;3:562–4.
Moolgavkar SH, Meza R, Turim J. Pleural and peritoneal mesotheliomas in SEER: age effects and temporal trends, 1973–2005. Cancer Causes Control. 2009;20:935–44.
Guo G, Chmielecki J, Goparaju C, Heguy A, Dolgalev I, Carbone M, Seepo S, Meyerson M, Pass HI. Whole-exome sequencing reveals frequent genetic alterations in BAP1, NF2, CDKN2A, and CUL1 in malignant pleural mesothelioma. Cancer Res. 2015;75:264–9.
Nasu M, Emi M, Pastorino S, Tanji M, Powers A, Luk H, Baumann F, Y-AA Z, Gazdar A, Kanodia S, Tiirikainen M, Flores E, Gaudino G, Becich MJ, Pass HI, Yang H, Carbone M. High incidence of somatic BAP1 alterations in sporadic malignant mesothelioma. J Thorac Oncol. 2015;10:565–76.
Yoshikawa Y, Emi M, Hashimoto-Tamaoki T, Ohmuraya M, Sato A, Tsujimura T, Hasegawa S, Nakano T, Nasu M, Pastorino S, Szymiczek A, Bononi A, Tanji M, Pagano I, Gaudino G, Napolitano A, Goparaju C, Pass HI, Yang H, Carbone M. High-density array-CGH with targeted NGS unmask multiple noncontiguous minute deletions on chromosome 3p21 in mesothelioma. Proc Natl Acad Sci U S A. 2016;113:13432–7.
Bueno R, Stawiski EW, Goldstein LD, Durinck S, De Rienzo A, Modrusan Z, Gnad F, Nguyen TT, Jaiswal BS, Chirieac LR, Sciaranghella D, Dao N, Gustafson CE, Munir KJ, Hackney JA, Chaudhuri A, Gupta R, Guillory J, Toy K, Ha C, Chen Y-J, Stinson J, Chaudhuri S, Zhang N, Wu TD, Sugarbaker DJ, de Sauvage FJ, Richards WG, Seshagiri S. Comprehensive genomic analysis of malignant pleural mesothelioma identifies recurrent mutations, gene fusions and splicing alterations. Nat Genet. 2016;48:407–16.
Hmeljak J, Sanchez-Vega F, Hoadley KA, Shih J, Stewart C, Heiman D, Tarpey P, Danilova L, Drill E, Gibb EA, Bowlby R, Kanchi R, Osmanbeyoglu HU, Sekido Y, Takeshita J, Newton Y, Graim K, Gupta M, Gay CM, Diao L, Gibbs DL, Thorsson V, Iype L, Kantheti H, Severson DT, Ravegnini G, Desmeules P, Jungbluth AA, Travis WD, Dacic S, et al. Integrative molecular characterization of malignant pleural mesothelioma. Cancer Discov. 2018;8:1548–65.
Bott M, Brevet M, Taylor BS, Shimizu S, Ito T, Wang L, Creaney J, Lake RA, Zakowski MF, Reva B, Sander C, Delsite R, Powell S, Zhou Q, Shen R, Olshen A, Rusch V, Ladanyi M. The nuclear deubiquitinase BAP1 is commonly inactivated by somatic mutations and 3p21.1 losses in malignant pleural mesothelioma. Nat Genet. 2011;43:668–72.
Ismail IH, Davidson R, Gagné J-P, Xu ZZ, Poirier GG, Hendzel MJ. Germline mutations in BAP1 impair its function in DNA double-strand break repair. Cancer Res. 2014;74:4282–94.
Bononi A, Giorgi C, Patergnani S, Larson D, Verbruggen K, Tanji M, Pellegrini L, Signorato V, Olivetto F, Pastorino S, Nasu M, Napolitano A, Gaudino G, Morris P, Sakamoto G, Ferris LK, Danese A, Raimondi A, Tacchetti C, Kuchay S, Pass HI, Affar EB, Yang H, Pinton P, Carbone M. BAP1 regulates IP3R3-mediated Ca2+ flux to mitochondria suppressing cell transformation. Nature. 2017;546:549–53.
Alakus H, Yost SE, Woo B, French R, Lin GY, Jepsen K, Frazer KA, Lowy AM, Harismendy O. BAP1 mutation is a frequent somatic event in peritoneal malignant mesothelioma. J Transl Med. 2015;13:122.
Chirac P, Maillet D, Leprêtre F, Isaac S, Glehen O, Figeac M, Villeneuve L, Péron J, Gibson F, Galateau-Sallé F, Gilly FN, Brevet M. Genomic copy number alterations in 33 malignant peritoneal mesothelioma analyzed by comparative genomic hybridization array. Hum Pathol. 2016;55:72–82.
Joseph NM, Chen Y, Nasr A, Yeh I, Talevich E, Onodera C, Bastian BC, Rabban JT, Garg K, Zaloudek C, Solomon DA. Genomic profiling of malignant peritoneal mesothelioma reveals recurrent alterations in epigenetic regulatory genes BAP1, SETD2, and DDX3X. Mod Pathol. 2017;30:246–54.
Kato S, Tomson BN, Buys TPH, Elkin SK, Carter JL, Kurzrock R. Genomic landscape of malignant mesotheliomas. Mol Cancer Ther. 2016;15:2498–507.
Sheffield BS, Tinker AV, Shen Y, Hwang H, Li-Chang HH, Pleasance E, Ch’ng C, Lum A, Lorette J, McConnell YJ, Sun S, SJM J, Gown AM, Huntsman DG, Schaeffer DF, Churg A, Yip S, Laskin J, Marra MA. Personalized oncogenomics: clinical experience with malignant peritoneal mesothelioma using whole genome sequencing. PLoS One. 2015;10:e0119689.
Singhi AD, Krasinskas AM, Choudry HA, Bartlett DL, Pingpank JF, Zeh HJ, Luvison A, Fuhrer K, Bahary N, Seethala RR, Dacic S. The prognostic significance of BAP1, NF2, and CDKN2A in malignant peritoneal mesothelioma. Mod Pathol. 2016;29:14–24.
Ugurluer G, Chang K, Gamez ME, Arnett AL, Jayakrishnan R, Miller RC, Sio TT. Genome-based mutational analysis by next generation sequencing in patients with malignant pleural and peritoneal mesothelioma. Anticancer Res. 2016;36:2331–8.
Leblay N, Leprêtre F, Le Stang N, Gautier-Stein A, Villeneuve L, Isaac S, Maillet D, Galateau-Sallé F, Villenet C, Sebda S, Goracci A, Byrnes G, McKay JD, Figeac M, Glehen O, Gilly FN, Foll M, Fernandez-Cuesta L, Brevet M. BAP1 is altered by copy number loss, mutation, and/or loss of protein expression in more than 70% of malignant peritoneal mesotheliomas. J Thorac Oncol. 2017;12:724–33.
Sugarbaker PH, Chang D. Long-term regional chemotherapy for patients with epithelial malignant peritoneal mesothelioma results in improved survival. Eur J Surg Oncol. 2017;43:1228–35.
Lee H-S, Jang H, Choi JM, Zhang J, de Rosen VL, Wheeler TM, Lee J, Tu T, Jindra PT, Kerman RH, Jung SY, Kheradmand F, Sugarbaker DJ, Burt BM. Comprehensive immunoproteogenomic analyses of malignant pleural mesothelioma. JCI insight. 2018;3:e98575.
Thapa B, Salcedo A, Lin X, Walkiewicz M, Murone C, Ameratunga M, Asadi K, Deb S, Barnett SA, Knight S, Mitchell P, Watkins DN, Boutros PC, John T. The immune microenvironment, genome-wide copy number aberrations, and survival in mesothelioma. J Thorac Oncol. 2017;12:850–9.
Calabrò L, Morra A, Fonsatti E, Cutaia O, Amato G, Giannarelli D, Di Giacomo AM, Danielli R, Altomonte M, Mutti L, Maio M. Tremelimumab for patients with chemotherapy-resistant advanced malignant mesothelioma: an open-label, single-arm, phase 2 trial. Lancet Oncol. 2013;14:1104–11.
Calabrò L, Morra A, Fonsatti E, Cutaia O, Fazio C, Annesi D, Lenoci M, Amato G, Danielli R, Altomonte M, Giannarelli D, Di Giacomo AM, Maio M. Efficacy and safety of an intensified schedule of tremelimumab for chemotherapy-resistant malignant mesothelioma: an open-label, single-arm, phase 2 study. Lancet Respir Med. 2015;3:301–9.
Maio M, Scherpereel A, Calabrò L, Aerts J, Perez SC, Bearz A, Nackaerts K, Fennell DA, Kowalski D, Tsao AS, Taylor P, Grosso F, Antonia SJ, Nowak AK, Taboada M, Puglisi M, Stockman PK, Kindler HL. Tremelimumab as second-line or third-line treatment in relapsed malignant mesothelioma (DETERMINE): a multicentre, international, randomised, double-blind, placebo-controlled phase 2b trial. Lancet Oncol. 2017;18:1–13.
Fennell DA, Kirkpatrick E, Cozens K, Nye M, Lester J, Hanna G, Steele N, Szlosarek P, Danson S, Lord J, Ottensmeier C, Barnes D, Hill S, Kalevras M, Maishman T, Griffiths G. CONFIRM: a double-blind, placebo-controlled phase III clinical trial investigating the effect of nivolumab in patients with relapsed mesothelioma: study protocol for a randomised controlled trial. Trials. 2018;19:233.
Calabrò L, Morra A, Giannarelli D, Amato G, D’Incecco A, Covre A, Lewis A, Rebelatto MC, Danielli R, Altomonte M, Di Giacomo AM, Maio M. Tremelimumab combined with durvalumab in patients with mesothelioma (NIBIT-MESO-1): an open-label, non-randomised, phase 2 study. Lancet Respir Med. 2018;2600:1–10.
DePristo MA, Banks E, Poplin R, Garimella KV, Maguire JR, Hartl C, Philippakis AA, del Angel G, Rivas MA, Hanna M, McKenna A, Fennell TJ, Kernytsky AM, Sivachenko AY, Cibulskis K, Gabriel SB, Altshuler D, Daly MJ. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat Genet. 2011;43:491–8.
Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010;38:e164.
Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative genomics viewer (IGV): high-performance genomics data visualization and exploration. Brief Bioinform. 2013;14:178–92.
Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, Batut P, Chaisson M, Gingeras TR. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21.
Anders S, Pyl PT, Huber W. HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31:166–9.
Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11:R106.
Hughes CS, Foehr S, Garfield DA, Furlong EE, Steinmetz LM, Krijgsveld J. Ultrasensitive proteome analysis using paramagnetic bead technology. Mol Syst Biol. 2014;10:757.
Hughes CS, McConechy MK, Cochrane DR, Nazeran T, Karnezis AN, Huntsman DG, Morin GB. Quantitative profiling of single formalin fixed tumour sections: proteomics for translational research. Sci Rep. 2016;6:34949.
Szklarczyk D, Franceschini A, Wyder S, Forslund K, Heller D, Huerta-Cepas J, Simonovic M, Roth A, Santos A, Tsafou KP, Kuhn M, Bork P, Jensen LJ, von Mering C. STRING v10: protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2014;43:D447–52.
Shrestha R, Hodzic E, Sauerwald T, Dao P, Wang K, Yeung J, Anderson S, Vandin F, Haffari G, Collins CC, Sahinalp SC. HIT’nDRIVE: patient-specific multidriver gene prioritization for precision oncology. Genome Res. 2017;27:1573–88.
Yoshihara K, Shahmoradgoli M, Martínez E, Vegesna R, Kim H, Torres-Garcia W, Treviño V, Shen H, Laird PW, Levine DA, Carter SL, Getz G, Stemke-Hale K, Mills GB, RGW V. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.
Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, Hoang CD, Diehn M, Alizadeh AA. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12:453–7.
Mermel CH, Schumacher SE, Hill B, Meyerson ML, Beroukhim R, Getz G. GISTIC2.0 facilitates sensitive and confident localization of the targets of focal somatic copy-number alteration in human cancers. Genome Biol. 2011;12:R41.
Yu H, Pak H, Hammond-Martel I, Ghram M, Rodrigue A, Daou S, Barbour H, Corbeil L, Hebert J, Drobetsky E, Masson JY, Di Noia JM, Affar EB. Tumor suppressor and deubiquitinase BAP1 promotes DNA double-strand break repair. Proc Natl Acad Sci. 2014;111:285–90.
Mouw KW, Goldberg MS, Konstantinopoulos PA, D’Andrea AD. DNA damage and repair biomarkers of immunotherapy response. Cancer Discov. 2017;7:675–93.
Alley EW, Lopez J, Santoro A, Morosky A, Saraf S, Piperdi B, van Brummelen E. Clinical safety and activity of pembrolizumab in patients with malignant pleural mesothelioma (KEYNOTE-028): preliminary results from a non-randomised, open-label, phase 1b trial. Lancet Oncol. 2017;18:623–30.
Shrestha R, Le Bihan S, Collins C: Whole exome sequencing and whole transcriptome sequencing of malignant peritoneal mesothelioma. Data sets. figshare. 2019 doi: https://doi.org/10.6084/m9.figshare.7653500
Shrestha R, Le Bihan S, Collins C: Mass spectrometry of malignant peritoneal mesothelioma. Data sets. figshare. 2019 doi: https://doi.org/10.6084/m9.figshare.7653503
We acknowledge the contribution and support of Drs. Jessica McAlpine, Anna Tinker, and Jeff Simko as well as Emily Taylor (Mount Sinai Hospital), Donald Donaldson, Rose Schweigert, Matthew Sturgen, Sarah Padilla at Vancouver Coastal Health for supporting and coordinating the procurement of PeM tumors. We are thankful for the support of regulatory and ethics bodies: Margaret Luk, Wylo Kyle, Jacqueline Lee, Zahra Karim, Nandita Chowdhury, Jessica Gagliardi, Suzanne Richardson, and Sheila O’Donoghue at Vancouver General Hospital; Chantal Lackan at Mount Sinai Hospital; Ida Deichaite and Oudone Sisanachandeng at Moores Cancer Centre. Thanks to Shane Colborne and Dr. Christopher Hughes for help and suggestions regarding mass spectrometry experiments and data analyses. The authors thank all members of the Collins’, Wang’s, Hach’s, and Sahinalp’s labs for helpful suggestions. The results published here are in part based upon data generated by the TCGA Research Network: http://cancergenome.nih.gov/.
This study is funded by: BC Cancer Foundation, Mitacs, Canadian Institutes of Health Research (CIHR), and Terry Fox Research Institute (C.C.C. and Y.Z.W.). R.S. and N.N. is supported by Mitacs Accelerate Awards.
Availability of data and materials
The whole exome and whole transcriptome sequencing data from this study is available in the European Genome-phenome Archive (EGA) under accession number EGAS00001002820 (https://www.ebi.ac.uk/ega/studies/EGAS00001002820 and https://doi.org/10.6084/m9.figshare.7653500) . The proteome data from mass spectrometry is available in the PRIDE Archive under accession number PXD008873 (http://proteomecentral.proteomexchange.org/cgi/GetDataset?ID=PXD008873 and https://doi.org/10.6084/m9.figshare.7653503) .
Ethics approval and consent to participate
This study was approved by the Institutional Review Board of the University of British Columbia and the Vancouver Coastal Health (REB Number. H1500902 and V15-00902). All samples and information were collected with written and signed informed consent from the participating patients. This research study is conformed to the Declaration of Helsinki.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.