Molecular patterns of resistance to immune checkpoint blockade in melanoma

Immune checkpoint blockade (ICB) has improved outcome for patients with metastatic melanoma but not all benefit from treatment. Several immune- and tumor intrinsic features are associated with clinical response at baseline. However, we need to further understand the molecular changes occurring during development of ICB resistance. Here, we collect biopsies from a cohort of 44 patients with melanoma after progression on anti-CTLA4 or anti-PD1 monotherapy. Genetic alterations of antigen presentation and interferon gamma signaling pathways are observed in approximately 25% of ICB resistant cases. Anti-CTLA4 resistant lesions have a sustained immune response, including immune-regulatory features, as suggested by multiplex spatial and T cell receptor (TCR) clonality analyses. One anti-PD1 resistant lesion harbors a distinct immune cell niche, however, anti-PD1 resistant tumors are generally immune poor with non-expanded TCR clones. Such immune poor microenvironments are associated with melanoma cells having a de-differentiated phenotype lacking expression of MHC-I molecules. In addition, anti-PD1 resistant tumors have reduced fractions of PD1+ CD8+ T cells as compared to ICB naïve metastases. Collectively, these data show the complexity of ICB resistance and highlight differences between anti-CTLA4 and anti-PD1 resistance that may underlie differential clinical outcomes of therapy sequence and combination.


Introduction
Immune checkpoint blockade (ICB) therapy has had a major clinical success in advanced stage melanoma.Objective response rates for anti-CTLA4, anti-PD1 and combination therapy of anti-CTLA4 with anti-PD1 were 19%, 45% and 58%, respectively 1 .Despite the clinical progress a large fraction of melanoma patients will not benefit from ICB.The majority of nonresponders are primary resistant, and a smaller fraction acquires resistance to ICB during treatment.Primary resistance manifests shortly after treatment and is accompanied by progressive disease (PD), whereas acquired resistance is observed after a period of time with initial complete response (CR) or partial response (PR) 2 .In principle, factors leading to disease progression can be pre-existing at baseline, acquired genetically, or adapted non-genetically, with possible interplay between these resistance mechanisms 3,4 .
In baseline pre-treatment samples, a wide range of factors that predict ICB outcome has been reported [5][6][7][8] .Tumor-cell intrinsic factors include tumor mutational burden [8][9][10] , mutational subsets 11 , clonality 12 , aneuploidy 13 , immune evasion 14,15 , MHC antigen presentation 16,17 and interferon gamma signaling 10,18,19 .Most other predictive factors derive from T cell immunity, such as presence and infiltration of CD8 + T cells 20 , cytotoxicity 9 , expression of T cell checkpoints 21 , T cell receptor repertoire 22 and T cell sub-populations 23 , in particular naïve T cells expressing TCF7 24 .Yet, additional cells from the tumor microenvironment can modulate ICB outcome, such as B cells via the formation of tertiary lymphoid structures 25 , or fibroblasts via immune cell exclusion 26 .Notably, the constitution of the gut and tumor microbiome affects therapy outcome 27 .In contrast to baseline samples, few ICB resistant samples have been studied so far.Here, loss-of-function (LoF) mutations in B2M, JAK1 and JAK2 28 , as well as genomic loss of B2M 16 , alone or co-existing were reported in samples with acquired resistance.
In addition, there is evidence of neoantigen loss 29 and T cell re-exhaustion 30 in progressing tumor lesions.With regards to acquired resistance, it is informative that CRISPR screens of tumor cells evading either PD1 blockade or T cell co-culture converge on inactivation of two pathways: Interferon-gamma signaling and MHC-I antigen presentation [31][32][33] .Despite these efforts, a full picture of the molecular mechanisms explaining ICB resistance is lacking, due to a paucity of tumor samples available at or after ICB progression.In addition, it is unclear whether CTLA4-and PD1 blockade resistant samples are substantially different.We undertook a comprehensive molecular exploration of tumor intrinsic and immune microenvironmental features to further unravel resistance to PD1 or CTLA4 blockade.

Genetic analysis of melanoma metastases resistant to ICB
To dissect molecular alterations associated with tumors resistant to different ICB regimens, we have collected tissue samples at progression from ICB treatment at the national Center for Cancer Immune Therapy (CCIT-DK) in Copenhagen, Denmark.Collectively, 23 metastases were from patients progressing on anti-CTLA4 monotherapy and 21 metastases from patients progressing on anti-PD1 monotherapy.17 out of 21 patients resistant to anti-PD1 had received and progressed or relapsed on prior anti-CTLA4 treatment.All patients resistant to anti-CTLA4 were naïve to PD1 blockade and instead, the majority of anti-CTLA4 resistant patients had received prior IL-2 treatment.Seven patients had also relapsed on BRAF inhibitor (BRAFi) and twelve anti-PD1 resistant samples were taken at day 7 during treatment with the BRAFi vemurafenib according to a study protocol 34 .Most patients displayed primary resistance to ICB, except six cases that clinically had acquired resistance.(Table 1).Site of primary melanoma was cutaneous skin or unknown primary, except three cases from patients with primary mucosal melanoma that were resistant to anti-CTLA4 (Table 1).Using whole-exome sequencing data, tumor mutational burden, defined here by the total number of non-silent somatic mutations, was not different between anti-CTLA4 and anti-PD1 resistant tumors.Moreover, we compared our data to distant metastases from the cancer genome atlas (TCGA) cohort (n=68), where prior systemic therapy did not include ICB (Fig. 1A) 35 , and to two public datasets at ICB baseline (Fig. S1A) 36,37 , and did not observe any differences in mutational burden.Together, mutational frequencies were similar between metastases progressing or relapsing on either anti-CTLA4 or anti-PD1 and ICB naïve melanomas.
In summary, the frequencies of genetic alterations in melanoma driver genes were not different in anti-CTLA4 and anti-PD1 resistant and ICB naïve metastatic melanomas.

Genetic alterations in immune regulatory pathways
Next, we specifically analyzed genetic alterations occurring in immune regulatory pathways.
Previously, loss of B2M was reported to be associated with resistance to ICB in melanoma and is essential for HLA class I assembly and presentation on the cell surface 16 .In this study, one anti-CTLA4 resistant case had a frameshift deletion and two anti-PD1 resistant lesions had deep deletions of B2M.Moreover, HLA-A, HLA-B, or HLA-C lacked LoF mutations, however, the HLA locus had LOH in three ICB resistant cases (9%), thereof one in anti-CTLA4 resistant and two in anti-PD1 resistant patients (Fig. 1C).Together antigen presentation was impaired in 18% (n=6) of ICB resistant tumors of which 12% of the anti-CTLA4 resistant and 24% of the anti-PD1 resistant (Fig. 1D).In comparison, the treatment-naïve TCGA data had LOH at the HLA-A locus in 14% of cases (n=9 of 66 evaluable cases), and B2M and TAP2 had one LoF mutation each, resulting in a similar frequency of MHC-I inactivation (Fig. S1B).In addition, in the interferon-gamma pathway that has been implicated in ICB resistance 19 , we found two deep deletions of the JAK2 gene and a frameshift mutation in the IFNGR1 gene, all cases being anti-CTLA4 resistant.We also compared tumor biopsies from patients that clinically demonstrated intrinsic or acquired resistance and did not find any differences (Fig.

S1C).
Together, genetic alterations in genes belonging to antigen presentation-or interferon-gamma pathways occurred in different samples, however in total only accounting for 26% of ICB resistant cases suggesting that other still unkown immune evasive mechanisms exist.

Immune transcriptional programs are different in anti-CTLA4 and anti-PD1 resistant melanomas
As the genetic landscape in the cohort only explained a minority of ICB resistance we performed transcriptomic profiling using RNA sequencing.In total, 17 (non-mucosal) melanoma metastases from anti-CTLA4 resistant and 21 from anti-PD1 resistant cases were analyzed.Anti-PD1 resistant metastases taken during BRAFi treatment were treated as a separate group as previous studies have demonstrated an influx of immune cells in tumors during BRAFi treatment 38 .Indeed, in this cohort, immune cell signatures 39,40 were increased in BRAFi treated tumors and this was specifically pronounced when comparing within the anti-PD1 resistant tumors alone (Fig. 2A).Moreover, the cell cycle module was upregulated in anti-PD1 resistant melanomas not taken during BRAFi (P = 0.009) whereas the immune module was upregulated in anti-CTLA4 resistant melanomas (P = 0.09).In addition, a variety of immune cell type signatures 15,24,40 , such as T-cells, NK-cells, monocytes and dendritic cells, had higher scores in anti-CTLA4 resistant melanomas (Fig. 2A, Fig. S3A).Exhaustion signatures were also at higher levels in anti-CTLA4 resistant melanomas, probably due to a higher infiltration of immune cells.Further, we also observed single gene expression of MHC-I, MHC-II, interferon gamma signaling, Tumor/T-cell interaction, inflammation and cytotoxicity genes, all generally at higher levels in anti-CTLA4 resistant melanomas, albeit not always crossing significance thresholds (Fig. S3B).In contrast, immune exclusion genes, e.g., MYC, demonstrated an upregulation in anti-PD1 resistant melanomas.Gene set enrichment analysis of discriminating genes confirmed that cell cycle processes were elevated in anti-PD1 resistant tumors (FDR<0.001)and immune processes, specifically from the adaptive immune system, were elevated in the anti-CTLA4 resistant cases (FDR<0.001)(Fig. 2B).In addition, CD3 immunofluorescence staining confirmed an increased frequency of CD3 + T cells in the anti-CTLA4 resistant tumors (P = 0.12) (Fig. 2C).Finally, we used T cell receptor clonotype sequencing and found that anti-CTLA4 resistant tumors had a higher TCR clone abundance as well as a higher frequency of patients with dominant TCR clones compared to anti-PD1 resistant melanoma (P = 0.047, Fig. 2D), suggesting a higher number of tumor-reactive T cells in the anti-CTLA4 resistant cases.Moreover, TCR clonality results demonstrated that anti-PD1 resistant melanomas taken during BRAFi treatment had an increased clonotype count, which is also supported by abundance of CD3 + T cells in such tumors.In contrast, anti-PD1 resistant melanomas taken during BRAFi treatment did not show an increase in clonal expansion as compared to anti-PD1 resistant melanomas without BRAFi (Fig. 2D).Thus, TCR clonality analysis suggests that BRAFi induces an influx of T cells in melanomas but not expansion of specific TCR clones.
To understand differences between ICB regimens with respect to tumor intrinsic properties, we selected tumor cell regions using S100B/PMEL antibodies for digital spatial profiling.Most genes in the panel were upregulated in anti-CTLA4 resistant as compared to anti-PD1 resistant melanomas with PMEL being one of the most differentially expressed genes (P = 0.1, Fig. 2E).
This suggests that anti-PD1 resistance may lead to downregulation of melanocytic antigens that subsequently leads to immune escape.
In summary, the observations indicate a sustained immune response in some tumors despite progressing on anti-CTLA4; in contrast, a particularly immune-poor microenvironment was observed in anti-PD1 resistant melanoma.

Distinct tumor cell states exist in melanoma metastases
Several intrinsic tumor cell states have been suggested for melanoma, which generally align on a gradient of MITF-low to MITF-high expression levels [41][42][43] .These cell states can switch to adapt to external cues.The association of such melanoma cell states to the immune microenvironment and immunotherapy resistance has not been thoroughly investigated.We therefore performed Visium sequencing on six melanoma metastases (four anti-CTLA4 resistant, one anti-PD1 resistant and one ICB naïve), and defined melanoma cell states as characterized by differential expression of MITF, MKI67, NGFR, AXL and TAP1.Consensus clustering of 2,766 tumor cell-enriched spots resulted in five groups where four were characterized mainly by differential expression of MITF and TAP1 (Fig. 3A).The fifth group had decreased levels of MITF and increased levels of NGFR gene expression.By morphologically mapping such melanoma states on H&E stainings we found an extensive heterogeneity of melanoma cell states in all ICB resistant tumors (Fig. 3B).Whereas the ICB naïve metastasis harbored predominantly MITF high /TAP1 high melanoma cells suggesting an immunogenic state, which was supported by numerous spots with increased gene expression of immune cell markers across that tumor section (Fig. S4).To expand on these findings, we used mIF staining of 10 anti-CTLA4 and 15 anti-PD1 resistant metastatic lesions, of which seven were taken under BRAFi treatment, and added metastases from 53 stage IV ICB naïve melanoma patients.First, we could confirm the results from the digital spatial profiling analysis showing a decreased expression of MITF in anti-PD1 resistant melanomas, however not reaching statistical significance when compared to anti-CTLA4 resistant or ICB naïve melanomas (Fig. 3C).Notably, one of the anti-PD1 resistant MITF high melanomas harbored a B2M genetic alteration.The percentage of NGFR + melanoma cells was not associated to ICB resistance (P = 0.96) as very few melanoma tumors harbored notable fractions of NGFR + melanoma cells.We then mimicked the melanoma cell states identified by the Visium sequencing using mIF antibodies for SOX10, MITF, B2M and NGFR.As expected, we found a correlation between presence of B2M + /MITF high melanoma cell populations and CD8 + (Spearman cor.0.75, P <0.001) and CD3 + T cell abundance (Spearman cor.0.66, P<0.001).
In conclusion, melanoma metastases harbor multiple tumor cell populations that are correlated to T cell infiltration.

B cells and tertiary lymphoid structures are rare in anti-PD1 resistant metastatic melanomas
Tertiary lymphoid structures (TLS) are ectopic immune cell niches with resemblances to secondary lymphoid organs 44 .We and others recently observed tertiary lymphoid structures in ICB naïve melanoma metastases and found that such structures correlated to patient survival and ICB clinical response 25,45 .Here, in line with CD3 + T cell abundance, we observed a higher frequency of CD20 + B cells in anti-CTLA4 resistant lesions as compared to anti-PD1 resistant lesions that contained extremely few CD20 + B cells (P = 0.06).When compared to melanoma metastases from ICB naïve patients no significant difference was found most likely due to the vast heterogeneity of CD20 + B cell presence observed in ICB naïve cases (Fig. 4A).Two anti-CTLA4 resistant metastases harbored multiple immature TLSs, while six ICB naïve metastases had at least one TLS which appeared as more mature than TLSs found in anti-CTLA4 resistant metastases (Fig. 4B-C).None of the anti-PD1 resistant metastatic melanomas had TLSs.However, one anti-PD1 resistant melanoma harbored an increased frequency of CD20 + B cells, which were not organized in TLS and instead were scattered and localized at the tumor margin (Fig. 4D).Intriguingly, this anti-PD1 resistant melanoma had a massive infiltration of CD3 + T cells and had predominantly B2M positive melanoma cells in contrast to most other anti-PD1 resistant melanomas.To further understand lymphocyte phenotypes in ICB resistant and naïve metastatic melanomas we performed single cell RNA sequencing of four tumors with TLS or B cells (1 anti-CTLA4 resistant, 1 anti-PD1 resistant and 2 ICB naïve).The 26,053 sequenced cells stemmed from a wide range of cell types including B cells (Fig. 5A).Using a set of B cell and TLS specific genes (Table S1), clustering revealed eight distinct B cell clusters, four plasma cell, two naïve B cell, one germinal centerlike B cell and one memory B cell-like group (Fig. 5B).ICB naïve 1 contained predominantly naïve B cells and germinal center-like B cells, suggesting presence of highly mature TLSs, whereas ICB naïve 2 contained mainly plasma cells together with additional B cell phenotypes (Fig. 5C).The anti-CTLA4 resistant metastatic melanoma consisted of naïve and memory-like B cells and only smaller fractions of plasma cells.Finally, the anti-PD1 resistant metastatic melanoma, with B cells and massive T cell infiltration, had predominantly (>80%) memory-like B cells.Memory B cells from the anti-PD1 resistant melanoma were mainly IgA + cells, in contrast to the samples with conventional TLS, which had large fractions of IgG + memory B cells (Fig. 5D).Indeed, IgA expressing B cells were reported to have immunosuppressive consequences by inducing distinct T cell phenotypes 46 .With this in mind, we observed 11 T cell clusters in the single cell RNA sequencing data, including naïve T cells, T follicular helper cells, T regulatory cells and effector/exhausted T cells (Fig. 5E).Strikingly, T cells from the anti-PD1 resistant melanoma consisted almost exclusively of effector/exhausted CD8 + T cells (Fig. S5A).Some of the CD8 + T cells expressed PD1, they however lacked TCF7 expression (Fig. 5F), suggesting a lack of a replenishing T cell reservoir.
Altogether, scRNAseq data reveal distinct B cell phenotypes in ICB resistant metastatic melanomas that may be linked to T cell phenotype.

Distinct T cell phenotypes are enriched in ICB resistant metastatic melanomas
Consequently, we went on to investigate T cell presence and different phenotypes in ICB resistant compared to ICB naïve cases using mIF.A wide range of CD3 + and CD8 + T cell abundance was observed in ICB naïve metastatic melanomas (Fig. 6A), and no significant difference was observed between ICB resistant and naïve cases.However, the majority of anti-PD1 resistant melanomas had very few infiltrating T cells, whereas abundance of CD3 + and CD8 + T cells in anti-CTLA4 resistant tumors was indistinguishable from ICB naïve melanomas.Immunosuppressive T regulatory cells are specifically characterized by expression of the transcription factor FOXP3. Intriguingly, we found an increase of FOXP3 + T cells in anti-CTLA4 resistant tumors as compared to anti-PD1 resistant and ICB naïve melanomas (P = 0.047, Fig. 6B).Such FOXP3 + T cells were closely co-localized with CD8 + T cells (Fig. 6C).
Moreover, increased frequencies of TCF7 + CD8 + T naïve/stem cells have recently been associated with an improved clinical response to ICB 24 .Moreover, tumor-associated T cells lacking PD1 and TCF7 expression are suggested to be bystander T cells, specific to tumorunrelated targets 47 .Therefore, we classified CD8 + T cells using a combination of TCF7 and PD1 expression, using mIF (Fig. 6D-E).As expected, PD1 + cells were more proliferating (Ki67 + ) than double positive and TCF7 + CD8 + T cells 48 .Moreover, the double negative CD8 + T cells also had elevated proliferation rates (Fig. S5B).We found anti-PD1 resistant melanomas to have a significantly lower frequency of PD1 + CD8 + T cells (P = 0.03) and consequently an increased frequency of double negative (PD1 -/TCF7 -) CD8 + T cells (Fig. 6F).No difference of TCF7 + or double positive CD8 + T cells was identified between anti-CTLA4 resistant and naïve melanomas.
In summary, the majority of CD8 + T cells infiltrating anti-PD1 resistant samples do not express PD1 and presumably are not tumor reactive.In contrast, CD8 + T cell phenotype in anti-CTLA4 resistant melanoma were indistinguishable from ICB naïve melanoma however instead an increase of FOXP3 + T cells was observed.These results indicate that ICB immune microenvironmental resistance mechanisms are different in anti-PD1 and anti-CTLA4 resistant tumor lesions.

Discussion
A complete picture of genetic and molecular effector mechanisms explaining ICB resistance has so far not been identified.In this study, we combined analyses of the tumor immune microenvironment and tumor intrinsic features on human melanoma specimens taken at progression from patients receiving PD1 or CTLA4 blockade monotherapy.Previous reports have converged on genetic alterations in two major pathways explaining ICB resistance: the interferon-gamma and antigen presentation pathways 16,28 .In this study, we found only 18% harboring genetic alterations in the major genes within the antigen presentation pathway, and 9% had genetic alterations in genes belonging to the interferon-gamma pathway.Importantly, genetic alteration in either pathway can be sufficient to develop resistance to ICB based on CTLA4 or PD1 targeting 49 , and we could not definitely determine differences between patients relapsing on CTLA4 or PD1 blockade.However, immune response transcriptional signatures demonstrated increased expression of such genes in anti-CTLA4 resistant samples.Indeed, CTLA4 blockade has demonstrated an increased influx of T cells in post-treatment samples from melanoma patients 50 .In this study, tumor infiltrating T cells in anti-CTLA4 resistant samples had a significantly increased number of expanded TCR clones as compared to anti-PD1 resistant samples suggesting that such T cells have an increased tumor-reactivity.
Intriguingly, we found an increased FOXP3 + T cell abundance in anti-CTLA4 resistant tumor lesions and there are reports describing that the immunosuppressive properties of FOXP3 + T cells are dependent on TCR signaling 51 suggesting that the increased TCR clonality reflects an increased immunosuppressive environment mediated by FOXP3 + T regulatory cells.
Interestingly, patient samples taken during BRAFi treatment also had a high fraction of CD3 + T cells but no indication of expansion of distinct TCR clones.Several studies have described that BRAFi induces recruitment of immune cells 38 and our data further add evidence to this.
The PD1 blockade resistant samples, instead, contained very few T cells.Only a small fraction of the T cell infiltrate is considered to be tumor-reactive, and expression of PD1 is a relevant marker to distinguish tumor-reactive from bystander T cells 47 .In addition, T cells expressing TCF7 were found to be more effectively reinvigorated by ICB and have been associated with improved response to ICB in human melanoma 24 .Here, we found that an inferior number of CD8 + /PD1 + T cells are present in anti-PD1 resistant samples.Instead, an increased number of bystander (PD1 -/TCF7 -) CD8 + T cells were found in anti-PD1 resistant melanomas.Overall, this suggests that the tumor specific immune response is considerably hampered when anti-PD1 resistance is developed.
One of the strongest predictive biomarkers to ICB across cancer diagnoses is the formation of TLS 25,45 .In this study, we found immature TLSs in anti-CTLA4 resistant tumors, however no TLS was identified in anti-PD1 resistant melanoma.One anti-PD1 resistant melanoma had a massive infiltration of effector/exhausted T cells that was accompanied by spatially scattered IgA + memory B cells.Indeed, IgA + B cells have been described to also confer regulatory functions 46 , and this tumor may sustain a suppressive immune cell niche.These results demonstrate that we need to know more about the different functional subsets and contexts of tumor-associated lymphocytes.
T cells are activated by exposure to specific antigens.Melanoma tumors can potentially present many different neoantigens as melanoma harbors extensive tumor mutational burden.
Further, antigen presentation in tumor cells can be expanded to highly expressed melanocyte differentiation self-antigens, which is regulated by immune tolerance 52 .In melanoma, T cell recognition can be reduced using down-regulation of such melanocyte differentiation antigens 53 .Moreover, several melanoma cell states exist that have different molecular and functional properties 41 .We used spatial transcript sequencing and found five different melanoma cell states.Interestingly, we found a vast heterogeneity of the spatial distribution of the different melanoma states that is similar to previous findings 54 .Our work describes that dedifferentiated melanoma cells lacking B2M expression are frequently observed in anti-PD1 resistant melanomas.This is in line with previous reports and indicates that a de-differentiated melanoma state represents a pan-therapy resistance feature 55 .We further demonstrate that such melanoma state is anti-correlated to CD8 + T cell presence suggesting that they escape the recognition of the immune system.
In conclusion, our work provides a comprehensive view of molecular and immune cell changes occurring at relapse on ICB.Our data further highlight that development of anti-CTLA4 and anti-PD1 resistance occur through different molecular mechanisms.This study highlights the molecular complexity existing in development of ICB resistance.

Patients
ICB resistant melanoma patients were referred to CCIT, Denmark for enrollment in three different clinical trials on adoptive cell therapy 34,56,57 1.

Whole exome sequencing
Whole exome sequencing data were generated as described previously 58 on a NextSeq500 instrument (Illumina) using patient tumor and blood samples.Alignment to the human reference genome (hg38) and post-alignment processing were performed using SAREK pipeline 59 , as described previously 60 .Median target coverage of tumor samples ranged from 38 to 165 (median = 81) and that of patient-matched normal samples from 37 to 106 (median = 73).Mutations were called using the intersection of VarScan 2.4.2 61 and Strelka2 62 SNV variant calls, with default settings for Strelka2 and filtering of VarScan variants as described previously 60 .Exonic and splice site mutations 63 with a variant allele-frequency > 10% were retained.Indels were called using VarScan 2.4.2 as described previously 58 .The final variant set is available as Supplementary Data 1.Loss-of-function mutations were defined as frameshift, nonsense, or splice site mutations or multiple gene events from different categories.
Maftools oncoplot 64 was used to screen the data for recurrently mutated genes.Mutational contexts were retrieved by deconstructSigs 65 with signatures.cosmic 66as reference.Loss-of-Heterozygosity of the HLA locus was called visually from plots of B-allele frequencies under the condition of heterozygous germline background.B-allele frequencies of common germline SNPs (dbSNP version 151) were obtained using samtools mpileup and bcftools 67 .Copy number data were generated using CONTRA 2.0.3 68 and segmented by GLAD 69 , and were merged with previously obtained copy number data 58 , on hg38 co-ordinates.Deep deletions and high amplifications were defined as values <(-1) and >1, respectively.Subclonal mutations were identified using ABSOLUTE 1.0.6 70 with settings as described previously 58 and defined as variants with a cancer cell fraction below 0.95 considering a 95% confidence interval.Public mutational data was downloaded from TCGA Pan-Cancer Atlas (gdc.cancer.gov/aboutdata/publications/pancanatlas),Liu et al 5 and Riaz et al 71 .For comparison of mutational burden, single nucleotide variants with VAF >= 10% and tumor depth >=7 reads were retained when such information was available, and each external cohort was combined with our cohort using the set of genes mutated in both datasets.

T cell receptor sequencing
T cell receptor (TCR) sequencing was performed as previously described 34 .Briefly, DNAse I (Thermo Scientific) treated tumor RNA samples were subjected to library preparation using AmpliSeq Immune Repertoire Panel (Illumina) and sequenced using NextSeq500.Data were analyzed with MiXCR 72,73 and then VDJtools 74 as before 34 , with the difference that in the current analysis, low frequency clonotypes were not discarded.

Bulk RNA sequencing
RNA sequencing data of bulk tumor samples were processed to FPKM values using Tophat2 75 and Cufflinks 76 as described previously 77 .Isoforms of the same gene were summed up, the data were limited to protein-coding genes and log-transformed as log2(data+1).RNAseq data were quantile-normalized together using limma 78 .Samples without ICB treatment were removed.Principal Components (PC) PC2 and PC6 values were differentially expressed between MM909 and the current MM1413/1414 data.However, PC2 and PC6 were not associated with immune-, pigmentation-or proliferation-related GO terms, instead with GO terms involving "metabolic process", "gene expression" and "translation", which frequently indicate batch effects in in-house datasets.We therefore removed PC2 and PC6 from the data using swamp 79 .Gene signatures 10,15,24,39,40,80,81 were used as mean expression scores of available genes.The processed gene expression dataset will be deposited at GEO. Gene set enrichment analysis 82 was performed for biological process ontology terms using clusterProfiler 83 , with 11,690 genes having a standard deviation >0.4 being ranked by t-test statistics.

Single cell RNAseq
We generated scRNAseq data from four tumor samples available as finely chopped cryopreserved material.Samples were gently thawed and dissociated using Dri Tumor & Tissue Dissociation Reagent (BD Horizon) according to manufacturer's protocol, with digestion incubation times up to 1 hour.Dead cells were removed using Dead Cell Removal Kit (Miltenyi Biotech) prior to processing the remaining single cell suspension using Chromium Single Cell 3' Kit with Dual Index Kit TT Set A sample barcodes (10x Genomics) according to manufacturer's recommendations.Libraries were sequenced on NovaSeq6000 (Illumina) with read length settings 28-10-10-90 as per 10x Genomics User Guide.The h5 files were processed and merged using the R package Seurat 4.0.1 84 , the data was reduced to proteincoding genes, translational (RPS/RPL) and mitochondrial genes (MT-), and genes which a maximum count <= 4 were removed.Cell with less than 500 expressed genes were removed.
The data were normalized using SCTransform 85 , counts that were zero before transformation were set back to zero, and data was log-transformed as log2(data+1).This resulted in a dataset of 11,606 genes and 26,053 cells.The data were visualized using UMAP on the top 30 principal components and clusters were identified using FindNeighbors (using top 30 PCs, k=15) and FindClusters (Louvain algorithm, resolution=0.3)functions of Seurat.Biological identities were assigned to the clusters after manual inspection.B cells were defined as cluster 8 and neighboring cluster 13 cells, without expression of CD3D, CD3E, MITF or SOX10 (n=559).T cells were defined as cluster 5 and 6 cells, without expression of CD79A, MITF and SOX10 (n=2,921).B and T cell data were separately re-normalized using SCTransform, zeros were restored, and data were log-transformed as log2(data+1).The data were then reduced to curated markers for B and T cells, respectively.The data were visualized using UMAP (without PC reduction) and clusters were identified using FindNeighbors (k=10, cosine distance) and FindClusters (Leiden algorithm, resolution=0.6).Biological identities were assigned to the clusters after manual inspection.CD8 T cells were defined as either expressing CD8A or CD8B, i.e., having non-zero expression.Similarly, presence of TCF7 or PD1 expression was defined by non-zero expression.

Nanostring GeoMx
Region of interests (ROI) from tissue microarray cores were selected using Immunofluorescence with CD3 (T-cell), CD20 (B-cell), PMEL/S100B (tumor cell) and DAPI antibodies.The Cancer Transcriptome Atlas assay was performed on the ROIs according to the manufacturer's instructions (Nanostring, Seattle, WA).Each sample was scaled by a factor to obtain the same 75% quantile, and the data were quantile-normalized and log-transformed as log2 (data+1).ROIs with limited tumor material in matched IHC cores as well as from patients not included in the study cohort were dismissed.The data were reduced to tumor cell specific genes, using two public single cell RNAseq datasets, GSE115978 15 and GSE120575 24 , which were processed as described previously 25 , and values >2 being considered as "expressed".Tumor cell-specific genes were defined as expressed in <20% combined CD4/ CD8 T-cells for both datasets and >10% malignant cells, respectively.To account for varying tumor cell content, tumor ROIs were divided by SOX10 expression.Multiple tumor ROIs of the same patient were merged by mean expression values.

Visium Spatial Transcriptomics
Fresh frozen tumor tissue was sectioned onto Visium Spatial Transcriptomics slides containing 4,992 barcoded spots with 55um diameter.After hematoxylin & eosin staining the sections were permeabilized according to the manufacturer's recommendations, with optimal permeabilization time previously determined to be 24 minutes following Tissue Optimization Guidelines (10x Genomics).Sequencing libraries were prepared in accordance with Visium Spatial Gene Expression User Guide with Dual Index Kit TT Set A sample barcodes (10x Genomics) according to manufacturer's recommendations.Libraries were sequenced on NovaSeq6000 (Illumina) with read length settings 28-10-10-90 as per 10x Genomics User Guide.Reads were processed together with histology images using SpaceRanger.Data was further processed using Seurat 84 .Specifically for sample MM909_37 a lymph node area was discarded.The samples contained a median of 9,997 median UMI per spot (range 3,987-19,121) and 3,427 median genes per spot (range 1,887-5,497).Protein coding genes were retained and translational genes (RPS, RPL) and mitochondrial genes (MT-) were further removed.Spots with less than 500 expressed genes were removed.The data was merged and normalized using SCTransform 85 (assay=Spatial), counts that were zero before transformation were set back to zero, and data was log-transformed as log2(data+1).Tumor cell-enriched spots were defined as SOX10 expression >=2, and CD3D, CD3E and CD79A expression <=1 (n=2,766 spots).The curated tumor marker genes were centered and clustered with ConsensusClusterPlus 86 , using Euclidean distance and 50 iterations of 80% of spots.

Multiplex immuno-Fluorescence
Paraffin-embedded TMA core tissue sections (3 μm) were baked for 1 hour at 65 °C and were subjected to deparaffinization and immuno-fluorescence staining in Roche's automatic samples preparation system (Ventana Discovery Ultra) in the following steps.1.
Antibodies and fluorophores used in the described staining steps are summarized in (Supplementary Tables 1-4).

Image Acquisition
All multiplex IF-stained (mIF) TMAs were scanned using the PhenoImager HT (Akoya Biosciences) at 20x magnification.Images were obtained through tile scanning using 7-color whole-slide unmixing filters.These filters included DAPI + Opal 570/690, Opal 480/620/780, and Opal 520.To ensure accurate signal specificity of the obtained images the synthetic Opal library in the image processing software InForm version 2.4.11(Akoya Biosciences) was used for spectral unmixing.Obtained tiles were subsequently stitched together with QuPath version 0.3.2 87using a QuPath script available in GitHub.

Digital Image Analysis
QuPath version 0.3.2 was used for all digital image analysis.Visual inspections of the tissue cores were performed to exclude samples of poor-quality including samples with high fluorescent background, insufficient amount of tumor cells and degraded tissue cores.Cell segmentation was performed using the StarDist extension running on the pre-trained model dsb2018_heavy_augment.pb in QuPath.The fluorescently labeled markers were analyzed using a machine learning classifier, random forest, trained on multiple measurements in QuPath.After establishing classifiers for each biomarker, they are subsequently combined and applied in a sequential manner.To avoid unexpected classes (e.g.CD20/SOX10, CD20/CD8…), marker calls were assigned to the cells using a pre-specified calling order.The analysis grouped together cells expressing MITF and/or SOX10 as melanoma cells.
Percentage of cells with a given call were calculated for each core, as a fraction of all cells of a core, including cells without a call; or as fraction within the relevant cell type (e.g., B2M in SOX10 + cells, or PD1 in CD8 + cells).Additionally, mean MITF intensity was calculated within SOX10 + cells, with MITF intensity defined as log2 (nucleus signal + 1).Multiple cores of the same sample were merged by mean percentage, and MITF intensity was merged by mean intensity.Panels were combined using NGFR in SOX10 + cells, CD3 and CD20 from panel 1, MITF/B2M in SOX10 + cells and MITF intensity in SOX10 + cells from panel 2, and CD8 and TCF7/PD1/Ki67 in CD8 + cells from panel 3.

Statistical Analyses
Bioinformatical analyses were performed using R version 4.0.5.For group comparisons, T-test and Wilcoxon test were used for two groups, Anova and Kruskal-Wallis test for more than two groups.Pearson correlation was used to compare numerical variables, and Fisher's exact test for categorical variables.All tests were two-sided.False discovery rate was calculated using Benjamini-Hochberg adjustment.Boxplots are displayed with the center-line as median, the box limits as lower and upper quartiles, and with whiskers covering the most extreme values within 1.5 x Interquartile-Range.

Figure 1 .
Figure 1.Mutational landscape of ICB resistant melanoma.A. Tumor mutational burden (TMB) calculated as total number of somatic non-silent mutations in anti-CTLA4-or anti-PD1 resistant as compared to the Cancer Genome Atlas (TCGA).B. Genetic aberrations of selected genes in anti-CTLA4-and anti-PD1 resistant melanomas.Tumor mutational burden and mutational signatures are indicated on top.Frequency plots of activating events for potential oncogenes and loss-of-function events for potential tumor suppressor genes are depicted on the right for ICB resistant tumors excluding mucosal samples, and for the ICB naïve control cohort, respectively.C. Genetic alterations of genes in the interferon gamma and MHC-I pathways.D. Frequency plot of immune regulatory pathways in anti-CTLA4 and anti-PD1 resistant melanomas, considering only loss-of-function events.

Figure 2 .
Figure 2. The immune transcriptomic landscape of ICB resistant melanoma.A. Heatmap of melanoma module-and immune pathway transcriptomic scores.The anit-PD1* group consists of anti-PD1 resistant samples taken at day 7 under BRAFi treatment.* P < 0.05 between anti-PD1 without and under BRAFi treatment.^ P < 0.05 between anti-PD1 without BRAFi and anti-CTLA4 resistant lesions.P-values from t-test.B. Top ten pathways from gene set enrichment analysis of genes differentiating anti-CTLA4 from anti-PD1 resistant melanomas.Red = enriched in anti-CTLA4, blue = enriched in anti-PD1.C. Fraction of CD3 + cells as determined by multiplex immunofluorescence.Representative images are shown.Pvalue from Wilcoxon test.D. Boxplots of T cell receptor (TCR) clonality data.Left boxplot shows the clonotype count (Efron Thisted).P-value from Kruskal-Wallis text.Right boxplot shows the evenness according to normalized Shannon Wiener index.P-value from Wilcoxon test F. Spatial gene expression data of tumor cell regions using GeoMx, normalized for SOX10 expression.Volcano plot showing differentially expressed genes between anti-CTLA4 and anti-PD1 resistant tumors.Boxplot of PMEL expression between the two groups.P-values from t-test.

Figure 3 .FractionFigure 4 .
Figure 3. Tumor cells states in ICB resistant melanoma.A. Heatmap displaying expression of melanoma state specific genes across 2766 tumor cell enriched spots in six melanoma tumors.Spots were divided in five distinct clusters based on consensus clustering.B. Mapping back melanoma cell states as defined in A) to the respective histological images.C. MITF multiplex immunofluorescence intensity of ICB naïve, anti-CTLA4 resistant, anti-PD1 resistant and anti-PD1 resistant under BRAFi treatment (PD1res*).MITF intensity was measured in SOX10 positive melanoma cells.P-value from Wilcoxon test D. Frequency plot of different melanoma cell states using multiplex immunofluorescence of MITF/B2M fractions and NGFR fractions.NGFR fractions are within SOX10 positive melanoma cells.Samples are sorted by CD8 fraction and grouped according to treatment.Green -ICB naïve, orange -anti-CTLA4 resistant, red -anti-PD1 resistant, purple -anti-PD1 resistant under BRAFi.

Figure 5 .
Figure 5. Functional B and T cell phenotype composition using single cell RNA sequencing.A. UMAP of 26,053 single cells from one anti-PD1 resistant, one anti CTLA4 resistant and 2 ICB naïve melanomas, that contained B cells.Cell types derived from manual annotation are indicated in the plot.B. UMAP of 559 B cells visualizing eight distinct clusters.B cell subsets are indicated in the plot after manual annotation.C. Pie charts describing the fraction of each B cell subset in the four melanomas.D. Violin plots showing the expression of IGHA1 and IGHG1 in memory-like B cells in the four melanomas.E. UMAP of 2,921 T cells visualizing 11 clusters, in the four melanomas.T cell subsets are indicated in the plot after manual annotation.F. CD8 T cell phenotype fractions based on non-zero expression of PD1 and TCF7 in CD8A or CD8B expressing T cells of all four melanomas.

Figure 6 .
Figure 6.T cell phenotypes in ICB resistant melanomas.A. CD3 + and CD8 + T cell fraction shown as boxplots.P values from Kruskal-Wallis test.B. FOXP3 + T cell fraction shown as boxplot.P value from Wilcoxon test.C. Multiplex immunofluorescence images from an anti-CTLA4 resistant melanoma that has infiltration of FOXP3 + T cells.FOXP3 + T cells are colocalized with CD8 + T cells and marked by an arrowhead.D. Multiplex immunofluorescence images from an anti-PD1 resistant melanoma that has a strong infiltration of CD8 + T cells.PD1 - /TCF7 -double negative CD8 + T cells areas are marked.E. Multiplex immunofluorescence images from an ICB naive melanoma that has a strong infiltration of CD8 + T cells.PD1 + CD8 + T cells area is marked.F. Boxplots of the fraction of PD1 + TCF7 -of CD8 + T cells and the fraction of double negative T cells (PD1 -/TCF7 -).P values from Wilcoxon test.
. All three trials are listed in clinicaltrials.gov,and all procedures were conducted in accordance with the Declaration of Helsinki and following approval from the Scientific Ethics Committee of the Capital Region of Denmark.Metastatic melanoma lesions from ICB naïve patients were collected at Skåne University Hospital in Sweden prior clinical introduction of immune checkpoint blockade under the ethical permit Dnr.101/2013 and 191/2007.Clinical data on ICB resistant cases are summarized on Table