Association of MUC16 Mutation With Response to Immune Checkpoint Inhibitors in Solid Tumors

Key Points Question What is the association between MUC16 mutation and response to immune checkpoint inhibitors (ICIs) in solid tumors? Findings In this cohort study of 3 groups of patients, including The Cancer Genome Atlas cohort with 10 195 patients across 30 solid tumor types, 56 patients with non–small cell lung cancer, and 145 patients with melanoma, MUC16 mutation was associated with greater response rates, prolonged overall survival, and genomic factors associated with ICI response, including tumor mutational burden and programmed cell death ligand–1 expression. Meaning These findings suggest that MUC16 mutation may be associated with superior response to ICIs in patients with solid tumors.

As the third most frequently mutated gene according to The Cancer Genome Atlas (TCGA), MUC16 (OMIM 606154) encodes a membrane-spanning mucin, comprising a tandem repeat region sandwiched between N-terminal and C-terminal domains. 3,4 Tumor antigen 125, residing in the tandem repeat region, is a common clinical biomarker to monitor epithelial ovarian cancer. 3 A recent study 5 reported that MUC16 mutation may be associated with elevated tumor mutational burden and superior overall survival (OS) in patients with gastric adenocarcinoma. However, a comprehensive analysis of the relationship of MUC16 mutation with ICI response across a large set of solid tumors is lacking. Here, we explored the associations of MUC16 mutation with ICI response based on multidimensional data from multiple solid tumors.

Multidimensional Data Collection
This study was exempted from approval by an institutional review board and from the need for informed consent because its data were gathered from publicly available data sets that have received approval, in accordance with 45 CFR §46. This study follows the Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) reporting guideline.
Somatic mutations, neoantigens, and gene expression data for TCGA samples from multiple solid tumor types were gathered. Specifically, somatic mutations of 10 195 samples and mRNA expression profiles of 9850 samples for 30 solid tumor types were downloaded from the TCGA Pan-Cancer Atlas project deposited in the cBioPortal database. 6,7 These somatic mutations were uniformly called by the Multi-Center Mutation Calling in Multiple Cancers working group to adjust for variance and batch effects introduced by the differences in DNA extraction, hybridization-capture, sequencing, and analytical methods over time, which enables robust cross-tumor-type analyses. 8 The mRNA expression levels were quantified using RNA-Seq by Expectation-Maximization algorithm 9 ; batch-corrected for sequencers (Illumina GAII and HiSeq), sequencing centers (University of North Carolina and British Columbia Cancer Agency), and a plate effect identified in prostate adenocarcinoma; and normalized using the upper quartile method by Hoadley et al. 10 Neoantigens of 5935 samples across 19 solid tumor types were obtained from The Cancer Immunome Atlas. 11 Abundance data of 64 immune and stromal cell types in the tumor microenvironment for 9358 TCGA tumor samples were inferred using xCell 12 and downloaded from the xCell website. xCell is a gene signatures-based method that uses a compensation technique to reduce spillover effects between closely related cell types. Overlapping associations of somatic mutations data with mRNA expression, neoantigens, and cell abundance data are shown in eFigure 1 in the Supplement.

Clinical Cohorts Treated With ICIs
Uniformly analyzed clinical cohorts would help to more robustly evaluate biomarkers of response to ICI therapy. 13 Thus, we obtained uniformly processed somatic mutations and clinical information of pretreatment tumors from patients with metastatic lung cancer or melanoma receiving ICI therapies.
Briefly, all these patients were gathered from 5 published cohorts. [13][14][15][16][17] To avoid the bias introduced by different studies, all raw whole-exome sequencing data (BAM files) and clinical annotations were processed through standardized quality control and mutation calling pipelines as described by Miao et al. 13 We excluded 1 patient with small cell lung cancer from the lung cancer cohort with 57 patients and thus obtained 56 patients with non-small cell lung cancer (NSCLC) treated with anti-PD-1/PD-L1 as our NSCLC cohort. The melanoma cohort initially consisted of 151 patients treated with anti-CTLA-4, anti-PD-1, anti-PD-L1, or a combination of these therapies. Given the substantial differences in response rates and survival among patients treated with anti-CTLA-4, anti-PD-1, and anti-CTLA-4 plus anti-PD-1, 18 we focused on 145 patients treated with anti-CTLA-4 alone to avoid potential bias introduced by different therapy classes.
With regard to the evaluation of clinical response to ICI therapy, several response definitions, 15,19,20

MUC16 Mutation and Tumor Mutational Burden
All nonsynonymous somatic mutations in MUC16, including missense, nonsense, nonstop, frameshift insertion and deletion, in-frame insertion and deletion, and splice site mutations, were considered. Tumor mutational burden was calculated as the total count of nonsynonymous mutations in coding sequence using maftools. 23

Association of MUC16 Mutation With Immune-Related Genes Expression Profile and Tumor Immune Microenvironment
First, we collected 40 immune-related genes, which were classified into 3 categories: immune checkpoint, T-effector and interferon-γ gene signature, and T cell receptor (eTable 1 in the Supplement). These genes were retrieved from previous reports. [24][25][26] We next extracted normalized expression levels of these genes from 9580 samples with both somatic mutations and expression data available. The expression levels were log 2 -transformed after adding an offset of 1 to avoid taking log of 0 before analysis. In addition, we stratified all samples into 4 tumor immune microenvironment (TIME) types on the basis of median expression levels of CD8A (OMIM 186910) and PD-L1 (OMIM 605402) according to the previous report, 27 which defined the positive as gene expression above the median level. On the basis of 8755 samples with both somatic mutations and cell abundances data, we further evaluated the extent of immune cells infiltration into TIME by generating a composite immune score as the sum of abundances for B cells, CD4 + T cells, CD8 + T cells, dendritic cells, eosinophils, macrophages, monocytes, mast cells, neutrophils, and natural killer cells, as described elsewhere. 12

Survival Analysis
A Kaplan-Meier estimator of OS was used to construct the survival curves. The log-rank test was performed to compare OS of patients stratified by MUC16 mutational status, and hazard ratios (HRs) with 95% CIs are provided. For the advanced melanoma cohort, complete clinical characteristics available included age at start of treatment and sex; 5 dominant COSMIC mutational signatures 28 (signature 7, exposure to UV light; signature 5, unknown environmental exposure; signature 11, prior chemotherapeutic treatment with alkylating agents; signature 1, age-related spontaneous deamination of 5-methyl-cytosine; and signature 6, defective DNA mismatch repair) were also derived from the previous study, 13 where mutational burden was not associated with ICI response after adjusting these dominant signatures, suggesting that it is a surrogate for an underlying mutagenic biological process. In melanoma, a large proportion of somatic mutations are known to arise from exposure to UV light. 29 Finally, a Cox proportional hazards model was fitted for OS, controlling for potential confounding factors (age, sex, and mutational signatures, except signature 1) and ruling out multicollinearity. Multicollinearity among variables was determined by variance inflation factor calculated using the rms package in R statistical software version 3.5.2 (R Project for Statistical Computing). All survival analyses were performed using the survival package in R statistical software.

Gene Set Enrichment Analysis
Gene set enrichment analysis (GSEA) was performed using the GSEA Desktop Application (version 3.0). 30,31 Specifically, genes with 0 count in more than 40% of samples were first removed, and then log 2 -transformed fold changes for expression of qualified genes in MUC16-mutated vs wild-type states were calculated as a ranking metric for GSEAPreranked to perform GSEA against 50 hallmark gene sets annotated by the Molecular Signatures Database (version 6.2). 32 This analysis involved 100 000 random permutations for gene set and weighted enrichment statistic. The normalized enrichment score was used as the magnitude of enrichment. The proportion of false-positives was controlled by calculating the false discovery rate (FDR). The FDR estimates the probability that a gene set with a given normalized enrichment score represents a false-positive finding. A significantly enriched gene set was expected at FDR < .001.

Statistical Analysis
The Mann-Whitney U test was used to determine the association between MUC16 mutation and a continuous variable. Fisher exact test and χ 2 test were used to compare the frequencies of clinical response to ICIs and TIME types in MUC16-mutated vs wild-type tumors, respectively. The difference in mean mRNA expression of each immune-related gene in MUC16-mutated vs wild-type patients for different cancer types was hierarchically clustered with the pheatmap package in R statistical software, using Euclidean distance as a distance metric for samples, and ward.D2 clustering method.
Statistical significance was expected at 2-tailed P < .05 unless otherwise specified. Statistical analyses were done in R statistical software version 3.5.2 (R Project for Statistical Computing). Data were obtained from October 1 through October 10, 2019, and were analyzed from October 11 through December 31, 2019. To investigate the association of MUC16 mutation with established genomic factors associated with ICI response, we began with tumor mutational burden, because it has been the most widely reproduced biomarker for ICI therapy. 33 We found that, in this pan-cancer data set, patients with MUC16 mutations exhibited significantly higher tumor mutational burden than those without them   (Figure 2A). To confirm the association between MUC16 mutation and TIME, we     Figure 2C, 2D, and 2E). The differences in mean expression levels of these genes in MUC16-mutated vs wild-type tumors from the respective cancer types are presented in eFigure 4 in the Supplement.

MUC16 Mutation Associated With Favorable Outcomes for ICI Treatment
We further examined the association of MUC16 mutation with outcomes in 2 ICI-treated cohorts, hazard ratio, 0.34; 95% CI, 0.12-0.99; P = .04, log-rank test) ( Figure 3A). There were 18 deaths in the NSCLC cohort.
We next confirmed the association between MUC16 mutation and OS in the melanoma cohort.
As shown in Figure 3B, significant OS improvement was also observed in patients with melanoma with MUC16 mutation (median, 27.03 vs 9.73 months; hazard ratio, 0.57; 95% CI, 0.36-0.90; P = .02, log-rank test). We used linear regression model to analyze the association of tumor mutational burden with neoantigen load and observed that each nonsynonymous mutation produced, on average, 2.79 estimated neoantigens, with a strong correlation between them (R 2 = 0.98;   (Figure 5A and 5B).

Enriched Biological Processes in MUC16-Mutated Tumors
To identify biological processes associated with MUC16 mutational status, we conducted GSEA on 50 hallmark gene sets, representing major biological processes, for 9580 patients with and without MUC16 mutations. Our analysis yielded 8 gene sets whose expression was significantly upregulated in patients with MUC16 mutations (FDR < .001) (eFigure 10 in the Supplement). Upon inspecting these gene sets, we were able to assign them to 3 biological themes involved in cell proliferation (E2F targets, G2/M checkpoint, MYC targets variant 1, and MYC targets variant 2), immune response (interferon-α response, interferon-γ response, and allograft rejection), and mTORC1 signaling.

Discussion
Our analysis of the TCGA pan-cancer data set across multiple solid tumor types demonstrated that suggesting that there exists potential adaptive immune resistance to anti-PD-1/PD-L1 therapies and that additional inhibitory pathways beyond the PD-1/PD-L1 axis might be targetable. 37 Taken together, these findings support the hypothesis that MUC16 mutation is associated with high immunogenicity and a responsive TIME with PD-L1-dependent or PD-L1-independent adaptive immune resistance. 38 Consistent with the hypothesis, to our knowledge, we first confirmed the association of MUC16 mutation with superior outcomes in cohorts receiving ICI treatment. MUC16 mutation was associated with improved OS and response rates in anti-PD-1/PD-L1-treated patients with NSCLC, and this finding was verified in the independent melanoma cohort treated with anti-

CTLA-4.
MUC16 has previously been implicated in tumor cell proliferation, metastasis, and modulation of the innate immune response by direct suppression of natural killer cell function. 3,4 Also, our GSEA analysis revealed that 8 biological processes regarding cell proliferation, immune response, and mTORC1 signaling were significantly upregulated in MUC16-mutated tumors, providing biological insights concerning the link between MUC16 mutation and ICI response.

Limitations
Our study has several limitations. First, although MUC16-mutated solid tumors were characterized by alteration of pathways involved in cell proliferation, immune response, and mTORC1 signaling, mechanistic underpinnings of MUC16 mutation relevant to ICI response remain elusive and merit further experimental work. Second, the sample size of both ICI-treated cohorts was limited. The limited number of deaths (18 deaths) in the NSCLC cohort restricted the ability to controlling