Transcriptional Mechanisms of Resistance to Anti–PD-1 Therapy

This article is featured in Highlights of This Issue, p. 2919

Complex interactions between the immune system and melanoma influence tumor invasion and metastasis, and ultimately define patient survival. Immunologic and genetic factors may contribute independently to prognosis and treatment resistance in melanoma; however, a deeper knowledge of their interactions is central to understanding disease evolution and developing synergistic combination treatment regimens. Melanoma has the capacity to metastasize to virtually any organ in the body; often, the full extent of metastasis is discovered only on postmortem examination. To date, studies of melanoma heterogeneity and metastatic evolution have been hampered by the limited availability of tissues from readily accessible anatomic sites in living patients. This limitation can be overcome by a rapid autopsy approach (1). This approach has been applied to pancreatic and other cancer types, revealing diverse patterns of clonal evolution from the primary lesion to metastases, in which intratumoral subclones and metastatic tumors accumulate private mutations superimposed on shared or ubiquitous mutations (2–4). In some cases, an unexpectedly long chronology from primary lesion to the development of metastasis has been observed.

Immune checkpoint blockade with PD-1/PD-L1 antibodies, recently approved for treating metastatic melanoma, Hodgkin lymphoma, and cancers of the lung, kidney, bladder, and head and neck, unleashes endogenous antitumor immunity to eliminate cancer cells (5). To date, most research into anti–PD-1 treatment resistance has focused on immunologic factors such as intratumoral PD-L1 expression and infiltration of cytolytic CD8⁺ T cells (6–8). However, the extent to which genetic and other factors may play equally important roles is unknown. We therefore combined genetic and immunologic analyses of tumor evolution within one patient's cutaneous melanoma over a 5-year clinical course from diagnosis to metastasis and death. This patient received anti–PD-1 therapy up to the time of death, at which time some metastases were responding to therapy while others were progressing. Surprisingly, genetic and immunologic factors were relatively homogenous among responding and progressing lesions, whereas a nongenomic program characterized the progressing metastases.

At autopsy, entire cutaneous metastases, or subsections in the case of very large lesions, were collected from the skin prior to opening the body. The body cavity was subsequently opened with standard autopsy techniques, and each grossly identified metastasis was removed. All tumor and normal tissues were immediately flash-frozen in liquid nitrogen and stored at minus 80°C. For each tumor, one half of the tissue was fixed in 10% buffered formalin, and the remaining tissue was stored at −80°C. Each metastasis was embedded and frozen in Tissue-Tek OCT for sectioning using a Leica Cryostat. For the formalin-fixed, paraffin-embedded (FFPE) samples, slides were cut for tissue macrodissection using a Leica Microtome. Once visualized, each tumor was macrodissected to remove surrounding nonneoplastic tissue. A 5-μm–thick section was used to create an H&E slide for estimating neoplastic cellularity with a standard microscope. We estimated that the neoplastic cellularity was >80% for all samples collected.

The established melanoma cell lines 397-, 537-, 624-, 1011-, 1088-, 1558-, and 2048-mel were generated from metastatic melanoma lesions, each from a unique patient, in our laboratory and maintained as previously described (9); they were characterized by morphology, expression of melanoma antigens, and HLA typing to confirm identity with the donor of origin. They were confirmed to be mycoplasma free. Growing cell cultures were harvested from flasks with trypsin and washed with centrifugation, and cell pellets were snap frozen on dry ice for subsequent RNA isolation.

Genomic DNA (gDNA) was extracted from each tissue piece using phenol and chloroform, followed by ethanol precipitation in ethanol, then quantified by LINE assay [i.e., long interspersed elements (LINE)] using real-time PCR. This method is highly sensitive for calculating the concentration of gDNA. The LINE forward primer used was 5′-AAAGCCGCTCAACTACATGG-3′, and the reverse primer used was 5′-TGCTTTGAATGCGTCCCAGAG-3′. The real-time PCR protocol was 50°C for 2 minutes, 95°C for 2 minutes, 40 cycles of 94°C for 10 seconds, 58°C for 15 seconds, and 70°C for 30 seconds, 95°C for 15 seconds, and 60°C for 30 seconds. All PCR reactions utilized Platinum SYBR Green qPCR master mix (Invitrogen). Only high-quality tissue samples with adequate concentrations were used for further study. For FFPE samples, gDNA was extracted using standard protocols and quantified with a Qbit prior to the generation of sequencing libraries.

Frequently mutated regions of BRAF, KIT, NRAS, and PIK3CA were assessed in a diagnostic tumor biopsy (lesion M26, Table 1) by next-generation sequencing (Ion AmpliSeq Cancer HotSpot Panel v2; Life Technologies).

Exome sequencing was performed on four premortem FFPE tumor specimens, 22 autopsy tumor specimens (FFPE or fresh frozen), and two matched normal samples acquired at autopsy (FFPE and fresh frozen normal kidney) using between 500 ng and 1 μg of gDNA captured by hybridization using the SureSelect XT Human All Exon V4 (Agilent Technologies; Table 1). Samples were prepared according to the manufacturer's instructions. PCR amplification of the libraries was carried out for 6 cycles in the precapture step, and either 8 or 10 to 12 cycles post-capture for frozen and FFPE materials, respectively. Samples were barcoded and run on a Hiseq 2500 in a 75 bp paired-end run using the TruSeq SBS Kit v3 (Illumina). The average depth of coverage in target regions was 180× and 86×, respectively, for tumor and matched normal samples, with an average duplication rate of ∼14% across all specimens, and 96.7% of the targeted regions were sequenced to 30-fold coverage or higher. Raw sequencing reads were aligned and mutations detected as previously described. Briefly, paired-end sequencing data from exome capture libraries were aligned to the reference human genome (hg19) using BWA (10). Read de-duplication, base quality recalibration, and multiple sequence realignment were performed using the Picard Suite and GATK (11, 12). Coordinate sorted BAM files were further processed for quality control, and point mutations and small insertions and deletions were detected with MuTect and HaplotypeCaller, respectively (11, 13). In total, 21,550 candidate somatic mutations were detected, of which 2207 were unique variants across the tumor samples.

Total, allele-specific, and integer DNA copy number genome wide were determined with the FACETS algorithm from exome data in all cases (14). The FACETS framework segments total and allele-specific DNA copy number simultaneously from the coverage and genotypes of polymorphic SNPs genome wide. Allele-specific segmentation is based on the log odds ratio of allele fractions at SNPs identified as heterozygous in the matched normal sample. A fit is applied to the resulting segments, identifying in each sample the (i) log ratio corresponding to diploidy, (ii) purity, and (iii) average ploidy. Major and minor integer copy number is then assigned to each segment by maximum likelihood. The median purity and ploidy across all pre- and postmortem specimens was estimated from sequencing data to be 0.75 and 2.88, respectively. The presence of whole genome duplication was determined from a major integer copy number of 4 for the majority of the autosomal genome not affected by copy number alterations (CNA) with the presence of a large-scale (chromosome/arm-length) heterozygous loss in total copy number that appeared to be in allelic balance with a major and minor copy number of 4 and 2. Finally, the purity and integer copy number results from FACETS analysis, along with coverage levels and allele frequencies, were used to estimate the fraction of cancer cells harboring each mutation (cancer cell fraction, CCF) in all specimens. For mutations in regions of genomic gain, two CCFs were calculated, assuming the minimum and maximum possible number of copies (15).

To identify driver gene mutations, we utilized multiple methods to cross-reference the somatic variants relying on two pan-cancer studies (16, 17) that leveraged ratiometric and statistical approaches to identify driver genes and hotspot mutations. We also analyzed significantly altered genes from three recent sequencing studies of cutaneous melanoma (18–20) as well as anti–PD-1 therapy sequencing studies (21–24). The somatic mutations were queried and manually reviewed for quality and biochemical effects. Sequence data have been deposited at the European Genome-phenome Archive (EGA), which is hosted by the European Bioinformatics Institute (EBI) and the Centre for Genomic Regulation (CRG), under accession number EGAS00001002195. Further information about EGA can be found at https://ega-archive.org and “The European Genome-phenome Archive of human data consented for biomedical research” (http://www.nature.com/ng/journal/v47/n7/full/ng.3312.html).

Evolutionary analysis was performed from a subset of the highest confidence somatic mutations. Phylogeny construction was based on the subset of somatic mutations (n = 795) with greater than 20-fold coverage in every tumor and normal specimen, with a minimum allelic fraction of 10% in at least one tumor specimen, and no more than a single variant read or 2% allelic fraction in both the frozen and FFPE matched normal samples. The phylogenetic tree was rooted using the patient's matched normal sample, and the leaves comprised the primary tumor, local and distant premortem recurrences, and postmortem metastatic samples. We generated a binary classification to describe the presence or absence of each somatic mutation in each sample based on an allelic fraction of >5% in a given sample. We then utilized the discrete character parsimony algorithm to construct the most parsimonious tree of this patient's disease utilizing PHYLIP (25). The resulting phylogeny indicated that some mutations were either lost during evolution or undetected during sequencing, possibly due to sample cellularity and sequencing sensitivity of chromosome-level genetic mutations.

The patient samples were HLA typed and analyzed as described previously (26). Neoantigen analysis discovered 812 unique nonamers across the tumors predicted to have ≤500 nm HLA binding affinity.

Ten posttreatment cutaneous metastases harvested at autopsy were submitted for routine H&E staining as well as immunohistochemistry (IHC). Cases were reviewed by a board-certified dermatopathologist (J.M. Taube) for the presence of post-autopsy autolysis. Percentage tumor necrosis and the presence of lymphoid aggregates were scored for each specimen on H&E. In addition, the architectural pattern of immune cell infiltration was scored as mild, moderate, or severe, as previously described (27). IHC for CD3, CD4, CD8, CD68, CD20, and FoxP3 was performed according to standard automated methods. IHC stains for PD-L1, PD-1, and LAG-3 were performed as previously described (27–29). TIM-3 IHC was performed using clone F38-2E2 (BioLegend) at a concentration of 1.5 μg/mL. Antigen retrieval was performed in citrate buffer, pH 6.0, using a decloaking chamber (Biocare Medical) set to 120°C for 10 minutes. The primary antibody was incubated at 4°C overnight. An anti-mouse IgG1-biotin secondary antibody (BD Biosciences) was used at 1.0 μg/mL. The signal was amplified (Dako K1500) followed by addition of streptavidin-HRP, and visualization with DAB. Membranous PD-L1 expression on tumor or infiltrating immune cells was scored separately as none, <5%, 5% to 10%, and at increasing 10% intervals. The geographic association of PD-L1 expression with the presence or absence of immune infiltrates was also noted. Lymphocyte subset markers (PD-1, CD20, FoxP3, LAG-3, TIM-3) were scored on a 0 to 3 scale, representing the proportion of lymphocytes present expressing the marker of interest (0, none; 1, <5%; 2, 5–50%; 3, >50%; ref. 29). The CD4:CD8 ratio of cells was also determined for each lesion.

Quantitative analysis was performed for CD8⁺ T-cell infiltration. Slides stained for CD8 were scanned using the Aperio ImageScope (Leica Biosystems). The slide images were then annotated and quantitatively analyzed using Halo Image Analysis Software (Indica Labs). After the tumor border was delineated and necrotic regions were excluded, 25 μm zones on each sides of the tumor border were generated. The intratumoral zones were designated as “T” zones, and the extratumoral (i.e., peripheral) zones were designated as “P” zones (Supplementary Fig. S2B), with zones T1 and P1 representing the 25 μm zone closest to the tumor border on the respective sides. The density of CD8⁺ T cells in each zone was quantified using the Immune Cell Module of the Halo software. CD8⁺ T-cell densities were compared between progressing and regressing lesions using the Wilcoxon rank-sum test.

LAMA3 protein expression was assessed by IHC per standard methods using Dako Target Retrieval Solution, a mouse mAb raised against human LAMA3 protein (Atlas Antibodies, catalog number AMAb91123), and a Dako EnVision+ Dual Link Kit for detection of immunolabeled protein. Sections were incubated with a 1:100 dilution of anti-LAMA3 overnight at room temperature. Both the intensity of labeling (0–3) and the percent of melanoma cells labeling (0%–100%) were recorded and used to generate a histoscore (intensity × percent positive cells) that ranged from 0 to 300. Cytokeratin 1 (KRT-1) was detected with the mAb LHK1 (ThermoFisher/Invitrogen, catalog number MA5-16032). Low and high molecular weight cytokeratins were detected with the mAb cocktail AE1/AE3 (Agilent/Dako, catalog number M351529-2).

For FFPE tumor specimens, areas of tumor were identified with H&E staining on neighboring tissue sections, and tissue was manually dissected by scraping with a scalpel for RNA isolation with the High Pure RNA Paraffin Kit (Roche) as previously described (27). Alternatively, areas of tumor that were free from necrotic or normal tissues were isolated with laser capture microdissection as described (27). For cultured melanoma cells, total RNA was isolated with the RNeasy Mini Kit (Qiagen) according to manufacturer's instructions. One hundred fifty nanograms of total RNA was reverse transcribed in a 10 μL reaction volume using qScript cDNA SuperMix (Quanta Biosciences) per protocol. From each RT reaction, 7.5 μL was preamplified in a total volume of 30 μL using a 14-cycle PCR reaction per PreAmp protocol (Applied Biosystems). Fourteen microliters of each preamplification reaction was expanded into a 440 μL total volume reaction mix and added to custom-made TaqMan Low-Density Array Micro Fluidic Cards per protocol (TLDA; Applied Biosystems). One TLDA Card contained 60 candidate immune genes of interest and four endogenous control genes as previously described (28). Other TLDA Cards contained genes derived from whole genome microarray results (Supplementary Table S3). Results from FFPE tissue specimens were analyzed using the 2^−ΔΔCt method (30) with the Welch t test and visualized with TIBCO Spotfire software. Results from melanoma cell lines were visualized by using GraphPad software.

Global gene expression in tumor specimens from melanoma cutaneous metastases [regressing on therapy (R), n = 6; or progressing on anti–PD-1 therapy (P), n = 4] was measured by DASL (cDNA-mediated Annealing, Selection, extension, and Ligation) assays arrayed on the Illumina Human HT-12 WG-DASL V4.0 R2 expression BeadChip, per the manufacturer's specifications (Illumina). This platform detects 29,377 annotated transcripts and is designed to detect partially degraded mRNAs such as typically found in FFPE tissue specimens. Details are provided in Supplementary Methods. For each Illumina probe on the DASL array, the expression levels in R and P lesions were compared. Principal component analysis based on a list of 866 Illumina probes corresponding to 827 transcripts with R versus P FC magnitude ≥1.5 was performed using Partek Genomics Suite. For functional analysis, two lists of genes (i.e., Illumina Probe IDs) for which the FC was ≥1.5 (up in R), or at or below −1.5 (up in P),were examined for significant enrichment in gene ontology (GO) categories and in functional categories including KEGG pathways using the DAVID 6.7 Web tool (http://david.abcc.ncifcrf.gov/; refs. 31, 32). Specifically, for each list of Illumina Probe IDs submitted to DAVID, the DAVID Functional Annotation Tool was run on the resulting Homo sapiens genes recognized in the DAVID database, followed by Functional Annotation Clustering. The Illumina Probe IDs from the submitted list that were in each Term (gene set) in the DAVID output file were then replaced by the corresponding, currently available (as of January 2016) NCBI/HUGO gene symbol.

To generate a heat map, genes derived from the R versus P class comparison using the criteria of Welch t test P ≤ 0.05 and FC magnitude ≥1.7 (calculated using the geometric mean of the expression levels in each group) were visualized using the Stanford Cluster 3.0 and TreeView software as elsewhere described (33). The gene expression data have been deposited in NCBI's Gene Expression Omnibus (34) and will be accessible through GEO Series accession number GSE79691 (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE79691). The GEO platform for this dataset is GPL14951.

Briefly, 140 probes were mapped to genes according to IlluminaHumanWGDASLv4.db (ref. 35; Supplementary Table S3), resulting in 133 unique gene annotations among which 124 were also represented in The Cancer Genome Atlas (TCGA) RNAseq dataset. An additional seven genes were mapped between platforms manually, through updated probe annotations, and/or gene aliases, to reach a total of 131 genes. The same mapping was used for the 13 genes confirmed as differentially expressed by RT-PCR. Logarithmic signal values and Manhattan distance were used for the clustering procedure. Kaplan–Meier survival curves were performed, in R, as elsewhere described, by applying the Peto & Peto modification of the Gehan–Wilcoxon test (36, 37). We also ran the same analysis on GSE78220 (21).

A 60-year-old man presented with a primary cutaneous melanoma of the left upper extremity (7 mm in Breslow thickness). Clinical testing indicated that no mutations were present in the most frequently mutated regions of BRAF, KIT, NRAS, and PIK3CA. Over the ensuing 2 years, the patient underwent surgical resection of locoregional recurrences, including lymph node and in-transit cutaneous metastases. Distant metastases then appeared, for which he received two 4-dose courses of ipilimumab (anti–CTLA-4) over 15 months. The patient's disease was initially stable for several months but then progressed. The patient subsequently began a 6-month course of nivolumab (anti–PD-1), during which he died suddenly of an unrelated cause. A rapid autopsy was performed within 6 hours of death, under consent on an institutional review board (IRB)–approved protocol.

Twenty-six tumor specimens were selected for study, including four premortem specimens (the primary melanoma, biopsies of two locoregional recurrences, and a distant cutaneous metastasis; Table 1). Twenty-two metastases were sampled at autopsy. Postmortem specimens were derived from six organ sites including skin, lymph nodes, stomach, small bowel, lung, and brain. Ten cutaneous metastases were specifically selected, as they showed unequivocally divergent clinical behavior on anti–PD-1 therapy, including six responding and four progressing lesions (Fig. 1), thereby facilitating the study of anti–PD-1 response/resistance in multiple lesions derived from a single organ in an individual patient.

Overall, at the time of autopsy, this patient had a classical high mutational burden cutaneous melanoma. Mutational signature decomposition analysis revealed the expected abundance of C>T/G>A mutations associated with prior ultraviolet light exposure (38). Whole exome sequencing of all 26 specimens confirmed that this tumor was wild type for canonical activating mutations in RAS, RAF, and KIT and was NF1-driven with clonal biallelic inactivation of the gene present in all primary, recurrent, and metastatic specimens collected pre- and postmortem (39, 40). All tumors harbored a truncating NF1 Q519* mutation that was accompanied by a focal heterozygous loss of the single remaining wild-type allele (Fig. 2).

We sought to explore the evolutionary trajectory of this melanoma from diagnosis until the patient's death and throughout immunotherapy with anti–CTLA-4 followed by anti–PD-1. At the genomic level, the tumors did not display a significant number of genetic subclones with respect to cancer cell fraction (Fig. 2A). Indeed, the majority of mutations were either present in all samples or private to one. The homogeneity of the genomic profile observed at the somatic mutational level was also found in DNA CNA analysis. All tumors possessed a high burden, but very stable pattern, of copy number changes, with few changes acquired after diagnosis (Supplementary Fig. S1). Consequently, phylogenetic analysis indicated that although all tumors harbored a small number of private mutations, there was only modest ongoing genomic evolution after the initial branch point defined by 586 shared truncal somatic mutations (Fig. 2B). Review of the somatic mutations that defined the first major branch point did not reveal evidence of additional driver genes or CNAs as the underlying basis of these two distinct subclonal populations.

As most late-arising branching in the phylogeny was supported by few mutations, we can infer only major patterns in this patient. The local recurrence in the skin of the left arm after primary tumor resection (lesion M13) was the closest in evolutionary relationship to the primary tumor (M26), although the local recurrence had evidence of ongoing genomic evolution in the form of additional private mutations. Moreover, the distant cutaneous metastasis resected 4 years after diagnosis, as well as those present postmortem, were highly similar to each other and appear to have emerged through linear clonal evolution from the primary tumor. The cutaneous metastases were relatively distant evolutionarily from most, but not all, visceral metastases, some of which branched off at an evolutionary point prior to diagnosis (Fig. 2B). Notably, no somatic mutations cleanly divided the two groups of cutaneous metastases that were either responding or progressing on anti–PD-1 therapy. Furthermore, none of the mutated genes were differentially expressed at the transcriptional level between the two groups (see gene expression profiling). Collectively, these findings indicate this patient's melanoma underwent a singular branching event prior to diagnosis, and the two major branches each proceeded by linear evolution with nearly all samples acquiring private mutations. Importantly, both the responding and progressing cutaneous metastases were derived from a single subclone, whereas visceral metastases arose from one of two distinct subclones that branched early in the phylogenetic history of this neoplasm. Similar to the mutational data, neoantigen prediction analysis did not indicate any neoepitopes common to the regressing or progressing cutaneous metastases, although many putative neoantigens were detected across the tumor samples (Supplementary Table S1).

The occurrence of 10 metastases in a single organ site (skin) that had clearly dichotomous clinical behavior in response to anti–PD-1 therapy but highly similar genomic profiles led us to explore candidate immunologic factors as the basis for differential response. Tissue integrity in all postmortem specimens was intact without significant histologic evidence of autolysis. We found no significant differences between regressing and progressing cutaneous lesions with regard to the extent of tumor necrosis, presence of lymphoid aggregates, or geographic pattern of immune cell infiltration. We then assessed candidate immune markers with IHC labeling on FFPE specimens. When present, the major PD-1 ligand, PD-L1, was displayed in association with immune cell infiltrates, suggesting “adaptive” expression in response to intratumoral inflammation (27). There were no significant differences in the prevalence of PD-L1 expression on either tumor or infiltrating immune cells, or on immune cell subsets, when comparing cutaneous metastases that were responding or progressing during anti–PD-1 therapy (Supplementary Fig. S2A). Furthermore, although others have correlated the density and geography of tumor infiltrating CD8 T cells with response to anti–PD-1 (7), in a quantitative analysis of these parameters we found no significant differences in CD8 T-cell densities in intratumoral or peritumoral locations in progressing vs. regressing cutaneous metastases (Fig. 3 and Supplementary Fig. S2B). The tumor microenvironments of both the progressing and regressing lesions showed maximum CD8 T-cell densities in the peritumoral stroma, particularly within 100 μm of the tumor border (Fig. 3B and Supplementary Fig. S2B). Similar geographic findings were observed in the primary melanoma and three metastatic lesions biopsied before death, although the density of CD8 cells was much lower in these lesions compared to those collected on anti–PD-1 therapy (Fig. 3C). Moreover, there were no significant differences between responding and progressing lesions among additional immune regulatory molecules and cellular subsets, including LAG-3 and TIM-3 (Supplementary Fig. S2A), and PD-1, CD20, and the CD4:CD8 ratio (not shown).

In melanoma, coordinate gene expression has been observed for PD-L1 (CD274), an immunosuppressive ligand found on tumor and stromal cells in the melanoma microenvironment, and genes whose expression is associated with CD8 T-cell activation (CD8A, IFNG, PRF1, CCL5), antigen presentation (CD163, TLR3, CXCL1, LYZ), and immunosuppression [PD1(PDCD1), LAG3, IL10; ref. 28]. To assess whether these or other candidate immune-related genes were differentially expressed between regressing and progressing cutaneous melanoma metastases, we used the same 60-gene multiplex qRT-PCR array as in our prior melanoma study. Gene expression was normalized to the pan-immune cell marker PTPRC (CD45). We found that none of the screened molecules were significantly differentially expressed according to clinical response characteristics (data not shown), consistent with the IHC results described above for selected immune markers. Thus, neither major genomic events, differences in immune cell infiltration, or candidate immune-related gene expression correlated with the divergent clinical behavior of progressing versus regressing cutaneous metastases on anti–PD-1 therapy.

Given the lack of a genomic or immunologic profile to explain the differential clinical response patterns in the 10 cutaneous metastases on anti–PD-1 therapy, we next turned to transcriptome-wide profiling. RNA was isolated from tumor-enriched areas in six regressing and four progressing lesions by manually scraping unstained tissue sections guided by neighboring sections stained with hematoxylin and eosin (H&E). By comparing the profiles of regressing versus progressing metastases, we identified 827 differentially expressed genes between the two groups, based on fold change (FC) magnitude ≥1.5. Principal component analysis illustrated a segregation of gene expression profiles between regressing and progressing lesions (Supplementary Fig. S3A). Enrichment analyses of the overexpressed genes in responding lesions yielded only four DAVID gene sets (terms), including Signal and Immunoglobulin domain (DAVID Benjamini FDR ≤0.05). In contrast, a similar analysis of genes overexpressed in progressing lesions yielded 42 gene sets involving epithelial, adhesion, and cell differentiation functions (DAVID Benjamini FDR ≤0.010; Supplementary Table S2).

Next, to better differentiate regressing versus progressing metastases, more stringent criteria were applied for gene selection (P ≤ 0.05 and expression FC magnitude ≥1.7), resulting in 140 probes corresponding to 138 distinct genes differentially expressed between the two groups (Supplementary Table S3). Inspection showed upregulation in progressing lesions of keratin genes such as KRT1, KRT5, KRT10, KRT15, KRT78; cellular adhesion genes such as loricrin (LOR), filaggrin (FLG2), desmocollin (DSC1, DSC3), galectin 7 (LGALS7), laminin (LAMA3), and kallikrein (KLK7); and the Wnt pathway genes WNT3 and WNT5A. Among 138 unique genes that were differentially expressed between regressing and progressing cutaneous melanoma metastases after anti–PD-1 therapy, 13 were also found in TCGA “keratin” subclass associated with poor prognosis in patients with melanoma, representing a substantial enrichment: KRT1, KRT5, KRT15, KRT78, KRTDAP, EVPL, KLK7, KLK11, LGALS7, LOR, SCEL, SPINK5, and SPRR1A (18).

Following global gene expression profiling, PCR-based custom arrays were designed to further characterize the differential expression of unique gene targets in regressing versus progressing cutaneous metastases. First, a multiplex array was assembled with 61 genes selected from among the 138 genes that were differentially expressed in microarray analysis (Supplementary Table S3). Criteria used for gene selection included the following: expression FC magnitude ≥2 and P ≤ 0.01, comparing regressing versus progressing tumors, biological associations, and overlap with genes included in the TCGA “keratin” profile (18). Twenty-eight among the 61 genes tested were confirmed to be significantly differentially expressed (FC magnitude ≥2 and P ≤ 0.05) in cutaneous metastases with divergent clinical outcomes (not shown). In particular, this assay confirmed upregulation of genes associated with keratin/epithelial pathways in progressing cutaneous metastases. Among them, KRT1 showed the greatest differential expression, being upregulated >200,000-fold in progressing lesions. To further investigate this finding, IHC for the cytokeratin 1 protein was performed on the 10 cutaneous metastases with clinically divergent behavior. Notably, we failed to observe tumor cell expression of KRT1 protein but instead saw robust expression in normal overlying squamous epithelium and eccrine sweat gland ducts (not shown). Similarly, IHC with an antibody cocktail directed against low and high molecular weight keratins, including the protein products of KRT1, KRT10, and KRT15, failed to detect cytokeratin expression in tumor cells among 21 different melanoma specimens collected from this patient. Retrospective H&E staining of FFPE slides that were manually scraped to isolate RNAs revealed that normal overlying skin had been inadvertently and selectively captured from the progressing skin lesions, generating an artifactual “keratin” tumor gene expression profile.

Next, to insure RNA isolation from pure tumor areas free from normal epithelium or necrosis, laser capture microdissection under direct visualization was performed on tissue sections from the 10 cutaneous metastases (27). Multiplex qRT-PCR arrays were used to assess expression of 122 among 138 genes that were previously found to be differentially expressed in whole genome microarray analysis, excluding noncoding RNAs and those for which commercial probes were unavailable (Supplementary Table S3). Thirteen among 122 genes were confirmed to be differentially expressed in these highly tumor-enriched RNA specimens from progressing versus regressing cutaneous metastases (FC magnitude ≥2 and P ≤ 0.1) when normalized to either 18S or GUSB (Table 2). Genes upregulated in the progressing metastases were associated with cellular adhesion and extracellular matrix formation (LAMA3, CCM2L, CST2, DACT1) and with neutrophil function (FAM183B, PTPRC, CXCR2). Among them, LAMA3 encoding the laminin subunit alpha 3 found in laminins-5, -6, and -7 showed the greatest degree of differential expression, being overexpressed ∼2,000-fold in progressing subcutaneous metastases (P ≤ 0.005). This gene was also well expressed in cultured melanoma cell lines (Fig. 4), supporting tumor cell–intrinsic expression.

IHC was utilized to assess LAMA3 protein expression in the melanoma primary lesion and metastases. As shown in Fig. 5, tumor cells in progressing cutaneous metastases expressed higher levels of LAMA3 protein than did those in regressing metastases. Furthermore, progressing metastases in other organs such as lung, lymph node, and bowel also expressed the LAMA3 protein. Interestingly, LAMA3 was not expressed in the primary melanoma (Fig. 5). Mutations and CNAs in LAMA3 were not observed in the primary or metastatic tumors from this patient, suggesting that LAMA3 expression was acquired during tumor evolution through nongenomic, transcriptional, or epigenetic events.

Melanoma gene expression classification in TCGA, based on mRNA sequencing of treatment-naïve tumors from 281 patients, identified a “keratin” subclass associated with decreased survival, an “immune” subclass associated with prolonged survival, and an intermediate “MITF-low” subclass (18). In our study of multiple melanoma metastases from one patient, a 13-gene profile was identified that differentiated lesional clinical behavior (progressing vs. regressing) after anti–PD-1 therapy. To investigate whether this gene expression profile might correlate with overall survival in the melanoma patient population, we downloaded and analyzed TCGA melanoma gene expression data. An unsupervised cluster analysis based on these 13 genes did not provide any interpretable picture on TCGA expression data, and there was no significant separation of overall survival curves based on this analysis when evaluated by the Wilcoxon test (Supplementary Fig. S4A). This suggests that our 13-gene expression profile, which specifically discriminates melanoma lesional behavior following anti–PD-1 therapy, is not a strong prognostic indicator for patients with melanoma in general. Furthermore, a similar in silico analysis was applied to gene expression data from pretreatment melanoma biopsies in 27 patients who did or did not respond to anti–PD-1 therapy (21). Although there was an apparent association of gene expression with response status that was seen in an unsupervised clustering analysis of these samples, this was not robust to small changes in the gene signature (Supplementary Fig. S4B). Therefore, our 13-gene profile appeared to be specifically relevant to the on-treatment behavior of individual metastatic melanoma lesions during anti–PD-1 therapy.

Melanoma is a disease often characterized by a prolonged “disease-free” interval between diagnosis and clinically apparent metastasis, which can exceed 5 years. Melanoma is also one of the most immunogenic human cancers, and anti-melanoma immune responses are readily detected in vitro (9). Thus, primary melanoma lesions likely metastasize early in the course of tumor evolution, and a dynamic balance between melanoma and the immune system in some patients may durably maintain the tumor in a microscopic state.

The findings in the patient described here are consistent with these ideas. This patient's melanoma was apparently driven by a functional loss of NF1 (neurofibromin), previously shown to induce RAS activation in melanoma cells (39) and to confer resistance to therapy with selective BRAF inhibitors via MAPK pathway reactivation (40). NF1-altered tumors were also recently defined by TCGA as a major genomic class in melanoma (18). Although the clinical history of this patient's melanoma spanned 5 years, and distant metastases were not detected until 3 years after initial diagnosis, relatively little genetic heterogeneity was observed in the multiple metastases sequenced. Instead, the phylogeny of this tumor's evolution revealed a singular early branching event, with each branch then proceeding by linear evolution. Under the selection pressure of anti–PD-1 therapy, several distant cutaneous metastases regressed, whereas others progressed. In a similar scenario with molecularly targeted drugs such as MAPK inhibitors, most studies have focused on genetic alterations as the basis for treatment response/resistance, although more recently, epigenetic, tumor microenvironmental, and immunologic events have also been invoked (41, 42). In this study, we found remarkable genetic similarities among regressing/progressing cutaneous metastases, indicating that all were derived from a single subclone that had not diverged significantly. There were also immunologic similarities among these metastases, and no defining characteristics in candidate immune molecules were associated with treatment response or resistance. However, unbiased global gene expression profiling revealed a distinct signature differentiating progressing from regressing tumors. In particular, the most highly differentially expressed gene, LAMA3, was also found to be overexpressed at the protein level in progressing metastases. LAMA3 is a subunit of laminin-5, which has been associated with epithelial-to-mesenchymal transition and poor prognosis in lung cancer (43), and which promotes the aggressive behavior of melanoma cells in 3-dimensional culture (44). Laminin-5 is a secreted molecule that plays an important role in tissue architecture and adhesion of normal epithelial cells to the basement membrane (45). The normal adhesive role of laminin presents one possible explanation for the selective contamination of progressing tumor specimens with overlying normal epidermis in the manually scraped tissue sections used for RNA isolation in our preliminary experiments. It also raises the possibility that overexpression of laminins might create a barrier excluding immune cells from penetrating the tumor and thus impairing immunotherapy, similar to a recent report of aberrant β-catenin overexpression in melanoma (46). Although we did not observe differential exclusion of CD8 effector cells from progressing versus regressing metastases in the current study, it is possible that dynamic exclusion effects were not captured at the time of tissue procurement, or that other cell types were involved. Interestingly, PTPRC, a pan-leukocyte marker, was overexpressed in progressing metastases. Because we did not detect differences in T and B lymphocyte and macrophage markers with IHC or qRT-PCR in progressing versus regressing cutaneous metastases, this raises the possibility that PTPRC reflects neutrophil infiltration or activation in progressing tumors. This is consistent with the coordinate overexpression of FAM183B (neutrophil-associated acyloxyacyl hydrolase) and CXCR2 (IL8 receptor mediating neutrophil migration) in these specimens. CXCR2-mediated accumulation of tumor-infiltrating neutrophils has been associated with intratumoral immunosuppression in other studies (47).

A recent report of a group of 27 patients with melanoma receiving anti–PD-1 therapy correlated innate treatment resistance with a complex transcriptional signature marked by upregulation of genes regulating mesenchymal transition, cell adhesion, extracellular matrix remodeling, T-cell suppressive inflammation, angiogenesis, and wound healing (21). Our report of a single patient receiving anti–PD-1, which focuses on the posttreatment dichotomy of progressing versus responding metastases, illustrates a more highly focused and potentially adaptive transcriptional profile that characterizes the progressing lesions. Transcriptional tumor profiling offers a powerful dimension in our understanding of response and resistance to cancer therapies. However, it also poses serious challenges stemming from the instability and fragmentation of RNA, especially in formalin-fixed specimens, and the possible contamination of tumor with normal tissues. Especially in melanoma, with many immediate pretreatment and on-treatment biopsies obtained from biopsies of cutaneous and superficial lymph node metastases, contamination with normal epithelium is a concern and indeed contributed to an artifactual “keratin” profile found in our preliminary analysis. This pitfall can be avoided only by tumor microdissection under direct visualization.

New therapies designed to address the complexities of cancer should be based on an understanding of how genetic and immunologic diversity intersect. Although studies of tumor biology in living patients with melanoma are facilitated by easy access to cutaneous and lymph nodal metastases which often occur, they are also limited by the typically small volumes of these disease sites and the inaccessibility of most visceral metastases. This study based on the autopsy of an individual with widespread melanoma offers a window on how cancer evolution can be characterized concordantly at the genetic, transcriptional, and immunologic levels and how we can begin to understand the basis of resistance to anti–PD-1 therapy.

S.L. Topalian reports receiving commercial research grants from Bristol-Myers Squibb, holds ownership interest (including patents) in Five Prime Therapeutics, is a consultant/advisory board member for Five Prime Therapeutics and GlaxoSmithKline, and is related to the coinventor of patents licensed to Bristol-Myers Squibb and Potenza Therapeutics; the family member is also a consultant/advisory board member for MedImmune, Merck, Pfizer, Potenza Therapeutics, and Sanofi. E.J. Lipson reports receiving commercial research grants from Genentech and Merck and is a consultant/advisory board member for Amgen, Bristol-Myers Squibb, and Merck. D.M. Pardoll reports receiving commercial research grants from and holds ownership interest (including patents) in Bristol-Myers Squibb. J.M. Taube reports receiving commercial research grants from Bristol-Myers Squibb and is a consultant/advisory board member for AstraZeneca, Bristol-Myers Squibb, and Merck. No potential conflicts of interest were disclosed by the other authors.

Conception and design: J.M. Taube, R.H. Hruban, D.B. Solit, C.A. Iacobuzio-Donahue, S.L. Topalian

Development of methodology: M.L. Ascierto, J. Fan, G. Kaunitz, T. Cottrell, T.A. Chan, B.S. Taylor, D.B. Solit, S.L. Topalian

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): A. Makohon-Moore, E.J. Lipson, J.M. Taube, T.L. McMiller, T.A. Chan, B.S. Taylor, D.B. Solit, C.A. Iacobuzio-Donahue

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M.L. Ascierto, A. Makohon-Moore, E.J. Lipson, J.M. Taube, A.E. Berger, G. Kaunitz, T. Cottrell, A. Favorov, V. Makarov, N. Riaz, T.A. Chan, L. Cope, B.S. Taylor, D.B. Solit, C.A. Iacobuzio-Donahue, S.L. Topalian

Writing, review, and/or revision of the manuscript: M.L. Ascierto, A. Makohon-Moore, E.J. Lipson, J.M. Taube, T.L. McMiller, A.E. Berger, G. Kaunitz, T. Cottrell, L. Cope, R.H. Hruban, D.M. Pardoll, D.B. Solit, C.A. Iacobuzio-Donahue, S.L. Topalian

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M.L. Ascierto, A. Makohon-Moore, E.J. Lipson, Z. Kohutek, T.A. Chan, D.B. Solit, C.A. Iacobuzio-Donahue

Study supervision: D.B. Solit, C.A. Iacobuzio-Donahue, S.L. Topalian

The authors thank Haiying Xu, Aleksandra Ogurtsova, and George A. Crabill (Johns Hopkins University School of Medicine) and Rajya Kappagantula (Memorial Sloan Kettering Cancer Center) for technical assistance, Steven Hashagen (Indica Labs) for assistance with CD8 T-cell mapping, Annamalai Selvakumar (Memorial Sloan Kettering Cancer Center) for assisting with HLA genotyping, and Ludmila Danilova (Johns Hopkins) for advice on statistical analyses.

This study was supported by the Melanoma Research Alliance (to C.A. Iacobuzio-Donahue and S.L. Topalian), the Bloomberg∼Kimmel Institute for Cancer Immunotherapy (to J.M. Taube, D.M. Pardoll, and S.L. Topalian), the Barney Family Foundation (to S.L. Topalian), Moving for Melanoma of Delaware (to S.L. Topalian), the Laverna Hahn Charitable Trust (to S.L. Topalian), the National Cancer Institute (R01 CA142779 to J.M. Taube, D.M. Pardoll, and S.L. Topalian; 5T32 CA193145 to T.R. Cottrell; R01 CA179991 to C.A. Iacobuzio-Donahue; P30-CA008748 to Memorial Sloan Kettering Cancer Center), the MSKCC TROT Fellowship (to A. Makohon-Moore), the Starr Cancer Consortium (to T.A. Chan), the Pershing Square Sohn Cancer Research Foundation (to T.A. Chan), the Immunogenomics and Precision Oncology Platform (to T.A. Chan), and the Marie-Josée and Henry R. Kravis Center for Molecular Oncology (to D.B. Solit).

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.