Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma (original) (raw)

Cell. Author manuscript; available in PMC 2017 Mar 24.

Published in final edited form as:

PMCID: PMC4808437

NIHMSID: NIHMS765463

Willy Hugo,1,6,9 Jesse M. Zaretsky,2,6,9 Lu Sun,1,6 Chunying Song,1,6 Blanca Homet Moreno,3 Siwen Hu-Lieskovan,3 Beata Berent-Maoz,3 Jia Pang,3 Bartosz Chmielowski,3 Grace Cherry,3 Elizabeth Seja,3 Shirley Lomeli,1,6 Xiangju Kong,1,6 Mark C. Kelley,7 Jeffrey A. Sosman,8 Douglas B. Johnson,8 Antoni Ribas,2,3,4,5,6 and Roger S. Lo1,2,5,6,*

Blanca Homet Moreno

3

Siwen Hu-Lieskovan

3

Beata Berent-Maoz

3

Bartosz Chmielowski

3

Grace Cherry

3

Elizabeth Seja

3

Mark C. Kelley

7

Jeffrey A. Sosman

8

Douglas B. Johnson

8

1

2

3

4

5

6

7

8

9Co-first authors

*Correspondence to: Dr. Roger S. Lo at 52-121 CHS Dept. of Medicine/Dermatology, 10833 Le Conte Ave, Los Angeles, CA 90095-1750 or at ude.alcu.tendem@olr.

Supplementary Materials

1.

GUID: F9E11DFC-6E62-4B8E-B39E-877E1581111E

2.

GUID: 8C385867-B26D-44E4-82D6-E13A9D1E6F40

3.

GUID: E791AE88-FB7B-4AAD-BDB2-52FE3B495468

4.

GUID: 3D993A5D-8FD4-4B5C-9496-4005AD8FB1CE

5.

GUID: E56B23B2-8668-43E5-A7DA-F9ABEAED7CA1

6.

GUID: 422A88FF-4598-4850-A7D4-9B2375B95EA6

SUMMARY

PD-1 immune checkpoint blockade provides significant clinical benefits for melanoma patients. We analyzed the somatic mutanomes and transcriptomes of pretreatment melanoma biopsies to identify factors that may influence innate sensitivity or resistance to anti-PD-1 therapy. We find that, while overall high mutational loads associate with improved survival both in responding and non-responding patients, responding tumors are specifically enriched for mutations in the DNA repair gene BRCA2. Innately resistant tumors display a transcriptional signature (referred to as the IPRES or Innate anti-PD-1 Resistance) indicating concurrent upexpression of genes involved in the regulation of mesenchymal transition, cell adhesion, ECM remodeling, angiogenesis and wound-healing. Notably, MAPK-targeted therapy (MAPKi) induces similar signatures in melanoma, suggesting that a non-genomic form of MAPKi resistance mediates cross-resistance to anti-PD-1 therapy. Validation of the IPRES in other independent tumor cohorts defines a transcriptomic subset across distinct types of advanced cancer. These findings suggest that attenuating the biological processes that underlie IPRES may improve anti-PD1 response in melanoma and other cancer types.

Graphical Abstract

An external file that holds a picture, illustration, etc. Object name is nihms-765463-f0001.jpg

INTRODUCTION

PD-1 immune checkpoint blockade therapy induces a high rate of anti-melanoma response and provides unprecedented clinical benefits (Hamid et al., 2013; Topalian et al., 2012). This therapeutic approach has also been shown to be active against a growing list of human malignancies, and clinical testing of combinations of PD-1 (or PD-L1) with other treatment targets has already begun (Sharma and Allison, 2015). However, effective clinical use of anti-PD-1 agents is encumbered by a high rate of innate resistance (60-70%) in advanced metastatic melanoma. The mechanistic basis for the variation in response patterns or long-term clinical benefits (i.e., survival) remains poorly explained.

In melanoma, the extent of pretreatment and especially treatment-induced intra-tumoral T cell infiltration correlates with clinical responses (Tumeh et al., 2014), supporting unleashing of tumor-specific T cells as the primary mechanistic basis of anti-PD-1 therapy. Preliminary retrospective analyses of clinical data hinted at prior failure of MAPK-targeted therapy being a negative factor for subsequent response to immune checkpoint blockade in melanoma (Puzanov et al., 2015; Ramanujam et al., 2015; Simeone et al., 2015). In this context, acquired resistance to MAPK-targeted therapy has been correlated with depletion of intra-tumoral T cells, exhaustion of CD8 T cells and loss of antigen presentation (Hugo et al., 2015).

At the genomic level, the overall mutation load has been correlated with clinical responses to anti-PD-1 therapy and linked to smoking in non-small cell lung cancer or mismatch repair deficiency, in colon cancer. (Le et al., 2015; Rizvi et al., 2015). Whether it was non-small cell lung tumors on anti-PD-1 therapy or melanoma tumors on anti-CTLA-4 therapy (Snyder et al., 2014; Van Allen et al., 2015), the range of somatic mutation and neoepitope loads of the responding pretreatment tumors overlaps significantly with that of the non-responding tumors, despite statistically significant differences in their medians. As such, the lack of high mutational loads did not necessarily preclude clinical responses, and, conversely, the presence of high mutational loads not always correlates with responses. Thus, additional genomic or non-genomic features seem also to contribute to anti-PD-1 response patterns. Here, we sought to assess omic-scale features related to clinical response and survival patterns in order to gain insights into potential strategies for patient stratification and identification of anti-PD-1 combinatorial therapies.

RESULTS AND DISCUSSION

High Mutational Load Does Not Associate with Tumor Response But Correlates with Improved Patient Survival

We analyzed the whole exome sequences (WES) of 38 pretreatment (pembrolizumab, nivolumab) melanoma tumors (responding, n=21; non-responding, n=17; total 34 of 38 pretreatment; 4 of 38 early on-treatment; 14 of 38 patients with prior MAPKi treatment; Table S1A) and patient-matched normal tissues for germline references. Responding pretreatment tumors were derived from patients who went on to have complete or partial responses or stable disease control (with mixed responses excluded) in response to anti-PD-1 therapy. Non-responding tumors were derived from patients who had progressive disease. These response patterns were based on irRECIST (Hoos et al., 2015; Wolchok et al., 2009). We also analyzed the transcriptomes through RNASeq of responding (n=15) and non-responding (n=13) pretreatment tumors (total 27 of 28 pretreatment; 1 of 28 early on-treatment) with available high-quality RNA. WES achieved a median of 140X coverage in both tumor and normal tissues (Table S1B). We detected a median of 489 non-synonymous somatic mutations in the 38 tumors (range 73 to 3,985, which is similar to that in a different set of melanoma tissues (Van Allen et al., 2015)).

We found that responding pretreatment tumors on anti-PD-1 therapy harbored harbors more non-synonymous single nucleotide variants (nsSNVs) compared to the non-responding tumors, albeit the statistical significance cutoff was not met (median nsSNVs responding=495 and non-responding=281, P=0.30, Mann-Whitney; Table S1B). Increased predicted HLA class I and class II neoepitope loads were also detected in the responding pretreatment tumors, although these differences were not statistically significant either (median HLA class I neoepitopes responding=231 and non-responding=156, P=0.41; median HLA class II neoepitopes responding=130 and non-responding=95, P=0.36, Mann-Whitney; Table S1B). Even when we considered only expressed nsSNV and neoepitope loads, the statistical significance of the differences between the responding versus non-responding tumors was not augmented. The comparison of these two groups of tumors was not likely biased by small differences in mean tumor purities or depth of sequencing (Figure S1A and S1B). The numbers of predicted HLA class I and II neoepitopes were strongly correlated with the number of nsSNVs (Figure S1C). We did not identify any recurrent predicted neoepitope or experimentally validated neoantigens (Table S1C). Previous work analyzing melanoma tumors sampled prior to anti-CTLA-4 antibody therapy had associated responses with a tetrapeptide signature (Snyder et al., 2014). However, we did not observed enrichment of this peptide motif in the pretreatment tumors that responded to anti-PD-1 therapy (Figure S1D). Likewise, analysis of an independent cohort of 110 melanoma tumors pre-anti-CTLA-4 therapy also did not yield enrichment of this tetrapeptide motif among responding tumors (Van Allen et al., 2015).

In addition to examining the relationship between non-synonymous somatic mutational loads in pretreatment tumors and anti-tumor responses (and lack thereof) elicited by anti-PD-1 antibodies, we also examined their relative potential influences on clinical benefits of anti-PD-1 immunotherapy as reflected by patient survival. Notably, a mutational load in the top third (compared to the bottom third) was significantly associated with improved survival (Figure 1A). We also observed a trend toward higher mutational load being associated with better survival among melanoma patients not treated with anti-PD-1 antibodies (TCGA, 2015), although this association did not reach statistical significance (Figure S1E), suggesting that the prognostic power of a high mutational load is augmented in the setting of anti-PD-1 therapy. As expected, a positive association between objective tumor responses and survival was highly statistically significant (Figure 1B). However, when we divided each non-responding and responding tumor group into sub-groups with low or high mutational loads (i.e., below or above the median total somatic nsSNVs of each response group) (Figure 1C), patients with responding tumors of low mutation loads significantly outlived patients with non-responding tumors of high mutation loads (Figure 1D). This is despite the fact that mutational loads of these two groups were significantly different, with no overlap across the two distributions (Figure 1C). Hence, factors beyond the mutational load also influence shorter-term tumor response patterns and longer-term patient survival.

An external file that holds a picture, illustration, etc. Object name is nihms-765463-f0002.jpg

Mutational Correlates of Innate Sensitivity to Anti-PD-1 Therapy

(A) Overall survival of anti-PD-1-treated patients whose melanoma tumors harbored high (top third) versus low (bottom third) mutational (somatic nsSNVs) loads. P values, log-rank test.

(B) Overall survival of anti-PD-1-treated melanoma patients whose pretreatment tumors responded (n=20) or did not respond (n=17). P value, log-rank test.

(C) Total number of nsSNVs detected in anti-PD-1 responding and non-responding melanoma tumors harboring high (above the respective group's median) or low (below the group median) mutational loads. P value, log-rank test.

(D) Overall survival of anti-PD-1-treated melanoma patients whose pretreatment tumors responded or did not respond and harboring high (above the group median) or low (below the group median) mutational loads. P value, log-rank test.

(E) Recurrent exomic alterations (nsSNVs and small insertion/deletions or INDELs) in pretreatment tumors of responding versus non-responding patients on anti-PD-1 therapy. Copy number alterations were annotated for the same gene as a reference. Top, mutations of melanoma signature genes. Middle, mutations recurrent in responding versus non-responding tumors (recurrence in 25% in one group and at most one occurrence in the opposite group, Fisher exact test, FDR-corrected P≤0.05 on enrichment against the background mutation frequency). Bottom, the total nsSNV load of each melanoma tumor.

(F) Schematics of impact of non-synonymous missense and nonsense mutations in the BRCA2 protein and its domains.

(G) Total number of nsSNVs detected in melanomas with or without BRCA2 non-synonymous mutations. P value, Mann Whitney test.

See also Table S1 and Figure S1.

Enrichment for BRCA2 Mutations in Anti-PD-1 Responsive Melanoma

We then sought to identify mutations (nsSNVs and small insertion-and-deletions or INDELs; Table S1D) that (i) were recurrently and selectively associated with either responding or non-responding tumors (recurrence ≥ 25% in one group and at most one hit in the other group) and (ii) occurred in genes at rates higher than background rates (Fisher exact test, FDR-corrected p≤0.05) (Figure 1E; Table S1E). The background mutation rate of each gene was calculated from the WES data of 469 melanoma tumors (Hodis et al., 2012; TCGA, 2015). Analysis of copy number variation (CNVs) did not identify any recurrent alterations exclusive to either group. BRCA2 harbored nsSNVs in six of 21 responding tumors (28%) but only one of 17 non-responding tumors (6%) (Figure 1E). With a background mutational rate estimated at 6% (28 of 469 melanoma tumors), BRCA2 was significantly more frequently mutated in the responding tumors than expected (Fisher P=0.002, odds ratio=6.2). The pattern of mutations in disparate BRCA2 protein domains suggested loss-of-function mutations (Figure 1F): one in the N-terminal NPM1-interacting region; one in the POLH-interacting domain; and four in the helical domain critical for FANCD2 interaction. Intriguingly, the somatic mutational load of the tumors with BRCA2 nsSNVs was significantly higher than those with wild type BRCA2 in this cohort of tumors (Figure 1G) as well as two additional cohorts of melanoma tumors (Figure S1F). Thus, BRCA2 LOF mutations, which are expected to produce defects in homologous recombination and double-stranded DNA break repair (Holloman, 2011), may produce specific mutational signatures or unknown effects (e.g., induction of cell death) which contribute to anti-PD-1 responsiveness.

Co-enriched Transcriptomic Signatures in a Major Subset of Anti-PD-1 Resistant Melanoma

We then addressed whether transcriptomic features would differentiate between responding (n=15) versus non-responding (n=13) tumors sampled prior to anti-PD-1 therapy (total 27 of 28 pretreatment tumors and 1 of 28 early on-treatment). We compared the transcriptomes of the two tumor groups using two approaches: (i) analysis of differentially expressed genes (DEGs) (Figure 2A top and Figure 2B) across the two aggregate groups (Table S2A) coupled with GO term enrichment analysis of DEGs (Figure 2C) and (ii) differential signature enrichment based on single-sample Gene Set Variance Analysis or GSVA scores using publicly available (C2 chemical and genetic perturbation C6 oncogenic, and C7 immunologic subsets of the Molecular Signature Database, Broad Institute) and self-curated (see below), perturbation-induced gene signatures (Table S2B; Figure 2D).

An external file that holds a picture, illustration, etc. Object name is nihms-765463-f0003.jpg

Transcriptomic Signatures of Innate Resistance to Anti-PD-1 Therapy

(A) (Top) Heatmap showing differentially expressed genes in the pretreatment tumors derived from patients who responded versus who did not respond to anti-PD-1 treatment (gene expression with inter-quartile range (IQR) ≥ 2; median fold-change (FC) difference ≥ 2; Mann-Whitney P ≤ 0.05). (Middle) mRNA expression levels of genes with hypothetical roles in modulating response patterns to anti-PD-1 therapy. (Bottom) Overall number of nsSNVs, HLA class 1 and 2 neoepitopes (predicted).

(B) mRNA levels of genes (which control tumor cell mesenchymal transition, tumor angiogenesis and macrophage and monocyte chemotaxis) that were differentially expressed between the responding versus non-responding pretreatment tumors. P values, Mann Whitney test.

(C) GO enrichment of genes that were expressed higher in the responding tumors.

(D) Heatmap showing the Gene Set Variance Analysis (GSVA) scores of gene signatures differentially enriched in the responding versus non-responding pre-anti-PD-1 tumors (absolute median GSVA score difference ≥ 10%, FDR-corrected Welch t-test p≤0.25 or nominal Welch t-test p≤0.1). For comparison, enrichment scores of interferon signatures are also displayed.

(E) Overall survival of anti-PD-1-treated melanoma patients with presence (n=10) or absence (n=16) of co-enriched Innate Anti-PD-1 RESistance (IPRES) signatures. P value, log-rank test.

See also Table S2 and Figure S2.

From analysis of DEGs (cutoff, two-fold difference between the absolute medians of normalized expressions in the two groups; nominal Mann-Whitney p≤0.1), we made observations suggesting that mesenchymal and inflammatory tumor phenotypes may be associated with innate anti-PD-1 resistance. First, 693 genes were differentially expressed between the responding versus non-responding pretreatment tumors, and the transcriptomes of non-responding tumors were dominated by relative gene up-expression events compared with the transcriptomes of responding tumors (Table S2A; Figure 2A top, showing only genes whose differential expression met nominal Mann-Whitney p≤0.05). Second, DEGs that were expressed higher in non-responding pretreatment tumors included mesenchymal transition genes (AXL, ROR2, WNT5A, LOXL2, TWIST2, TAGLN, FAP), immunosuppressive genes (IL10, VEGFA, VEGFC), and monocyte and macrophage chemotactic genes (CCL2, CCL7, CCL8 and CCL13) (Figure 2A and 2B). In addition to mesenchymal genes, genes associated with wound healing and angiogenesis, which are considered T cell-suppressive (Motz and Coukos, 2011; Schafer and Werner, 2008; Voron et al., 2014), were expressed higher among non-responding relative to responding pretreatment tumors. Interestingly, a recent study using a mouse melanoma model showed that VEGFA and CCL2 expression was associated with innate anti-PD-1 resistance (Peng et al., 2015). CDH1, which is typically down-expressed by mesenchymal cancer cells, was also down-expressed by non-responding (versus responding) pretreatment tumors. Third, genes with putative roles in modulating immune checkpoint sensitivity were not differentially expressed between responding versus non-responding tumor groups (Figure 2A bottom; Figure S2). GZMA, PRF1 (CD8 T cell cytolytic score), PDCD1LG2 (PD-L2) and CTLA4 were expressed higher in the pretreatment melanoma tumors of patients who derived benefit from CTLA-4 antibodies (Van Allen et al., 2015). However, these genes, along with other T cell-related genes such as CD8A/B, PD-L1, LAG3 (T cell checkpoint genes) and IFNG, did not present higher expression in anti-PD-1-responsive tumors (Figure 2A bottom; Figure S2A). Similarly, we did not observe higher enrichment of multiple interferon signatures in the anti-PD-1-responsive group (Figure 2C bottom). Previously, an interferon gamma signature was found to be differentially up-expressed in the pretreatment tumor biopsies from responding patients when a restricted set of immune genes were analyzed (Ribas et al., 2015). However, the technical approach may not be comparable to our whole tumor transcriptomic approach. We did note that the expression levels of HLA class I genes (HLA-A, -B, -C) trended higher among the responding tumors, although the differences were not statistically significant. Lastly, the complete loss of PTEN was reported to promote resistance to immune checkpoint blockade (Peng et al., 2015), but there was only one case of homozygous PTEN deletion (with nearly undetectable PTEN mRNA expression; Figure S2A) in our cohort (in the non-responsive sub-group), limiting our ability to draw meaningful associations in this dataset. Generally, we did not observe a statistically significant difference in PTEN expression between anti-PD-1 responding versus non-responding tumors. Thus, individual gene-based expression analysis suggested mesenchymal and T cell-suppressive inflammatory or angiogenic tumor phenotypes as being associated with innate anti-PD-1 resistance.

We then queried biological processes represented by DEGs. While gene ontology (GO) enrichment analysis of genes up-expressed among responding tumors produced no significantly enriched terms, genes up-expressed among non-responding tumors were enriched for cell adhesion, ECM organization, wound healing and angiogenesis (FDR-adjusted p-values of GO gene sets shown in Figure 2C). Using independently derived perturbation-based transcriptomic signatures (Molecular Signature Database; Table S2C), we tested for differentially enriched processes in the responding versus non-responding pretreatment tumors (cutoff, 10% difference between the absolute medians of GSVA scores in the two groups; FDR-corrected Welch t-test p≤0.25). Gene sets meeting these standard cutoffs formed the core sets (Figure 2D upper, in bold) from which we compiled additional concurrently enriched (nominal Welch t-test p≤0.1) and functionally related gene sets (Figure 2D upper, Table S2B). We considered these statistically weaker gene set enrichments biologically meaningful given the functional coherence of these gene signatures with the core signatures (Subramanian et al., 2005).

Importantly, a group of 26 transcriptomic signatures were co-enriched en bloc in 9 of 13 non-responding versus 1 of 15 responding pre-anti-PD-1 tumors (see Experimental Procedures). Co-enrichment of these signatures, collectively referred to as the Innate anti-PD-1 Resistance or IPRES signature, again indicated heightened mesenchymal transition, angiogenesis, hypoxia and wound healing. The concurrence of a tumor cell mesenchymal phenotype with an angiogenesis- and wound healing-related inflammatory microenvironment has been documented in the literature (Chen et al., 2015a; Chen et al., 2015b; Mak et al., 2015). Interestingly, this set of 26 IPRES signatures included signatures induced by MAPK inhibitor (MAPKi) treatment of melanoma tumors and cell lines (Table S2C). We have shown recently that MAPKi treatment of melanoma cells induces transcriptome-wide re-programming leading to concurrent phenotype switches (Song et al., 2015). Notably, MAPKi-induced signatures of mesenchymal-invasive transition, angiogenesis, and wound healing signatures were detected in the residual melanoma tumors from patients on MAPKi therapy, suggesting that induction of these signatures may negatively impact responsiveness to combinatorial anti-PD-1/L1 therapy.

IPRES (Innate anti-PD-1 Resistance) Signatures Define a Transcriptomic Subset Across Cancers

The observations that IPRES content signatures were co-enriched in the same tumors (Figure 2D) and that MAPKi induced these signatures concurrently (Table S2C) implied co-regulated tumor phenotypes that together define a transcriptomic subset. To evaluate whether co-enrichment of IPRES content signatures was an exclusive feature of our cohort, we queried three additional cohorts of metastatic melanoma-derived RNASeq (Hugo et al., 2015; TCGA, 2015; Van Allen et al., 2015), including a cohort consisting of only V600_BRAF_ mutant melanomas (cohort 3) (Hugo et al., 2015). We found that IPRES content signatures co-enriched not only in the same tumors but also in about a third of total samples in each of the four independent transcriptomic data sets (cohort 1 from this study, 10 IPRES-enriched tumors of 28 total tumors; cohort 2, 15 of 42; cohort 3, 11 of 32; cohort 4, 90 of 282) (Figure 3A). Considering 126 among 384 total tumors as the background prevalence for co-enrichment of IPRES content signatures in metastatic melanoma, we determined that this IPRES-enriched transcriptomic subset was over-represented among the anti-PD-1 non-responding pretreatment tumors (Fisher P=0.013, odds ratio=4.6) and under-represented among the responding pretreatment tumors (Fisher P=0.04, odds ratio=0.15) within cohort 1. In contrast, co-enrichment of IPRES signatures was neither over- nor -under-represented among the responding or non-responding pre-anti-CTLA-4 melanoma tumors in cohort 2 (Figure S2B) (Van Allen et al., 2015), which suggests that mechanisms of innate resistance to anti-PD-1 and anti-CTLA-4 are not necessarily similar.

An external file that holds a picture, illustration, etc. Object name is nihms-765463-f0004.jpg

Co-enrichment of Innate Anti-PD-1 Resistance-associated Signatures Defines a Transcriptomic Subset in Melanoma and Multiple Cancers

(A) Heatmap showing GSVA scores of IPRES signatures across four independent RNASeq data sets derived from metastatic melanoma. Cohort 1, pretreatment (anti-PD-1) tumors; cohort 2, pretreatment (anti-CTLA-4) tumors; cohort 3, pretreatment (MAPKi) tumors; cohort 4, TCGA cutaneous melanoma (metastatic only).

(B) Heatmap showing GSVA scores of IPRIM signatures across TCGA RNASeq data sets (metastatic melanoma or SKCM, lung adenocarcinoma or LUAD, colon adenocarcinoma or COAD, kidney clear cell carcinoma or KIRC, and pancreatic adenocarcinoma or PAAD).

See also Figure S3.

Furthermore, co-enrichment of the IPRES signatures defined a transcriptomic subset within not only melanoma but also all major common human malignancies analyzed (Figure 3B). The IPRES-enriched transcriptomic subset of certain cancers such as pancreatic adenocarcinoma made up the majority of tumors. Within a side-by-side comparison, only six of 69 primary cutaneous melanomas showed co-enrichment of IPRES signatures, in contrast to 90 of 282 metastatic (TCGA) melanomas (P=3.9e-5, odds ratio=0.2) (Figure S3), consistent with mesenchymal transition and metastasis gene sets among IPRES signatures. Thus, co-enrichment of IPRES signatures defines a distinct transcriptomic program that exists across cancers of distinct histology.

This study highlights the utility of both exome and transcriptome sequencing data generated from pretreatment tumor samples for the identification of potential determinants of response to anti-PD-1. Although the overall somatic mutational loads of anti-PD-1-responsive melanoma tumors were not significantly higher than those of non-responsive tumors, higher mutational loads associated significantly with better survival after anti-PD-1 therapy. This finding is still consistent with the notion that neoepitopes derived from somatic non-synonymous mutations are critical for deriving clinical benefits from anti-PD-1 therapy in melanoma. However, objective tumor responses, although strongly associated with survival benefits, did not appear to be driven overwhelmingly by the overall somatic mutational loads. That is to say, a relatively low mutational load did not preclude a tumor response. This is consistent with findings from gastrointestinal cancers where low mutational loads did not preclude tumor infiltration by mutation-reactive, class I and II-restricted T cells (Tran et al., 2015). Thus, overall somatic or predicted neoepitope loads of pretreatment melanoma tumors are not enough to predict response patterns to anti-PD-1 therapy.

In our cohort, responsive tumors were significantly enriched for (likely) loss-of-function mutations in BRCA2. As one would predict from the known function of BRCA2 in DNA repair, _BRCA2_-mutated melanomas harbored higher mutational loads than _BRCA2_-wildtype melanomas. Although it is conceivable that defective BRCA2-DNA repair results in specific mutational motifs (as opposed to the general increase in mutational load) that enhance responsiveness, it is also possible that cellular stress resulting from defective DNA repair could lead to increased cell death and anti-tumor immunity. Moreover, these data support the notion that tumor cell phenotypic plasticity (i.e., mesenchymal transition) and the resultant impacts on the microenvironment (e.g., ECM remodeling, cell adhesion, angiogenesis- features of immune suppressive wound healing) are critical barriers to anti-PD-1 responses. The limited number of patients in our melanoma cohort posed certain challenges to our analysis. For example, we relaxed the statistical stringency in single gene-based differential expression analysis (bypassing multiple hypothesis correction) to derive enough genes for gene ontology enrichment analysis. However, converging findings from alternative analysis (i.e., GSVA) of the transcriptome data helped to mitigate potential caveats. Finally, in separate work, we found that mutation-targeted therapy (i.e., MAPKi) induces tumor cell-autonomous changes (e.g., mesenchymal transition) (Song et al., 2015) and upregulates anti-PD-1 resistance-associated processes in residual tumors that have regressed in response to MAPKi treatment. Thus, while our findings in this study necessitate confirmation in independent tissue cohorts, the identification of transcriptomic features associated with anti-PD-1 resistance suggests that mitigation of IPRES-related biological processes may enhance response rates to anti-PD-1 (and anti-PD-1 plus MAPKi) therapy.

EXPERIMENTAL PROCEDURES

Tumor Specimens and Profiling

All tissues in this study were obtained with the approval of Institutional Review Boards and patients’ consents. All patients received either pembrolizumab or nivolumab as the anti-PD-1 therapy for their metastatic melanoma. Thirty-eight melanoma specimens (thirty-two pre-treatment tumors, two pretreatment tumor-derived cultures, three early on-treatment tumors without response, and one early on-treatment tumor with response) and their patient-matched normal tissues were analyzed by whole exome sequencing (WES). Among these thirty-eight samples with WES data, twenty-eight with sufficient RNA quality were also analyzed by RNASeq. This set include another RNASeq dataset derived from a second-site pre-treatment tumor biopsy from patient #27. However, this second-site, pre-treatment tumor-derived WES dataset was excluded in our aggregate mutation analysis to avoid double-counting two tumor exomes from the same patient (Table S1A).

Thirty eight tumor specimens and their respective normal tissues were subjected to whole exome sequencing (WES) (Table S1A). WES was performed using pair-end sequencing with read length of 2×100 bps based on the Illumina HiSeq2000 platform. RNA from a subset of twenty eight tumors were pair-end sequenced with read length of 2 × 100 bps (Illumina HiSeq2000). We included two tumors from Pt27 for transcriptomic analyses but not for mutation and neoepitope analyses since the tumors may not share the same transcriptomic profile but they essentially contain the same set of non-synonymous somatic mutations.

Whole Exome Sequencing

We called single nucleotide variant (SNV) and small insertion-deletion (INDEL) as reported (Shi et al., 2014) using a stand-alone version of Oncotator (Ramos et al., 2015). Copy numbers were called using the intersection of the copy number calls derived from Sequenza (Favero et al., 2015) and VarScan2 (Koboldt et al., 2012). Tumor purities and ploidies (Table S1B) were calculated based on the calls of Sequenza using WES data with default parameters. The impact of BRCA2 nsSNVs was visualized using the domain information in the INTERPRO protein domain database (Mitchell et al., 2015).

HLA Types and Neoepitopes

The 4-digit HLA class 1 and 2 types of each patient were called using ATHLATES (Liu et al., 2013) using the WES sequencing reads from the normal tissue. To ensure concordance, we manually compared ATHLATES’ calls of the normal versus tumor samples and ascertained there was at least no two-digit HLA typing discrepancy between any normal-tumor pair. For each non-synonymous coding mutation from a tumor, we predicted its impact on the patient's HLA class I and II binding using the stand-alone version of the programs NetMHCpan v2.8 (Hoof et al., 2009; Nielsen et al., 2007) and NetMHCIIpan v3.0 (Karosiene et al., 2013), respectively. Specifically, for HLA class I binding prediction using netMHCpan v2.8, we tested all 9-11-mer peptides containing the mutated amino acids for binding to the patient's HLA-A, -B and -C. A peptide was defined as a neoepitope based on two criteria : i) predicted binding affinity ≤ 500nM, and ii) rank percentage ≤ 2% (default cutoff). For HLA class II binding prediction using netMHCIIpan v3.0, we tested the 9-19-mers containing the mutated amino acids for binding to the patient-specific, ATHLATES-predicted DPA-DPB, DQA-DQB and DRA-DRB allele pairs. We also applied the same predicted binding affinity and rank percentage cutoff as we did for HLA class I to nominate the HLA class II-binding neoepitopes. Expressed non-synonymous mutations and neoepitopes were defined based on corresponding genes with normalized expression levels ≥ 1 (in FPKM). Statistical differences of nsSNV, HLA class I and II neoepitopes, WES coverages and tumor purities were computed using two-sided Mann-Whitney test.

Mutation Recurrence

To estimate the statistical significance of the recurrence of gene mutations in the responding or non-responding tumors, we used an independent batch of 469 melanomas’ whole exome sequence datasets (Hodis et al., 2012; TCGA, 2015) to estimate each gene's background mutation frequency. Significance was computed by Fisher exact test followed by FDR adjustment for multiple hypothesis testing. We listed genes that fulfilled the criteria: i) recurrence in at least 25% of the responder/non-responder, ii) occurrence of at most once in the opposite group and iii) Fisher exact test FDR adjusted p-value ≤ 0.05. These genes were illustrated in Figure 1A and all genes that fulfilled i) and ii) and tested for multiple hypotheses were listed in Table S1E. The association between BRCA2 nsSNVs and overall nsSNV counts were tested using two-sided Mann-Whitney test and validated in independent WES datasets (Hodis et al., 2012; TCGA, 2015).

RNASeq and Gene Set Enrichment

Paired-end transcriptome reads were mapped to the UCSC hg19 reference genome using Tophat2 (Kim et al., 2013). Normalized expression levels of genes were expressed in FPKM values as generated by cuffquant and cuffnorm. The program were run with the option “--frag-bias-correct” and “--multi-read-correct” to improve sensitivity (Roberts et al., 2011). A gene was defined as differentially expressed between the responding and non-responding tumor groups when its median expression differed by at least two-fold between the groups with a nominal two-sided Mann-Whitney p-value ≤ 0.1 (Table S2A). Applying multiple hypothesis correction of FDR p ≤ 0.25 only yielded 3 differentially expressed genes; ALDH1L2 and MFAP2 in the non-responding and CDH1 (E-cadherin) in the responding group. As such, the genes meeting the uncorrected, nominal Mann-Whitney p-value ≤ 0.1 that were expressed higher either in the responding or non-responding group were separately analyzed for GO term enrichments using the online functional annotation tools DAVID (Huang et al., 2008). Enriched GO terms were selected from the GO biological process terms in DAVID's Fat database (Huang et al., 2009). GO terms which were highly overlapping, as defined by functional clustering in DAVID's website, were represented by the terms with the best FDR-adjusted p-values.

To calculate single-sample gene set enrichment, we used the GSVA program (Hanzelmann et al., 2013) to derive the absolute enrichment scores of previously experimentally validated gene signatures as follow: i) the C2 CGP (chemical and genetic perturbation sets), ii) the C6 and C7 subset of the Molecular Signature Database (Subramanian et al., 2005), iii) self-curated MAPK inhibitor-induced gene signatures using cell lines and patient-derived tumors (Song et al., 2015), iv) post-operation wound signature (Inkeles et al., 2015), and v) melanoma invasive/proliferative signatures (Hoek et al., 2008). To derive the GSVA score of each signature in each tumor sample, we computed from raw RNASeq read counts by HTSEQ COUNT program and then normalized them to log2 CPM values using EdgeR (McCarthy et al., 2012). We removed batch effects using the edgeR function RemoveBatchEffect when we combined RNAseq data from multiple experiments (Figure 3A). The normalized log2 CPM values were then passed on as input for GSVA in the RNASeq mode. Differentially enriched core gene sets between the responding and non-responding tumor groups were defined by GSVA score differences of ≥ 10% and FDR-corrected, two-sided Welch T-test p-value ≤ 0.25 (we used T-test because the GSVA scores were normally distributed around 0). Two gene sets, INGRAM_SHH_TARGETS_DN and WONG_ENDMETRIUM_CANCER_DN, were not included in the core set because they did not specifically point to a cellular process and/or relate to the other six gene sets in the core set (Table S2B, top 8). We also collected gene sets that met the GSVA score differences of ≥ 10% and nominal Welch T-test p-value ≤ 0.1 (Table S2B) and included those which were concordantly enriched and functionally related to the core gene sets to make up the full list of IPRES signatures (Figure 2D).

To compare co-enrichment of IPRES signatures across multiple melanoma cohorts, we combined and batch-corrected the log2 CPM values of four melanoma transcriptome cohorts: i) our current pre-anti-PD-1 tumors (n=28), ii) pre-anti-CTLA-4 tumors (n=42), iii) pre-MAPKi tumors (n=32) and iv) the metastatic subset of TCGA melanoma (n=282). We row-normalized the GSVA scores of each gene set in the IPRES signature across the samples from the four cohorts. For this comparative study, we excluded the gene sets “JAEGER_METASTASIS_UP,” “YE_METASTATIC_LIVER_CANCER,” “KARAKAS_ TGFB1_SIGNALING,” and “JEON_SMAD6_ TARGETS_ DN” from the IPRES set because they showed weaker co-enrichment with rest of the gene sets (see Figure 2D upper panel). The IPRES (enrichment) score was defined as the average Z-score across all gene sets in the IPRES signature, and we applied an average Z-score of 0.35 as the cutoff for IPRES signature enrichment in a tumor sample. This resulted in IPRES co-enrichment in 9 non-responding tumors and 1 responding tumor in our anti-PD-1 cohort (this cutoff was chosen because it provided the largest average Z-score separation between the samples with and without IPRES co-enrichment). Since the average Z-score was not comparable between different cohorts, we used the 90th highest IPRES score in the TCGA metastatic melanoma cohort as the IPRES score cutoff (since there were 90 of 282 tumors showing IPRES co-enrichment in this TCGA metastatic cohort; Figure 3A) for analyses performed to yield Figures 3B and S3. This allowed for a non-parametric comparison across multiple TCGA datasets at the IPRES co-enrichment level established in our anti-PD-1 cohort.

Source Data

Analysis of differential non-synonymous mutational hits in responders versus non-responders to ipilimumab was based on the mutation calls as reported (Van Allen et al., 2015). We curated published CD8 T cell exhaustion genes (Wherry, 2011) to minimize those likely to be expressed by melanoma cells by excluding genes whose maximum log2 FPKM was 1 in an in-house melanoma cell line-derived RNASeq database (n = 26 cell lines). This resulted in the inclusion of genes for surface receptors PDCD1 (PD-1), LAG3, HAVCR2 (Tim-3), CD160, and CD244 as well as transcription factors EOMES, PRDM1 (Blimp-1) and TBX21 (T-bet). We assessed co-enrichment of IPRES content signatures in the i) anti-CTLA-4 pretreatment cohort (Van Allen et al., 2015), ii) MAPKi pretreatment cohort (Hugo et al., 2015; Song et al., 2015), iii) TCGA melanoma (metastatic and primary subsets separately analyzed) (TCGA, 2015), iv) TCGA pancreatic ductal adenocarcinoma (TCGA, provisional, see https://tcga-data.nci.nih.gov/tcga/tcgaHome2.jsp), v) TCGA lung adenocarcinoma (TCGA, 2014), and vi) TCGA colorectal adenocarcinoma (TCGA, 2012), and vii) TCGA kidney clear cell carcinoma (TCGA, 2013).

​

Article Highlights

High mutational loads may predict better survival, but not response to anti-PD-1

BRCA2 mutations are enriched melanomas responsive to anti-PD-1

A transcriptional signature (IPRES) is related to innate anti-PD-1 resistance

IPRES defines a transcriptomic subset across distinct types of advanced cancer

Supplementary Material

1

2

3

4

5

6

ACKNOWLEDGEMENTS

This work has been funded by the National Institutes of Health (R35 CA197633 to A.R., P01 CA168585 to A.R. and R.S.L., 1R01CA176111 to R.S.L; K12CA0906525 to D.B.J.), the Ressler Family Foundation (to R.S.L. and A.R.), Grimaldi Family Fund (to A.R. and R.S.L.), Burroughs Wellcome Fund (to R.S.L), the Ian Copeland Melanoma Fund (to R.S.L.), Wade F.B. Thompson/Cancer Research Institute CLIP Grant (to R.S.L.), Steven C. Gordon Family Foundation (to R.S.L.), American Skin Association (to H.W.), Dr. Robert Vigen Memorial Fund the Garcia-Corsini Family Fund (to A.R.), the ASCO Conquer Cancer Career Development Award (to K.B.D), and American Cancer Society Research Professorship (to J.A.S.). Patient-informed consents were obtained for the research performed in this study.

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

AUTHOR CONTRIBUTION

W.H., J.M.Z., L.S., C.S., A.R. and R.S.L. generated, analyzed and interpreted data. S. H-L., B.H.M., B.C., G.C., E.S., M.C.K., J.A.S., D.B.J., A.R. and R.S.L. evaluated patients and provided tissue reagents. B. H-M., J.P., S.L. and X.K. processed tissues for analysis. All authors contributed to manuscript preparation. W.H. and R.S.L developed the concepts, supervised the project, and wrote the paper.

ACCESSION NUMBER

The GEO accession number for the transcriptome data is GSE78220.

REFERENCES