Transcriptional downregulation of MHC class I and melanoma de- differentiation in resistance to PD-1 inhibition

Transcriptomic signatures designed to predict melanoma patient responses to PD-1 blockade have been reported but rarely validated. We now show that intra-patient heterogeneity of tumor responses to PD-1 inhibition limit the predictive performance of these signatures. We reasoned that resistance mechanisms will reflect the tumor microenvironment, and thus we examined PD-1 inhibitor resistance relative to T-cell activity in 94 melanoma tumors collected at baseline and at time of PD-1 inhibitor progression. Tumors were analyzed using RNA sequencing and flow cytometry, and validated functionally. These analyses confirm that major histocompatibility complex (MHC) class I downregulation is a hallmark of resistance to PD-1 inhibitors and is associated with the MITFlow/AXLhigh de-differentiated phenotype and cancer-associated fibroblast signatures. We demonstrate that TGFß drives the treatment resistant phenotype (MITFlow/AXLhigh) and contributes to MHC class I downregulation in melanoma. Combinations of anti-PD-1 with drugs that target the TGFß signaling pathway and/or which reverse melanoma de-differentiation may be effective future therapeutic strategies. A significant proportion of patients develop innate or acquired resistance to immune checkpoint inhibitors. Here, the authors show that resistance to anti-PD-1 blockade is associated with TGF-beta driven major histocompatibility complex I (MHCI) down-regulation and a de-differentiated phenotype in melanoma patients.

I mmunotherapy has transformed the treatment of melanoma patients, with monoclonal antibodies blocking the immune regulator programmed cell death protein 1 (PD-1; pembrolizumab and nivolumab) receiving rapid FDA approval following demonstration of a survival benefit compared with chemotherapy and the cytotoxic T-lymphocyte-associated protein 4 (CTLA-4) inhibitor, ipilimumab [1][2][3] . Despite durability of response, innate resistance to PD-1 inhibition occurs in 30% of melanoma patients [4][5][6] and approximately 25% of responding patients will develop acquired resistance, defined as disease progression following initial objective response, within two years of PD-1 inhibitor treatment 4 .
The mechanisms responsible for failure of PD-1 inhibition are diverse and incompletely understood, although resistance effectors have been identified in a small subset of patients. These include the expression of alternate immune checkpoint inhibitors (TIM-3 and LAG-3) 7,8 , loss of major histocompatibility complex (MHC) class I expression, abnormalities in the interferon-γ (IFNγ) immune effector signaling pathway (JAK1, JAK2, IFNGR1, STAT1 mutations) [9][10][11][12][13] , oncogenic signaling (elevated ß-catenin/WNT) that leads to immune exclusion 14 , T-cell induced secretion of immunosuppressive colony-stimulating factor 1 15 and an hypoxic tumor micro-environment that may impair T-cell function 16 . Furthermore, several immune and gene-expression signatures predictive of PD-1 inhibitor response have been reported, but few have been validated in independent patient cohorts 11,[17][18][19] . For example, the innate PD-1 inhibitor resistance (IPRES) signature, which includes 26 gene signatures associated with dedifferentiation and BRAF/MEK inhibitor resistance, was associated with lack of PD-1 inhibitor response in pre-treatment melanoma biopsies in one study 17 , but was not associated with PD-1 inhibitor response in other melanoma cohorts 11,19 .
In this study, we perform transcriptome and flow cytometric analysis on 94 longitudinal melanoma biopsies in a large cohort of melanoma patients receiving PD-1 inhibitors. Analysis of pretreatment and on-treatment tumors, including those responding to therapy (RES) and those that progressed (PROG) due to innate or acquired resistance. We provide insights into the complex and heterogeneous response of individual metastases to PD-1 inhibition and the heterogeneous immune transcriptome profile observed in synchronous and longitudinal biopsies. In addition, we demonstrate that down-regulation of MHC class I expression, rather than complete loss of MHC class I molecules, is common in melanoma and potently driven by TGFß signaling and dedifferentiation.

Results
Patient and tumor characteristics. Transcriptome analysis was performed on RNA sequence data (n = 79 tumors; 55 patients) and flow cytometric analysis on single cell suspensions (n = 31; 24 patients) from a total of 94 melanoma tumors derived from 68 patients treated with PD-1 inhibitor monotherapy (Supplementary Data 1). In total 53 tumor biopsies were pre-treatment (PRE) and 41 were taken while on-treatment (Fig. 1A). PRE tumors from patients who subsequently underwent complete response (CR) or partial response (PR) by irRC criteria 20 were termed responding-PREs (n = 31), and PRE biopsies from patients who had stable disease (SD) or progressive disease (PD) by irRC (n = 22) were termed non-responding-PREs. Of the 41 on-treatment tumor specimens, six biopsies were taken from clinically responding lesions (RES) and 35 were biopsies taken at time of progression (PROG). All PROG lesions were characterized as either innate PROG tumors (n = 22) or acquired PROG tumors (n = 13). Innate PROGs were defined as pre-existing metastases that did not undergo tumor shrinkage or new metastases identified within 6 months of starting treatment, and acquired PROGs were defined as pre-existing tumors that initially responded but subsequently progressed on PD-1 inhibitor or new metastases identified after 6 months of starting treatment (Fig. 1A).
Heterogeneous tumor responses and predicting anti-PD-1 response. We initially examined the predictive accuracy of seven transcriptome signatures associated with clinical response to PD-1 inhibition 11,15,[17][18][19]21,22 in the 44 PRE tumors with available RNA sequence data (Supplementary Data 1). None of the published immune-predictive signatures, including signatures indicative of inflammation (i.e., CD8+ T-cell, CYT score and the 18immune gene set) accurately defined responding (CR/PR) or nonresponding patients (SD/PD) in our cohort (Fig. 1C, Supplementary Data 2). Additionally, we found that none of these predictive transcriptomic signatures were consistently and significantly associated with irRC response in three separate immune-checkpoint inhibitor treated melanoma patient RNAseq datasets (Supplementary Figure 1). Further, we did not detect any significant differentially expressed genes (FDR-adjusted pvalues <0.05) or gene signatures (ssGSEA score differences between two groups; FDR-adjusted p-values <0.05) between responders and non-responders in our cohort.
We hypothesized that accurately predicting melanoma patient response to PD-1 inhibitor therapy based solely on the characteristics of a single pre-treatment biopsy may be confounded by intra-patient heterogeneity 23 . An examination of lesion-specific responses to PD-1 inhibition revealed that 16/68 (24%) patients had heterogeneous tumor responses (Supplementary Data 1). Of these 16 patients, five had PRE core biopsies, allowing lesion-specific assessment of the response. Interestingly, four patients who had irRC PD underwent CR, PR, PR and SD in the lesion biopsied at baseline, and the one patient who had irRC PR underwent PD in the lesion biopsied ( Table 2, Fig. 1D/E). The predictive accuracy of the seven anti-PD-1 predictive transcriptome signatures did not improve, however, when these lesionspecific responses were included in the ROC analyses (Supplementary Fig. 2A).
Transcriptome evolution during PD-1 inhibitor progression. We extended the transcriptome analysis to all patients (n = 55) and tumors (n = 79) with RNA sequencing data, including 44 PRE, 6 RES and 29 PROG melanomas (11 acquired, 18 innate). We initially examined the intra-tumoral cytolytic activity (CYT) score of each tumor, a quantitative measure of immune cytolytic activity based on transcript levels of two cytolytic effectors, granzyme A (GZMA) and perforin (PRF1) 24 . As expected, the CYT score correlated with computational estimates of immune cell fractions and immune activation signatures ( Supplementary  Fig. 2B, 2C). The CYT score was highest, but not exclusively elevated, in the RES biopsies ( Fig. 2A). Fifteen of 44 (34%) PRE and 15/29 (52%) PROG tumors (6/11 acquired and 9/18 innate) showed evidence of an active inflamed transcriptome (i.e CYT score within the CYT score range of the 6 RES tumors) ( Fig. 2A).
The patterns and evolution of tumor inflammation were explored in eight patients who had paired PRE (Fig. 2B). Reanalyses of a separate series of longitudinal melanoma biopsies collected pre-treatment, early on therapy and on progression in patients undergoing treatment with sequential CTLA-4 and PD-1 inhibition confirmed intra-patient heterogeneity of tumor inflammation ( Supplementary Fig. 2D) 19 .
HLA-ABC downregulation is associated with dedifferentiation. We reasoned that mechanisms of PD-1 inhibitor resistance reflect the degree of immune cell infiltration and thus, we explored melanoma transcriptome signaling relative to the level of T-cell activity as defined by the CYT score 24 . As expected, the CYT score was positively correlated with the expression of the MHC class I genes in the 79 melanoma biopsies with RNA sequence data (Pearson correlations of 0.61, 0.77, 0.74 and 0.74 for HLA-A, -B, -C and B2M). However, the correlation of HLA-A transcript expression with the CYT score was diminished in PD-1 progressing biopsies compared to the PRE and RES tumors (Fig. 3A, Supplementary Fig. 3A). We did not identify any expressed alterations in the B2M, HLA-A, -B or -C genes in the transcriptome data from the 79 melanoma biopsies, and the mutations identified in other genes associated with antigen presentation and IFNγ signaling (i.e., JAK1, JAK2, STAT1, STAT2, IFNGR1, IFNGR2, TAP1, TAP2) were not exclusively identified in poor responders or PROG biopsies (Supplementary Data 3). Thus, in our patient cohort, T-cell activity as measured by the CYT score was not strongly associated with HLA-A transcript expression in PD-1 progressing tumors. In order to examine mechanisms that may influence the expression of HLA genes independently of CYT score, we stratified tumors according to CYT score, and identified 19 CYT The downregulation of HLA-A in these CYT-score matched tumors was not associated with diminished CD8+ T-cell content ( Supplementary Fig. 4A), or alterations in the frequency of other leukocyte subsets (based on CIBERSORT profiling; Supplementary Fig. 4B and Supplementary Data 6). HLA-A downregulation was also not related to IFN-γ gene sets and although we detected downregulation of B2M in a number of HLA-A low tumors, this did not reach statistical significance ( Supplementary Fig. 4C). It is also worth noting that PRE tumors with downregulated HLA-A transcript expression were not enriched in patients who had received MAPK inhibitor therapy (Supplementary Data 1) prior to anti PD-1 treatment (Fisher's exact test, p > 0.99). Interestingly, in the PROG tumor derived from patient 53 with a pre-existing, dysfunctional STAT1 S316L mutation 27 ( Supplementary Fig. 5), low transcript expression of HLA-A occurred in the absence of   (46) AJCC, American Joint Committee on Cancer; LDH, lactate dehydrogenase; ULN, upper limit of normal; CR, complete response; PR, partial response; SD, stable disease; PD, progressive disease. a One patient had both a BRAF (G469E) and an NRAS (Q61L/R/P) mutation b Including BRAF non-V600 mutations c Patients were stratified into response groups based on immune-related response criteria. Patients with CR or PR were classified as responders, while patients with SD and PD were classified as non-responders elevated SNAI1 expression (Fig. 3). The second PROG tumor derived from patient 53 (also with the STAT1 316L mutation, Supplementary

HLA-ABC downregulation is common in melanoma biopsies.
To explore the frequency and contribution of HLA-A downregulation to immunotherapy resistance we analyzed 31 tumor dissociates derived from PD-1 PRE (n = 15, including 6 with RNA sequence data; Supplementary Data 1) and PD-1 PROG tumors (n = 16, including 10 from with RNA sequence data). Using flow cytometry, melanoma HLA-ABC expression was categorized as either normal or downregulated (ratio < 0.65) relative to the matching tumor infiltrating lymphocytes from the same biopsy sample. The immune cells served as an internal control that standardized the modulatory effects within the microenvironment. A representative gating strategy is shown in Supplementary Fig. 6. We observed varying degrees of cell surface HLA-ABC downregulation in 11/31 (35%) melanoma tumors (PRE 6/15 (40%) and PROG 5/16 (31%) samples) (Fig. 4A, B). Eight of these eleven tumors had an activating BRAF or NRAS mutation, and although BRAF V600E has been associated with the internalization of HLA-ABC from the cell surface 28 , we did not detect any genotype-associated differences in the cell surface expression of HLA-ABC in 31 melanoma tumors ( Supplementary  Fig. 7A). It is worth noting that 16 melanoma PRE and PROG biopsies had matching flow cytometry and RNA sequence data and the HLA-ABC cell surface expression and HLA-A transcript expression were concordant in these samples (Spearman correlation 0.67, p < 0.01; Supplementary Fig. 7B). We also examined the cell surface expression of HLA-ABC and HLA-DR, which was recently shown to correlate with response to PD-1 inhibition 29 , in the 15 PRE tumors. Although the tumor numbers were small, the cell surface expression of HLA-ABC or HLA-DR at PRE did not accurately reflect patient response ( Supplementary Fig. 7C & D). The importance of HLA-ABC downregulation was confirmed using a short-term co-culture model of melanoma cells and autologous tumor infiltrating lymphocytes. Diminished expression of HLA-ABC (~80% reduction in expression; Fig. 4C) was driven by two independent B2M-targeting silencing molecules and this significantly reduced the immune recognition of melanoma cells, leading to approximately a 70% reduction in the levels of IFNγ production (Fig. 4C). As a comparison, pretreatment of melanoma cells with an HLA-ABC blocking antibody reduced IFNγ secretion by immune cells by over 95% (Fig. 4D).

Discussion
This study confirms that MHC class I downregulation associated with the de-differentiation phenotype is a hallmark of both innate and acquired resistance to PD-1 inhibitors. Despite recent advances in the treatment of melanoma with immune checkpoint inhibitors, innate and acquired resistance remains a major challenge, with only modest improvements observed with second-line salvage therapies 35,36 . Genetic alterations in antigen presentation (B2M, HLA-A) and IFNγ signaling (IFNGR, STAT1 and JAK1/2) in non-responding melanoma patients 9,11-13 are uncommon, and  Data are means ± s.d. and individual data points represent the average of technical triplicates. Data were compared using one-way ANOVA with the Geisser-Greenhouse correction, *p < 0.05. d IFNγ production 72 h after co-culture of WMD-084 melanoma cells with the patient-matched TILs expanded from the same tumor biopsy. Lymphocytes were pre-treated for 1 h with 10 µg/ml HLA-ABC blocking or an isotype-matched antibody prior to co-culturing. Isotype antibody-treated control cells were compared to HLA-ABC blocking antibody-treated data using a paired t-test, *p < 0.05.  [39][40][41] and is associated with intrinsic resistance to PD-1 inhibitor monotherapy in melanoma patients 17 . AXL high melanoma are enriched during BRAF/ MEK and PD-1 inhibitor therapy 41 and are induced via microenvironmental cues including T-cell-induced inflammatory stimuli (e.g., TNFα) 42,43 and cancer-associated fibroblast activity (e.g., TGFß signaling) 41,44,45 . The MAPK inhibitor and immune resistance effectors in de-differentiated melanoma may be distinct, however. In particular, the upregulation of multiple receptor tyrosine kinases, including AXL, PDGFRß and EGFR drive MAPK inhibitor resistance 39 whereas the global downregulation of melanocytic antigens 42,43 and the diminished upregulation of MHC class I expression (this report) mediate immune evasion. Consistent with our data, AXL expression diminished STAT1 phosphorylation and MHC class I expression in a mouse mammary tumor model 46 and induced mesenchymal transition via TGFß/SNAI1-signaling, leading to the downregulation of MHC class I expression in prostate cancer 47 . Although, tumor cells lacking normal expression of MHC class I molecules should activate and be cleared by natural killer cells, the activation and function of these immune cells are likely impaired by the presence of stromal TGFß 48,49 .
Recent data have also confirmed that fibroblast-derived TGFß restrained anti-urothelial cancer cell immunity and PD-1 response by restricting T-cell movement within the microenvironment 25 and elevated TGFβ1 diminished the positive effect of T-cell infiltration on melanoma patient survival outcomes 22 . Further, PD-1/PD-L1 blockade invoked limited response in preclinical colon cancer and mammary carcinoma mouse models with TGFß-activated stroma 25,50 . The combined inhibition of PD-1/PD-L1 with TGFß receptor kinase inhibitors has also been shown to enhance tumor regression in preclinical tumor models 51,52 and a bifunctional fusion protein targeting PD-L1 and TGFß in solid tumors showed encouraging efficacy in a recent Phase I trial 53 . Collectively, these data confirm that environmental TGFß cues contribute to immune evasion and PD-1 inhibitor escape by limiting T-cell infiltration and downregulating MHC class I expression.
In our large patient cohort of carefully annotated PROG samples, biopsied at time of progression, from a confirmed progressing metastases, we found that previously described alterations in antigen presentation (B2M, HLA-A loss) and IFNγ signaling (IFNGR, STAT1 and JAK1/2 alterations) were rare mechanisms of treatment failure 9,11-13 . None of our PROG tumors expressed B2M, IFNGR1, JAK1 or JAK2 mutations, and although we identified 5 PROG tumors, derived from three patents, with STAT1 mutations (patient 25-STAT1 S735F , patient 40-STAT1 S316L and patient 50-STAT1 P633L ), only patient 40 displayed concurrent STAT1 mutation and HLA-A downregulation.
Finally, we were unable to define a baseline transcriptome signature that accurately classified clinical response to PD-1 inhibitors nor could we confirm the predictive value of previously reported gene expression signatures 11,15,17,18 . Several recent reports have also been unable to validate published predictors of immunotherapy response. For instance, the IPRES gene expression signature did not accurately predict PD-1 inhibitor response in several independent melanoma cohorts 11,19,21 . CIBERSORT, cytolytic activity and IFN-γ gene expression signatures were poor predictors of PD-1 and CTLA-4 inhibitor response across several publicly available patient datasets 21,54 and the validity and reproducibility of the IMPRES gene expression signature in predicting immune checkpoint inhibitor response remains contentious 55,56 .
Although baseline predictive information is eagerly sought to guide the selection of first-line immune therapies, our findings suggest that a single-point, robust pre-treatment biomarker is unlikely to represent the heterogeneous nature of cancer, or predict the rapidly evolving profile of tumors under the selective pressure of immunotherapies. The longitudinal analyses of multiple pre-treatment biopsies derived from various sites will be required to validate the influence of tumor heterogeneity on clinical outcomes and on the accuracy of baseline predictive signatures. In this study we observed significant variation in the response of individual lesions to immune checkpoint inhibition within melanoma patients, with heterogeneity of tumor responses observed in 85% and 69% of melanoma patients who did not achieve an objective response to pembrolizumab monotherapy or combination ipilimumab plus nivolumab, respectively 23,57 . This heterogeneity presumably reflects genetic, molecular and cellular variables, which may include the presence of neoantigens in only a subset of tumor clones, the expression of PD-L1 in a small proportion of tumor cells, variable expression of PD-1 on immune cells and intra-tumoral differences in T-cell density and clonality [58][59][60] . It will be interesting to explore whether these limitations influence the predictive value of tumor mutation load, which has been shown to predict clinical benefit of immune checkpoint blockade in multiple cancers, including melanoma [61][62][63] . Recent data suggest that tumor mutation burden may not change significantly during treatment 64 , although intra-patient mutation burden heterogeneity has been observed 65 . Further, when tumor samples were stratified according to melanoma subtype (i.e., cutaneous, occult, acral or mucosal), tumor mutation burden did not predict response to anti-PD-1 based immune checkpoint therapies 66 . Unfortunately, germline variant data was not available in this study and it was not possible to accurately estimate tumor mutation burden in this study, We conclude that MHC class I downregulation associated with the de-differentiation phenotype is a common mechanism of resistance to PD-1 inhibitors. With the availability of synchronous and longitudinally collected biopsy samples from melanoma patients, we were able to address the issue of immune and intrapatient heterogeneity as a significant limiting factor in identifying predictive signatures and in devising strategies to overcome immune checkpoint inhibitor resistance. Nevertheless, combination immunotherapy strategies such as the addition of CTLA-4 inhibitors that utilize additional anti-immunity pathways without solely depending on T-cell mediated anti-tumor immunity, and mechanisms that restore antigen presentation by inhibiting TGFß signaling, may improve the outcomes of patients with metastatic melanoma who are progressing on PD-1 inhibitors.

Methods
Patient, response assessment and tumor biopsies. This study included 68 metastatic melanoma patients who were treated with PD-1 inhibitors (pembrolizumab or nivolumab) at Melanoma Institute Australia (MIA) and affiliated hospitals. Written consent was obtained from all patients (Human Research ethics committee protocols from Royal Prince Alfred Hospital; Protocol X15-0454 & HREC/11/RPAH/444). Tumor response was assessed using the immune-related response criteria (irRC 20 ), with heterogeneous response defined as the presence of progressing and/or new metastases in conjunction with at least one responding metastasis on first restaging imaging. Clinicopathological characteristics including American Joint Committee on Cancer (AJCC) stage, lactate dehydrogenase (LDH) and mutation status were recorded (Supplementary Data 1), and follow-up duration was calculated from the date of first dose of systemic therapy to the following three dates: date of death, loss to follow-up or 30 th November 2018.
RNA isolation and sequencing. Total RNA was isolated from 79 fresh frozen tissue sections using the AllPrep DNA/RNA/miRNA Universal kit (Qiagen, Hilden, Germany) 16 . cDNA synthesis and library construction were performed using the TruSeq RNA Library Prep Kit (Illumina) and paired-end 100 bp sequencing, with each sample yielding 40-50 million read. Sequencing was performed on the Illumina Hiseq 2500 platforms at the Australian Genome Research Facility in Melbourne.
RNA-sequence data processing. Trimming of Illumina TruSeq paired-end sequencing adapter sequences and bases with a quality of <20 from each end was done using cutadapt 67 . Reads less than 50 bases long after trimming were discarded from subsequent analysis. The filtered reads were mapped to reference genome hg38 using STAR with 10% mismatches allowed. Reads that mapped equally to more than one genomic location were discarded. Reads were imported into R with GenomicAligments read GAlignmentPairs function. strandMode was set to 2. GENCODE Genes version 26 was used as the gene reference database. Counting of reads overlapping with exonic regions of each gene was done with the countOverlaps function from GenomicAligments.
Differentially expressed gene analysis. RNA count data were normalized using trimmed mean of M-values (TMM) and transformed with voom to log2-counts per million with associated precision weights 68 . To identify differentially expressed genes associated with variable levels of HLA-A transcript expression we compared all PRE and PROG tumor pairs for CYT score and selected tumor pairs with similar CYT scores (ratio between 0.9-1.1). These were ranked according to HLA-A transcript differences and tumor pairs with HLA-A expression differences greater than two-fold were selected. Nineteen tumors pairs were categorized as low HLA-A transcript relative to CYT score and high HLA-A transcript relative to CYT score. Differentially expressed genes between these two groups were determined using the moderated t-tests (implemented in R package version 3.6.0 limma version 3.40.2) based on an empirical Bayesian approach to estimate gene expression changes 69 . Similarly, we grouped patients into responders (complete and partial response) and non-responders (stable and progressive disease) and also applied moderated t-test (implemented in limma) to identify gene expression associated with response.
Gene set and cell type enrichment analysis. Rank ordering of TMM-voom transformed gene expression data was carried out using the linear model for microarray module (limma package in R/Bioconductor) 69 and analyses was performed using gene set enrichment analysis in pre-ranked mode provided by GenePattern 70 . The Hallmark gene sets 71 of the Molecular Signature Database version 6.2 72 , with stromal and IPRES cell gene signatures 17,25,26 were considered. A false discovery (FDR) corrected p-value <0.05 was used for comparisons between CYT-matched melanoma tumors.
To obtain abundance values corrected for transcript lengths as required by the single-sample gene set enrichment analysis (ssGSEA 73 ), RSEM was used to derive the FPKM estimates using GENCODE Genes version 26 as the reference transcript database. Absolute signature enrichment scores were determined using the ssGSEA (version 9.1.1) implementation provided by GenePattern 70 with the gene sets described above. Subsequently, differential expression analyses on ssGSEA enrichment scores was performed using the moderated t-test (limma package in R/ Bioconductor) 69 . A false discovery (FDR) adjusted p-value <0.05 was used for comparisons between tumor groups. The same FPKM values were also used to infer the relative proportions of 22 types of infiltrating immune cells using the CIBERSORT web portal (http://cibersort.stanford.edu/). These 22 cell subsets were further grouped into 11 major leukocyte subtypes 74 and a moderated t-test (implemented in limma) was applied to identify cell subsets associated with HLA-A transcript expression.
The correlation between transcript expression and ssGSEA enrichment scores was calculated using the Spearman's Rank correlation coefficient in the nearest neighbor algorithm within the Morpheus web based tool (https://software. broadinstitute.org/morpheus/).
The PD-1 predictive signatures applied to our dataset were as follows: IPRES signature (average Z scores of ssGSEA scores of the gene sets in the IPRES signatures 17 ), IMPRES signature 21 , CD8A/CSF1R ratio (numeric difference between TMM-voom transformed CD8A and CSF1R expression data 15 ), 18immune gene signature (ssGSEA score of 18-genes in the immune signature 18 ) TIDE (calculated using TMM-voom sequencing counts normalized for each gene by subtracting the average gene values among all samples at http://tide.dfci.harvard. edu 22 , CYT score (average of TMM-voom transformed GZMA and PRF1 expression data 24 ) and CIBERSORT estimated relative proportion of CD8+ T cells using FPKM expression estimates (http://cibersort.stanford.edu 74 ). These predictive scores were derived for each pre-treatment melanoma biopsy (n = 44) and used with patient response data (complete (CR) and partial response (PR) versus stable (SD) and progressive disease (PD)) to generate receiver operator characteristic (ROC) curves in order to measure the performance of each indicated signature in predicting PD-1 inhibitor responses in our patient cohort. The performance of these seven PD-1 predictive signatures in predicting responding (irRC: CR and PR) were also evaluated in three publicly-available pre-treatment melanoma RNA-seq datasets with response data: (1) 49 patients treated with the PD-1 inhibitor, nivolumab 11 , 26 patients treated PD-1 inhibitors nivolumab or pembrolizumab 17 and 41 patients treated with anti-CTLA4 75 . The area under the ROC curve was calculated using GraphPad version 8.2.1 using non-parametric estimates and 95% confidence interval based on the hybrid method of Wilson and Brown 76 .
Single nucleotide variant (SNV) analysis. SNVs were called against the reference genome using VarScan2. Minimum variant frequency was set to 20% and other parameters were left at their default values. Briefly, the SAMtools mpileup utility provided a summary of the read coverage, and the mpileup output was processed using VarScan2 to call variants and produce a VCF format file with variants that passed the minimum read and allele frequency thresholds. Insertion and deletion calls were not included due to positional ambiguity and low alignment accuracy 77 . Visualization of the resulting VCF files and analysis was performed through the use of Ingenuity Variant Analysis software (https://www.qiagenbioinformatics.com/ products/ingenuity-variant-analysis) from Qiagen.
Tissue processing and cell isolation. Tumor biopsies were manually minced and enzymatically processed, then dissociated into single-cell suspensions using the human Tumor Dissociation Kit and gentleMACS Dissociator (Miltenyi Biotec), according to the manufacturer's instructions.
Co-culture melanoma:TIL assays. IFNγ production in melanoma:TIL co-cultures was analyzed by ELISA or flow cytometry. To measure IFNγ release, 1 × 10 4 melanoma cells were cultured with 1 × 10 4 TILs (1:1 effector to tumor ratio) in a 96-well plate in a total volume of 100 µl TIL medium, and each experimental setup was performed in triplicate. After two days culture, supernatant was collected, spun down to remove cell debris, and stored at −20°C for IFNγ analysis using the Human IFNγ DuoSet ELISA (R&D Systems). ELISA was performed according to manufacturer's protocol. Alternatively, frequency of IFNγ-producing CD8+ T cells was measured by flow cytometry, For flow cytometry, 1×10 5 melanoma cells were cultured with 1×10 5 autologous TILs in a 24-well plate in 0.5 ml of TIL medium. Four hours post co-culture, 5 µg/ml of Brefeldin A and 5 µg/ml monensin (both from Sigma) were added to stop cytokine release and the co-cultures were incubated overnight prior to staining (see below). For TGFβ treatment, melanoma cells were pre-treated with 10 ng/ml TGFβ or vehicle control for 72 h prior to co-culture with TILs.
For B2M silencing experiments, melanoma cells were first transduced with negative control shRNA or B2M shRNA before co-culturing with autologous TILs, as described above. HLA-ABC blocking was performed using a monoclonal mouse anti-HLA-ABC antibody (clone W6/32, Cat No. 311409, Biolegend). Melanoma cells were pre-treated with 20 µg/ml anti-HLA-ABC antibody or mouse IgG2a isotype control antibody (Biolegend) for 1 before co-culture with autologous TILs.
Constructs and lentivirus transductions. The B2M shRNA constructs correspond to nucleotides 144-162 and 402-420 (Genebank accession number NM_004048.2) 78 . The non-silencing negative control shRNA did not show complete homology to any known human transcript and had the following sequence: 5-TTAGAGGCGAGCAAGACTA-3. The shRNA were cloned into pSIH-H1-puro (System Biosciences) lentiviral vector. Lentiviruses were produced in HEK293T cells as described previously 79 . Cells were infected using a multiplicity of infection of 5 to provide an efficiency of infection above 90%. Cells were used 72 h post transduction.
Immunoblotting. Total cellular proteins were extracted at 4°C using RIPA lysis buffer containing protease inhibitors and phosphatase inhibitors (Roche). Proteins Where indicated, membranes were incubated with REVERT 700 total protein stain (LI-COR, Lincoln, NE), imaged using Odyssey CLx imaging system, washed and blocked using LI-COR Odyssey blocking buffer.
All samples were acquired on BD LSRFortessa X20 flow cytometer (BD Biosciences) and analyzed using FlowJo software v10.4 or later (TreeStar, Ashland, OR). At least 10,000 live events were acquired, while all events were collected for the dissociated tumor samples. Relative marker expression levels were calculated by dividing the median fluorescence intensity (MFI) of the antibody-stained sample by the unstained control, unless otherwise specified. For dissociated tumor samples, melanoma HLA-ABC expression was calculated relative to TILs (geometric mean fluorescence intensity, MFI HLA-ABC melanoma/MFI HLA-ABC TILs).
Statistical analyses. For statistical analysis, we used GraphPad Prism software v8.1.1. Figure legends specify the statistical analysis used and define error bars.

Data availability
TCGA gene expression data was downloaded from the The National Cancer Institute (NCI) Genomic Data Commons (GDC) database using the R/Bioconductor package 'TCGAbiolinks'63. BROAD javaGSEA standalone version can be downloaded from http://www.broadinstitute.org/ gsea/downloads.jsp. The RNAseq data is deposited in the European Genome-phenome Archive (EGA) under study accession number EGAS00001001552, and dataset accession EGAD00001005738. All other data is available within the Article, Supplementary Information or available from the authors upon reasonable request.