Abstract
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.
Similar content being viewed by others
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%, respectively1. Despite the clinical progress a large fraction of patients with melanoma will not benefit from ICB. The majority of non-responders 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 mechanisms3,4.
In baseline pre-treatment samples, a wide range of factors that predict ICB outcome has been reported5,6,7,8. Tumor-cell intrinsic factors include tumor mutational burden8,9,10, mutational subsets11, clonality12, aneuploidy13, immune evasion14,15, antigen presentation16,17 and interferon gamma signaling10,18,19,20. Most other predictive factors derive from T cell immunity, such as presence and infiltration of CD8+ T cells21, cytotoxicity9, expression of T cell checkpoints22, T cell receptor repertoire23 and T cell sub-populations24, in particular naïve T cells expressing TCF725. Yet, additional cells from the tumor microenvironment can modulate ICB outcome, such as B cells via the formation of tertiary lymphoid structures26, or fibroblasts via immune cell exclusion27. Notably, the constitution of the gut and tumor microbiome affects therapy outcome28. In contrast to baseline samples, few ICB resistant samples have been studied so far. Here, loss-of-function (LoF) mutations in B2M, JAK1 and JAK229, as well as genomic loss of B2M16, alone or co-existing were reported in samples with acquired resistance. In addition, there is evidence of neoantigen loss30 and T cell re-exhaustion31 in progressing tumor lesions. With regard 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 presentation32,33,34.
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. In this work, we undertake a comprehensive molecular exploration of tumor intrinsic and immune microenvironmental features to further unravel resistance to PD1 and CTLA4 blockade.
Results
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 on ICB treatment at the national Center for Cancer Immune Therapy (CCIT-DK) in Copenhagen, Denmark. Collectively, 23 metastases were from patients resistant to anti-CTLA4 monotherapy (CTLA4res) and 21 metastases from patients resistant to anti-PD1 monotherapy (PD1res). Importantly, all biopsies from metastases were taken after progression on either CTLA4 or PD1 blockade. Moreover, 17 out of the 21 PD1res patients had received and progressed or relapsed on prior anti-CTLA4 treatment. All CTLA4res patients were naïve to PD1 blockade and instead, the majority of CTLA4res patients had received prior IL-2 treatment. Seven patients had also relapsed on BRAF inhibitor (BRAFi). In addition, biopsies from 11 PD1res patients were taken at day 7 during treatment with BRAFi according to a study protocol (PD1res*). Moreover, biopsy from one CTLA4res patient was taken at day 7 during BRAFi treatment35. 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). Metastatic lesions from patients with mucosal melanoma were excluded from downstream statistical analyses due to their distinct biological characteristics36. Hence, we believe that that resistance mechanisms in such melanomas may be different from cutaneous melanomas. Using whole-exome sequencing data, tumor mutational burden, defined here by the total number of non-silent somatic mutations, was not different between CTLA4res and PD1res 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)37, and to two public datasets at ICB baseline (Supplementary Fig. 1A)5,38, 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.
As expected, samples from ICB resistant patients contained BRAF V600 (n = 24, 71%, not including mucosal samples) and NRAS Q61 (n = 7, 21%) mutations in a mutual exclusive way (Fig. 1B). Moreover, CDKN2A had the highest frequency of inactivation with nine deep deletions and four loss-of-function (LoF) mutations in ICB resistant cases (n = 13, 34%) however, this frequency was not different as compared to the treatment-naïve distant metastases from TCGA (Fig. 1B, Supplementary Fig. 1B). Other known melanoma driver genes including TP53 (n = 4, 12%) and PTEN (n = 3, 9%) also harbored inactivating events in ICB resistant cases. Driver mutations were predominantly clonal (Supplementary Fig. 2A). In a genome-wide analysis, we did not observe novel genes with high frequencies of alterations (n > = 4 hotspot mutation, LoF mutation, deletion or amplification) in ICB resistant cases (Supplementary Fig. 2B–D). Interestingly, APC and PLCB4 had three nonsense mutations each.
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 surface16. In this study, one CTLA4res case had a frameshift deletion and two PD1res lesions had deep deletions of B2M that consequently also had loss of the B2M protein in the tumor cells (Supplementary Fig. 3). Moreover, HLA-A, HLA-B, or HLA-C lacked LoF mutations, however, the HLA locus had LOH in three ICB resistant cases (9%), one in CTLA4res and two in PD1res patients (Fig. 1C). Together antigen presentation was impaired in 18% (n = 6) of ICB resistant tumors of which 12% of the CTLA4res and 24% of the PD1res (Fig. 1D). In comparison, the treatment-naïve TCGA data had LOH at the HLA locus in 13% of cases, and B2M and TAP2 had one LoF mutation each, resulting in a similar frequency of MHC-I inactivation (Supplementary Fig. 1B). In addition, in the interferon-gamma pathway that has been implicated in ICB resistance19, we found two deep deletions of the JAK2 gene and a frameshift mutation in the IFNGR1 gene, all cases being CTLA4res. We also compared tumor biopsies from patients that clinically demonstrated intrinsic or acquired resistance and did not find any differences (Supplementary Fig. 1C). Further, we compared PD1res patients with prior relapse to anti-CTLA4, to anti-PD1 naïve samples with prior relapse to anti-CTLA4 from two public datasets5,38 to isolate the effect of PD1 blockade. However, we found similar mutational landscapes across the datasets (Supplementary Fig. 1D).
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 unknown 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 (CTLA4res) and 21 from anti-PD1 resistant (PD1res) cases were analyzed. Anti-PD1 resistant metastases taken during BRAFi treatment were treated as a separate group (PD1res*) as previous studies have demonstrated an influx of immune cells in tumors during BRAFi treatment39. The single CTLA4res patient where biopsy was taken at day 7 during BRAFi treatment was excluded from downstream statistical analyses. Indeed, in this cohort, immune cell signatures40,41 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 PD1res tumors (P = 0.006) whereas the immune module was upregulated in CTLA4res melanomas (P = 0.02). In addition, a variety of immune cell type signatures15,25,41, such as T-cells, NK-cells, monocytes and dendritic cells, had higher scores in CTLA4res melanomas (Fig. 2A, Supplementary Fig. 4A). Exhaustion signatures were also at higher levels in CTLA4res 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 CTLA4res melanomas, albeit not always crossing significance thresholds (Supplementary Fig. 4B). In contrast, immune exclusion genes, e.g., MYC, demonstrated an upregulation in PD1res melanomas. Gene set enrichment analysis of discriminating genes confirmed that cell cycle processes, such as DNA replication, were elevated in PD1res tumors (FDR < 0.001) and immune processes, specifically from the adaptive immune system, were elevated in the CTLA4res cases (FDR < 0.001) (Fig. 2B). However, CD3 immunofluorescence staining could not confirm a statistically increased frequency of CD3+ T cells in the CTLA4res tumors (P = 0.12) (Fig. 2C). Finally, we used T cell receptor clonotype sequencing and found that CTLA4res tumors had a higher TCR clone abundance as well as a higher frequency of patients with dominant TCR clones compared to PD1res melanoma (P = 0.047, Fig. 2D), suggesting a higher number of tumor-reactive T cells in the CTLA4res cases. Moreover, TCR clonality results demonstrated that PD1res* had an increased clonotype count, which is also supported by abundance of CD3+ T cells in such tumors. In contrast, PD1res* did not show an increase in clonal expansion as compared to PD1res melanomas (Fig. 2D). Thus, TCR clonality analysis suggests that BRAFi is associated with a higher abundance 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 CTLA4res as compared to PD1res melanomas with PMEL being one of the most differentially expressed genes (P = 0.1, Fig. 2E). This suggests that anti-PD1 resistance correlates with decreased melanocytic antigens which may facilitate 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 levels42,43,44. 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 (three CTLA4res, one PD1res, one anti-CTLA4 resistant mucosal 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 MITFhigh/TAP1high melanoma cells suggesting an immunogenic state, which was supported by numerous spots with increased gene expression of immune cell markers across that tumor section (Supplementary Fig. 5). To expand on these findings, we used multiplex immunofluorescence (mIF) staining of 10 CTLA4res and 15 PD1res metastatic lesions, of which seven were PD1res*, and added staining on metastases from 53 patients with stage IV ICB naïve melanoma. First, we could confirm the results from the digital spatial profiling analysis showing a decreased expression of MITF in PD1res melanomas, however not reaching statistical significance when compared to CTLA4res or ICB naïve melanomas (Fig. 3C). Notably, one of the PD1res MITFhigh melanomas harbored a B2M genetic alteration. The percentage of NGFR+ melanoma cells was not associated with 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+/MITFhigh melanoma cell populations and CD8+ (Spearman cor. 0.75, P < 0.0001) and CD3+ T cell abundance (Spearman cor. 0.66, P < 0.0001). PD1res tumors were depleted of MITFhighB2M+ melanoma cell populations and enriched in MITFlowB2M− populations (Fig. 3D).
In conclusion, melanoma metastases harbor multiple tumor cell populations that are associated with 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 resemblance to secondary lymphoid organs45. We and others recently observed TLSs in ICB naïve melanoma metastases and found that such structures correlated to patient survival and ICB clinical response26,46. Here, in line with CD3+ T cell abundance, we observed a higher frequency of CD20+ B cells in CTLA4res lesions as compared to PD1res lesions that contained very 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 CTLA4res 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 CTLA4res metastases (Fig. 4B, C). None of the PD1res metastatic melanomas had TLSs. However, one PD1res 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 PD1res melanoma had a massive infiltration of CD3+ T cells and had predominantly B2M positive melanoma cells in contrast to most other PD1res 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 CTLA4res, 1 PD1res 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 (Supplementary Table 1), clustering revealed eight distinct B cell clusters: four plasma cell, two naïve B cell, one germinal center-like B cell and one memory B cell-like group (Fig. 5B). The ICB naïve 1 metastasis contained predominantly naïve B cells and germinal center-like B cells, suggesting presence of highly mature TLSs, whereas the ICB naïve 2 metastasis contained mainly plasma cells together with additional B cell phenotypes (Fig. 5C). The CTLA4res metastatic melanoma consisted of naïve and memory-like B cells and only smaller fractions of plasma cells. Finally, the PD1res metastatic melanoma, with B cells and massive T cell infiltration, had predominantly (> 80%) memory-like B cells. Memory B cells from the PD1res 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 have been reported to have immunosuppressive consequences by inducing distinct T cell phenotypes47. 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 PD1res melanoma consisted almost exclusively of effector/exhausted CD8+ T cells (Supplementary Fig. 6 A). 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 PD1res melanomas had very few infiltrating T cells, whereas abundance of CD3+ and CD8+ T cells in CTLA4res 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 CTLA4res tumors as compared to PD1res and ICB naïve melanomas (P = 0.047, Fig. 6B). Such FOXP3+ T cells were closely co-localized with CD8+ T cells (Fig. 6C).
Increased frequencies of TCF7+ CD8+ T naïve/stem cells have recently been associated with an improved clinical response to ICB25. Moreover, tumor-associated T cells lacking PD1 and TCF7 expression are suggested to be bystander T cells, specific to tumor-unrelated targets48. 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 cells49. Moreover, the double negative CD8+ T cells also had elevated proliferation rates (Supplementary Fig. 6B). We found PD1res 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 observed between CTLA4res and naïve melanomas.
In summary, the majority of CD8+ T cells infiltrating PD1res samples do not express PD1 and presumably are not tumor reactive. In contrast, CD8+ T cell phenotypes in CTLA4res 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. Notably, the majority of anti-PD1 resistant patients herein had before relapsed on CTLA4 blockade. Previous reports have converged on genetic alterations in two major pathways explaining ICB resistance: the interferon-gamma and antigen presentation pathways16,29. Specifically, a landmark study reported JAK1, JAK2 and B2M LoF mutations in three of four investigated ICB resistant samples, respectively29. Another study highlighted 5 of 12 patients progressing on ICB to harbor B2M LoF mutations or LOH16. In a recent study, 22 cell lines from 18 patients that had progressed on anti-PD1 or anti-PD1/anti-CTLA4 therapy, contained one JAK2 and two B2M inactivating events and two HLA LOH events, next to other potential ICB resistance mechanisms50. In the present 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 targeting51, 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 patients with melanoma52. 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 signaling53 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 cells39 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 cells48. 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 melanoma25. 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 in PD1 resistant cases, which may either have developed during resistance or has pre-existed and was selected for due to the high response rate of PD1 blockade.
One of the strongest predictive biomarkers to ICB across cancer diagnoses is the formation of TLS26,46. 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 functions47, 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 tolerance54. This renders melanoma to be a potentially highly immunogenic cancer type, however, in some cases immunogenicity can be low due to e.g., reduced neoantigen presentation, or down-regulation of melanocyte differentiation antigens50. Moreover, several melanoma cell states exist that have different molecular and functional properties42. 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 findings55. Our work describes that de-differentiated 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 feature56. We further demonstrate that such melanoma state is anti-correlated to CD8+ T cell presence suggesting that they escape the recognition by the immune system.
In conclusion, our work provides a comprehensive view of the molecular and immune cell landscape 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 in development of ICB resistance.
Methods
Patients
Patients with ICB resistant melanoma were treated with ICB regimens as per standard of care in Denmark until progression and were subsequently referred to CCIT, Denmark for enrollment in three different clinical trials on adoptive cell therapy35,57,58. Sample collection was done at the time of enrollment in the trials and informed consent was obtained. All three trials (NCT00937625, NCT02379195 and NCT02354690) 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 to clinical introduction of immune checkpoint blockade under the ethical permit Dnr. 101/2013 and 191/2007. Patients signed an informed consent before sample was collected. Clinical data on ICB resistant cases are summarized in Table 1.
Whole exome sequencing
Whole exome sequencing data were generated as described previously59 on a NextSeq500 instrument (Illumina) using patient tumor and blood samples. Alignment to the human reference genome (hg38) and post-alignment analyses were performed using SAREK pipeline60, as described previously61. 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.262 and Strelka263 single nucleotide variant (SNV) calls, with default settings for Strelka2 and filtering of VarScan variants as described previously61. Exonic and splice site mutations64 with a variant allele-frequency >10% were retained. Indels were called using VarScan 2.4.2 as described previously59. 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 oncoplot65 was used to screen the data for recurrently mutated genes. Mutational contexts were retrieved by deconstructSigs66 with signatures.cosmic67 as 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 bcftools68. Copy number data were generated using CONTRA 2.0.369 and segmented by GLAD70, and were merged with previously obtained copy number data59, 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.671 with settings as described previously59 and defined as variants with a cancer cell fraction below 0.95 considering a 95% confidence interval. Public mutational data were downloaded from TCGA Pan-Cancer Atlas (gdc.cancer.gov/about-data/publications/pancanatlas), Liu et al.5 and Riaz et al.38. For comparison of mutational burden, SNVs 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. For comparison of mutation landscapes with prior anti-CTLA4 relapse, for the Liu data, mucosal and acral samples and samples taken before/on anti-CTLA4 and on anti-PD1 treatment were removed, “HDEL” and “HIGH_AMP” were used as deep deletions and high amplifications, respectively, and LOH of the HLA locus was visually called from LOH plots of the region. For the Riaz data, single nucleotide variants with VAF > = 10% and tumor depth > = 7 reads from non-acral/mucosal/uveal anti-CTLA4 progressed samples were plotted.
T cell receptor sequencing
T cell receptor (TCR) sequencing was performed as previously described35. 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 MiXCR72,73 and then VDJtools74 as before35, with the difference that in the current analysis, low frequency clonotypes were not discarded.
Bulk RNA sequencing
RNA sequencing data of bulk tumor samples, in part used in a previous study of adoptive T cell therapy59, were processed to FPKM values using HISAT2 version 2.1.075 with hg38 reference and StringTie76. Transcripts with the same gene name 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 limma77. Six samples without ICB treatment were removed. Principal Component (PC) PC2 values were differentially expressed between a batch variable. PC2 was not associated with immune-, pigmentation- or proliferation-related GO terms, instead with GO terms involving “metabolic process” and “gene expression”, which frequently indicate batch effects in in-house datasets. We therefore removed PC2 from the data using swamp78. Gene signatures41 were used as mean expression scores of available genes. In addition, a raw count matrix was generated for the same set of samples and genes, where again transcripts with the same gene name were summed up. These data were utilized to test for differential expression of single genes, using DESeq279 and adjusting for PC2 values. Gene set enrichment analysis80 was performed for biological process ontology terms using clusterProfiler81, with 12,513 genes having a standard deviation >0.4 of transformed FPKM values and being ranked by DESeq2 test statistics from raw counts.
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 h. 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.182, the data were reduced to protein-coding 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 SCTransform83, 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. The processed bulk and single cell RNAseq gene expression datasets are deposited at GEO with accession number GSE244984.
Spatial RNA 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, GSE11597815 and GSE12057525, which were processed as described previously26, 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 were further processed using Seurat82. 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 were merged and normalized using SCTransform83 (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 ConsensusClusterPlus84, using Euclidean distance and 50 iterations of 80% of spots.
Multiplex immunofluorescence
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. Deparaffinization in EZ prep (70 °C 8 min). 2. Cell conditioning was applied (CC1, 95 °C 40 min). 3. Blocked with inhibitor CM (37 °C 4 min). 4. Primary antibody incubation. 5. HRP-conjugated antibody incubation (37 °C 16 min) (Roche Discovery OmniMap anti-Rb or anti-Ms HRP). 6.Tyramide-coupled fluorescent dye incubation (37 °C 16 min). 7. Antibody denaturation (CC2, 100 °C 8 min). Steps 4-7 were repeated until all intended markers had been fluorescently labeled. Counterstaining was performed using DAPI (0.75 mM, 37 °C 8 min). Antibodies and fluorophores used in the described staining steps are summarized in Supplementary Tables 2-5.
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.285 using a QuPath script available in GitHub.
Digital image analysis
QuPath version 0.3.2 was used for all digital image analysis85. 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 QuPath85. 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.
Three mIF panels consisting of 6 markers each were evaluated (panel 1: NGFR/MITF/SOX10/CD3/CD20/CD31, panel 2: B2M/MITF/SOX10/CD3/CD20/Ki67, panel 3: CD8/PD1/TCF7/Ki67/FOXP3/SOX10). Cores with insufficient amount of tumor cells or high background were removed after visual inspection. Using QuPath software, initial marker calls were screened, and a set of manually confirmed calls was used to train marker calling. Final marker calls were assigned to the cells using a pre-specified calling order. B2M, NGFR and MITF were called within SOX10+ cells. TCF7, PD1 and Ki67 were called within CD8+ 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 (Supplementary Tables 2-5).
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.
Reporting summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Data availability
Public mutational data was downloaded from TCGA Pan-Cancer Atlas (gdc.cancer.gov/about-data/publications/pancanatlas), Liu et al.5 and Riaz et al.38. Gene signatures were obtained from referenced publications, respectively. Publicly available data with accession numbers GSE115978 and GSE120575 were downloaded from Gene Expression Omnibus (GEO) and were used to identify tumor-specific genes.
Processed bulk RNA sequencing and single cell RNA sequencing data have been deposited at GEO with accession number GSE244982 and GSE244983. Raw data are not available for these GEO submissions, as due to Swedish and Danish laws, the patient consent, and the risk that the sequencing data contains personally-identifiable information and hereditary mutations, we cannot deposit the short sequencing read data in a public access repository.
Spatial transcriptomics data have been deposited under accession number GSE261347.
Whole exome sequencing- and T cell receptor sequencing data were deposited in European Genome Archive (EGA) under EGAD50000000380 and EGAD50000000379, respectively. These data are available under restricted access. Data access can be granted via the EGA under collaborative conditions and when aligned with current ethical approval, and data will be available for duration of the proposed project. Somatically called mutations are available as Supplementary Data 1.
All other remaining data are available within the Article, Supplementary Information or as Source data file. Source data are provided with this paper.
References
Larkin, J. et al. Five-Year Survival with Combined Nivolumab and Ipilimumab in Advanced Melanoma. N. Engl. J. Med. 381, 1535–1546 (2019).
Gide, T. N., Wilmott, J. S., Scolyer, R. A. & Long, G. V. Primary and Acquired Resistance to Immune Checkpoint Inhibitors in Metastatic Melanoma. Clin. Cancer Res. 24, 1260–1270 (2018).
Sharma, P., Hu-Lieskovan, S., Wargo, J. A. & Ribas, A. Primary, Adaptive, and Acquired Resistance to Cancer Immunotherapy. Cell 168, 707–723 (2017).
Marine, J. C., Dawson, S. J. & Dawson, M. A. Non-genetic mechanisms of therapeutic resistance in cancer. Nat. Rev. Cancer 20, 743–756 (2020).
Liu, D. et al. Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma. Nat. Med. 25, 1916–1927 (2019).
Miao, D. et al. Genomic correlates of response to immune checkpoint blockade in microsatellite-stable solid tumors. Nat. Genet. 50, 1271–1281 (2018).
Jenkins, R. W., Barbie, D. A. & Flaherty, K. T. Mechanisms of resistance to immune checkpoint inhibitors. Br. J. Cancer 118, 9–16 (2018).
Litchfield, K. et al. Meta-analysis of tumor- and T cell-intrinsic mechanisms of sensitization to checkpoint inhibition. Cell 184, 596–614 e514 (2021).
Van Allen, E. M. et al. Genomic correlates of response to CTLA-4 blockade in metastatic melanoma. Science 350, 207–211 (2015).
Cristescu, R. et al. Pan-tumor genomic biomarkers for PD-1 checkpoint blockade-based immunotherapy. Science 362, eaar3593 (2018).
Litchfield, K. et al. Escape from nonsense-mediated decay associates with anti-tumor immunogenicity. Nat. Commun. 11, 3800 (2020).
McGranahan, N. et al. Clonal neoantigens elicit T cell immunoreactivity and sensitivity to immune checkpoint blockade. Science 351, 1463–1469 (2016).
Davoli, T., Uno, H., Wooten, E. C. & Elledge, S. J. Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy. Science 355, eaaf8399 (2017).
Spranger, S., Bao, R. & Gajewski, T. F. Melanoma-intrinsic beta-catenin signalling prevents anti-tumour immunity. Nature 523, 231–235 (2015).
Jerby-Arnon, L. et al. A Cancer Cell Program Promotes T Cell Exclusion and Resistance to Checkpoint Blockade. Cell 175, 984–997 e924 (2018).
Sade-Feldman, M. et al. Resistance to checkpoint blockade therapy through inactivation of antigen presentation. Nat. Commun. 8, 1136 (2017).
Rodig, S. J. et al. MHC proteins confer differential sensitivity to CTLA-4 and PD-1 blockade in untreated metastatic melanoma. Sci. Transl. Med. 10, eaar3342 (2018).
Ayers, M. et al. IFN-gamma-related mRNA profile predicts clinical response to PD-1 blockade. J. Clin. Invest. 127, 2930–2940 (2017).
Shin, D. S. et al. Primary Resistance to PD-1 Blockade Mediated by JAK1/2 Mutations. Cancer Discov. 7, 188–201 (2017).
Newell, F. et al. Multiomic profiling of checkpoint inhibitor-treated melanoma: Identifying predictors of response and resistance, and markers of biological discordance. Cancer Cell 40, 88–102 e107 (2022).
Tumeh, P. C. et al. PD-1 blockade induces responses by inhibiting adaptive immune resistance. Nature 515, 568–571 (2014).
Auslander, N. et al. Robust prediction of response to immune checkpoint blockade therapy in metastatic melanoma. Nat. Med. 24, 1545–1549 (2018).
Roh, W. et al. Integrated molecular analysis of tumor biopsies on sequential CTLA-4 and PD−1 blockade reveals markers of response and resistance. Sci. Transl. Med. 9, eaan3788 (2017).
Gide, T. N. et al. Distinct Immune Cell Populations Define Response to Anti-PD-1 Monotherapy and Anti-PD-1/Anti-CTLA-4 Combined Therapy. Cancer Cell 35, 238–255 e236 (2019).
Sade-Feldman, M. et al. Defining T Cell States Associated with Response to Checkpoint Immunotherapy in Melanoma. Cell 175, 998–1013 e1020 (2018).
Cabrita, R. et al. Tertiary lymphoid structures improve immunotherapy and survival in melanoma. Nature 577, 561–565 (2020).
Mariathasan, S. et al. TGFbeta attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature 554, 544–548 (2018).
Gopalakrishnan, V. et al. Gut microbiome modulates response to anti-PD-1 immunotherapy in melanoma patients. Science 359, 97–103 (2018).
Zaretsky, J. M. et al. Mutations Associated with Acquired Resistance to PD-1 Blockade in Melanoma. N. Engl. J. Med. 375, 819–829 (2016).
Anagnostou, V. et al. Evolution of Neoantigen Landscape during Immune Checkpoint Blockade in Non-Small Cell Lung Cancer. Cancer Discov. 7, 264–276 (2017).
Koyama, S. et al. Adaptive resistance to therapeutic PD-1 blockade is associated with upregulation of alternative immune checkpoints. Nat. Commun. 7, 10501 (2016).
Lawson, K. A. et al. Functional genomic landscape of cancer-intrinsic evasion of killing by T cells. Nature 586, 120–126 (2020).
Patel, S. J. et al. Identification of essential genes for cancer immunotherapy. Nature 548, 537–542 (2017).
Manguso, R. T. et al. In vivo CRISPR screening identifies Ptpn2 as a cancer immunotherapy target. Nature 547, 413–418 (2017).
Borch, T. H. et al. Clinical efficacy of T-cell therapy after short-term BRAF-inhibitor priming in patients with checkpoint inhibitor-resistant metastatic melanoma. J. Immunother. Cancer 9, e002703 (2021).
Newell, F. et al. Comparative Genomics Provides Etiologic and Biological Insight into Melanoma Subtypes. Cancer Discov. 12, 2856–2879 (2022).
Cancer Genome Atlas, N. Genomic Classification of Cutaneous Melanoma. Cell 161, 1681–1696 (2015).
Riaz, N. et al. Tumor and Microenvironment Evolution during Immunotherapy with Nivolumab. Cell 171, 934–949 e915 (2017).
Cooper, Z. A. et al. Distinct clinical patterns and immune infiltrates are observed at time of progression on targeted therapy versus immune checkpoint blockade for melanoma. Oncoimmunology 5, e1136044 (2016).
Cirenajwis, H. et al. Molecular stratification of metastatic melanoma using gene expression profiling: Prediction of survival outcome and benefit from molecular targeted therapy. Oncotarget 6, 12297–12309 (2015).
Becht, E. et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 17, 218 (2016).
Rambow, F., Marine, J. C. & Goding, C. R. Melanoma plasticity and phenotypic diversity: therapeutic barriers and opportunities. Genes Dev. 33, 1295–1318 (2019).
Tsoi, J. et al. Multi-stage Differentiation Defines Melanoma Subtypes with Differential Vulnerability to Drug-Induced Iron-Dependent Oxidative Stress. Cancer Cell 33, 890–904 e895 (2018).
Hoek, K. S. et al. In vivo switching of human melanoma cells between proliferative and invasive states. Cancer Res 68, 650–656 (2008). [pii] 10.1158/0008-5472.CAN-07-2491.
Lauss, M., Donia, M., Svane, I. M. & Jonsson, G. B Cells and Tertiary Lymphoid Structures: Friends or Foes in Cancer Immunotherapy? Clin. Cancer Res 28, 1751–1758 (2022).
Helmink, B. A. et al. B cells and tertiary lymphoid structures promote immunotherapy response. Nature 577, 549–555 (2020).
Shalapour, S. et al. Inflammation-induced IgA+ cells dismantle anti-liver cancer immunity. Nature 551, 340–345 (2017).
van der Leun, A. M. & Schumacher, T. N. An atlas of intratumoral T cells. Science 374, 1446–1447 (2021).
Li, H. et al. Dysfunctional CD8 T Cells Form a Proliferative, Dynamically Regulated Compartment within Human Melanoma. Cell 176, 775–789 e718 (2019).
Lim, S. Y. et al. The molecular and functional landscape of resistance to immune checkpoint blockade in melanoma. Nat. Commun. 14, 1516 (2023).
Nielsen, M. et al. Coexisting Alterations of MHC Class I Antigen Presentation and IFNgamma Signaling Mediate Acquired Resistance of Melanoma to Post-PD-1 Immunotherapy. Cancer Immunol. Res 10, 1254–1262 (2022).
Sharma, A. et al. Anti-CTLA-4 Immunotherapy Does Not Deplete FOXP3(+) Regulatory T Cells (Tregs) in Human Cancers-Response. Clin. Cancer Res. 25, 3469–3470 (2019).
Levine, A. G., Arvey, A., Jin, W. & Rudensky, A. Y. Continuous requirement for the TCR in regulatory T cell function. Nat. Immunol. 15, 1070–1078 (2014).
Lo, J. A. et al. Epitope spreading toward wild-type melanocyte-lineage antigens rescues suboptimal immune checkpoint blockade responses. Sci. Transl. Med. 13, eabd8636 (2021).
Karras, P. et al. A cellular hierarchy in melanoma uncouples growth and metastasis. Nature 610, 190–198 (2022).
Haas, L. et al. Acquired resistance to anti-MAPK targeted therapy confers an immune-evasive tumor microenvironment and cross-resistance to immunotherapy in melanoma. Nat. Cancer 2, 693–708 (2021).
Andersen, R. et al. T cells isolated from patients with checkpoint inhibitor-resistant melanoma are functional and can mediate tumor regression. Ann. Oncol. 29, 1575–1581 (2018).
Andersen, R. et al. Long-lasting complete responses in patients with metastatic melanoma after adoptive cell therapy with tumor-infiltrating lymphocytes and an attenuated IL-2 regimen. Clin. Cancer Res 22, 3734–3745 (2016).
Lauss, M. et al. Mutational and putative neoantigen load predict clinical benefit of adoptive T cell therapy in melanoma. Nat. Commun. 8, 1738 (2017).
Garcia, M. et al. Sarek: A portable workflow for whole-genome sequencing analysis of germline and somatic variants. F1000Res. 9, 63 (2020).
Sanna, A. et al. Tumor genetic heterogeneity analysis of chronic sun-damaged melanoma. Pigment Cell Melanoma Res. 33, 480–489 (2020).
Koboldt, D. C. et al. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 22, 568–576 (2012).
Kim, S. et al. Strelka2: fast and accurate calling of germline and somatic variants. Nat. Methods 15, 591–594 (2018).
Wang, K., Li, M. & Hakonarson, H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res 38, e164 (2010).
Mayakonda, A., Lin, D. C., Assenov, Y., Plass, C. & Koeffler, H. P. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res 28, 1747–1756 (2018).
Rosenthal, R., McGranahan, N., Herrero, J., Taylor, B. S. & Swanton, C. DeconstructSigs: delineating mutational processes in single tumors distinguishes DNA repair deficiencies and patterns of carcinoma evolution. Genome Biol. 17, 31 (2016).
Alexandrov, L. B. et al. Signatures of mutational processes in human cancer. Nature 500, 415–421 (2013).
Li, H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics 27, 2987–2993 (2011).
Li, J. et al. CONTRA: copy number analysis for targeted resequencing. Bioinformatics 28, 1307–1313 (2012).
Hupe, P., Stransky, N., Thiery, J. P., Radvanyi, F. & Barillot, E. Analysis of array CGH data: from signal ratio to gain and loss of DNA regions. Bioinformatics 20, 3413–3422 (2004).
Carter, S. L. et al. Absolute quantification of somatic DNA alterations in human cancer. Nat. Biotechnol. 30, 413–421 (2012).
Bolotin, D. A. et al. MiXCR: software for comprehensive adaptive immunity profiling. Nat. Methods 12, 380–381 (2015).
Bolotin, D. A. et al. Antigen receptor repertoire profiling from RNA-seq data. Nat. Biotechnol. 35, 908–911 (2017).
Shugay, M. et al. VDJtools: Unifying Post-analysis of T Cell Receptor Repertoires. PLoS Comput Biol. 11, e1004503 (2015).
Kim, D., Paggi, J. M., Park, C., Bennett, C. & Salzberg, S. L. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat. Biotechnol. 37, 907–915 (2019).
Pertea, M., Kim, D., Pertea, G. M., Leek, J. T. & Salzberg, S. L. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat. Protoc. 11, 1650–1667 (2016).
Ritchie, M. E. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43, e47 (2015).
Lauss, M. et al. Monitoring of technical variation in quantitative high-throughput datasets. Cancer Inf. 12, 193–201 (2013).
Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550 (2014).
Subramanian, A. et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl Acad. Sci. USA 102, 15545–15550 (2005).
Wu, T. et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innov. (Camb.) 2, 100141 (2021).
Hao, Y. et al. Integrated analysis of multimodal single-cell data. Cell 184, 3573–3587 e3529 (2021).
Hafemeister, C. & Satija, R. Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression. Genome Biol. 20, 296 (2019).
Wilkerson, M. D. & Hayes, D. N. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics 26, 1572–1573 (2010).
Bankhead, P. et al. QuPath: Open source software for digital pathology image analysis. Sci. Rep. 7, 16878 (2017).
Jiang, P. et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med 24, 1550–1558 (2018).
Tirosh, I. et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science 352, 189–196 (2016).
Acknowledgements
This work was supported by Swedish Research Council (Vetenskapsrådet Dnr 2018-02786, Dnr 2022-00871, GJ), Swedish Cancer Society (19 0458 Pj, 22 2105 Pj, GJ), Berta Kamprad Foundation (GJ and BP) and the governmental funding for healthcare research (ALF and GJ), Knut and Alice Wallenberg Foundation (KAW 2022.0066, GJ) and Göran Gustafsson Foundation (GJ). The authors would like to acknowledge Clinical Genomics Lund, SciLifeLab and Center for Translational Genomics (CTG), Lund University, for providing expertise and service with sequencing and analysis. The computations and data handling were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) partially funded by the Swedish Research Council through grant agreement no. 2018-05973. We also thank Nanostring for GeoMx experimental analyses.
Funding
Open access funding provided by Lund University.
Author information
Authors and Affiliations
Contributions
M.L. conducted all data analysis, design of study and writing the manuscript. B.P. conducted all immunostaining and analysis of such data. T.H.B. retrieved and interpreted clinical data from immunotherapy resistant cases. K.H. conducted whole-exome-, RNA-, TCR and single cell RNA sequencing. K.K. performed single cell RNA sequencing. A.E. and I.H. performed sectioning for Visium sequencing. J.Y. conducted and analysed single cell RNA sequencing data from B cells. K.N., C.I., A.C. and K.I. collected clinical information on immune checkpoint blockade naïve patients and set up protocols for collection of viable tumor tissue from patients with melanoma. K.P. performed and analysed multiplex immunofluorescence data. I.M.S. designed study, collected and interpreted clinical data from immunotherapy resistant cases. M.D. supervised and designed study, collected and interpreted clinical data from immunotherapy resistant cases. G.J. supervised and designed study and wrote the manuscript. All authors have read and approved the manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors have no competing interests.
Peer review
Peer review information
Nature Communications thanks Antoni Ribas and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Source data
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Lauss, M., Phung, B., Borch, T.H. et al. Molecular patterns of resistance to immune checkpoint blockade in melanoma. Nat Commun 15, 3075 (2024). https://doi.org/10.1038/s41467-024-47425-y
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41467-024-47425-y
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate.