Multi-organ landscape of therapy-resistant melanoma

Metastasis and failure of present-day therapies represent the most common causes of mortality in patients with cutaneous melanoma. To identify the underlying genetic and transcriptomic landscapes, in this study we analyzed multi-organ metastases and tumor-adjacent tissues from 11 rapid autopsies after treatment with MAPK inhibitor (MAPKi) and/or immune checkpoint blockade (ICB) and death due to acquired resistance. Either treatment elicits shared genetic alterations that suggest immune-evasive, cross-therapy resistance mechanisms. Large, non-clustered deletions, inversions and inter-chromosomal translocations dominate rearrangements. Analyzing data from separate melanoma cohorts including 345 therapy-naive patients and 35 patients with patient-matched pre-treatment and post-acquired resistance tumor samples, we performed cross-cohort analyses to identify MAPKi and ICB as respective contributors to gene amplifications and deletions enriched in autopsy versus therapy-naive tumors. In the autopsy cohort, private/late mutations and structural variants display shifted mutational and rearrangement signatures, with MAPKi specifically selecting for signatures of defective homologous-recombination, mismatch and base-excision repair. Transcriptomic signatures and crosstalks with tumor-adjacent macroenvironments nominated organ-specific adaptive pathways. An immune-desert, CD8+-macrophage-biased archetype, T-cell exhaustion and type-2 immunity characterized the immune contexture. This multi-organ analysis of therapy-resistant melanoma presents preliminary insights with potential to improve therapeutic strategies.

Metastasis and failure of present-day therapies represent the most common causes of mortality in patients with cutaneous melanoma. To identify the underlying genetic and transcriptomic landscapes, in this study we analyzed multi-organ metastases and tumor-adjacent tissues from 11 rapid autopsies after treatment with MAPK inhibitor (MAPKi) and/or immune checkpoint blockade (ICB) and death due to acquired resistance. Either treatment elicits shared genetic alterations that suggest immune-evasive, cross-therapy resistance mechanisms. Large, non-clustered deletions, inversions and inter-chromosomal translocations dominate rearrangements. Analyzing data from separate melanoma cohorts including 345 therapy-naive patients and 35 patients with patient-matched pre-treatment and post-acquired resistance tumor samples, we performed cross-cohort analyses to identify MAPKi and ICB as respective contributors to gene amplifications and deletions enriched in autopsy versus therapy-naive tumors. In the autopsy cohort, private/late mutations and structural variants display shifted mutational and rearrangement signatures, with MAPKi specifically selecting for signatures of defective homologous-recombination, mismatch and base-excision repair. Transcriptomic signatures and crosstalks with t um or -a dj acent m ac ro en vi ro nments nominated organ-specific adaptive pathways. An immune-desert, CD8 + -macrophage-biased archetype, T-cell exhaustion and type-2 immunity characterized the immune contexture. This multi-organ analysis of therapy-resistant melanoma presents preliminary insights with potential to improve therapeutic strategies.
Cutaneous melanoma (CM) exhibits a UV-related high mutational burden 1,2 . Mutually exclusive BRAF and NRAS mutations drive MAPK addiction in ~70% of metastatic CM 3 . CM genomes also harbor a high burden of structural variants (SVs) 4 and chromothripsis 5 . Current knowledge of this mutational landscape is derived from tumors naive to highly active treatments developed recently and inclusive of earlier-stage disease. How MAPK inhibitor (MAPKi)/immune checkpoint blockade (ICB) therapies alter the mutational landscape and, thereby, cause death remains largely unknown.  Table 2). These TMBs are higher than or similar to previous estimates (Supplementary Table 7). We detected by Genomic Identification of Significant Targets In Cancer (GISTIC2.0) 48 significantly amplified and 75 significantly deleted regions (Extended Data Fig. 2a and Supplementary Table 8). Amplified regions harbored BRAF (80% of tumors), ACTA1 (66%) and TERT (42%); deleted regions harbored IFN cluster genes and CDKN2A/B (51%), B2M and SPRED1 (45%), BRCA1 (41%) and JAK2 and CD274 (35%) (Extended Data Fig. 2b). Tumors from six of six MAPKi-only and MAPKi+ICB cases harbored BRAF amplification, whereas tumors from two of four ICB-only cases harbored BRAF amplification. Copy number loss of JAK2, previously associated with acquired ICB resistance 21,22 , was observed in MAPKi-only cases. The findings of mechanisms of acquired MAPKi resistance in tumors with acquired ICB resistance and vice versa suggest cross-resistant mechanisms that converge on immune evasion.
To match RAM tumors, we selected TCGA-SKCM tumors driven by BRAF or NRAS mutations. By comparing RAM versus TCGA-SKCM copy number alteration (CNA) frequencies with Fisher's exact test, we detected 571 significantly amplified and 132 significantly deleted genes in the RAM cohort (Supplementary Table 9). We observed significant overlaps of amplified (but not deleted) genes between findings from GISTIC and RAM-enriched (versus TCGA-SKCM) CNAs (P = 0.0482198, hypergeometric test) (Extended Data Fig. 2c). We tested the hypothesis that a subset of RAM-enriched (versus TCGA-SKCM) CNAs is due to acquired resistance to MAPKi/ICB therapy. We identified significant CNAs in post-treatment (versus pre-treatment) melanoma from patients with either MAPKi-only or ICB-only treatment (between pre and post biopsies) histories. Notably, gene amplifications overlapped significantly between the RAM (versus TCGA-SKCM) and the MAPKi-only post (versus pre) frequency-enriched CNAs (P = 1.205907 × 10 −18 , hypergeometric test) (Fig. 1b). In contrast, gene deletions overlapped significantly between the RAM (versus TCGA-SKCM) and the ICB-only post (versus pre) frequency-enriched CNAs (P = 6.52371 × 10 −33 , hypergeometric test) (Fig. 1b). Thus, acquired resistance to ICB versus MAPKi distinctly contributes to the RAM CNA landscape.
We also nominated significantly mutated genes (SMGs) (Fig. 1c  and Supplementary Table 10). To circumvent a limited cohort size, we inflated type I error by not performing multiple testing but reduced false positives by enforcing an expression cutoff. MutSig2CV (at a raw P < 0.05 cutoff) called 110 SMGs. Among these, we nominated 62 SMGs (Extended Data Fig. 2d) with RNA (Supplementary Table 10) or high RNA (Fig. 1c and Supplementary Table 10) expression, using cohort-matched transcriptomes. Among high-expression SMGs, we identified resistance driver mutations in BRAF and NRAS and predicted loss-of-function B2M mutations (Fig. 1c). Consistent with resistance causation, B2M was significantly mutated in ICB-only post (but not pre) tumors (Extended Data Fig. 2e). Moreover, the 62 nominated SMGs enriched for immune, cell death and senescence regulation (Extended Data Fig. 2f). We cross-referenced the 62 RAM-nominated SMGs to SMGs reported by large-scale, melanoma-specific or pan-cancer studies and observed highly significant overlaps (Extended Data Fig. 2g). As context, among these publications, the consensus SMGs comprised small fractions (1.49% and 11.5% among melanoma and pan-cancer cohorts, respectively) (Extended Data Fig. 2h). Seven high-expression RAM SMGs (BRAF, NRAS, B2M, PTPRC, RAC1, TP53 and CDKN2A) were previously reported as SMGs (Fig. 1c).
Among genes affected by overlapping CNAs (Fig. 1b and Extended Data Fig. 2c), eight amplified genes harbored mutations in at least one copy (Fig. 1d,e). BRAF displayed a mean ratio of variant-to-normal allelic frequency of 1.87 (Fig. 1f), indicating selective amplification of the mutant allele. Among affected tumors, deleted genes predominantly showed single-copy loss (Fig. 1e). Several genes (BRCA1, RSPO3 Although MAPKi/ICB have become standard-of-care therapies for patients with metastatic CM in developed countries, clinical relapse occurs commonly, with multi-therapy resistance being highly lethal. Acquired MAPKi resistance in CM has been evaluated at an omics scale in a few cohorts, but knowledge of clinically acquired ICB resistance is limited [10][11][12][13][14][15][16][17][18] . Metastases to accessible anatomic sites overrepresent current datasets on acquired resistance. Monitoring (for example, liquid biopsy) and therapeutic strategies to counter resistance require insights into multi-organ mechanisms and heterogeneity. Whole-exome sequencing (WES) and whole-genome sequencing (WGS) can characterize patterns termed 'mutational signatures' that reflect imprints of DNA mutagenic processes and defective DNA damage repair processes 19,20 . It is unknown whether rare signatures of a particular malignancy, such as UV-related CM, are common with respect to late mutations, potentially due to the influence of a particular therapy. Such signatures that emerge later during tumor evolution might represent targetable pathway defects or synthetic lethalities.
Hence, we assembled a rapid autopsy melanoma (RAM) cohort from patients with BRAF MUT or NRAS MUT CM who were treated with and responded initially to MAPKi/ICB therapies but later died because of disease progression. This cohort includes 71 distinct metastatic tumors, 41 tumor-adjacent 'normal' (AN) tissues representing organ-specific tumor macroenvironments and 38 tumor-non-adjacent normal (NAN) tissues. We generated and analyzed WES from tumors and patient-matched NANs as well as WGS from a subset. To dissect the contribution of distinct therapies, we comparatively analyzed WES data from longitudinal pre-and post-tumors from patients with CM who had progressed on either MAPKi-only or ICB-only therapy. Moreover, we developed organ-specific metastatic signatures based on tumorcell-enriched transcriptomes and analyzed ligand-receptor signaling between tumor and AN tissues. Finally, we deconvolved tumor, AN and NAN transcriptomes to decipher organ-specific immune contextures.

Altered mutational spectra and signatures
We previously identified altered mutational spectra associated with somatic mutations unique to acquired MAPKi resistance 14 . Here, we divided somatic mutations into early, intermediate and late mutations and analyzed mutant allelic frequencies (Extended Data Fig. 3d), mutational spectra (Extended Data Fig. 3e) and single-base substitution (SBS) signatures based on COSMIC version 3.3 (Fig. 2b). The mean mutant allelic frequencies of early, intermediate and late somatic single-nucleotide variants (SNVs) were, respectively, 0.39, 0.36 and 0.27. Early and intermediate mutational spectra displayed a case-specific but no site-specific or treatment-specific pattern. However, late mutational spectra clustered independently of case, site or treatment history (Extended Data Fig. 3f), enriching for C>A, T>C and T>G (Extended Data Fig. 3e). Among 10 SBS signatures, UV signatures (SBS7a and SBS7b) dominated early and intermediate somatic mutations (Fig. 2b) in a case-specific pattern (Fig. 2c). Notably, non-UV-related signatures dominated late mutations ( Fig. 2b) with extensive intra-patient and inter-patient heterogeneity but treatment-elicited convergence (for example, signatures of defective homologous recombination repair (HHR) and MMR clustering with MAPKi) (Fig. 2c). Notably, SBS3 (defective HRR) was detected among late mutations in nine of 10 patients (Fig. 2b). We also detected SBS5 (clock-like), SBS9 (polymerase eta somatic hypermutation), SBS26 (defective MMR), SBS30 (defective DNA base excision repair (BER) due to NTHL1 mutations), SBS31 (platinum treatment) and SBS87 (thiopurine treatment) among late mutations in most patients. SBS11 (temozolomide) was detected among late somatic mutations in association with ICB treatment (Fig. 2b). We sought to validate an association between non-UV-related signatures and MAPKi/ICB therapies. We assembled two MAPKi validation cohorts of patient-matched melanoma with WES data: (1) a clinical BRAF MUT cohort (88 tumors from 28 patients) consisting of pre-treatment and acquired-resistant tumors (along with normal genomic DNA (gDNA)) (Supplementary Table 5) and (2) a PDX BRAF MUT or NRAS MUT cohort (29 tumors from eight models) consisting of model-matched vehicle-treated and acquired-resistant tumors derived in patient sex-matched NSG mice (information on gender of source patients in Supplementary Table 6) (along with normal gDNAs). All pre-treatment melanomas in both MAPKi cohorts were ICB naive and MAPKi sensitive (Supplementary Tables 5 and 6). We extracted somatic mutations unique to acquired MAPKi resistance and identified the frequencies of SBS signatures, including 12 non-UV signatures of defective HRR, MMR and BER (Fig. 2d,e). In contrast, UV signatures dominated the SNVs in pre/sensitive tumors. We also assembled WES data from patient-matched pre/sensitive and post/acquired ICB-resistant melanoma (14 tumors from seven patients) and normal gDNAs. All pre/ sensitive melanoma in this ICB validation cohort, except one (V52), were MAPKi naive (Supplementary Table 5). After extracting SBS signatures from somatic SNVs in pre-treatment tumors and unique to acquired resistance, we detected enrichment of UV signatures in four of seven acquired ICB-resistant tumors (Fig. 2f). Thus, MAPKi (versus ICB) therapy selects for non-UV-related mutational signatures.

Somatic whole-genomic alterations
We generated WGS data (median coverage of 23×; range, 10×-43×) from a subset of RAM tumors (one sensitive and 21 acquired-resistant; eight cases; six sites) (Supplementary Table 3). We observed no significant differences in the median SNVs or IDs based on treatment histories. The mean WGS-estimated TMB was 97 mutations per Mb (Fig. 3a), which is higher than WES-estimated TMB (Fig. 1a) and WGS-estimated TMB (39.6 mutations per Mb) in earlier-stage, MAPKi/ICB-naive CM (Supplementary Table 7). We identified a mean of 493 SVs per tumor (range, 254-1,708) ( Fig. 3a), in contrast to a mean of 342 or 106 in mucosal or CM cohorts, respectively 4 . We first classified SV/rearrangements as clustered or non-clustered 24,25 and found 67% as non-clustered translocations, 17% as non-clustered deletions and 13% as non-clustered inversions. Rearrangement signature (RS) 2 (ref. 24), defined by large (>100 kb) non-clustered deletions, inversions and inter-chromosomal translocations, was most frequent, regardless of case, site or treatment (Fig. 3b). To dissect the temporality of RSs, we focused on two RAM cases with multiple tumors to reconstruct the phylogeny and determine early ( Fig. 3c and Extended Data Fig. 4a), intermediate (Fig. 3c) and late SVs ( Fig. 3c and Extended Data Fig. 4a). Analysis of RS1 to RS6     (10 pre and 17 post from 10 patients), we observed that RS2 enrichment dominated late SVs in both pre and post tumors (Extended Data Fig. 4b).
To assess DNA double-strand break (DSB) repair mechanisms, we analyzed the breakpoint-junctional sequences of early ( Fig. 3d and Extended Data Fig. 4c), intermediate (Fig. 3d) and late SVs ( Fig. 3d and Extended Data Fig. 4c). Among early SVs, 76% of breakpoints displayed a homologous sequence (HS) size of 0-1 base pairs (bp), supporting non-homologous end joining (NHEJ) as a key DSB repair mechanism. Among intermediate and late SVs, we inferred either NHEJ or NHEJ + alternative NHEJ. We validated the importance of NHEJ across temporal SVs using WGS data from the clinical pre and post MAPKi cohort (Extended Data Fig. 4d). Thus, NHEJ and alternative NHEJ represent potential targets to blunt SV-driven melanoma progression.

Organ site-specific transcriptomic features
Consistent with 'contamination' of bulk tumors by organ-specific cell types, we observed an inverse correlation between WES-based tumor cell purities and enrichment of NAN gene expression (Extended Data Fig. 5a). However, bulk tumor transcriptomes did not segregate by cases, treatment histories or sites, possibly because of wide-ranging tumor purities (<25% to 90%) ( Fig. 4a and Supplementary Table 2). We then devised a strategy to identify tumor-cell-enriched signatures by detecting differential gene sets (DGSs) and differentially expressed genes (DEGs), where DGSs and DEGs of organ-specific metastasis are (1) depleted bioinformatically of NAN DGSs and DEGs across case-matched, cross-organ, pair-wise comparisons and (2) consistent across multiple such comparisons in ≥3 cases ( Fig. 4b and Extended Data Fig. 5b). We observed upregulated and downregulated DGSs and DEGs specific to brain, cardiac, liver, splenic, lung and ST metastases ( Fig. 4c and Extended Data Fig. 5c). We also performed Gene Ontology (GO) enrichment analysis of recurrent DEGs (Extended Data Fig. 5d). Brain metastases organ-specifically upregulated IFN response signatures and associated with oxidative phosphorylation and PI3K-AKT-mTOR signaling 27,28 . Cardiac metastases upregulated oxidative phosphorylation and response to reactive oxygen species and downregulated neurotransmitter and anabolic hypoxia genes. Liver and splenic metastases both upregulated neural genes/pathways. Liver metastases upregulated the complement pathway but downregulated IFN response genes. Lung metastases upregulated snoR-NAs and pigment biosynthesis, whereas ST metastases upregulated epidermal differentiation genes. We confirmed that organ-specific, tumor-cell-enriched transcripts were not differentially expressed by corresponding NANs (Extended Data Fig. 5e). At the protein level, we validated the glutaminergic versus the GABAnergic phenotypes of splenic (Fig. 4d) versus liver (Fig. 4e) metastases.

Organ-specific tumor immune contextures
We next evaluated organ-specific (intra)tumor immune microenvironments (TIMEs) and immune macroenvironments (within AN and NAN tissues). Analysis of 12 pan-cancer immune archetypes 29 via unsupervised clustering revealed immune macroenvironments, more than TIMEs, as organ-specific or patient-specific (Fig. 6a). We then assigned each tissue sample to the highest enrichment-scoring immune archetype and calculated the distribution of archetypes for each organ (Fig. 6b and Extended Data Fig. 6a). AN and NAN immune archetypes were similar in each organ, and TIMEs uniformly (~100%) enriched for the immune-desert, CD8 + -macrophage-biased archetype linked to T-cell exhaustion and the worst survival in patients from the TCGA-SKCM cohort 29 . This contrasted with a lower frequency (~60%) of the same archetype in TCGA-SKCM tumors 29 , 30% of which belong to two immune archetypes: T cell-centric macrophage biased and T cell-centric dendritic cell biased 29 . These and the immune-rich CD8 + or CD4 + archetypes were all but absent in our RAM cohort. Interestingly, ~5% of liver TIMEs comprised the myeloid-centric, cDC2-biased archetype (associated with tumor fibrosis). Using CIBERSORTx, we evaluated RAM versus TCGA-SKCM immune contextures ( Fig. 6c and Extended Data Fig. 6b). Across 22 immune cell types, we observed consistent patterns between AN and NAN tissues, except for plasma cells, regulatory T cells and naive/ memory B cells. Relative to RAM tumors, TCGA melanomas displayed a higher proportion of anti-tumorigenic macrophages (Fig. 6c). Across RAM organs, brain metastases enriched for pro-tumorigenic macrophages, eosinophils and resting mast cells. Using multi-spectral immunofluorescence, we confirmed a higher pro-to anti-tumorigenic tumor-associated macrophage (TAM) ratio in brain metastases (Fig. 6d and Extended Data Fig. 6c). In addition, lung metastases preferentially comprised neutrophils, potentially related to tumor-derived CCL8 (Fig. 5). Finally, RAM (versus TCGA-SKCM) enriched for T-cell exhaustion and type-2 immunity (Fig. 6e).

Discussion
This RAM study begins to build foundational insights into highly evolved and lethal CMs that resist MAPKi and/or ICB therapies. By comparative analysis of acquired-resistant CM (preceded by only one of the two types of therapies) with patient-matched pre-treatment tumors, we resolved further how each therapy distinctly and convergently shapes the high mutational, CNA and SV burdens of acquired-resistant CM. SMGs and genes altered by CNAs and SVs enrich in immune-evasive processes (for example, BRAF MUT amplification and loss-of-function alterations in B2M, JAK2, CD274/PD-L1 and PTEN; Supplementary Table 13) that may confer cross-therapy resistance, accelerating lethal disease progression. Notably, evolution of MAPKi (versus ICB) resistance shifts the mutational signatures, implicating therapy-elicited DNA damage and/or deficiency in repair pathways (for example, MMR, BER and HRR) as culprits. The evolution of late/private SVs, regardless of treatment history, enriches for RS2. Analysis of breakpoint-junctional sequences of SVs suggests NHEJ as a MAPKi or ICB co-target. Overall, multiple forms of genomic instability may cause and/or result from resistance evolution, with therapeutic implications that warrant mechanistic studies.
Acquired therapy resistance co-evolves with, and may also promote, metastatic progression. However, comparative analysis of RAM versus earlier-stage and MAPKi/ICB-naive CM cohorts is limited by cross-study technical variables (sequencing depth, read lengths, tumor purities and library preparation) and by RAM's relatively small sample size, which may contribute to false-positive SMGs. We mitigated false positives by requiring SMGs to display gene expression and by analyzing validation cohorts. Future RAM studies should expand the current cohort size and increase representation of CM subtypes, ethnic and ancestral diversities and the treatment-naive landscape.
Our analysis, by computationally depleting bulk metastatic tumor transcriptomes of patient-matched and organ-matched normal tissue-derived transcriptomes, sheds light on organ-specific metastatic signatures. Liver and spleen metastases display neural differentiation, suggesting therapeutic targets 30,31 . Melanoma brain metastasis (MBM) displays signatures of IFN signaling, oxidative phosphorylation and PI3K-AKT signaling 27,28,32 . The brain-specific macroenvironment appears to be a predominant source of IFN ligands. Overall, RAM tumors, including MBM, strongly display an immune-desert but CD8 + -macrophage-biased archetype with enrichment of T-cell exhaustion. For MBM, loss of antigen presentation and enrichment of type-2 immunity suggest TGFβ blockade and upregulation of cytotoxic natural killer (NK) cell-mediated or CD4 + T cell-mediated anti-tumor immunity 33 as potential therapeutic strategies. The pro-tumorigenic TAM phenotype of MBM also suggests therapeutic co-targets. Thus, we have uncovered a preliminary set of organ-specific metastatic signatures, tumor macroenvironment crosstalks and immune contextures that characterize therapy-resistant CM, justifying expanded RAM-based and functional analyses.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41591-023-02304-9.

Rapid autopsies and samples
We performed warm autopsies with informed consent at the University of North Carolina at Chapel Hill and the University of California, Los Angeles (UCLA). In brief, patients or persons holding the healthcare power of attorney signed the Autopsy Authorization Form and the institution-specific tumor tissue procurement and banking consent form. Research in this study involving autopsy specimens does not meet the regulatory definition of human subject research. The last tissue sample of each case was excised and stored no longer than 6 hours from death. We collected metastatic tumor, AN (≤1 cm away from tumor border) and NAN (>1 cm) tissues from 11 RAM cases. We collected formalin-fixed, paraffin-embedded (FFPE) tissues from all available sites. If tissue was sufficient, we also collected snap-frozen tissues and stored them at −80 °C. An autopsy pathologist (Leigh B. Thorne) reviewed hematoxylin and eosin tissue sections. Histopathologic analysis selected for high tumor content and against necrosis, AN with ≤10% tumor cell contamination and NAN without tumor cell contamination. We further selected for tumor purity >30% based on Sanger-sequencing-estimated BRAF or NRAS mutant allele frequencies (BRAF forward primer, 5′-GACTCTAAGAGGAAAGATGAAGTAC-3′; BRAF reverse primer, 5′-GTTGAGACCTTCAATGACTTTCTAG-3′; NRAS forward primer, 5′-GGCTTGAATAGTTAGATGCTTATTTAACCTTGGC-3′; and NRAS reverse primer, 5′-GCTCTATCTTCCCTAGTGTGGTAACCTC-3′).

Clinical samples
Tumor tissues from living patients with CM and patient-matched normal tissues were collected with the approval of institutional review boards at UCLA and Vanderbilt Ingram Cancer Center and with informed consent of each patient or the patient's legal representative. We analyzed by WES 88 tumor samples from 28 patients at UCLA with BRAF MUT CM obtained before treatments with and responses to MAPKi and then at the time of disease progression (that is, acquired resistance), along with patient-matched normal tissues (Supplementary Table 5). From PDX models collected from patients at UCLA, we analyzed by WES 29 tumors and patient-matched normal tissues (from eight patients with BRAF MUT or NRAS MUT CM) treated with vehicle or trametinib (at sufficient in vivo doses to induce tumor regression) in NSG mice (Supplementary Table 6). Also, we analyzed by WES 14 tumor samples from seven patients with CM at Vanderbilt obtained before treatments with and responses to ICB and then at the time of disease progression (Supplementary Table 5). Sex/gender was self-reported and not considered in the study design given that each cohort size was small. Participation in research was not compensated.

PDXs and treatments
To develop PDX models, tumor fragments derived from metastatic CM, with approval by UCLA institutional review boards and the UCLA Animal Research Committee, were transplanted subcutaneously in sex-matched NSG mice (4-6 weeks old) from the UCLA vivarium or Jackson Laboratory. We conducted all animal experiments in accordance with approved protocol and regulations (ARC 2016-086). We adhered to the maximal tumor size for experimental endpoints, which was ~1,500-2,000 mm 3 without mobility impairment and with body condition score >2. We implanted one tumor fragment in each mouse on the flank and measured tumors with a calliper every 2 days. Tumor volumes were calculated by (length × width 2 )/2, and we excluded data from occasional animals that died before final analysis. We assigned tumors with volumes ~500 mm 3 randomly into experimental groups. A special mouse chow (Test Diet) incorporated trametinib (LC Laboratories) to achieve daily dosing at 5 mg/kg/day.

Somatic SNVs and CNAs
We used BWA for mapping and Picard for removal of duplications. We identified somatic SNVs and IDs of tumors 34,35 by using patient-matched normal tissues for germline reference. We called SNVs using the Unified Genotyper tool of GATK, MuTect and VarScan2 and IDs using GATK-UGF, SomaticIndelDetector of GATK (IndelLocator) and VarScan2 (calls made by at least two of three algorithms). SNV/IDs were supported by Article https://doi.org/10.1038/s41591-023-02304-9 at least five reads in the tumor samples and none in the patient-matched normal tissues. Somatic SNV/IDs were then annotated by Oncotator 36 . Finally, we used Sequenza 37 to detect tumor purity, ploidy, somatic CNAs and loss-of-heterozygosity regions.

Significant CNA genes in RAM tumors
We applied GISTIC2.0 (ref. 38) to identify the significantly deleted and amplified regions using each RAM tumor's copy number segmentation file as the input. We generated circos plots with both q values and G-scores representing the amplitude of the aberration. Only regions with q values less than 0.05 were defined as significantly altered regions. We performed Fisher's exact test to identify differentially amplified or deleted genes in the RAM tumor cohort versus the BRAF or NRAS mutated TCGA-SKCM cohort (n = 345 tumors: seven stage 0, 65 stage 1, 100 stage 2, 130 stage 3, 17 stage 4 and 26 unknown). We downloaded TCGA-SKCM CNA data from cBioPortal and compiled the frequencies of CNA genes (Supplementary Table 16). We applied Fisher's exact test to individual-level amplified or deleted events, which were counted in a given RAM case if they were identified from at least one tumor of that case. We identified RAM frequency-enriched (versus TCGA-SKCM) CNA genes as amplified or deleted genes with a false discovery rate (FDR)-adjusted P < 0.05 and a higher frequency in RAM (versus TCGA-SKCM) cohort.

CNA genes enriched in post (versus pre) acquired MAPKi resistance
We compared the patient-level frequency of each amplified or deleted gene in post-treatment (acquired-resistant) versus pre-treatment melanoma by using Fisher's exact test. When multiple post tumors were available from a given patient, they were considered as an entirety, and an amplified or deleted gene was counted in this patient if it was identified from at least one of the post tumors. Multiple pre-treatment tumors from one patient were also counted as an entirety. We nominated a post-enriched amplified or deleted gene if its frequency is higher in post (versus pre) tumors and the FDR-adjusted P value is less than 0.05.

CNA genes enriched in post (versus pre) acquired ICB resistance
An amplified or deleted gene enriched in acquired ICB resistance was defined by its amplification or deletion frequency in the post tumors being >2× that in the pre tumors and by its detection in ≥2 patients' post tumors.

SMGs
We applied MutSig2CV to identify SMGs using each RAM case's mutational profile as the input and each RAM case (not each tumor sample) as an identifier. A mutation exists in a given RAM case if it was identified in ≥1 tumor. To circumvent the limitation of a small sample size, we inflated type I error by not performing multiple testing 39 and identified by MutSig2CV genes at P < 0.05. To reduce false positives, we nominated SMGs by filtering for RNA expression. Based on the mean values of normalized expression levels (log 2 counts per million (CPM)) of RAM tumors, we annotated MutSig2CV SMGs as no expression (mean log 2 CPM < 0), expressed (0 ≤ mean log 2 CPM < 4) or highly expressed (mean log 2 CPM > 4) and excluded those with no expression. We performed Fisher's exact test to compare the patient-level frequency of each SMG in the RAM versus TCGA-SKCM cohorts and used the FDR approach to adjust the P values. We then performed GO enrichment analysis with the clusterProfiler package to detect the significant biological processes for expressed SMGs. Moreover, we identified SMGs with MutSig2CV for post tumors from patients with MAPKi-only or ICB-only treatments by using each patient's mutational profile of post tumor(s) as input and each patient (not each post tumor) as an identifier. For patients with multiple post tumors, a mutation was considered to exist in this patient if it was detected in ≥1 tumor. We defined post/ acquired-resistant SMGs as genes identified by MutSig2CV at P < 0.05 in the post tumors and not identified in the patient-matched pre tumor(s).

Phylogeny and mutational signatures
We performed phylogenetic analysis using the PHYLIP program with the parsimony algorithm 35 and annotated each tree with potential drivers of tumorigenesis and/or resistance.

ITH and preferentially mutated genes
We conducted a subclonal analysis for each RAM case using PyClone-VI 41 and assessed each mutation's cancer cell fraction (CCF). Mutations were clonal if the CCF approaches 1; otherwise, mutations were subclonal. The ratio of subclonal mutations to all mutations determined ITH. We determined preferentially mutated genes for each organ as follows: (1) selected for non-synonymous mutations; (2) calculated ΔCCF (see formula in Extended Data Fig. 3b) of each mutation in the tumor of one specific organ (that is, T A ) versus tumors of other organs (that is, T B …T N ); (3) identified organ-specific enriched mutations with a ΔCCF of >0.2 for each RAM case; and (4) defined a gene as 'preferentially mutated genes' of organ A when organ-A-enriched mutations in this gene occurred in ≥3 patients.

WGS-based SV analysis
We mapped WGS reads to GRCh38/hg38 human reference genome using BWA-MEM 42

Reconstruction of focal amplifications
Focal amplicon identification and elucidation of circular extrachromosomal DNAs (ecDNAs) and complex genomic rearrangements (CGRs) using WGS data were carried out by AmpliconArchitect 50 . In brief, we determined the list of potential intervals for each amplicon, for which copy numbers and SVs were estimated using read depth and discordant read signatures. It then constructed a breakpoint graph. Simple cycles were then decomposed from the breakpoint graphs and amplicons classified into ecDNAs, CGRs and linear amplicons. We used CNVKit

Nature Medicine
Article https://doi.org/10.1038/s41591-023-02304-9 to infer the initial set of CNV seed regions. SV view of amplicons was generated using functions available in AmpliconArchitect.

WGS-based mutation analysis
Somatic SNVs and IDs were identified using Strelka2 (ref. 51) with default parameters and then subjected to mutational and ID signature analyses using the R package deconstructSigs 40 and MutationalPatterns 52 , respectively, with COSMIC SBS/ID signatures version 3.3 as reference. Classification of early, intermediate and late somatic SNVs was in accordance with that for WES-based mutations.

RNA-seq analysis
We analyzed single-end and paired-end RNA-seq data 3,35 by mapping transcriptome reads to the GRCh38/hg38 human reference genome using HISAT2. Gene-level counts were estimated by the htseqcount program. The normalized expression level of each gene, log 2 CPM, was calculated by the R package edgeR 53 and batch corrected with the 'removeBatchEffect' function in the limma R package 54 . We performed principal component analysis (PCA) using the prcomp function in R package stats to visualize the clustering of samples.

Organ tissue expression in RAM tumors
We first used the DESeq2 package to detect DEGs between NANs from one specific organ versus other organs. We defined organ-specific normal gene signatures as the top five significantly upregulated genes, all of which were confirmed to display organ-specific expression in the human protein atlas 55 . For RAM tumors, we performed single-sample gene set enrichment analysis (GSEA) to generate the enrichment scores of the organ-matched normal organ-specific gene signature. CPM values were input into the gene set variation analysis (GSVA) program using the default 'kcdf=Gaussian' option.

DGS analysis
GSEA via the fgsea package used the human hallmark (H) gene sets (MSigDB). Genes were pre-ordered by the log 2 -transformed expression fold change metrics (log 2 FC). We calculated the enrichment nominal P values by permutation test (100,000 permutations), with Benjamini-Hochberg FDR correction for multiple testing. We then performed GSEA to identify the DGSs for both tumors and NANs of one organ versus other organs. We defined normal-corrected DGSs (NC-DGSs) in tumors of organ A as those DGSs detected in the tumors of organ A (versus other organs) but not in the NAN compartment of organ A (versus other organs).

DEG analysis
To filter out normal tissue-specific DEGs from bulk tumor transcriptomes, we used RNA-seq derived from each tumor's organ-matched NANs and public datasets of normal organ/tissue gene expressions from the Illumina Human Body Map and GTEx. DEGs between a tumor pair from two organs (for example, T A versus T B ) were corrected for expression of the same genes between the normal tissue pair of the same two organs (for example, N A versus N B ). We calculated normal-corrected fold change (NC-FC) of each gene between T A versus T B as the FC of this gene between T A versus T B divided by the FC of this gene between N A versus N B . We defined the genes with NC-FC > 2 or NC-FC < 0.5 as the upregulated or downregulated NC-DEGs between these two organs. We then computed the recurrence (≥30%) of NC-DEGs by collecting NC-DEGs of tumors from one specific organ against all other tumors from other organs of the same RAM case across all RAM cases. GO enrichment analysis via the 'clusterProfiler' package 56 identified the top five significant GO biological processes for recurrently upregulated or downregulated NC-DEGs of each organ-specific metastasis.

Analysis of tumor-macroenvironment interactions
Only organ sites with ≥4 tumor-AN pairs and four tumor-NAN pairs available from ≥3 RAM cases were included for this analysis using the CellChat R package 57 . We culled significant ligand-receptor signaling that was detected from tumor-AN pairs of one specific organ site but not from tumor-NAN pairs of the same organ site. The netAnalysis_ signalingRole_heatmap function was used for visualization. The circle plots depicting tumor-AN interactions of each organ were generated by applying netVisual_aggregate with the options layout = 'circle'.

Analysis of immune contextures
We performed single-sample GSEA to generate the absolute enrichment scores of the 12 immune archetypes for each tissue, using CPM values of all expressed genes as input for GSVA 58 in the default 'kcdf=Gaussian' option. We then assigned each tissue to the highest enrichment scoring immune archetype. CIBERSORTx 59 was used in the 'absolute mode' to estimate infiltration levels of 22 immune cell types with CPM values as input. We downloaded the gene-level normalized read counts (RSEM, file name: Batch normalized from Illumina HiSeq_RNASeqV2) of TCGA-SKCM RNA-seq (inclusive of only BRAF or NRAS mutant tumors) from cBioPortal. CIBERSORTx estimated the absolute abundance of 22 immune cell types with the normalized expression level as an input. We calculated the enrichment scores of two signatures, 'T cell exhaustion' and 'type 2 immunity' 33 , by GSVA using the default 'kcdf=Gaussian' option.

Immunofluorescence
FFPE RAM tissue sections were heated at 90 °C for 25 minutes and immersed in xylene and gradient ethanol to achieve deparaffinization and re-hydration. Then, tissue sections were subjected to heat at 95 °C 10 mM citrate buffer (pH 6.0) for 15 minutes to retrieve antigens. After permeabilization and blocking with 0.1% Triton X-100/10% normal goat serum in PBS for 1 hour, tissue sections were incubated with primary antibodies, including anti-GRIK4 (Invitrogen, MA5

Multi-spectral immunofluorescence
Using Ventana Discovery Ultra (Roche) and Opal fluorophores (Akoya Biosciences), we deparaffinized 5-μm-thick tissue sections using EZ-Prep reagent (Roche) and retrieved antigens in CC1 buffer (pH 9, 95 °C; Roche). Discovery Inhibitor (Roche) was applied to inhibit enzymatic activities, followed by six sequential rounds of staining. Each round included the addition of a primary antibody followed by secondary antibody detection using either OmniMap anti-Ms HRP (Roche, 760-4310, ready-to-use) for mouse or OmniMap anti-Rb-HRP (Roche, 760-4311, ready-to-use) for rabbit following the manufacturer's specifications. We amplified signals by using Opal fluorophores at 1:400. Between rounds of staining, the tissue sections underwent heat-induced epitope retrieval to remove the primary/secondary-HRP antibody complexes before staining with the subsequent antibody. the Vectra Polaris imaging system (Akoya Biosciences). After image capture, unmixing of the spectral libraries was performed with inForm software (Akoya Biosciences). Unmixed images were then imported into HALO (Indica Labs) for stitching, cell segmentation and cell phenotyping. We analyzed whole tumor regions from each slide. Data were exported and graphed with Prism (GraphPad). Representative images were exported from HALO after spectral unmixing.

Statistical methods
We conducted statistical analyses in R 4.02, Python 3.8.0, Python 2.7.17 and Prism.

Reporting Summary
Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Data availability
The BAM files of WES, WGS and RNA-seq data are deposited in the European Genome-phenome Archive (https://www.ebi.ac.uk/ega/) with accession number EGAS00001006644. Access requires a data sharing agreement. Source data are provided with this paper. with two tissue/data sets. RAM cohort (tumor, n = 22 for WGS, n = 74 for WES, n = 93 for RNA-seq; adjacent 'normal', n = 68 for RNA-seq; non-adjacent normal, n = 67 for RNA-seq, n = 10 for WES). Pre-and-post cohort consists of patientmatched cutaneous melanoma tumors biopsied pretreatment and post acquired resistance (ICB, n = 7 pairs; MAPKi, n = 59 pairs; both subgroups with patient-matched normal tissues for all patients; n = 102 tumor WES datasets). TCGA-SKCM cohort, only BRAF MUT (n = 233) or NRAS MUT (n = 125) tumors included (13 tumors with both BRAF and NRAS mutations; total n = 345 tumors). Bolded arrows, analyses compared to derive insights into temporality or causality.