Predictive Immune-Checkpoint Blockade Classifiers Identify Tumors Responding to Inhibition of PD-1 and/or CTLA-4

Although innovations in immune-checkpoint blockade (ICB) of metastatic melanoma have led to improved responses in a considerable number of patients, ICB resistance before start of treatment or acquired later remain formidable challenges (1). Although resistance is observed in the majority of melanoma patients treated with either anti–PD-1 or anti–CTLA-4 monotherapy, combining anti–PD-1 with anti–CTLA-4 has significantly improved these response rates (2–5). However, this comes at a cost, as combination ICB is associated with an increase in the number of patients suffering from immune-related adverse events (irAE), while also the cost of combining these treatments is exceedingly high. Therefore, there is a dire need to stratify individual patients before the start of treatment for their likelihood of responding to either anti–PD-1 or anti–CTLA-4 monotherapy, or to anti–PD-1 + anti–CTLA-4 combination therapy. In particular, clinicians would be aided if they could make an informed decision on whether a patient would need anti–CTLA-4 in addition to anti–PD-1 antibody.

Although it is established that anti–PD-1 and anti–CTLA-4 monotherapy–induced immune responses are driven by distinct mechanisms (6), this has not yet been captured faithfully in predictive gene-expression signatures. Whereas currently available signatures predict general ICB response, they fail to distinguish between response to either anti–CTLA-4 or anti–PD-1 in the combination setting (7, 8). This would be an important feature, given that independent response prediction would allow for the stratification of patients likely responding to either anti–PD-1 or anti–CTLA-4, as well as identify those patients who would benefit from anti–CTLA-4 in addition to anti–PD-1 treatment. Furthermore, these classifiers lack the ability to distinguish between tumor cell–intrinsic and tumor cell–extrinsic [tumor microenvironment (TME) versus stromal] signals (7–9). However, there is ample evidence that both tumor cell–extrinsic signals (e.g., the abundance of immune cells and specific immune cell types; refs. 7, 10) and tumor cell–intrinsic signals (e.g., specific cellular states like phenotype switching; refs. 9, 11) influence response to ICB.

Therefore, we hypothesized that a framework of predictive signatures ought to comprise separate tumor cell–intrinsic and –extrinsic elements. To disentangle these signals, we combined XenofilteR (12), a software tool we developed previously that distinguishes sequence reads from mouse and human origin, with our stage IV metastatic melanoma patient-derived xenograft (PDX) cohort (13). This allowed us to develop two gene-expression signatures (“InTumor” and “ExTumor”), which when used independently predict clinical response to either anti–PD-1 or anti–CTLA-4 monotherapy. Most importantly, we show that when used in conjunction, these signatures can be used to stratify patients before start of treatment for their likelihood of responding to either anti–PD-1 or anti–CTLA-4 monotherapy, or to anti–PD-1 + anti–CTLA-4 combination therapy.

The collection and use of human tissue were approved by the Medical Ethical Review Board [institutional review board (IRB)] of the Antoni van Leeuwenhoek. All patients provided written consent. The study was conducted in accordance with ethical guidelines (1964 Helsinki declaration and its later amendments or comparable ethical standards). Animal experiments were approved by the animal experimental committee of the institute and performed according to Dutch law. PDX were generated as described in Kemper and colleagues (13). PDX-derived cell lines were generated after enzymatic dissociation (described in ref. 1) by growing the cells in 2D in 10% fetal bovine serum (FCS, Sigma), 2 mmol/L glutamine, 100 U/mL penicillin, and 0.1 mg/mL streptomycin (all Gibco) under standard conditions. A full report of PDX relevant information following the guidelines described in Minimal Information for Patient-Derived Tumor Xenograft Models (PDX-MI) is available (Supplementary Table S1; ref. 14).

STR profiling was performed on gDNA from all PDX and PDX-derived cell lines (Supplementary Table S2). gDNA was isolated by using the DNA Easy Blood and Tissue Kit (Qiagen) according to the manufacturer's protocol. STR profiling was performed using the Powerplex 16 HS kit (DC2101, Promega) according to the manufacturer's protocol. The samples were run on an ABI3700 and analyzed using GeneMapper. In addition to the STR profiling, we performed an SNP analysis on both the RNA sequence data from the PDX cell lines as well as on RNA sequence data of the first passage PDX. The percentage of overlap in SNPs was tested between all cell lines, between all first passage PDX samples and between PDX and PDX cell lines, which confirmed PDX cell line origin (Supplementary Fig. S9).

RNA was isolated from PDX and PDX-derived cell lines using TRIzol according to the manufacturer's protocol. Next, samples were cleaned twice by using the RNeasy MinElute Cleanup Kit (74204, Qiagen) according to the manufacturer's protocol, before they were submitted for quality control by using the Agilent 2100 Bioanalyzer system according to the manufacturer's protocol. A total of 95 (PDX) and 22 (PDX-derived cell lines) samples were sequenced in two batches. PDX samples were uniquely barcoded and pooled into a single-stranded library and sent for sequencing. To minimize the possible lane effect, we adopted the strategy of pooling all uniquely labeled PDX samples and running the pool on multiple lanes on Illumina HiSeq 2000 sequencers. An average 56 million unique read pairs were sequenced per sample ranging from 33–78 million. The same pooled strategy was used for the 22 PDX cell lines, but these samples were sequenced to an average of 39 million unique reads pairs, ranging from 27–54 million. All samples were mapped with Tophat2 v2.1.0 (15) with the following parameters: –library-type fr-firststrand -g 1 -p 8 -G ENSEMBL_Annotation_v82.gtf. All samples were mapped to human reference GRCh38 (ENSEMBL v82) and mouse reference genome GRCm38 (ENSEMBL v82).

Sequence data from PDX samples are contaminated with sequence reads of mouse origin (12). Therefore, computational removal of these sequence reads is pivotal. Previously we developed XenofilteR, which showed accurate filtering, both for RNA sequencing (RNA-seq) and DNA-seq data. To further test the performance of XenofilteR in our RNA-seq samples, we mixed the sequence reads from a mouse cell line (∼25%) with the sequence reads from a human cell line (∼75%). Both the original mouse, human and mixed sequence reads were mapped to the human and mouse reference genome. The resulting bam file from the mixed sample was subsequently filtered using XenofilteR (version 1.4, github.com/PeeperLab/XenofilteR) with default settings. Furthermore, the RNA isolation protocols and sequencing for these two cell lines were identical to that of the PDX RNA-seq described here. Filtering with XenofilteR retained 99.14% of human sequence reads while removing 99.7% of mouse reads (Supplementary Table S3). The removal of mouse sequence reads in the mixed sample resulted in read count data similar to the read counts when only the human sequence reads were mapped. Thus, XenofilteR shows accurate filtering of mouse sequence reads while retaining sequence reads of human origin.

Next, we applied XenofilteR (12) to the 95 PDX samples. The percentage of reads from mouse origin in each PDX RNA-seq sample was measured by counting the unique read pairs in the bam files with mouse-specific reads as generated with XenofilteR divided by the total number of unique read pairs (sum of unique read pairs mapped to either mouse or human after filtering with XenofilteR).

To generate read count data, the filtered bam files were name sorted with picard followed by counting reads with HTseq-count (HTSeq-0.6.1p1; ref. 16) with settings: -m intersection-nonempty -a 10 -i gene_id -s reverse -f bam. Count data generated with HTseq-count were analyzed with DESeq2 (17). Centering of the normalized gene-expression data is performed by subtracting the row means and scaling by dividing the columns by the standard deviation to generate a z-score.

The curated tumor cell–intrinsic read count data of 95 PDX samples were used to identify the gene-expression signals from the human TME in total patient tumors (tumor + TME). For this, we used a data set of 133 patient melanoma samples (18) for which preprocessing was performed identical as for the PDX samples. To identify the human TME-specific signals from the patient samples, we selected the genes with a minimal read count of 20 and with at least a log₂-fold difference of 8 in expression between curated tumor cell–intrinsic signals (tumor) from PDX and total patient tumors (tumor + TME). This yielded 787 TME-specific genes.

To verify these TME-specific signals, we tested the expression in single-cell RNA sequence data (scSeq) downloaded from the NCBI GEO database (GSE72056; ref. 19). Count data from each cell were normalized to 1 million reads per cell. Classification into somatic and germline and specific cell types was taken from the original publication (6). The expression of the 787 TME-specific genes was tested in the scSeq data from the 3 patients (sample ID 79, 80, and 88; Fig. 2C; Supplementary Fig. S1B) with the highest number of cells sequenced for tumor cells as well as TME cells [T cells, endothelial cells, B cells, cancer-associated fibroblasts (CAF), and natural killer (NK) cells].

To further specify the 787 TME-specific signals, we performed cluster analysis using the The Cancer Genome Atlas (TCGA) skin cutaneous melanoma (SKCM) data set (Supplementary Fig. S2A). For this, raw read count data were downloaded using the “TCGAbiolinks” package (20) and preprocessed and normalized similar to the PDX data set. Cluster analysis was performed using the metastatic melanoma samples from the TCGA database (21) and the 767 genes (Supplementary Fig. S2B). This analysis revealed two clusters with highly coregulated genes. The first cluster was dominated by keratin genes and signals that originate from keratinocytes; the second cluster (361 genes) was highly enriched for immune cell–related signaling. The latter was confirmed by Gene Ontology analysis (Supplementary Fig. S2C) showing the GO pathways “inflammatory response” and “Immune response” as top hits. These 361 genes comprise the ExTumor signature (Supplementary Table S7).

The curated tumor cell–intrinsic read count data of 95 PDX samples (13) were used to perform a principal component analysis (PCA) to identify the signals that explain the largest variation in the data set (Supplementary Fig. S3A). The first three principal components (PC) together explained 22.5% of the variance in the data set (11.3%, 6.9%, and 4.4%, respectively). These three PCs were used to generate three gene-expression signatures by selecting the genes for each of the three principle components with the highest rotation value and comprise 39, 89, and 76 genes respectively (Supplementary Table S4; Supplementary Fig. S3B). Next, we performed gene set enrichment analysis (GSEA) to correlate the expression of the three gene sets with response to immune-checkpoint blockade in two data sets (monotherapy anti–PD-1, ref. 9; monotherapy anti–CTLA-4, ref. 22). The only significant association between the PCs and response to anti–PD-1 was the second PC with response to anti–PD-1. The gene set generated using PC2 comprises 89 genes and is named the InTumor signature.

The InTumor signature was developed on PDX samples without signals from the TME. However, this does not preclude the possibility that these genes are expressed also in stromal cells. To rule out any possible influence from signals of TME origin in patient samples, we checked the average expression of the signature genes in both tumor and stromal cells using scSEQ data. For each gene in the InTumor signature, we calculated the average expression in tumor cells and the stromal cells using scSEQ data (19). Samples were normalized to 1M reads per sample, and the three samples with high number of tumor as well as stromal cells were used for the analysis (samples 79, 80, and 88). Only genes with an arbitrary 10-fold higher expression in tumor cells compared with stromal cells comprise the InTumor exclusive signature (Supplementary Fig. S5A). Out of the 86 genes, we identified 14 genes with a 10-fold higher expression in tumor compared with stromal cells (Supplementary Table S5). These 14 genes were used to measure the InTumor signature in patient samples.

Signature values (ExTumor and InTumor) were generated by calculating the average gene expression of the signature genes per sample. This was performed on the z-scores as calculated per data set. Because the z-score is a mean-centered expression value over the data set, this will result in an equal distribution of samples in a data set that are divided into signature^HI and signature^LO groups. These cutoffs are not trained for and are used identically in all data sets including the ones with patients treated with the combination of anti–CTLA-4 + anti–PD-1 and the separation of patients in the ICB response quadrant.

Gene-expression data for the TCGA (21) melanoma samples were downloaded as raw count data using TCGAbiolinks (20) and analyzed with DESeq2 (17). Centering of the normalized gene-expression data per data set was performed by subtracting the row means and scaling by dividing the columns by the standard deviation. The generated z-score was used to calculate the signature values. For the analysis, only metastatic samples were selected.

The collection and use of human tissue were approved by the Medical Ethical Review Board (IRB) of the Antoni van Leeuwenhoek. All patients provided written consent. The study was conducted in accordance with ethical guidelines (1964 Helsinki declaration and its later amendments or comparable ethical standards). Gene expression of 51 samples was assessed using the nCounter PanCancer Immune Profiling assay by NanoString Technologies. RNA was isolated from formalin-fixed, paraffin-embedded tumor samples using the QIAgen AllPrep DNA/RNA kit on the QIAcube according to standard manufacturer's protocol. Preprocessing and data normalization were performed using the nSolver (Nanostring) software using default settings. The ExTumor signature value was calculated using the 151 genes available on the platform. All samples were taken before start of treatment with anti–CTLA-4. Available clinical variables included survival data for all patients and response to anti–CTLA-4 for 42 out of the 51 patients. Patients with progressive disease (PD) were indicated as nonresponders. Patients with stable disease (SD), partial response (PR), and complete response (CR) were indicated as responders.

For all patient data sets that were used to validate the ExTumor and InTumor signatures, raw sequence data (fastq files) were downloaded. All samples were mapped to the human reference genome (Homo.sapiens.GRCh38.v82) using STAR (2.6.0c; ref. 23) using 2-pass mapping with default settings. Count data generated with HTseq-count were analyzed with DESeq2 (17). Centering of the normalized gene-expression data per data set was performed by subtracting the row means and scaling by dividing the columns by the standard deviation. The generated z-score was used to calculate the signature values. For all data sets the identical settings were used except for the OpACIN sequence data as these were sequenced single-end, whereas all other data sets were sequence paired-end. For the validation analyses, only cutaneous melanomas were used, omitting the uveal and mucosal melanoma samples. These samples were identified based on the clinical data as reported in each publication.

To assess tumor mutation burden (TMB) for the melanoma patient data sets [van Allen (22), Riaz (24) and OpACIN (25)], raw sequence reads were downloaded (fastq) and mapped to the human reference genome using BWA (26). Subsequently, the GATK4 pipeline was used to mark duplicate sequence reads and recalibrate base quality scores followed by mutation calling using Mutect2 (27). TMB was calculated by summarizing the number of nonsynonymous mutations with a minimum of 15 sequence reads and a variant allele frequency of 5%. For the Hugo data sets, raw sequence reads were not available; therefore, the TMB was taken as reported in the supplementary data of the Hugo and colleagues paper (9).

Three previously published gene-expression signatures that were designed to predict response to ICB were tested and its predictive power compared with the predictive power of the InTumor signature. The cytolytic score (10) was defined as the average expression (z-score) of the genes PRF1 and GZMA for each sample. The IFNg (8) score was defined as the average expression (based on the z-score) of 10 genes (CXCL11, IDO1, CXCL9, IFNG, PRF1, STAT1, HLA-DRA, CCR5, GZMA, and CXCL10). The EMT score was based on a pan-cancer EMT signature (28) with 76 genes (24 upregulated in epithelial cells and 52 upregulated in mesenchymal cells). The EMT ratio represents the ratio between the two gene sets. To calculate the IMPRES (7) score, we generated FPKM values based on the read count data for each data set. For each sample, the ratio between 15 gene pairs was assessed (“CD274/VSIR,” “CD28/CD276,” “CD86/TNFRSF4,” “CD86/CD200,” “CTLA4/TNFRSF4,” “PDCD1/TNFRSF4,” “CD80/TNFSF9,” “CD86/HAVCR2,” “CD28/CD86,” “CD27/PDCD1,” “CD40/CD274,” “CD40/CD80,” “CD40/CD28,” “CD40/CD274,” and “TNFRSF14/CD86”). The IMPRES score is the sum of the number of pairs where the expression of the first gene is higher than the second. This value will be between 0 and 15.

Identical preprocessing was used for each data set to acquire read count data using STAR and HTseq-count as described above. Based on these read counts, we calculated the FPKM values, and from those the IMPRES scores. These results deviate from the original scores as reported in Auslander and colleagues (7). In this publication, the presented IMPRES scores were based either on normalized read count data (Riaz and colleagues) or FPKM data (Hugo and colleagues) as provided in the original publication. We chose to use one single method to calculate the IMPRES scores for all data sets based on FPKM values. FPKM was chosen because the IMPRES signature was trained on FPKM data in the original publication.

GSEA based on the principal components was performed using the BROAD javaGSEA standalone version (http://www.broadinstitute.org/gsea/downloads.jsp) and the curated “hallmark genesets” (http://software.broadinstitute.org/gsea/msigdb/collections.jsp). Analyses were run using 10,000 permutations. Genes were ranked based on the Signal2Noise metric. Gene sets with an FDR <0.1 were considered significant.

cDNA was generated using either Superscript III Reverse Transcriptase or Maxima First-Strand cDNA synthesis kit for RT-PCR according to manufacturer's protocol. qPCR was performed as described previously (29).

CD8, PD-L1, and vimentin stainings were performed manually. FFPE slides were deparaffinized and subjected to antigen retrieval using TRIS/EDTA buffer. Slides were blocked in 4% BSA and 4% normal goat serum (NGS) in PBS for 30 minutes, and subsequently incubated overnight at 4°C with CD8 (M7103, DakoCytomation), PD-L1 (1795–1-AP, Proteintech), or vimentin antibody (M0725, DAKO) at a dilution of 1:2,000 in 1% BSA and 1.25% NGS in PBS. After rinsing with 0.05% Tween-20 in PBS, the slides were incubated with secondary antibody (Goat, Southern Biotech) for 30 minutes at a dilution of 1:100, and then incubated with streptavidin/HRP (DakoCytomation; P0397) for 30 minutes at a dilution of 1:200, both in 1% BSA and 1.25% NGS in PBS. Stainings were visualized using 3,3-diaminobenzidine (DAB) chromophore (Sigma; D-5905) and counterstained with hematoxylin. Sections were reviewed with a Zeiss Axioskop2 Plus microscope (Carl Zeiss Microscopy) and images were captured with a Zeiss AxioCam HRc digital camera and processed with AxioVision 4 software (both from Carl Zeiss Vision).

All melanoma cell lines used for in vivo adoptive T-cell transfer (ACT) models were obtained from the Peeper laboratory cell line stock. PDX were generated as described previously (13). Cell lines were regularly confirmed to be Mycoplasma-free by PCR. Melanoma cell lines were cultured in DMEM (Gibco), with fetal bovine serum (Sigma), and 100 U/mL penicillin and 0.1 mg/mL streptomycin (both Gibco) under standard conditions. MART-126-35 and HLA-A2 were introduced in the cell lines (D10, A375, BLM, and SkMel-147) using lentiviral constructs. Constructs for lentivirus were packaged using two helper plasmids (psPax and MS2G, Addgene) in HEK293T cells. MART-126-35-Katushka–positive cells were sorted by flow cytometry.

MART-1 (1D3) TCR retrovirus was produced in a packaging cell line as described previously (30). Peripheral blood mononuclear cells were isolated from healthy donor buffycoats (Sanquin) by density gradient centrifugation using Lymphoprep (Stem Cell Technologies). CD8⁺ T cells were purified using CD8 Dynabeads (Thermo Fisher Scientific), activated for 48 hours on a nontissue culture-treated 24-well plate that was precoated overnight with αCD3 and αCD28 antibodies (eBioscience, 16-0037-85 and 16-0289-85) at 2×10⁶ per well. Activated CD8 T cells were harvested and mixed with TCR retrovirus and spinfected on a Retronectin-coated (Takara, 25 μg per well) nontissue culture-treated 24-well plate for 2 hours at 2,000 × g. After 24 hours, T cells were harvested and maintained in RPMI (Gibco) containing 10% human serum (One Lamda), 100 units/mL of penicillin, 100 μg/mL of streptomycin, 100 units/mL IL2 (Proleukin, Novartis), 10 ng/mL IL7 (ImmunoTools), and 10 ng/mL IL15 (ImmunoTools).

All animal studies were approved by the animal ethics committee of the Netherlands Cancer Institute (NKI) and performed in concordance with ethical and procedural guidelines established by the NKI and Dutch legislation. Xenograft studies were performed as previously described (30); 1 × 10⁶ tumor cells were injected subcutaneously into NSG-β2Mnull mice (Jax). Growth was monitored three times per week with calipers, and tumor size was calculated using the following formula: ½ × length (mm) × width (mm). When tumors reached indicated sizes, mice were randomized over different treatment groups in a blinded fashion and were administered 5×10⁶ human CD8 T cells, corrected for transduction efficiency, intravenously in the tail vein. In vivo persistence of T cells was stimulated by administering 100.000 U IL2 (Proleukin, Novartis) intraperitoneally daily for three consecutive days. All experiments ended for individual mice either when the tumor volume exceeded 1,000 mm³, when the tumor showed ulceration, in case of serious clinical illness, when the tumor growth blocked the movement of the mouse, or when tumor growth assessment had been completed.

RNA sequence data of the 95 PDX samples and 22 PDX-derived cell lines have been deposited in the NCBI GEO database under accession number: GSE129127.

TCGA gene-expression data were downloaded using the R/Bioconductor package TCGAbiolinks (20) as read count values. RNA-seq data from the Hugo and colleagues papers (9, 18) were downloaded from NCBI's GEO database (GSE78220 and GSE65186). RNA sequence data from the OpACIN trial (25) are available from the EGA database: EGAS00001003099. Read count data from the CA209–038 trial data are available as Supplementary Table S7. Single-cell RNA sequence data were downloaded from the NCBI GEO database: GSE72056.

Anti–CTLA-4 and anti–PD-1 treatments differentially affect tumor and TME (6, 31), and hence we argued that both effects should be accounted for when building gene-expression signatures to predict clinical response. To distinguish between tumor cell–intrinsic and –extrinsic gene expression, we exploited our previously established stage IV metastatic melanoma PDX cohort (n = 95; ref. 13). Computational analysis comparing melanoma RNA expression of both tumor and TME genes from patients (from TCGA, n = 367) with genes from these PDX demonstrated that mouse TME had almost completely replaced human TME upon PDX outgrowth after the first passage (Fig. 1A). This was confirmed by the gene expression–based TME deconvolution tool MCPcounter and by IHC staining with species-specific antibodies (Fig. 1B–D; Supplementary Fig. S1A).

This observation allowed us to use an algorithm we recently developed (XenofilteR; ref. 12) to distinguish RNA-seq reads of mouse (TME) and human (tumor cell) origin in PDX (Supplementary Table S2). By removing the sequence reads of mouse origin, we obtained human-specific tumor cell–intrinsic gene-expression values from PDX. These curated signals allowed us to subsequently computationally subtract tumor cell–intrinsic gene-expression signals from bulk tumor (that is, tumor cell + TME) gene-expression signals from patients in order to establish a TME gene-expression signature (Fig. 2A). This yielded a set of 767 curated genes expressed in patients' bulk melanomas, which were (almost) lacking from PDX (Fig. 2B; Supplementary Table S6). Single-cell sequencing (scSeq) data (19) confirmed that these genes are predominantly expressed by nontumor cells such as T, B, NK, and endothelial cells (Fig. 2C; Supplementary Fig. S1B). These data show that by applying the XenofilteR algorithm to PDX data, the RNA signals from tumor cells and TME can be computationally dissected.

In the TCGA, patient melanomas have been classified based on three genetic signatures (“MITF-low,” “Keratin,” and “Immune”; ref. 21). To determine the contribution to this classification of the 767 TME genes that we identified above, we performed a cluster analysis on the TCGA melanoma samples using these TME and immune-related genes only. This revealed a remarkably high overlap with the three melanoma subgroups that had previously been classified using bulk tumor signals (Fig. 2D). Cluster analysis with the TME genes in a cohort of 57 samples from stage IV melanoma patients (32) recapitulated these findings (Supplementary Fig. S1C). Thus, by computationally subtracting tumor cell–intrinsic gene-expression signals from those in bulk tumor, we uncovered hundreds of TME-specific gene signals, which were sufficient to recapitulate the classification of melanoma subgroups in two independent tumor data sets. These results illustrate the importance of the signals from the TME, demonstrating that the TCGA classification of melanomas is largely based on gene-expression signals that are derived from the tumor signals in the TME, rather than the tumor cell–intrinsic ones.

The main goal of this study was to generate a predictive classifier for combination ICB treatment. Previously, ICB responses were shown to be affected by both tumor cell–intrinsic and TME factors (6, 31), arguing that a predictive classifier should recapitulate both distinct features. When inspecting the 767 TME genes, we noted that 361 genes showed a high degree of coregulation (Supplementary Fig. S2A and S2B) and were immune related (Supplementary Fig. S2C). The expression of these 361 genes was highly correlated with that of multiple immune cell types, including T cells (R = 0.9) and B-cells (R = 0.7; Fig. 2E). With these immune-related genes, which is a measure for the abundance of immune cells in the sample (33), we built a so-called ExTumor signature (Supplementary Fig. S2D; Supplementary Table S7) and determined whether it is capable of predicting clinical outcome of melanoma patients treated with ICB, whether anti–PD-1 or anti–CTLA-4.

We applied GSEA on gene-expression data from a cohort of patient samples obtained before the start of treatment with anti–PD-1 monotherapy (9). The ExTumor signature failed to show a correlation with the clinical outcome of this patient cohort (Supplementary Fig. S2E). In contrast, when tested on tumors prior to receiving anti–CTLA-4 monotherapy, we observed a significant association between high expression of ExTumor signature genes and clinical response (FDR < 0.02; Fig. 2F and G), but not with the TMB (Supplementary Fig. S2F). This association was validated in an independent cohort of patient samples before start of treatment with anti–CTLA-4 (Fig. 2H; ref. 22). Again, the ExTumor signature was associated with clinical outcome, showing significantly higher expression in the patients with clinical benefit (P < 0.01) and a better overall survival for those patients with the highest expression (highest quartile) of the ExTumor signature (log-rank P = 0.028; Fig. 2I). Thus, the ExTumor signature is strongly correlated with response to anti–CTLA-4, but not anti–PD-1 treatment.

Given that the ExTumor signature failed to predict response to anti–PD-1 therapy, we investigated whether, instead, this response could be predicted by tumor cell–intrinsic factors. We again used the curated PDX signals (Fig. 3A) and applied a PCA (which identifies gene sets that explain the largest variance) of the tumor cell–intrinsic RNA expression data (Supplementary Fig. S3A). One principal component (PC2) was associated with immune-related signaling pathways (TGFb, IL2/STAT5, and IFNg; Supplementary Fig. S3B) in the PDX set and showed an association with response to anti–PD-1 in patient samples before start of treatment (FDR = 0.097; Fig. 3B). The genes from this PC were used to build a so-called InTumor gene-expression signature (Supplementary Fig. S3C). They were enriched for genes encoding cell membrane proteins and include genes involved in “cell–cell signaling” (NRP1, IL1B), “inflammatory response” [neuronal growth factor receptor (NGFR), CCR1], and “signal transduction” (ITPKA; Supplementary Table S8). We confirmed that the InTumor signature was specific for tumor cell (rather than TME) genes, as shown by (i) its high concordance between PDX and matched PDX cell lines (Supplementary Fig. S4A and S4B), (ii) IHC with species-specific antibodies (Supplementary Fig. S3D–S3F), and (iii) lack of a significant correlation between the InTumor signature expression and immune cells (Supplementary Fig. S3G) when using the TCGA database.

We next determined whether the InTumor signature predicts response to anti–PD-1. To preclude any effect from nontumor cells (which are abundantly present in bulk patient tumor sequence data), we first reduced the number of genes of this signature to the 14 showing at least 10 times higher expression in tumor compared with TME (Supplementary Fig. S5A). We then performed GSEA on gene-expression data from a cohort of patient samples obtained prior to anti–PD-1 monotherapy (9). This analysis revealed a significantly larger number of nonresponding patients in the InTumor^HI group (Χ² P < 0.05; Fig. 3C). Furthermore, the InTumor signature outperformed three previous signatures [IMPRES (7), IFNg (8), and cytolytic score (Fig. 3D; ref. 10)] and previously described markers [CD8A and CD274 (PD-L1); Supplementary Fig. S5B]. Only an EMT signature outperformed the InTumor signature in this data set. Splitting the patient samples into two groups based on the median expression resulted in a response rate of 75% in the InTumor^LO patients and 40% in the InTumor^HI patients (Fig. 3E). The predictions were significant, yet imperfect. This could be partly explained by the notion that misclassified predicted responders showed a significantly lower TMB (Fig. 3F), in keeping with the association between TMB and ICB response (34). These results were reproduced in a second cohort (24) of cutaneous melanoma patients treated with anti–PD-1 monotherapy, demonstrating the predictive power of the InTumor signature in patient samples before start of treatment (Supplementary Fig. S5C–S5F), as well as on treatment (Supplementary Fig. S5G–S5I). No differences in TMB were observed between the InTumor^HI and InTumor^LO groups for either pretreatment data set (Supplementary Fig. S5J), and no link was observed between the InTumor signature and response to anti–CTLA-4 (Supplementary Fig. S5K). Thus, the InTumor signature shows a strong inverse correlation with response to anti–PD-1, but not anti–CTLA-4.

To biologically validate the inverse correlation between the InTumor signature and immune sensitivity, we determined whether this signature predicts T-cell killing in a well-defined in vivo setting that we developed recently (30). We selected two melanoma PDX derived from biopsies from the same patient taken before and after progression on vemurafenib (M009 and M009R; ref. 13; Supplementary Fig. S6A and S6B), which showed a low and high InTumor signature, respectively (Fig. 3G). These PDX were transplanted into mice and subsequently challenged with ACT of TCR-matched cytotoxic T cells (30). This setup allowed us to challenge these InTumor^LO and InTumor^HI tumors under identical TME and immune-infiltrate conditions. Complete T-cell–mediated tumor control was observed for the InTumor^LO PDX (M009; Fig. 3H). In contrast, the InTumor^HI PDX (M009R) treated with matched T cells expanded exponentially and indistinguishably from the control ACT treatment (Fig. 3H). These results were confirmed in four melanoma cell lines with different InTumor signature scores (Supplementary Fig. S6C–S6E). No correlation between response and the infiltration of T cells was observed (Supplementary Fig. S6D and S6E). These results confirm that an InTumor^HI signature is associated with resistance to T-cell killing in vivo.

The establishment of the ExTumor signature predicting anti–CTLA-4 response and the InTumor signature predicting anti–PD-1 response allowed us to ask next whether these signatures, alone, had predictive value for the clinical response to either anti–CTLA-4 or anti–PD-1 in the combination treatment setting. This was determined for stage IV cutaneous melanoma patients (CA209-038 trial, NCT0162149023; refs. 24, 35; Supplementary Table S9). ROC analysis revealed that both individual signatures had predictive value, but underestimated the number of patients responding to combination ICB (Supplementary Fig. S7A–S7D).

Because the ExTumor and InTumor signatures independently underestimated the number of responders to combination ICB, we next determined whether they, when used together, have better predictive power for combination ICB. Indeed, the signatures when used in conjunction correctly predicted outcome for combined ICB for 5/6 patients who had relapsed and 14/15 patients who remained relapse-free (Χ² P = 0.003, AUC = 0.92; Fig. 4A). The sensitivity and specificity of this ICB response quadrant model outperformed previous signatures [IMPRES (7), IFNg (8), and cytolytic score (10); Supplementary Fig. S7E]. Furthermore, the combined signature prediction showed a significant difference in progression-free survival (PFS) between predicted combination ICB responders and nonresponders (log-rank, P = 0.0001; Fig. 4B). Importantly, we confirmed these results in the independent OpACIN trial (ref. 25; NCT02437279) in which stage III melanoma patients were treated also with anti–PD-1 + anti–CTLA-4. Again, the ExTumor and InTumor signatures used in conjunction correctly predicted outcome for combined ICB for 4/6 nonresponders and 11/12 responders (Χ² P = 0.04, AUC = 0.792; Supplementary Fig. S7F and S7G), which was independent of TMB (Supplementary Fig. S7H). Furthermore, the combined signature prediction showed a significant difference in PFS (log-rank, P = 0.0007) between predicted responders and nonresponders (Supplementary Fig. S7I).

Lastly, we combined the InTumor with the ExTumor signatures to determine the percentage of responders to monotherapy [either anti–PD-1 (9, 24) or anti–CTLA-4 (22)] and combination ICB therapy (anti–PD-1 + anti–CTLA-4; refs. 25, 35) by classifying them into so-called response quadrants (5 cohorts, n = 131 patients; Fig. 4C). Both for anti–PD-1 monotherapy and the combination therapy, the classification of patients in these cohorts into the quadrants was significantly associated with clinical benefit (Χ², P < 0.001 and P < 0.05, respectively; Fig. 4C). More importantly, this analysis revealed that patients with an InTumor^HI/ExTumor^HI classification respond significantly better to combination (anti–CTLA-4 + anti–PD-1) treatment (80% of patients had clinical benefit) than to either anti–CTLA-4 or anti–PD-1 alone (36% and 33% respectively, Χ² P < 0.05; Fig. 4D; Supplementary Fig. S8). Furthermore, the patients classified as InTumor^LO/ExTumor^LO did not significantly benefit from anti–CTLA-4 in addition to anti–PD-1 (Χ² P < 0.05; Fig. 4D; Supplementary Fig. S8), demonstrating that this ICB response quadrant model has the potential to perform personalized response prediction and thereby guide treatment decisions.

CTLA-4 and PD-1 blockade have become standard checkpoint inhibitor therapies for various malignancies including melanoma, bladder cancer, and renal cancer. These antibodies are administered either as mono- or combination therapy, the latter of which typically shows more durable clinical benefit (2, 4). However, the toxicity profile of the combination treatment is significantly worse; therefore, clinicians need tools to decide whether PD-1 monotherapy is likely to give sufficient clinical benefit, or whether combination with CTLA-4 blockade is required for improved clinical outcome. Although there are signatures available that predict general response to combination ICB, they cannot inform clinicians on which of the two individual treatment arms will be most effective (7, 8).

Given that the ICB response is likely governed by both tumor cell–intrinsic and –extrinsic signals (6), we set out to develop such individual signatures. We took advantage of an algorithm we developed recently (XenofilteR; ref. 12) and which allows for the computational curation of gene-expression signals from mouse and human origins. This was applied on a melanoma PDX platform we developed previously, for which we had noted that soon after transplantation, human TME cells were swiftly replaced by murine cells. In silico subtraction of the curated human PDX tumor cell RNA reads from the bulk tumor from patients (that is, tumor cell + TME) allowed us to obtain a signature comprising patient nontumor cell signals (that is, TME).

This approach yielded tumor cell–derived “InTumor” and tumor cell–extrinsic “ExTumor” gene-expression profiles. Biologically, we demonstrate that melanomas classified with a high InTumor signature score are refractory to T-cell elimination in mice. Of note, similar to several other types of gene-expression signatures (7, 9, 36), whether the classifiers contain one or more genes that mechanistically contribute to these biological and immunologic traits will need further study. A possible candidate is NGFR, one of the 89 genes in the InTumor signature. We recently showed that it contributes to resistance to multiple therapies including ICB (11). For AXL, also part of the InTumor signature, similar results were demonstrated with high expression in samples with a lack of response to therapy (37, 38). With the validation of these classifiers both in animal models and in several independent patient cohorts, we demonstrate their clinical relevance as critical biomarkers for single and combination ICB treatment response.

With the two signatures we were able to develop a classification for individual melanoma patients using four treatment response quadrants. Each quadrant covers a distinct response to both individual ICB monotherapies and to the ICB combination therapy. Probably the most applicable outcome is the possibility to predict whether patients need anti–CTLA-4 in addition to anti–PD-1 for durable clinical benefit. This is a clear benefit over previously described signatures, which only generate a general prediction, regardless of whether this relates to aCTLA-4, aPD1, or the combination. Our results will need to be validated in additional larger patient cohorts before they may become a decision tool in daily clinical practice. In conclusion, these signatures may enable personalized treatment advice based on the predicted clinical benefit from either anti–CTLA-4 or anti–PD-1 monotherapy or from their combination, or point to the need for alternative checkpoint inhibitor combination therapy.

O. Krijgsman reports grants from Oncode Institute during the conduct of the study, as well as grants from Bristol-Myers Squibb outside the submitted work; in addition, O. Krijgsman has a patent for WO2020005068A8 pending. K. Kemper reports other support from and is currently an employee of GenMab outside the submitted work. M.A. Ligtenberg reports other support from Immagene BV outside the submitted work. P. Ross-Macdonald reports being an employee of Bristol-Myers Squibb at the time of this work. J.B.A.G. Haanen reports grants and personal fees from BioNTech, BMS, Novartis, MSD, and Neogene Tx; grants from Amgen; and personal fees from Immunocore, Merck Serono, Molecular Partners, Pfizer, Sanofi, Roche, PokeAcel, T-Knife, Third Rock Ventures, Achilles Tx, and Iotech outside the submitted work. D.J. Adams reports grants from Netherlands Cancer Institute during the conduct of the study. C.U. Blank reports grants from advisory roles at BMS, MSD, Roche, Novartis, GSK, AstraZeneca, Pfizer, Lilly, GenMab, Pierre Fabre, and Third Rock Ventures, as well as research funding from BMS, Novartis, and NanoString during the conduct of the study; C.U. Blank also reports stock ownership in Uniti Cars and is co-founder of Immagene BV. D.S. Peeper reports to be co-founder, shareholder, and advisor of Immagene, which is unrelated to the current study. No disclosures were reported by the other authors.

O. Krijgsman: Conceptualization, data curation, software, formal analysis, validation, investigation, visualization, writing–original draft. K. Kemper: Conceptualization, validation, investigation, writing–review and editing. J. Boshuizen: Investigation, writing–review and editing. D.W. Vredevoogd: Investigation. E.A. Rozeman: Resources, investigation. S. Ibanez Molero: Investigation. B. de Bruijn: Investigation. P. Cornelissen-Steijger: Investigation. A. Shahrabi: Investigation. M. Del Castillo Velasco-Herrera: Resources, data curation. J.-Y. Song: Investigation. M.A. Ligtenberg: Investigation. R.J.C. Kluin: Software. T. Kuilman: Software, investigation. P. Ross-Macdonald: Data curation, validation, investigation. J.B.A.G. Haanen: Resources. D.J. Adams: Resources, investigation. C.U. Blank: Resources, supervision. D.S. Peeper: Conceptualization, supervision, project administration, writing–review and editing.

We would like to express our gratitude to the patients and their relatives for their cooperation in this study. We acknowledge NKI-AVL Animal Pathology facility for performing the IHC stainings and the NKI Research-HPC facility for computational resources. We would like to thank the Peeper and Blank research groups, Lodewyk Wessels and Karin de Visser for valuable input, Rachid Manumur for bioinformatical analysis, and Steven de Jong for the vimentin IHC protocol. The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007–2013)/ERC synergy grant agreement 319661 COMBATCANCER. This work was financially supported by a grant from the Dutch Cancer Society (NKI-2013-5799; K. Kemper, P. Cornelissen-Steijger and D.S. Peeper), Cancer Research UK and The Wellcome Trust (WT098051, D.J. Adams), and a Queen Wilhelmina Award by the Dutch Cancer Society and by Oncode Institute (D.S. Peeper).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.