Isolated BAP1 Genomic Alteration in Malignant Pleural Mesothelioma Predicts Distinct Immunogenicity with Implications for Immunotherapeutic Response

Simple Summary Malignant Pleural Mesothelioma (MPM) is a highly aggressive, therapy-resistant cancer with a well-established inflammatory etiology and has no cure. Immune checkpoint inhibition (ICI) therapy has shifted treatment paradigms in many types of cancers, and recent clinical data have shown promise for improving MPM treatment. However, response to ICI therapy has been neither uniform nor predictable. The genomic landscape of MPM is primarily characterized by genomic alterations in tumor suppressor genes (TSGs) (~70%), particularly BAP1, CDKN2A/B and NF2. The impact of an isolated TSG genomic alteration versus multiple concurrent TSG alterations on clinical outcome, treatment response and MPM biology and the immune tumor microenvironment are unclear. Here, we showed the effect of TSG alteration combinations on clinical outcome, therapeutic response, and molecular pathways in MPM. For example, tumors with alterations in BAP1 alone were (a) associated with a longer overall patient survival rate compared to tumors with CDKN2A/B and/or NF2 alterations with or without BAP1 and (b) comprised a distinct immunogenic subtype with altered transcription factor and pathway activity patterns. Abstract Malignant pleural mesothelioma (MPM), an aggressive cancer of the mesothelial cells lining the pleural cavity, lacks effective treatments. Multiple somatic mutations and copy number losses in tumor suppressor genes (TSGs) BAP1, CDKN2A/B, and NF2 are frequently associated with MPM. The impact of single versus multiple genomic alterations of TSG on MPM biology, the immune tumor microenvironment, clinical outcomes, and treatment responses are unknown. Tumors with genomic alterations in BAP1 alone were associated with a longer overall patient survival rate compared to tumors with CDKN2A/B and/or NF2 alterations with or without BAP1 and formed a distinct immunogenic subtype with altered transcription factor and pathway activity patterns. CDKN2A/B genomic alterations consistently contributed to an adverse clinical outcome. Since the genomic alterations of only BAP1 was associated with the PD-1 therapy response signature and higher LAG3 and VISTA gene expression, it might be a candidate marker for immune checkpoint blockade therapy. Our results on the impact of TSG genotypes on MPM and the correlations between TSG alterations and molecular pathways provide a foundation for developing individualized MPM therapies.


Introduction
Malignant pleural mesothelioma (MPM) is a rare, devastating cancer of the lining of the lung and thoracic cavity with a 5-year survival rate of <10%. About 70% of MPM cases

Association of Gene Signature Predictive of Response to Therapy with TSG Genotypic Groups in MPM
We performed gene set variation analysis (GSVA) [19] using published gene expressionbased treatment-response signatures to determine differences in patient responses to standard pemetrexed chemotherapy and palbociclib based on TSG genotypes for the TCGA MPM dataset since parallel RNA-seq (RNA-sequencing) data was not available for theMSK-IMPACT cohort. We found an expression-based signature derived from non-small cell lung cancer that predicted resistance to pemetrexed [20] (Supplementary Materials Table S5). We collected a signature of resistance to palbociclib derived from breast cancer [21] (Supplementary Materials Table S6) to determine differences in sensitivity to cyclin-dependent kinase inhibitors among the eight TSG genotypic groups. Figure 3A shows GSVA scores for the two signatures across the MPM TSG tumors grouped by TSG genotypic types (see Materials and Methods). On average, MPM tumors with CDKN2A/B alteration, with or without BAP1 or NF2 alterations, were more resistant to pemetrexed and palbociclib, whereas tumors with only BAP1 loss or with BAP1 and NF2 losses were more sensitive (a lower GSVA score).

Groups in MPM
We performed gene set variation analysis (GSVA) [19] using published gene expression-based treatment-response signatures to determine differences in patient responses to standard pemetrexed chemotherapy and palbociclib based on TSG genotypes for the TCGA MPM dataset since parallel RNA-seq (RNA-sequencing) data was not available for theMSK-IMPACT cohort. We found an expression-based signature derived from nonsmall cell lung cancer that predicted resistance to pemetrexed [20] (Supplementary Materials Table S5). We collected a signature of resistance to palbociclib derived from breast cancer [21] (Supplementary Materials Table S6) to determine differences in sensitivity to cyclin-dependent kinase inhibitors among the eight TSG genotypic groups. Figure 3A shows GSVA scores for the two signatures across the MPM TSG tumors grouped by TSG genotypic types (see Materials and Methods). On average, MPM tumors with CDKN2A/B alteration, with or without BAP1 or NF2 alterations, were more resistant to pemetrexed and palbociclib, whereas tumors with only BAP1 loss or with BAP1 and NF2 losses were more sensitive (a lower GSVA score).  We next compared immune checkpoint mRNA expression levels as a function of TSG genotypes and found that expression of LAG3 (lymphocyte activation gene-3) and VISTA (V-domain Ig suppressor of T cell activation; also known C10orf54) was higher in B∆ tumors ( Figure 3B), but mRNA levels for other immune checkpoint genes, including PD-1, PD-L1, TIGIT, and CTLA4, were not associated with TSG genotypes (Supplementary Materials Figure S1).
To evaluate the connection between TSG genotypes and resistance to PD-1 therapy, we also tested a gene expression signature predictive of benefit from immune checkpoint inhibitor (ICI) treatment [22] (Supplementary Materials Table S7). Briefly, Jang et al. [22] generated the signature based on MPM patients treated with nivolumab (4 responders and 4 non-responders) and then they used it to subgroup TCGA samples into the anti-PD-1responsive and anti-PD-1-resistant. Strikingly, we found that B∆ tumors were represented only in the anti-PD-1-responsive subgroup (p = 0.0032), whereas N∆C∆ tumors were observed only in the anti-PD-1 resistant subgroup (p = 0.0243) regardless of their histology ( Figure 3C, Supplementary Materials Table S8). We also applied CIBERSORTx to 86 TCGA tumors and inferred the infiltration of immune and tumor stromal cells using gene expression data. There was a lower proportion of B cells (p-value = 0.016, one-way ANOVA) and a higher proportion of natural killer (NK) (p-value = 0.031, one-way ANOVA) in B∆ tumors compared to other TSG genotype groups (Supplementary Materials Figure S2). However, we did not observe differences for other cell types including CD8+ T cells. In summary, tumors with only BAP1 alteration (B∆) without other frequent TSG alterations formed a distinct MPM subtype that was associated with significantly longer overall patient survival, a better therapeutic response signature (e.g., anti-PD-1 and pemetrexed), and higher expression of LAG3 and VISTA.

Transcription Factors and Pathways Associated with TSG Genotypic Groups in MPM
To determine whether the eight MPM TSG genotypes impacted similar or distinct molecular pathways, we performed sample-specific transcription factor activity analysis via the Integrated System for Motif Activity Response Analysis (ISMARA) [23] and sample-specific pathway enrichment analysis via GSVA [19] using gene expression data based on TCGA RNA-seq data. Figure 4A summarizes mean transcription factor (TFs) activities (false discovery rate [FDR] < 0.05) and mean pathway scores (FDR < 0.05) significantly associated with TSG genotypes. Tumors with CDKN2A alteration including B∆N∆, B∆N∆C∆ and N∆C∆ were associated with increased activity of E2F Targets (related to cell cycle) ( Figure 4A,B). Further, B∆N∆C∆, B∆N∆, C∆, and N∆ tumors were associated with increased activity of epithelial-mesenchymal transition (EMT) ( Figure 4A,B). Whereas B∆ tumors were associated with increased activity of interferon regulatory factors (IRFs) that have a role in immunity and decreased activity of BCL6B (also known as B-cell CLL/lymphoma 6 member B), which is a transcriptional repressor that interacts with the Notch, STAT, p53 and PI3K/AKT signaling pathways, all of which may be involved in inflammatory response regulation in cancer cells [24][25][26] (Supplementary  Materials Table S9). Consistent with the TF activity patterns, B∆ tumors were associated with increased interferon-alpha and interferon-beta signaling activity. Genes that showed increased expression that correlated with IRF TF activities and pathway score included interferon response genes, such as IFIT3, OAS2, OAS3, IFIH1, STAT1, DDX60, and DDX58; genes regulating the CGAS-induced type I interferon signaling pathway, such as TRIM14; and transporters associated with antigen processing such as TAP1 and TAP2 ( Figure 4C). In summary, we identified a group of MPM patients with loss-of-function mutations and/or copy number loss only in B∆ tumors, which define a distinct molecular subtype associated with a high interferon response. loss only in B∆ tumors, which define a distinct molecular subtype associated with a high interferon response.

Discussion
Precision oncology, which has been successful in treating a variety of cancers, is still in its infancy in MPM. While there are ongoing clinical trials for MPM, there are no defined molecular markers (e.g., genomic markers) that can be used to predict the efficacy of treatments. Genomic alterations of BAP1, NF2, and CDKN2A/B TSGs are thought to play a critical role in MPM pathogenesis. In this study, we provide a comprehensive analysis of data on MPM tumors with genomic alteration in one or more of these specific TSGs, which are thought to be critical drivers of MPM pathogenesis.
Tumors with alterations in only BAP1 showed a distinct pattern of expression of inflammatory tumor microenvironment genes, including activation of interferon signaling and IRF TFs and high LAG3 and VISTA expression. Interferon production is a defense response that recruits and activates immune cells and has been studied extensively in cancer. Activation of the interferon response in cancer cells [27], possibly by genomic instability through the cGAS-STING pathway [28], may affect immune cells in the tumor microenvironment and the tumor response to immunotherapies. In tumors, type I interferons are secreted by cancer cells and dendritic cells (DCs) in response to DNA fragments that activate the cGAS/STING pathway and result in T cell priming and antitumor activity. Thus, isolated BAP1 alteration may serve as a predictive and prognostic candidate biomarker for MPM to improve disease stratification and therapy. Interestingly, BAP1 alterations have recently been shown to be correlated with perturbed immune signaling in malignant peritoneal mesothelioma [29]. In addition, in vitro and in vivo studies have shown type I interferon (IFN-I) activation in mesothelioma cells [30,31]. Further, Hmeljak et al. [12] reported an IFN-I signature in pleural mesothelioma tumors with an inactivated BAP1 gene. Further, Yang et al. showed a negative correlation between BAP1 expression and a constitutively activated IFN-I response [32]. Regarding possible future validation experiments to confirm the tumor-promoting roles of IRF TFs and interferon signaling in B∆ MPM, expression levels of IRF TFs can be manipulated in cultured cells through overexpression or silencing, and the effects of IRF TFs on cancer cell proliferation, survival, motility, and invasion potential can then be evaluated. Ultimately, in vivo validation of the roles of IRF TFs expression on MPM progression can be studied using mouse models of mesothelioma crossed with mice harboring knockout of specific IRF TF genes [33].
Harnessing the combination of different immunotherapy approaches to improve outcomes of patients with MPM is an area of clinical interest [34]. We observed that expression of immune checkpoint genes LAG3 and VISTA were higher in B∆ tumors. VISTA is a member of the B7 family of B7-CD28 family of ligands and receptors that is expressed primarily on myeloid cells and T-lymphocytes [35,36]. When overexpressed, VISTA suppresses early T-cell activation and proliferation and reduces cytokine production. Hmeljak et al. [12] reported strong expression of VISTA in benign mesothelium by immunohistochemistry (IHC) and increased mRNA expression of VISTA in epithelioid MPM compared to other tumor types. Muller et al. [37] also confirmed VISTA expression by IHC in malignant pleural mesotheliomas in both tumor and infiltrating inflammatory cells. LAG3 is a coinhibitory receptor expressed on activated T cells and has now become part of the repertoire of combinatorial immunotherapeutics available for the treatment of metastatic melanoma [38]. Recently, Marcq et al. [39] reported monotherapy with PD-L1 and its combination with LAG-3 blockade, resulted in delayed tumor growth and significant survival benefit in malignant mesothelioma mouse models. Our study also provides a candidate biomarker for preclinical studies of PD-1 combination therapy with a VISTA or LAG3 inhibitor in B∆ MPM tumors.
A major limitation of our study is our small cohort size. Further, there are several differences in percentages for the different TSG status between the TCGA and the MSK-IMPACT datasets. One reason might be due to the fact that the TCGA data was based on whole exon sequencing whereas MSK-IMPACT was based on a targeted screen, which may have identified a higher percentage of whole exon deletions of BAP1. Notably, in early reports, performed with Sanger sequencing, revealed point mutations in 20-25% of Cancers 2022, 14, 5626 9 of 12 sporadic MPMs [40,41]. Subsequent deletion mapping and multiplex ligation-dependent probe amplification studies have identified alterations of BAP1 in 50-60% of MPMs, with the increase due to inactivating deletions of entire BAP1 exons [42,43]. Despite the limitations, our results indicate that isolated BAP1 could be candidate biomarker for selection of patient to immune checkpoint blockade therapies.

Data and Preprocessing
TCGA data were downloaded through the Broad Institute TCGA GDAC firehose tool. RNA-seq data were available for 86 samples. Genetic alteration data (copy number alteration and mutation) for CDKN2A/B, NF2, and BAP1 were retrieved from cBioPortal and clinical data from an online portal for data from the TCGA project and MSK-IMPACT [18].

Survival Analyses
Overall survival (OS) was defined as the time from diagnosis to death resulting from any cause. Survival curves were estimated using the Kaplan-Meier method. Hazard ratios (HR) were estimated using a Cox proportional hazards regression model. For each of the eight TSG genotypes, a separate univariate regression model was created with one binary variable indicating if a tumor had the exact genotype. To model the effect of genotypes on survival without interactions, a Cox proportional hazards regression model was created with three binary variables, one for each: CDKN2A/B, BAP1, and NF2. p-values were obtained for each coefficient in the Cox regression models with a Wald-test, with the null hypothesis that the HR is 1. Both the Cox regression and Kaplan-Meier analyses were done with the survival R package [44].

Motif Activity Analysis
To analyze activities of transcription factor binding motifs (TFBM) from TCGA MPM RNA-seq data, we used ISMARA [23].

Gene Set Enrichment Analysis
GSVA and single-sample gene set enrichment analysis (ssGSEA) were performed using the GSVA R package (version 1.40.1) [19] on the pemetrexed and palbociclib response signature (Supplementary Materials Tables S5 and S6) and pathway enrichment analysis. For pathway enrichment analysis, we obtained pathway annotations from the Molecular Signatures Database (MsigDB) [45], a collection of hallmarks of cancer and REACTOME pathways (c2.all.v7.1.symbols.gmt). The log-transformed TMM (trimmed mean of M values) normalized TPM (transcripts per million) counts based on TCGA MPM RNA-seq data were used as input to the GSVA package.

Analysis of the Tumor Immune Microenvironment (TIME)
The proportion of different immune, tumor and stromal cells in the tumor microenvironment was estimated from RNA-seq data using CIBERSORTx [46]. The CIBESORTx algorithm was run with default settings, excluding quantile normalization, for 100 permutations with our signature matrix based MPM CITE-seq data [47] from one tumor sample (GSE172155) to estimate the abundance of immune, stromal and malignant cells types (Supplementary Materials Figure S2). Then, we evaluated the association of cell types with each TSG genotype using one-way ANOVA.

Statistical Analysis and Visualization
All statistical tests in the exploratory analysis were performed using R version 4.1.1 and associated packages. The statistical analyses for differences in mRNA expression, GSVA score, and TF activity across TSG genotype groups were performed using one-way ANOVA with an FDR cut off of 0.05. Further, we used one-way ANOVA and post hoc Tukey's HSD (honestly significant difference), with an adjusted p-value cut off of 0.1 for pairwise comparison of genotypes for visualization with boxplots. One or more letters were assigned to each TSG genotype group using the multcompView package in R (version: 0.1-8). Assignments were made such that any two samples that had a statistically significant difference did not share any letters.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/cancers14225626/s1, Figure S1: Comparison of immune checkpoint gene mRNA expression levels as a function of TSG genotypes in 86 MPM samples from the TCGA cohort; Figure S2: CIBERSORTx analysis of TCGA MPM RNA-seq dataset with MPM scRNA-seq reference of ten major cell types; Table S1: Comparison of pleural mesothelioma histological subtypes as a function of TSG genotypes in 86 MPM samples from the TCGA cohort; Table S2: Comparison of pleural mesothelioma histological subtypes as a function of TSG genotypes in 61 MPM samples from the MSK-IMPACT cohort; Table S3: Multivariate cox regression analysis based on BAP1, NF2 and CDKN2A/B status for TCGA MPM dataset; Table S4: Multivariate cox regression analysis based on BAP1, NF2 and CDKN2A/B status for MSK-IMPACT dataset; Table S5: Genes for Pemetrexed response signature; Table S6: Genes for Palbociclib response signature; Table S7: Genes for anti-PD-1 resistance signature; Table S8: The anti-PD-1-resistant mRNA signature was used to predict the subgroups; Table S9: Candidate TF regulators (5% FDR) based on Figure 4A. References [20][21][22]47]   Data Availability Statement: Publicly available datasets were analyzed in this study. This data can be found here: (1) RNA-seq gene expression data from TCGA's Firehose data run (https:// confluence.broadinstitute.org/display/GDAC/Dashboard-Stddata, accessed on 28 September 2022); (2) Genetic alteration data (copy number alteration and mutation) and clinical data were retrieved from cBioPortal (http://www.cbioportal.org) [48] for both the TCGA project and MSK-IMPACT [18].