Highlights

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

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

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

  • IPRES defines a transcriptomic subset across distinct types of advanced cancer

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 overall high mutational loads associate with improved survival, and tumors from responding patients are 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 up-expression of genes involved in the regulation of mesenchymal transition, cell adhesion, extracellular matrix remodeling, angiogenesis, and wound healing. Notably, mitogen-activated protein kinase (MAPK)-targeted therapy (MAPK inhibitor) induces similar signatures in melanoma, suggesting that a non-genomic form of MAPK inhibitor 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-PD-1 response in melanoma and other cancer types.

Graphical Abstract

Figure thumbnail fx1

Graphical Abstract

Introduction

PD-1 immune checkpoint blockade therapy induces a high rate of anti-melanoma response and provides unprecedented clinical benefits (
Hamid et al., 2013
  • Hamid O.
  • Robert C.
  • Daud A.
  • Hodi F.S.
  • Hwu W.J.
  • Kefford R.
  • Wolchok J.D.
  • Hersey P.
  • Joseph R.W.
  • Weber J.S.
  • et al.

Safety and tumor responses with lambrolizumab (anti-PD-1) in melanoma.N. Engl. J. Med. 2013; 369: 134-144

,
Topalian et al., 2012
  • Topalian S.L.
  • Hodi F.S.
  • Brahmer J.R.
  • Gettinger S.N.
  • Smith D.C.
  • McDermott D.F.
  • Powderly J.D.
  • Carvajal R.D.
  • Sosman J.A.
  • Atkins M.B.
  • et al.

Safety, activity, and immune correlates of anti-PD-1 antibody in cancer.N. Engl. J. Med. 2012; 366: 2443-2454

). 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
  • Sharma P.
  • Allison J.P.

Immune checkpoint targeting in cancer therapy: toward combination strategies with curative potential.Cell. 2015; 161: 205-214

). 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 in 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
  • Tumeh P.C.
  • Harview C.L.
  • Yearley J.H.
  • Shintaku I.P.
  • Taylor E.J.
  • Robert L.
  • Chmielowski B.
  • Spasic M.
  • Henry G.
  • Ciobanu V.
  • et al.

PD-1 blockade induces responses by inhibiting adaptive immune resistance.Nature. 2014; 515: 568-571

), 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 mitogen-activated protein kinase (MAPK)-targeted therapy being a negative factor for subsequent response to immune checkpoint blockade in melanoma (Puzanov et al., 2015, Pigment Cell Melanoma Res., abstract; Ramanujam et al., 2015, Pigment Cell Melanoma Res., abstract; Simeone et al., 2015, Pigment Cell Melanoma Res., abstract). 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
  • Hugo W.
  • Shi H.
  • Sun L.
  • Piva M.
  • Song C.
  • Kong X.
  • Moriceau G.
  • Hong A.
  • Dahlman K.B.
  • Johnson D.B.
  • et al.

Non-genomic and Immune Evolution of Melanoma Acquiring MAPKi Resistance.Cell. 2015; 162: 1271-1285

).
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
  • Le D.T.
  • Uram J.N.
  • Wang H.
  • Bartlett B.R.
  • Kemberling H.
  • Eyring A.D.
  • Skora A.D.
  • Luber B.S.
  • Azad N.S.
  • Laheru D.
  • et al.

PD-1 Blockade in Tumors with Mismatch-Repair Deficiency.N. Engl. J. Med. 2015; 372: 2509-2520

,
Rizvi et al., 2015
  • Rizvi N.A.
  • Hellmann M.D.
  • Snyder A.
  • Kvistborg P.
  • Makarov V.
  • Havel J.J.
  • Lee W.
  • Yuan J.
  • Wong P.
  • Ho T.S.
  • et al.

Cancer immunology. Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer.Science. 2015; 348: 124-128

). 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
  • Snyder A.
  • Makarov V.
  • Merghoub T.
  • Yuan J.
  • Zaretsky J.M.
  • Desrichard A.
  • Walsh L.A.
  • Postow M.A.
  • Wong P.
  • Ho T.S.
  • et al.

Genetic basis for clinical response to CTLA-4 blockade in melanoma.N. Engl. J. Med. 2014; 371: 2189-2199

,
Van Allen et al., 2015
  • Van Allen E.M.
  • Miao D.
  • Schilling B.
  • Shukla S.A.
  • Blank C.
  • Zimmer L.
  • Sucker A.
  • Hillen U.
  • Foppen M.H.
  • Goldinger S.M.
  • et al.

Genomic correlates of response to CTLA-4 blockade in metastatic melanoma.Science. 2015; 350: 207-211

), the range of somatic mutation and neoepitope loads of the responding pretreatment tumors overlapped 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 does not always correlate with responses. Thus, additional genomic or non-genomic features seem likely 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 and 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 MAPK inhibitor (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
  • Hoos A.
  • Wolchok J.D.
  • Humphrey R.W.
  • Hodi F.S.

CCR 20th Anniversary Commentary: Immune-Related Response Criteria--Capturing Clinical Activity in Immuno-Oncology.Clin. Cancer Res. 2015; 21: 4989-4991

,
Wolchok et al., 2009
  • Wolchok J.D.
  • Hoos A.
  • O’Day S.
  • Weber J.S.
  • Hamid O.
  • Lebbé C.
  • Maio M.
  • Binder M.
  • Bohnsack O.
  • Nichol G.
  • et al.

Guidelines for the evaluation of immune therapy activity in solid tumors: immune-related response criteria.Clin. Cancer Res. 2009; 15: 7412-7420

). We also analyzed the transcriptomes through RNA-seq 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 140× 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
  • Van Allen E.M.
  • Miao D.
  • Schilling B.
  • Shukla S.A.
  • Blank C.
  • Zimmer L.
  • Sucker A.
  • Hillen U.
  • Foppen M.H.
  • Goldinger S.M.
  • et al.

Genomic correlates of response to CTLA-4 blockade in metastatic melanoma.Science. 2015; 350: 207-211

).
We found that responding pretreatment tumors on anti-PD-1 therapy harbored 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 (human leukocyte antigen) 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 (Figures 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
  • Snyder A.
  • Makarov V.
  • Merghoub T.
  • Yuan J.
  • Zaretsky J.M.
  • Desrichard A.
  • Walsh L.A.
  • Postow M.A.
  • Wong P.
  • Ho T.S.
  • et al.

Genetic basis for clinical response to CTLA-4 blockade in melanoma.N. Engl. J. Med. 2014; 371: 2189-2199

). However, we did not observe 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
  • Van Allen E.M.
  • Miao D.
  • Schilling B.
  • Shukla S.A.
  • Blank C.
  • Zimmer L.
  • Sucker A.
  • Hillen U.
  • Foppen M.H.
  • Goldinger S.M.
  • et al.

Genomic correlates of response to CTLA-4 blockade in metastatic melanoma.Science. 2015; 350: 207-211

).
Figure thumbnail figs1

Figure S1Genomic Features of Melanoma Tumors from Patients Treated with Anti-PD-1 Therapy, Related to Figure 1

Show full caption

(A) The difference of tumor purities between the responding versus non-responding tumors with WES; p value, Mann Whitney test.

(B) The difference of WES coverages between the responding versus non-responding tumors; p value, Mann Whitney test.

(C) Correlations between the number of somatic nsSNVs and the number of predicted HLA class 1 (left) or class 2 (right) neoepitopes. Correlation, Pearson, p values, Student’s t test.

(D) Recurrence of tetrapeptides (previously reported as enriched in responding pre-anti-CTLA-4 tumors) in non-responding and responding pre-anti-PD-1 tumors.

(E) Overall survival of TCGA melanoma patients whose tumors harbored high (top third) versus low (bottom third) mutational (somatic nsSNVs) loads; p value, log-rank test.

(F) Mutational loads (somatic nsSNVs) detected in melanoma with or without BRCA2 somatic nsSNVs in two datasets; p values, Mann Whitney test.

In addition to examining the relationship between non-synonymous somatic mutational loads in pretreatment tumors and anti-tumor responses elicited by anti-PD-1 antibodies, we also examined the potential influences of somatic mutational loads on clinical benefits derived from 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 (
Cancer Genome Atlas Network, 2015

Cancer Genome Atlas Network. Genomic Classification of Cutaneous Melanoma.Cell. 2015; 161: 1681-1696

), 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). Moreover, 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.
Figure thumbnail gr1

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

Show full caption

(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’s median) mutational loads; p value, Mann Whitney test.

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

(E) Recurrent exomic alterations (nsSNVs and small INDELs) in pretreatment tumors of responding versus non-responding patients on anti-PD-1 therapy. Copy-number alterations were annotated for the same genes 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’s 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 missense 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 [indels]; Table S1D) that (1) 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 (2) occurred in genes at rates higher than background rates (Fisher’s exact test, false discover rate [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
  • Hodis E.
  • Watson I.R.
  • Kryukov G.V.
  • Arold S.T.
  • Imielinski M.
  • Theurillat J.P.
  • Nickerson E.
  • Auclair D.
  • Li L.
  • Place C.
  • et al.

A landscape of driver mutations in melanoma.Cell. 2012; 150: 251-263

,
Cancer Genome Atlas Network, 2015

Cancer Genome Atlas Network. Genomic Classification of Cutaneous Melanoma.Cell. 2015; 161: 1681-1696

). Analysis of copy-number variations (CNVs) did not identify any recurrent alteration exclusive to either group. BRCA2 harbored nsSNVs in 6 of 21 responding tumors (28%), but only 1 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 that with wild-type BRCA2 in this cohort of tumors (Figure 1G), as well as two additional cohorts of melanoma tumors (Figure S1F). Thus, BRCA2 loss-of-function mutations, which are expected to produce defects in homologous recombination and double-stranded DNA break repair (
Holloman, 2011
  • Holloman W.K.

Unraveling the mechanism of BRCA2 in homologous recombination.Nat. Struct. Mol. Biol. 2011; 18: 748-754

), may produce specific mutational signatures or unknown effects (e.g., induction of cell death) that 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: (1) analysis of differentially expressed genes (DEGs; Figures 2A, top, and 2B) across the two aggregate groups (Table S2A) coupled with GO term enrichment analysis of DEGs (Figure 2C) and (2) 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 (Figure 2D; Table S2B).

Figure thumbnail gr2

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

Show full caption

(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 and 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 non-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, 2-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 (Figure 2A, top, showing only genes whose differential expression met nominal Mann-Whitney p ≤ 0.05; Table S2A). 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; Figures 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
  • Motz G.T.
  • Coukos G.

The parallel lives of angiogenesis and immunosuppression: cancer and other tales.Nat. Rev. Immunol. 2011; 11: 702-711

,
Schäfer and Werner, 2008
  • Schäfer M.
  • Werner S.

Cancer as an overhealing wound: an old hypothesis revisited.Nat. Rev. Mol. Cell Biol. 2008; 9: 628-638

,
Voron et al., 2014
  • Voron T.
  • Marcheteau E.
  • Pernot S.
  • Colussi O.
  • Tartour E.
  • Taieb J.
  • Terme M.

Control of the immune response by pro-angiogenic factors.Front. Oncol. 2014; 4: 70

), 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., 2016
  • Peng W.
  • Chen J.Q.
  • Liu C.
  • Malu S.
  • Creasy C.
  • Tetzlaff M.T.
  • Xu C.
  • McKenzie J.A.
  • Zhang C.
  • Liang X.
  • et al.

Loss of PTEN promotes resistance to T cell-mediated immunotherapy.Cancer Discov. 2016; 6: 202-216

). 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 (Figures 2A, bottom, and S2A). 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
  • Van Allen E.M.
  • Miao D.
  • Schilling B.
  • Shukla S.A.
  • Blank C.
  • Zimmer L.
  • Sucker A.
  • Hillen U.
  • Foppen M.H.
  • Goldinger S.M.
  • et al.

Genomic correlates of response to CTLA-4 blockade in metastatic melanoma.Science. 2015; 350: 207-211

). 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 (Figures 2A, bottom, and S2A). Similarly, we did not observe higher enrichment of multiple interferon signatures in the anti-PD-1-responsive group (Figure 2D, 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, J. Clin. Oncol., abstract). 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, or -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., 2016
  • Peng W.
  • Chen J.Q.
  • Liu C.
  • Malu S.
  • Creasy C.
  • Tetzlaff M.T.
  • Xu C.
  • McKenzie J.A.
  • Zhang C.
  • Liang X.
  • et al.

Loss of PTEN promotes resistance to T cell-mediated immunotherapy.Cancer Discov. 2016; 6: 202-216

), 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.
Figure thumbnail figs2

Figure S2Gene or Signature Expression Patterns in Pretreatment Melanoma Tumors on Anti-PD-1 or Anti-CTLA-4 Therapies, Related to Figure 2

Show full caption

(A) mRNA levels of genes (CD8 T cell markers, effectors, cytolytic scores; immune checkpoints, MHC class 1, and PTEN) between the responding versus non-responding pretreatment tumors; p values, Mann Whitney test.

(B) Heatmap showing GSVA scores of IPRES signatures across responding (n = 14) versus non-responding (n = 27) pre-anti-CTLA-4 tumors.

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, extracellular matrix (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 sought after 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 set (Figure 2D, top, in bold) from which we compiled additional concurrently enriched (nominal Welch t test p ≤ 0.1) and functionally related gene sets (Figure 2D, top; 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
  • Subramanian A.
  • Tamayo P.
  • Mootha V.K.
  • Mukherjee S.
  • Ebert B.L.
  • Gillette M.A.
  • Paulovich A.
  • Pomeroy S.L.
  • Golub T.R.
  • Lander E.S.
  • Mesirov J.P.

Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.Proc. Natl. Acad. Sci. USA. 2005; 102: 15545-15550

).
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 (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., 2014
  • Chen L.
  • Gibbons D.L.
  • Goswami S.
  • Cortez M.A.
  • Ahn Y.H.
  • Byers L.A.
  • Zhang X.
  • Yi X.
  • Dwyer D.
  • Lin W.
  • et al.

Metastasis is regulated via microRNA-200/ZEB1 axis control of tumour cell PD-L1 expression and intratumoral immunosuppression.Nat. Commun. 2014; 5: 5241

,
Chen et al., 2015
  • Chen L.
  • Heymach J.V.
  • Qin F.X.
  • Gibbons D.L.

The mutually regulatory loop of epithelial-mesenchymal transition and immunosuppression in cancer progression.OncoImmunology. 2015; 4: e1002731

,
Mak et al., 2016
  • Mak M.P.
  • Tong P.
  • Diao L.
  • Cardnell R.J.
  • Gibbons D.L.
  • William W.N.
  • Skoulidis F.
  • Parra E.R.
  • Rodriguez-Canales J.
  • Wistuba I.I.
  • et al.

A Patient-Derived, Pan-Cancer EMT Signature Identifies Global Molecular Alterations and Immune Target Enrichment Following Epithelial-to-Mesenchymal Transition.Clin. Cancer Res. 2016; 22: 609-620

). Interestingly, this set of 26 IPRES signatures included signatures induced by 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 (C.S., L.S., W.H., and R.S.L., unpublished data). 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 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 RNA-seq (
Hugo et al., 2015
  • Hugo W.
  • Shi H.
  • Sun L.
  • Piva M.
  • Song C.
  • Kong X.
  • Moriceau G.
  • Hong A.
  • Dahlman K.B.
  • Johnson D.B.
  • et al.

Non-genomic and Immune Evolution of Melanoma Acquiring MAPKi Resistance.Cell. 2015; 162: 1271-1285

,
Cancer Genome Atlas Network, 2015

Cancer Genome Atlas Network. Genomic Classification of Cutaneous Melanoma.Cell. 2015; 161: 1681-1696

,
Van Allen et al., 2015
  • Van Allen E.M.
  • Miao D.
  • Schilling B.
  • Shukla S.A.
  • Blank C.
  • Zimmer L.
  • Sucker A.
  • Hillen U.
  • Foppen M.H.
  • Goldinger S.M.
  • et al.

Genomic correlates of response to CTLA-4 blockade in metastatic melanoma.Science. 2015; 350: 207-211

), including a cohort consisting of only V600BRAF mutant melanomas (cohort 3;
Hugo et al., 2015
  • Hugo W.
  • Shi H.
  • Sun L.
  • Piva M.
  • Song C.
  • Kong X.
  • Moriceau G.
  • Hong A.
  • Dahlman K.B.
  • Johnson D.B.
  • et al.

Non-genomic and Immune Evolution of Melanoma Acquiring MAPKi Resistance.Cell. 2015; 162: 1271-1285

). 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 datasets (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
  • Van Allen E.M.
  • Miao D.
  • Schilling B.
  • Shukla S.A.
  • Blank C.
  • Zimmer L.
  • Sucker A.
  • Hillen U.
  • Foppen M.H.
  • Goldinger S.M.
  • et al.

Genomic correlates of response to CTLA-4 blockade in metastatic melanoma.Science. 2015; 350: 207-211

), which suggests that mechanisms of innate resistance to anti-PD-1 and anti-CTLA-4 are not necessarily similar.
Figure thumbnail gr3

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

Show full caption

(A) Heatmap showing GSVA scores of IPRES signatures across four independent RNA-seq datasets 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 RNA-seq datasets (metastatic melanoma, SKCM; lung adenocarcinoma, LUAD; colon adenocarcinoma, COAD; kidney clear cell carcinoma, KIRC; and pancreatic adenocarcinoma, 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 6 of 69 primary cutaneous melanomas showed co-enrichment of IPRES signatures, in contrast to 90 of 282 metastatic (The Cancer Genome Atlas [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 distinct cancer histologies.

Figure thumbnail figs3

Figure S3Co-enrichment of IPRES Signatures in Metastatic versus Primary Cutaneous Melanoma, Related to Figure 3

Show full caption

Heatmap showing GSVA scores of IPRES signatures across TCGA primary and metastatic melanoma tumors.

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
  • Tran E.
  • Ahmadzadeh M.
  • Lu Y.C.
  • Gros A.
  • Turcotte S.
  • Robbins P.F.
  • Gartner J.J.
  • Zheng Z.
  • Li Y.F.
  • Ray S.
  • et al.

Immunogenicity of somatic mutations in human gastrointestinal cancers.Science. 2015; 350: 1387-1390

). 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-wild-type melanomas. Although it is conceivable that defective BRCA2-DNA repair results in specific mutational motifs (rather than or in addition 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, and 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 GO 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; C.S., L.S., W.H., and R.S.L., unpublished data) 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. 38 melanoma specimens (32 pre-treatment tumors, 2 pretreatment tumor-derived cultures, 3 early on-treatment tumors without response, and 1 early on-treatment tumor with response) and their patient-matched normal tissues were analyzed by WES. Among these 38 samples with WES data, 28 with sufficient RNA quality were also analyzed by RNA sequencing (RNA-seq). This set included another RNA-seq 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).

38 tumor specimens and their respective normal tissues were subjected to 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 28 tumors were pair-end sequenced with read length of 2 × 100 bps (Illumina HiSeq2000). We included two tumors from patient #27 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 INDEL as reported (
Shi et al., 2014
  • Shi H.
  • Hugo W.
  • Kong X.
  • Hong A.
  • Koya R.C.
  • Moriceau G.
  • Chodon T.
  • Guo R.
  • Johnson D.B.
  • Dahlman K.B.
  • et al.

Acquired resistance and clonal evolution in melanoma during BRAF inhibitor therapy.Cancer Discov. 2014; 4: 80-93

) using a stand-alone version of Oncotator (
Ramos et al., 2015
  • Ramos A.H.
  • Lichtenstein L.
  • Gupta M.
  • Lawrence M.S.
  • Pugh T.J.
  • Saksena G.
  • Meyerson M.
  • Getz G.

Oncotator: cancer variant annotation tool.Hum. Mutat. 2015; 36: E2423-E2429

). Copy numbers were called using the intersection of the copy-number calls derived from Sequenza (
Favero et al., 2015
  • Favero F.
  • Joshi T.
  • Marquard A.M.
  • Birkbak N.J.
  • Krzystanek M.
  • Li Q.
  • Szallasi Z.
  • Eklund A.C.

Sequenza: allele-specific copy number and mutation profiles from tumor sequencing data.Ann. Oncol. 2015; 26: 64-70

) and VarScan2 (
Koboldt et al., 2012
  • Koboldt D.C.
  • Zhang Q.
  • Larson D.E.
  • Shen D.
  • McLellan M.D.
  • Lin L.
  • Miller C.A.
  • Mardis E.R.
  • Ding L.
  • Wilson R.K.

VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing.Genome Res. 2012; 22: 568-576

). 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
  • Mitchell A.
  • Chang H.Y.
  • Daugherty L.
  • Fraser M.
  • Hunter S.
  • Lopez R.
  • McAnulla C.
  • McMenamin C.
  • Nuka G.
  • Pesseat S.
  • et al.

The InterPro protein families database: the classification resource after 15 years.Nucleic Acids Res. 2015; 43: D213-D221

).

HLA Types and Neoepitopes

The four-digit HLA class 1 and 2 types of each patient were called using ATHLATES (
Liu et al., 2013
  • Liu C.
  • Yang X.
  • Duffy B.
  • Mohanakumar T.
  • Mitra R.D.
  • Zody M.C.
  • Pfeifer J.D.

ATHLATES: accurate typing of human leukocyte antigen through exome sequencing.Nucleic Acids Res. 2013; 41: e142

) 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 that 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 standalone version of the programs NetMHCpan v2.8 (
Hoof et al., 2009
  • Hoof I.
  • Peters B.
  • Sidney J.
  • Pedersen L.E.
  • Sette A.
  • Lund O.
  • Buus S.
  • Nielsen M.

NetMHCpan, a method for MHC class I binding prediction beyond humans.Immunogenetics. 2009; 61: 1-13

,
Nielsen et al., 2007
  • Nielsen M.
  • Lundegaard C.
  • Blicher T.
  • Lamberth K.
  • Harndahl M.
  • Justesen S.
  • Røder G.
  • Peters B.
  • Sette A.
  • Lund O.
  • Buus S.

NetMHCpan, a method for quantitative predictions of peptide binding to any HLA-A and -B locus protein of known sequence.PLoS ONE. 2007; 2: e796

) and NetMHCIIpan v3.0 (
Karosiene et al., 2013
  • Karosiene E.
  • Rasmussen M.
  • Blicher T.
  • Lund O.
  • Buus S.
  • Nielsen M.

NetMHCIIpan-3.0, a common pan-specific MHC class II prediction method including all three human MHC class II isotypes, HLA-DR, HLA-DP and HLA-DQ.Immunogenetics. 2013; 65: 711-724

), 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: (1) predicted binding affinity ≤ 500nM and (2) 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 fragments per kilobase of transcripts per million mapped reads [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
  • Hodis E.
  • Watson I.R.
  • Kryukov G.V.
  • Arold S.T.
  • Imielinski M.
  • Theurillat J.P.
  • Nickerson E.
  • Auclair D.
  • Li L.
  • Place C.
  • et al.

A landscape of driver mutations in melanoma.Cell. 2012; 150: 251-263

,
Cancer Genome Atlas Network, 2015

Cancer Genome Atlas Network. Genomic Classification of Cutaneous Melanoma.Cell. 2015; 161: 1681-1696

) to estimate each gene’s background mutation frequency. Significance was computed by Fisher’s exact test followed by FDR adjustment for multiple hypothesis testing. We listed genes that fulfilled the criteria: (1) recurrence in at least 25% of the responder/non-responder, (2) occurrence of at most once in the opposite group, and (3) Fisher’s exact test FDR-adjusted p value ≤ 0.05. These genes were illustrated in Figure 1A, and all genes that fulfilled the first and second criteria 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
  • Hodis E.
  • Watson I.R.
  • Kryukov G.V.
  • Arold S.T.
  • Imielinski M.
  • Theurillat J.P.
  • Nickerson E.
  • Auclair D.
  • Li L.
  • Place C.
  • et al.

A landscape of driver mutations in melanoma.Cell. 2012; 150: 251-263

,
Cancer Genome Atlas Network, 2015

Cancer Genome Atlas Network. Genomic Classification of Cutaneous Melanoma.Cell. 2015; 161: 1681-1696

).

RNA-Seq and Gene Set Enrichment

Paired-end transcriptome reads were mapped to the UCSC hg19 reference genome using Tophat2 (
Kim et al., 2013
  • Kim D.
  • Pertea G.
  • Trapnell C.
  • Pimentel H.
  • Kelley R.
  • Salzberg S.L.

TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions.Genome Biol. 2013; 14: R36

). Normalized expression levels of genes were expressed in FPKM values as generated by cuffquant and cuffnorm (
Trapnell et al., 2012
  • Trapnell C.
  • Roberts A.
  • Goff L.
  • Pertea G.
  • Kim D.
  • Kelley D.R.
  • Pimentel H.
  • Salzberg S.L.
  • Rinn J.L.
  • Pachter L.

Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks.Nat. Protoc. 2012; 7: 562-578

). The program was run with the option “–frag-bias-correct” and “–multi-read-correct” to improve sensitivity (
Roberts et al., 2011
  • Roberts A.
  • Trapnell C.
  • Donaghey J.
  • Rinn J.L.
  • Pachter L.

Improving RNA-Seq expression estimates by correcting for fragment bias.Genome Biol. 2011; 12: R22

). A gene was defined as differentially expressed between the responding and non-responding tumor groups when its median expression differed by at least 2-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 three 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., 2009a
  • Huang W.
  • Sherman B.T.
  • Lempicki R.A.

Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources.Nat. Protoc. 2009; 4: 44-57

). Enriched GO terms were selected from the GO biological process terms in DAVID's fat database (
Huang et al., 2009b
  • Huang W.
  • Sherman B.T.
  • Lempicki R.A.

Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists.Nucleic Acids Res. 2009; 37: 1-13

). 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 (
Hänzelmann et al., 2013
  • Hänzelmann S.
  • Castelo R.
  • Guinney J.

GSVA: gene set variation analysis for microarray and RNA-seq data.BMC Bioinformatics. 2013; 14: 7

) to derive the absolute enrichment scores of previously experimentally validated gene signatures as follows: (1) the C2 CGP (chemical and genetic perturbation sets), (2) the C6 and C7 subset of the Molecular Signature Database version 4.0 (
Subramanian et al., 2005
  • Subramanian A.
  • Tamayo P.
  • Mootha V.K.
  • Mukherjee S.
  • Ebert B.L.
  • Gillette M.A.
  • Paulovich A.
  • Pomeroy S.L.
  • Golub T.R.
  • Lander E.S.
  • Mesirov J.P.

Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.Proc. Natl. Acad. Sci. USA. 2005; 102: 15545-15550

), (3) self-curated MAPKi-induced gene signatures using cell lines and patient-derived tumors (C.S., L.S., W.H., and R.S.L., unpublished data), (4) post-operation wound signature (
Inkeles et al., 2015
  • Inkeles M.S.
  • Scumpia P.O.
  • Swindell W.R.
  • Lopez D.
  • Teles R.M.
  • Graeber T.G.
  • Meller S.
  • Homey B.
  • Elder J.T.
  • Gilliet M.
  • et al.

Comparison of molecular signatures from multiple skin diseases identifies mechanisms of immunopathogenesis.J. Invest. Dermatol. 2015; 135: 151-159

), and (5) melanoma invasive/proliferative signatures (
Hoek et al., 2008
  • Hoek K.S.
  • Eichhoff O.M.
  • Schlegel N.C.
  • Döbbeling U.
  • Kobert N.
  • Schaerer L.
  • Hemmi S.
  • Dummer R.

In vivo switching of human melanoma cells between proliferative and invasive states.Cancer Res. 2008; 68: 650-656

). To derive the GSVA score of each signature in each tumor sample, we computed from raw RNA-seq read counts by HTSEQ COUNT program and then normalized them to log2 counts per million (CPM) values using EdgeR (
McCarthy et al., 2012
  • McCarthy D.J.
  • Chen Y.
  • Smyth G.K.

Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation.Nucleic Acids Res. 2012; 40: 4288-4297

). We removed batch effects using the edgeR function RemoveBatchEffect when we combined RNA-seq data from multiple experiments (Figure 3A). The normalized log2 CPM values were then passed on as input for GSVA in the RNA-seq 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 zero). 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 eight). 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: (1) our current pre-anti-PD-1 tumors (n = 28), (2) pre-anti-CTLA-4 tumors (n = 42), (3) pre-MAPKi tumors (n = 32), and (4) 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, top). 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 nine non-responding tumors and one 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 IPRES score was not comparable across analyses in Figures 3A, 3B, and S3 (since the IPRES score is Z score-base), 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
  • Van Allen E.M.
  • Miao D.
  • Schilling B.
  • Shukla S.A.
  • Blank C.
  • Zimmer L.
  • Sucker A.
  • Hillen U.
  • Foppen M.H.
  • Goldinger S.M.
  • et al.

Genomic correlates of response to CTLA-4 blockade in metastatic melanoma.Science. 2015; 350: 207-211

). We curated published CD8 T cell exhaustion genes (
Wherry, 2011
  • Wherry E.J.

T cell exhaustion.Nat. Immunol. 2011; 12: 492-499

) 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 RNA-seq 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 (1) anti-CTLA-4 pretreatment cohort (
Van Allen et al., 2015
  • Van Allen E.M.
  • Miao D.
  • Schilling B.
  • Shukla S.A.
  • Blank C.
  • Zimmer L.
  • Sucker A.
  • Hillen U.
  • Foppen M.H.
  • Goldinger S.M.
  • et al.

Genomic correlates of response to CTLA-4 blockade in metastatic melanoma.Science. 2015; 350: 207-211

), (2) MAPKi pretreatment cohort (
Hugo et al., 2015
  • Hugo W.
  • Shi H.
  • Sun L.
  • Piva M.
  • Song C.
  • Kong X.
  • Moriceau G.
  • Hong A.
  • Dahlman K.B.
  • Johnson D.B.
  • et al.

Non-genomic and Immune Evolution of Melanoma Acquiring MAPKi Resistance.Cell. 2015; 162: 1271-1285

; C.S., L.S., W.H., and R.S.L., unpublished data), (3) TCGA melanoma (metastatic and primary subsets separately analyzed;
Cancer Genome Atlas Network, 2015

Cancer Genome Atlas Network. Genomic Classification of Cutaneous Melanoma.Cell. 2015; 161: 1681-1696

), (4) TCGA pancreatic ductal adenocarcinoma (TCGA, provisional, see https://tcga-data.nci.nih.gov/tcga/tcgaHome2.jsp), (5) TCGA lung adenocarcinoma ( ), (6) TCGA colorectal adenocarcinoma (
Cancer Genome Atlas Network, 2012

Cancer Genome Atlas Network. Comprehensive molecular characterization of human colon and rectal cancer.Nature. 2012; 487: 330-337

), and (7) TCGA kidney clear cell carcinoma (
Cancer Genome Atlas Research Network, 2013

Cancer Genome Atlas Research Network. Comprehensive molecular characterization of clear cell renal cell carcinoma.Nature. 2013; 499: 43-49

).

Author Contributions

W.H., J.M.Z., L.S., C.S., A.R., and R.S.L. generated, analyzed, and interpreted the 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.B.-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.

Acknowledgments

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

Accession Numbers

The accession number for the transcriptome data is GEO: GSE78220. The whole-exome sequencing data has been deposited to the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra), under the accession numbers SRA: SRP067938 (UCLA samples) and SRA: SRP090294 (Vanderbilt samples).

Supplemental Information

References

Cancer Genome Atlas Network

Comprehensive molecular characterization of human colon and rectal cancer.

Nature. 2012; 487: 330-337

Cancer Genome Atlas Network

Genomic Classification of Cutaneous Melanoma.

Cell. 2015; 161: 1681-1696

Cancer Genome Atlas Research Network

Comprehensive molecular characterization of clear cell renal cell carcinoma.

Nature. 2013; 499: 43-49

Cancer Genome Atlas Research Network

Comprehensive molecular profiling of lung adenocarcinoma.

Nature. 2014; 511: 543-550

Metastasis is regulated via microRNA-200/ZEB1 axis control of tumour cell PD-L1 expression and intratumoral immunosuppression.

Nat. Commun. 2014; 5: 5241

The mutually regulatory loop of epithelial-mesenchymal transition and immunosuppression in cancer progression.

OncoImmunology. 2015; 4: e1002731

Sequenza: allele-specific copy number and mutation profiles from tumor sequencing data.

Ann. Oncol. 2015; 26: 64-70

Safety and tumor responses with lambrolizumab (anti-PD-1) in melanoma.

N. Engl. J. Med. 2013; 369: 134-144

GSVA: gene set variation analysis for microarray and RNA-seq data.

BMC Bioinformatics. 2013; 14: 7

A landscape of driver mutations in melanoma.

Cell. 2012; 150: 251-263

In vivo switching of human melanoma cells between proliferative and invasive states.

Cancer Res. 2008; 68: 650-656

Unraveling the mechanism of BRCA2 in homologous recombination.

Nat. Struct. Mol. Biol. 2011; 18: 748-754

NetMHCpan, a method for MHC class I binding prediction beyond humans.

Immunogenetics. 2009; 61: 1-13

CCR 20th Anniversary Commentary: Immune-Related Response Criteria--Capturing Clinical Activity in Immuno-Oncology.

Clin. Cancer Res. 2015; 21: 4989-4991

Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources.

Nat. Protoc. 2009; 4: 44-57

Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists.

Nucleic Acids Res. 2009; 37: 1-13

Non-genomic and Immune Evolution of Melanoma Acquiring MAPKi Resistance.

Cell. 2015; 162: 1271-1285

Comparison of molecular signatures from multiple skin diseases identifies mechanisms of immunopathogenesis.

J. Invest. Dermatol. 2015; 135: 151-159

NetMHCIIpan-3.0, a common pan-specific MHC class II prediction method including all three human MHC class II isotypes, HLA-DR, HLA-DP and HLA-DQ.

Immunogenetics. 2013; 65: 711-724

TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions.

Genome Biol. 2013; 14: R36

VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing.

Genome Res. 2012; 22: 568-576

PD-1 Blockade in Tumors with Mismatch-Repair Deficiency.

N. Engl. J. Med. 2015; 372: 2509-2520

ATHLATES: accurate typing of human leukocyte antigen through exome sequencing.

Nucleic Acids Res. 2013; 41: e142

A Patient-Derived, Pan-Cancer EMT Signature Identifies Global Molecular Alterations and Immune Target Enrichment Following Epithelial-to-Mesenchymal Transition.

Clin. Cancer Res. 2016; 22: 609-620

Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation.

Nucleic Acids Res. 2012; 40: 4288-4297

The InterPro protein families database: the classification resource after 15 years.

Nucleic Acids Res. 2015; 43: D213-D221

The parallel lives of angiogenesis and immunosuppression: cancer and other tales.

Nat. Rev. Immunol. 2011; 11: 702-711

NetMHCpan, a method for quantitative predictions of peptide binding to any HLA-A and -B locus protein of known sequence.

PLoS ONE. 2007; 2: e796

Loss of PTEN promotes resistance to T cell-mediated immunotherapy.

Cancer Discov. 2016; 6: 202-216

Oncotator: cancer variant annotation tool.

Hum. Mutat. 2015; 36: E2423-E2429

Cancer immunology. Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer.

Science. 2015; 348: 124-128

Improving RNA-Seq expression estimates by correcting for fragment bias.

Genome Biol. 2011; 12: R22

Cancer as an overhealing wound: an old hypothesis revisited.

Nat. Rev. Mol. Cell Biol. 2008; 9: 628-638

Immune checkpoint targeting in cancer therapy: toward combination strategies with curative potential.

Cell. 2015; 161: 205-214

Acquired resistance and clonal evolution in melanoma during BRAF inhibitor therapy.

Cancer Discov. 2014; 4: 80-93

Genetic basis for clinical response to CTLA-4 blockade in melanoma.

N. Engl. J. Med. 2014; 371: 2189-2199

Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.

Proc. Natl. Acad. Sci. USA. 2005; 102: 15545-15550

Safety, activity, and immune correlates of anti-PD-1 antibody in cancer.

N. Engl. J. Med. 2012; 366: 2443-2454

Immunogenicity of somatic mutations in human gastrointestinal cancers.

Science. 2015; 350: 1387-1390

Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks.

Nat. Protoc. 2012; 7: 562-578

PD-1 blockade induces responses by inhibiting adaptive immune resistance.

Nature. 2014; 515: 568-571

Genomic correlates of response to CTLA-4 blockade in metastatic melanoma.

Science. 2015; 350: 207-211

Control of the immune response by pro-angiogenic factors.

Front. Oncol. 2014; 4: 70

T cell exhaustion.

Nat. Immunol. 2011; 12: 492-499

Guidelines for the evaluation of immune therapy activity in solid tumors: immune-related response criteria.

Clin. Cancer Res. 2009; 15: 7412-7420

Article Info

Publication History

Published: March 17, 2016; corrected online January 18, 2017

Accepted: February 29, 2016

Received in revised form: January 30, 2016

Received: December 25, 2015

IDENTIFICATION

DOI: 10.1016/j.cell.2016.02.065

Copyright

© 2016 Elsevier Inc. Published by Elsevier Inc.

User License

Elsevier user license |
How you can reuse

Elsevier user license

Permitted

For non-commercial purposes:

  • Read, print & download
  • Text & data mine
  • Translate the article

Not Permitted

  • Reuse portions or extracts from the article in other works
  • Redistribute or republish the final article
  • Sell or re-use for commercial purposes

Elsevier's open access license policy

ScienceDirect

Access this article on ScienceDirect

Linked Article

Cell Press Commenting Guidelines

To submit a comment for a journal article, please use the space above and note the following:

  • We will review submitted comments within 2 business days.
  • This forum is intended for constructive dialog. Comments that are commercial or promotional in nature, pertain to specific medical cases, are not relevant to the article for which they have been submitted, or are otherwise inappropriate will not be posted.
  • We recommend that commenters identify themselves with full names and affiliations.
  • Comments must be in compliance with our Terms & Conditions.
  • Comments will not be peer-reviewed.


Related news

Bluebottles woolgoolga accommodation
Placar para basquete de cadeira
Demarreur electrique scooter tires
Frasqueiras para viagem mercado livre roupas
Semantique synonyme de faire
Malacates de que sirve querer acordes de bajo
L'affirmation de soi chalvin dominique moceanu
Waipu motel accommodation in adelaide
Inis mor accommodation tripadvisor san francisco
Panoramique dentaire 3d movies