Recurrent Tumor Cell–Intrinsic and –Extrinsic Alterations during MAPKi-Induced Melanoma Regression and Early Adaptation

The only clinically validated approach (i.e., BRAF inhibitor + MEK inhibitor) to suppress acquired BRAF inhibitor (BRAFi) resistance in melanoma addresses only MAPK-reactivating mechanisms (1–7). Recent appreciation of the pervasiveness of nongenomic and immune evolution in melanoma with acquired or late MAPK inhibitor (MAPKi) resistance (2) suggests dynamic and continuous alterations early on therapy. This is consistent with the clinical experience that the majority of initial tumor responses to MAPKi are partial, with many acquired resistant tumors later growing in situ from these residual tumors. The equilibrium between the tumor cell and stromal/immune compartments is expected to be highly dynamic immediately after initiating therapy. Thus, there is a need to develop models that would recapitulate temporally either tumor cell–intrinsic or tumor cell–extrinsic adaptations in regressing/persisting melanoma early during MAPKi therapy.

It is known that melanoma tumors responding to MAPKi therapy display robust T-cell infiltration (8). However, the extent and nature of stromal/immune compartment remodeling and potential adaptations, including adaptive immune resistance, has not been systematically characterized in regressing tumors responding to targeted therapy. Compensatory survival signaling, through therapy-induced relief of feedback inhibition on the MAPK pathway (9) and rebound upregulation of the PI3K–AKT pathway (10, 11), has been shown to participate as adaptive mechanisms. However, it is not clear whether early adaptation of melanoma to MAPKi is limited to compensatory signaling or involves transcriptomic or epigenomic reprogramming. Recently, tumor cell mesenchymal transition has been linked to an immune-suppressive microenvironment (12, 13), wound-healing signatures, and innate anti–PD-1 resistance (14). Thus, in theory, large-scale phenotypic changes of tumor cells in response to targeted therapy can have profound effects on the tumor microenvironment and antitumor immunity. We postulated that macroscopic or clinically evident tumor shrinkage occurs in parallel with the development of therapy tolerance followed by microscopic foci of resistance. In this case, understanding subclinical resistance evolution may offer insights into combinatorial therapeutic targets to suppress such early resistance.

Despite the heterogeneous kinetics of melanoma regression and timing of tumor biopsies on MAPKi therapy, we sought insights into the early adaptive landscape by cataloging recurrent gene and gene set expression alterations across patient-matched tumors on treatment (On-Tx) versus pretreatment (Fig. 1A). Across distinct tumors or patients, transcriptome reprogramming, if any, may transit at variable rates along successive stages. Thus, we captured temporally dynamic transcriptomic states by profiling BRAF^V600-mutant human melanoma cell lines or immune-competent murine melanoma tumors after increasing durations of continuous MAPKi exposure (Fig. 1B). To nominate tumor cell–intrinsic and tumor cell–extrinsic (stromal/immune) transcriptomic alterations, we compared regressing clinical tumors against in vitro–treated human melanoma cell lines and in vivo–treated murine melanoma tumors, respectively. To analyze stromal/immune transcriptomic alterations in murine melanoma, we subtracted transcriptomic alterations in the CD45-negative subpopulation from those in the bulk regressing tumors. We also provided proof-of-concept functional studies to support adaptive, tumor cell–intrinsic and/or –extrinsic roles of RTK or PD-L2 (encoded by PDCD1LG2) upregulation.

We obtained patient-matched melanoma tumors (total n = 46) before MAPKi therapies (BRAFi, MEKi, or BRAFi + MEKi) at baseline (n = 20) and On-Tx (n = 26; 23 of 26 biopsies on day 22 or earlier) after objective clinical responses (Fig. 1A; Supplementary Table S1A; ref. 15). Residual melanoma cells in “responding” (vs. baseline) tumors showed marked decreases in Ki-67 staining or proliferation (Supplementary Fig. S1A). In vitro, BRAFi-sensitive BRAF^V600E melanoma cell lines on continuous BRAFi treatment reached a population nadir, consisting of slow-cycling, drug-tolerant persisters (DTP), prior to resuming growth as drug-tolerant proliferating persisters (DTPP; Fig. 1B; refs. 11, 16). In most cell lines examined, DTPs acquired a flattened morphology and displayed senescence-associated β-galactosidase (SA-βGal) activity, which was later lost in fibroblastoid DTPP cells (Supplementary Fig. S1B–S1D; refs. 3, 11).

To delineate dynamic transcriptomic states, we profiled parental (P) cell lines treated with DMSO/vehicle, temporal subpopulations (2 d, DTP, DTPP; days to weeks on BRAFi), and long-term sublines (months to years on BRAFi or BRAFi + MEKi, resulting in single drug–resistant (SDR) or double drug–resistant (DDR) lines (R-line), respectively; Fig. 1B; Supplementary Fig. S1E and S1F). We then performed unsupervised clustering based on the top 3,000 most variant genes expressed by 18 isogenic P-lines versus R-lines (SDR or DDR; Fig. 1C). Two major clusters emerged: one group of R-lines clustered with all the P-lines; another group of R-lines clustered away from their isogenic P-lines. The group of R-lines clustering with all the P-lines harbored specific drivers of MAPK-reactivation/addiction (i.e., BRAF^V600E ultra-amplification, M249 DDR4; BRAF^V600E amplification + MEK1^F129L, M249 DDR5; BRAF^V600E splicing, M397 SDR and M395 SDR1; BRAF^V600E amplification, M395 SDR2; NRAS^MUT, M249 SDR; refs. 3, 5, 6, 17). The second group of R-lines displayed transcriptomic reprogramming, relatively weak BRAF^V600E signature enrichment, and attenuated MAPK dependency (Supplementary Fig. S1G and S1H). Thus, there are two distinct functional and transcriptomic classes of MAPKi resistance: Resistance with MAPK addiction (Ra) versus Resistance with MAPK redundancy (Rr). Accordingly, we observed from principal component analysis (PCA) two temporal transcriptomic trajectories: (i) for Rr, successive transcriptomic reprogramming away from MAPKi-sensitive isogenic P-lines (positioned at the origin), and (ii) for Ra, transient transcriptomic reprogramming followed by reestablishment of a transcriptomic state closely matching the isogenic P-lines (Fig. 1D).

We then assessed how transcriptomic changes within On-Tx or disease-progressive (DP) tumors (vs. patient-matched baseline tumors) would compare relative to the PC trajectory observed for Rr lines (Fig. 1E). Notably, in contrast to the hetero-geneous patterns of DP tumors, nearly all On-Tx tumors shifted with the directionality of the DTP→DTPP→Rr trajectory, suggesting converging phenotypes. To understand the biological processes of Rr transcriptomic reprogramming, we performed gene ontology (GO) enrichment analysis of upregulated or downregulated genes unique to each stage, overlapping between two consecutive stages or among all three stages (Fig. 1F and G). For example, DTP was characterized by enhanced neural differentiation, DTPP by enhanced interferon signaling, and Rr by reduced melanocyte differentiation. Also, upregulated genes shared by DTPP and R or all three stages enriched for cytokine signaling/migration or extracellular matrix (ECM) remodeling/angiogenesis. In contrast, Ra cells displayed few unique or overlapping differentially expressed genes (DEG; Supplementary Fig. S1I). Thus, BRAFi can induce temporally distinct transcriptomic and phenotypic states leading to MAPK-redundant resistance.

Using gene signatures unique to or shared by DTP, DTPP, and Rr states, we detected recurrent and often concurrent positive enrichments in the On-Tx tumors (Fig. 2A). We assessed the significance of transcriptomic similarity between the On-Tx tumors and Rr lines and found 533 upregulated and 457 downregulated genes being shared between On-Tx tumors (≥50% of biopsies) and DTP, DTPP, Rr cells (>50% of cell samples in any of the states; Fig. 2B; Supplementary Tables S2A and S2B). This overlap of DEGs was significant across all pair-wise comparisons (Supplementary Table S2A). The 533 upregulated genes shared across tumor/cell line samples were enriched for the GO terms “extracellular matrix organization,” “response to IFNA,” “neuron projection,” and “mesenchyme development” (Fig. 2C). Although On-Tx (vs. baseline) tumors harbored higher immune/stromal contents (Supplementary Fig. S2A), the significant overlap of their DEGs with those in MAPKi-treated cell lines suggests that the mesenchymal phenotype could arise, in part, tumor cell autonomously.

We listed the gene sets (C2/6, Molecular Signature Database and Supplementary Table S2C) that were enriched in On-Tx tumors (≥50% of biopsies) and any of the in vitro temporal stages (Fig. 2D). We observed no gene set recurrently enriched (positive or negative) among the Ra transcriptomes. The top gene sets positively enriched were related to cancer-associated fibroblasts (CAF)/mesenchymal transition, ECM, IFN/IRF3, and JNK/TNF signaling. Consistently, network analysis of lineage-specific transcription factors (TF) and target genes (CellNET; ref. 18) showed that Rr lines expressed more fibroblast-specific TFs and target genes (Supplementary Fig. S2B). For recurrently (>50%) differentially expressed mesenchymal, neural crest, and melanocytic genes found during in vitro evolution (Supplementary Table S2C), we observed corresponding DEG patterns in On-Tx tumors (Supplementary Fig. S2C–S2E). The heterogeneity of DEG patterns in neural crest–melanocytic differentiation observed in On-Tx tumors may be related to heterogeneous intratumoral abundance of MITF-high versus MITF-low, AXL/NGFR-high melanoma cells in treatment-naïve tumors (19). We validated BRAFi-induced altered differentiation states based on protein levels of TFAP2A, MITF, FOXD3, CD271/NGFR, ZEB2, SNAI2, and TAGLN/SM22 (Supplementary Fig. S2F–S2I). Consistent with therapy-induced changes in chromatin regulators of melanoma cell lines (10, 16, 20, 21), we observed highly dynamic temporal expression patterns of chromatin regulator genes, including changing ratios of the atypical histone splice variants mH2A1.1 and mH2A1.2, during BRAFi resistance development (Supplementary Fig. S2J–S2N). Mutually exclusive splicing of mH2A1.1 versus mH2A1.2 is known to regulate cancer cell senescence (e.g., heterochromatin, secretory phenotype) and developmental genes (22, 23).

We then sought an epigenomic perspective of the invasive/mesenchymal switch in Rr lines versus On-Tx tumors. We first analyzed published RNA-sequencing (RNA-seq) and H3K27ac-sequencing data derived from proliferative (n = 9) and invasive (n = 2) short-term melanoma cultures (Supplementary Fig. S2O; ref. 24). Rr lines clustered with the invasive group, whereas MAPK-addicted P-lines and Ra lines clustered with the proliferative group. Because CpG methylation levels at H3K27ac peaks were generally anticorrelated with H3K27ac levels (24), we determined the differential CpG methylation levels of Rr lines at the sites of differential H3K27ac specific to invasive versus proliferative cultures. In both On-Tx tumors and Rr lines, we observed general hypermethylation at proliferative phenotype-associated H3K27ac peaks and an admixture of hypomethylation and hypermethylation at invasive phenotype-associated H3K27ac peaks (Fig. 2E and F). Temporally in cell lines (Supplementary Fig. S2P and S2Q), differential CpG methylation at invasive–proliferative switch-associated H3K27ac sites was stage-dependent, occurring first in DTPPs. Furthermore, we analyzed the correlation between differential CpG methylation and expression of proximal genes in Rr lines at the proliferative versus invasive H3K27ac peaks (Supplementary Table S2D). At proliferative H3K27ac peaks, most genes (e.g., MITF, DCT, SOX10, and MLANA) were hypermethylated and down-expressed (Fig. 2G) and enriched for “pigmentation during development” and “melanocyte differentiation.” Most of the genes in invasive H3K27ac peaks were up-expressed and enriched for angiogenic and promotility processes. Interestingly, the up-expressed invasive genes associated with either hypermethylation (e.g., FBN1/2, EGFR) or hypomethylation (e.g., AXL, EPHA2) at CpG sites or both hypermethylation and hypomethylation at distinct CpG sites within the same genes (e.g., LOXL2, WNT5A, VEGFC). Thus, Rr lines and On-Tx residual melanoma tumors underwent an epigenome-based mesenchymal–invasive–angiogenic switch.

To identify gene-level consequences of MAPKi-induced tumor cell–intrinsic transcriptomic reprogramming, we rank-ordered the recurrence of mRNA up-expression and down-expression in On-Tx tumors based on those genes displaying DEG patterns in >50% of any in vitro temporal stage (Fig. 3A; Supplementary Table S2B). DEGs were highly concordant among the DTP–DTPP–Rr stages but discordant between Rr and Ra groups. Importantly, recurrence of DEGs in On-Tx tumors was highly consistent with those observed during in vitro transcriptomic evolution of the Rr phenotype. Among the top 50 recurrently up-expressed genes were those involved in immunity, inflammation, wound healing, and cell motility, as well as growth and differentiation. Many up-expressed genes (e.g., GAS6, LGALS9, IL34, IL32, IL15) constituted a shared secretome response (Supplementary Fig. S3A and Table S3A). PDGFRB up-expression has been shown to be highly recurrent during disease progression (2) and was also highly recurrent early On-Tx, suggesting an adaptive role. Including only RTK genes that met minimal expression cutoff (FPKM ≥ 1) and dynamic range [interquartile range (IQR) ≥ 2-fold], we found that Rr lines formed a distinct cluster from all Ra lines and P lines (Supplementary Fig. S3B). Including only RTK genes meeting a cutoff ≥ 2 fold change (FC) in expression (≥1 pair-wise comparison; IQR ≥ 2), we identified differentially expressed RTKs in MAPK-redundant growth (Supplementary Fig. S3C). RTKs up-expressed by Rr lines and On-Tx tumors were supraphysiologic (i.e., above pretreatment levels) in range (Fig. 3B; Supplementary Fig. S3D and S3E). Among RTK genes recurrently up-expressed during DTP–DTPP–Rr evolution, PDGFRB (3) and other RTK genes (AXL, c-MET, EGFR) have been linked to MAPKi resistance (2, 11, 25, 26). Using validated antibodies (Supplementary Fig. S3F), we detected PDGFRB up-expression as tumor cell–intrinsic and concurrent with EGFR or c-MET up-expression (Fig. 3C).

As in On-Tx tumors, PDGFRB, EGFR, c-MET, and AXL in cell lines displayed supraphysiologic and concurrent up-expression in Rr (but not Ra) lines, along with higher pAKT but lower pERK levels (Supplementary Fig. S4A). We then tested whether concurrent PDGFRB and EGFR supraphysiologic up-expression could result in PDGF-BB–dependent EGFR transactivation. After serum starvation, PDGF-BB treatment strongly induced pPDGFRB, pEGFR, pAKT, and pERK levels in M229 and M238 SDR lines (versus P-lines; Supplementary Fig. S4B). This PDGF-BB→PDGFRB→EGFR signaling axis did not depend on the kinase activity of PDGFRB, as sunitinib (PDGFRBi) treatment inhibited PDGF-BB–induced pPDGFRB levels but not pEGFR accumulation. Interestingly, downstream pAKT and pERK induction was attenuated but not abolished (Supplementary Fig. S4C). However, cotreatment of the Rr lines with gefitinib (EGFRi) on top of PDGF-BB and sunitinib strongly reduced PDGF-BB–induced pAKT and pERK levels and their clonogenic growth (Supplementary Fig. S4D and S4E).

Among immune-regulomic genes (Supplementary Table S3A) with DEGs coherent between On-Tx tumors and DTP–DTPP–Rr temporal evolution were immune checkpoint genes such as PDCD1LG2/PD-L2, CD274/PD-L1, and LAG3. IFN signatures were enriched positively starting at DTPP in vitro, but, unlike in On-Tx tumors, there was no upregulation of IFNA/B/G (Supplementary Fig. S5A). This finding suggested that tumor cell–intrinsic IFN pathway activation in response to MAPKi could remodel the tumor immune microenvironment independently of an adaptive response to antitumor T cells. In On-Tx tumors, T-cell infiltration occurred along with a wholesale immune compartment expansion (Supplementary Fig. S5B and Table S3B), consistent with a tissue injury or wound-healing response, and with a lower ratio of T-cell exhaustion or checkpoint inhibitory genes to CD8A (Supplementary Fig. S5C). Because eventual disease progression is commonly associated with loss of CD8⁺ T-cell inflammation (2), the MAPKi-“responsive” phase may be a window of opportunity to augment antitumor immunity and suppress resistance development.

Among the top recurrently up-expressed immune genes, PD-L2 was expressed at a supraphysiologic range (Fig. 3E and F). Methylation levels (FC or absolute levels) at a specific CpG site in PD-L2 were significantly anticorrelated with mRNA levels in this melanoma tissue and cell line cohort and The Cancer Genome Atlas (TCGA) melanoma (Fig. 3G). The identified PD-L2 CpG site lost methylation with increasing duration of BRAFi selection and PD-L2 mRNA expression (Fig. 3H; Supplementary Fig. S5D and S5E) and located within H3K27ac peaks of short-term invasive melanoma cultures (Fig. 3I). Moreover, PD-L2 (vs. PD-L1) protein up-expression in Rr lines (versus P-lines) was more robust; neither PD-L1 nor PD-L2 was up-expressed in Ra lines (Fig. 3J; Supplementary Fig. S5F–S5H). In paired frozen tumors, we observed strong PD-L2 protein up-expression in residual tumor cells on MAPKi (Fig. 3K). Thus, PD-L2, which has not been shown as a preclinical or clinical target in cancer, may in fact be relevant therapeutically in the context of MAPKi-induced immune phenotypic switch.

Given a reported tumor cell–intrinsic, antiapoptotic role of PD-L1 (27), we first assessed a potential tumor cell–intrinsic, adaptive function of PD-L2′s supraphysiologic expression by testing whether PD-L2 knockdown by shRNAs in Rr lines would suppress clonogenic growth and induce apoptosis (Fig. 4A and B; Supplementary Fig. S6A–S6C). Indeed, partial knockdown of PD-L2 levels in Rr lines by five independent shRNAs consistently and strongly reduced clonogenic growth and induced apoptosis. Interestingly, coculture with HEK293T cells expressing PD-1 attenuated or rescued basal or PD-L2 knockdown–induced apoptosis in Rr lines (Supplementary Fig. S6D; Fig. 4C and D), which was reversed by coexpression of soluble PD-L2 in HEK293T cells (Supplementary Fig. S6E and S6F). Additionally, we induced apoptosis in Rr lines with staurosporine (STP) prior to cocultures (with STP washed away before coculture) and observed that PD-1 expression by neighboring cells protected Rr lines from basal and STP-induced apoptosis (Fig. 4E). This PD-1–mediated protection of apoptosis was not observed in the P-lines or Ra lines, which did not up-express PD-L2 (Supplementary Fig. S6G and S6H). Thus, PD-L2 up-expression promotes survival of MAPK-redundant resistance tumor cell–intrinsically and in a manner augmented by interaction with PD-1.

To gain insights into tumor cell–intrinsic antiapoptotic processes mediated by PD-L2–PD-1 interaction, we sorted two Rr lines (based on negative GFP expression) in the baseline growth condition or after brief STP treatments from cocultures with HEK293T cells, which did or did not overexpress PD-1. PCA of RNA-seq data generated from nonapoptotic, GFP-negative Rr cells showed that STP treatment was associated with a negative shift along the PC2 axis (Fig. 4F), where PC2 positive-driving genes were enriched for cell-cycle progression and PC2 negative-driving genes were enriched for cell death genes (Supplementary Fig. S6I). Although this reflected the effects of STP pretreatment on gene expression, we also observed that, regardless of STP pretreatment, coculture with PD-1–expressing cells invariably resulted in a positive shift along PC3 (Fig. 4F), where PC3 positive-driving genes were enriched for angiogenic processes and PC3 negative-driving genes were enriched for inflammation (Fig. 4G). Thus, binding of PD-1 with PD-L2 expressed by Rr lines may reduce the apoptotic propensity of melanoma cells under chronic MAPKi treatment stress by suppressing tumor cell–intrinsic inflammation, a potential mechanism corroborated by gene set enrichment analysis (Supplementary Fig. S6J) and further by in vivo experiments (see below).

We also assessed gene- and gene signature–level consequences of MAPKi-induced tumor stromal/immune compartment remodeling by comparing clinical On-Tx tumors to immune-competent murine melanoma temporally transcriptome profiled after BRAFi treatment. A syngeneic, transplantable murine melanoma model (YUMM1.7 Braf/Pten/Cdkn2a-mutant cells subcutaneously transplanted into C57BL/6 mice) was used for this comparative transcriptomic analysis, but the clinical relevance of this murine melanoma model was first assessed. YUMM1.7 tumors regressed soon after daily gavage with BRAFi (vemurafenib), and residual tumors persisted before resuming growth (Fig. 5A). Transcriptomes derived from vehicle-treated tumors [day 3 (d 3)] and BRAFi-treated tumors at d 6 (regressing), d 15 (residual), and d 30 (regrowing) were analyzed for expression levels of CD8A/B and CD8⁺ T-cell activation/cytolytic genes (Fig. 5B). This analysis showed heterogeneous but progressive induction of CD8⁺ T-cell activation/cytolytic genes early On-Tx, which peaked in the residual, maximally regressed tumors. Tumor regrowth, however, was associated with loss of CD8⁺ T-cell inflammation. We also observed progressive enrichment of a set of signatures related to innate anti–PD-1 resistance (IPRES; ref. 14), reflecting hypoxia, TGFβ1 signaling, and mesenchymal transition. Enrichment of inflammatory and IFN-induced signatures corresponded to upregulation of CD8A/B, peaked in maximally regressed YUMM1.7 tumors on BRAFi but ceased when the tumors regrew. Thus, YUMM1.7 tumors in response to BRAFi therapy recapitulated the transcriptomic evolution of clinical melanoma.

We first determined in patient-derived On-Tx tumors the differential expression patterns of immune genes recurrently and coherently up-expressed or down-expressed in any stage (regressing, residual, or progressing) of YUMM1.7 tumors on BRAFi (Fig. 5C; Supplementary Table S2B). By rank-ordering up-expressed immune genes in clinical On-Tx tumors based on the frequency of recurrence (Fig. 5C), we could appreciate the immune-regulomic similarity between On-Tx tumors and regressing/residual (but not progressing) YUMM1.7 tumors. To define tumor cell–intrinsic versus tumor cell–extrinsic immune-regulomic alterations, we annotated the recurrence profile (Fig. 5C) with DEG patterns in the CD45-negative fraction of residual YUMM1.7 tumors and in human melanoma cell samples (i.e., DTP, DTPP, Rr, or Ra). We first determined that DEG (in particular gene up-expression) overlapped significantly between MAPKi-selected human cell lines in vitro and the CD45-negative fraction of mouse tumors in vivo (Supplementary Fig. S7A). These overlapping up-expressed genes enriched for GO terms such as ECM organization and cell migration. Additionally, a gene set enrichment analysis that incorporated recurrent alterations in the CD45-negative fraction of residual YUMM1.7 tumors (Supplementary Fig. S7B) revealed concordant and recurrent enrichment of CAF, mesenchymal transition and ECM-remodeling signatures in On-Tx tumors, DTP–DTPP–Rr temporal populations, and CD45-negative residual YUMM1.7 tumors. By excluding the characteristic MAPKi-induced (tumor cell–intrinsic) transcriptomic alterations observed in human melanoma cell lines in vitro and in the CD45-negative fraction of residual YUMM1.7 tumors, we nominated the majority of recurrent gene up-expression in On-Tx tumors (Fig. 5C) as tumor cell–extrinsic alterations. Consistently, at the gene signature enrichment level (Fig. 5D), the top recurrently enriched gene signatures in the On-Tx tumors, with support in any of the temporally sampled YUMM1.7 tumors, were all related to immune cell function, indicating wholesale immune cell infiltration (including CD8⁺ T cells) as during wound healing. Interestingly, the loss of CD8⁺ T cells noted earlier (Fig. 5B) in progressing YUMM1.7 tumors was accompanied by loss of pan-immune subsets (Fig. 5D). Although the top recurrently up-expressed immune genes in On-Tx tumors were not commonly attributable to the tumor cell compartment, the top recurrently up-expressed genes in general (including the secretome; Supplementary Fig. S3A) were frequently attributable to the CD45-negative residual YUMM1.7 tumor cells or DTP–DTPP–R subpopulations in vitro (Supplementary Fig. S7C).

By comparing results in Figs. 3D versus 5C, we noted a few exceptions where recurrent immune gene (e.g., CD36, CXCL12, and CD274/PD-L1) up-expression in On-Tx tumors may be attributable, at least partially, to the tumor cell compartment. We also noticed that, although PDCD1LG2/PD-L2 was up-expressed robustly during the DTP–DTPP–Rr evolution in human melanoma cell lines (Fig. 3), it was not in the CD45-negative fraction of residual YUMM1.7 tumors (Fig. 5C). In fact, the full-length PD-L2 transcript was detected by RNA-seq only with sufficient intratumoral inflammation (e.g., when YUMM1.7 tumors were treated with BRAFi; Supplementary Fig. S7D). As a cell line in response to BRAFi treatment, YUMM1.7 expressed only the 3′ short transcript, which is missing the 5′ exons required to translate the extracellular domain (Supplementary Fig. S7E and S7F). Thus, this model presented an opportunity to assess the impact of PD-L2 induction specifically in the immune compartment. Furthermore, in YUMM1.7 tumors, we confirmed by Ki-67 and CD8 quantitative immunofluorescence highly dynamic alterations in the proliferative status and intratumoral CD8⁺ T-cell abundance (Fig. 5E and F). We also confirmed PD-L2 accumulation and observed its frequent colocalization with CD11c (a myeloid dendritic marker) in the immune compartment of regressing/residual YUMM1.7 tumors (Fig. 5G).

Importantly, cotreatment of BRAFi with a PD-L2–blocking antibody (aPD-L2) consistently delayed growth of resistant tumors, whereas aPD-L2 monotherapy did not lead to tumor growth inhibition or shrinkage (Fig. 6A and B; Supplementary Fig. S7G–S7I). We hypothesized that aPD-L2 acts by delaying the loss of CD8⁺ T-cell inflammation that occurs during the evolution of BRAFi resistance. Consistently, at d33, BRAFi + aPD-L2–responsive tumors, compared with tumors that had escaped either BRAFi alone or BRAFi + aPD-L2, displayed the highest number of CD8⁺ T cells (Fig. 6C). We also profiled via RNA-seq YUMM1.7 tumors that had escaped BRAFi (n = 3) or BRAFi + aPD-L2 (n = 2) or remained responsive to BRAFi + aPD-L2 (n = 2) at d 33 (Fig. 6D) and compared their PC positions relative to those from the BRAFi treatment time course (Fig. 5A). From this transcriptomic perspective, we found that d 33 YUMM1.7 tumors that remained responsive to BRAFi + aPD-L2 segregated with d 15 residual tumors on BRAFi alone, whereas d 33 tumors that escaped BRAFi or BRAFi + aPD-L2 shifted positively along the PC1 axis toward d 30 BRAFi-treated tumors (Fig. 6D). To understand this pattern further, we analyzed the GO term enrichments of PC1⁻ and PC1⁺ driving genes and found that the former genes enriched for T-cell immunity whereas the latter enriched for proliferation (Fig. 6E). Moreover, we harvested YUMM1.7 tumors before the tumor sizes differed significantly between the treatment groups (Supplementary Fig. S7J) and quantified the intratumoral CD8⁺ T-cell numbers. Consistent with the hypothesis that aPD-L2 acts by delaying the loss of CD8⁺ T-cell inflammation, we found that, in the residual tumors, aPD-L2 cotreatment already led to a significantly higher number of CD8⁺ T cells (Fig. 6F). Also, as MAPKi induced infiltration of PD-L2–positive immunocytes into regressing and residual tumors, we expected and observed that delayed aPD-L2 (versus BRAFi) treatments reduced its efficacy (Supplementary Fig. S7K). Moreover, aCD8 neutralization (Supplementary Fig. S7L) reduced the tumor-control efficacy of BRAFi + aPD-L2 but not of BRAFi alone (Supplementary Fig. S7M and S7N), suggesting that in this model CD8⁺ T cells contribute effective antitumor immunity only with PD-L2 blockade in the immune compartment.

Because BRAFi-induced PD-L2 expression in the immune compartment promoted melanoma adaptation in a T cell–dependent fashion, we next examined whether tumor cell–intrinsic PD-L2 expression could promote BRAFi resistance in an in vivo but immune deficient microenvironment [i.e., NOD/SCID gamma (NSG) mice]. Because YUMM1.7 melanoma cells lack full-length Pd-l2 expression, we engineered PD-L2 overexpression and, after subcutaneous transplantation, assessed tumor growth in NSG mice (n = 8 tumors per group). PD-L2 overexpression in YUMM1.7 melanoma cells did not significantly alter tumor growth in mice treated with the vehicle, but, with BRAFi treatment, accelerated resistance development (Fig. 7A). Histologically, BRAFi-treated YUMM1.7 tumors (at d 19; three tumors nearest the median weight were selected) with PD-L2 overexpression (versus no overexpression) displayed more densely packed sheets of spindled malignant-appearing cells (Fig. 7B). Blinded quantification of mitotic figures along with histiocytes and granulocytes revealed significantly higher levels of mitosis but lower levels of infiltration by histiocytes and granulocytes with PD-L2 overexpression (Fig. 7C), consistent with distinct tumor sizes and suggestive of tumor cell cross-talk with innate immune cells. Analysis of the RNA-seq profiles of tumors from all four groups (n = 3 per group) showed highest PD-L2 levels in the two groups with exogenous overexpression. We also observed BRAFi-inducible Pd-l2 up-expression in the control vector group, likely from the limited immune microenvironment. Regardless of BRAFi treatment, we also observed enhanced levels of the immune-suppressive factor VEGFA with PD-L2 overexpression (Fig. 7D), the latter of which was confirmed at the protein level on the cell surface (Fig. 7E). PCA of the RNA-seq data suggested that PC2 positive- and negative-driving genes and their functional pathways as critical to the tumor cell–intrinsic effects of PD-L2 in an immunodeficient in vivo environment (Fig. 7F). In this model, BRAFi negatively shifted the tumor transcriptomes along the PC2 axis, regardless of PD-L2 overexpression status, whereas PD-L2 overexpression positively shifted the tumor transcriptomes along PC2, “reversing” the effect of BRAFi. PC2 positive-driving genes were enriched for hypoxia and other processes, whereas PC2 negative-driving genes were enriched for inflammatory processes (Fig. 7G). Consistent with these results, gene set enrichment analysis showed that PD-L2 overexpression (regardless of BRAFi treatment) led to positive enrichment of hypoxia gene signatures and negative enrichment of multiple IRF-dependent and immunocyte signatures (including those of histiocytes and granulocytes; Fig. 7H). Thus, melanoma tumor cell–intrinsic overexpression of PD-L2 accelerates BRAFi resistance even in an immunodeficient in vivo microenvironment, potentially through upregulating hypoxia-driven oncogenic processes and downregulating inflammatory signaling.

The scale and nature of tumor cell–intrinsic and –extrinsic alterations, during the early course of MAPKi therapy for advanced melanoma, have not been defined systematically. This characterization is essential to start to address the potential adaptive nature of these changes and define therapeutic opportunities. We integrated transcriptomic analysis of patient-matched pretreatment and regressing/residual melanoma with analysis of MAPKi-induced temporal transcriptomic alterations in human melanoma cell lines in vitro and murine melanoma tumors in an immune-competent context. By focusing on recurrent transcriptomic alterations during each stage of response/adaptation in the experimental models, we uncovered similar highly recurrent gene and gene program expression alterations in patient-derived On-Tx tumors. These recurrent transcriptomic alterations originated in the tumor cell and/or immune compartment, as inferred by concordant changes in (i) On-Tx clinical tumors versus human cell lines (i.e., putative tumor cell–intrinsic alterations) or (ii) On-Tx versus mouse tumors (but absent in the CD45-negative fraction of mouse tumors, i.e., putative tumor cell–extrinsic alterations). This approach incorporating temporal dynamics in experimental systems enhanced our view of heterogeneous On-Tx clinical tumors, which were biopsied at varied times and responded to therapy with distinct kinetics. Importantly, this integrative analysis of patient-derived tumors and experimental systems enabled functional interrogation of clinically relevant and highly ranked alterations in patient tumors as potentially adaptive and targetable responses.

This study identified regressing/residual melanoma “responsive” to MAPKi therapy as a nidus of dynamic transcriptomic–CpG methylomic reprogramming and of tumor cell persistence and resistance through attenuated MAPK dependency. MAPK-redundant tumor growth control by multi-RTK upregulation associated with a mesenchymal/invasive/angiogenic switch, highlighting therapeutic challenges. Recent pan-cancer studies (12–14) have linked mesenchymal, angiogenic, and immune-suppressive transcriptomic features, which together are associated with clinical innate (14) and adaptive (28) anti–PD-1 resistance in melanoma. Thus, MAPKi induction of such microenvironmental alterations, including CD8⁺ T-cell infiltration, reflects an immune phenotypic switch that presages loss of T-cell inflammation and likely loss of anti–PD-1 responsiveness. In the immune-competent context, MAPKi should induce release by tumor cells of danger-associated molecular patterns (DAMP), which would promote noninfectious inflammation. Surviving malignant cells appear to undergo mesenchymal/angiogenic/fibroblastic transition and, through an expanded secretome and cell-surface expression of immune-modulatory molecules, actively participate in chronic wound healing or inflammation, which is immune suppressive. Hence, it may be critical to disrupt this immune-phenotypic transition early during “successful” MAPK-targeted therapy in order to prevent the eventual loss of T-cell inflammation.

As a proof of concept, we showed that, by inducing infiltration of PD-L2–expressing immunocytes and tumor cell–intrinsic PD-L2 upregulation, BRAFi unmasks the protumorigenic and resistance functions of PD-L2. These functions of PD-L2 expressed by nontumor cells in the immune microenvironment are at least partially dependent on the ability of PD-L2 to downregulate intratumoral CD8⁺ T-cell function and/or number. Tumor cell–intrinsic upregulation of PD-L2 and its interaction with PD-1 increased fitness of melanoma cells that have epigenomically adapted to MAPKi. Recent work showed that PD-L2 expression, independently of PD-L1 expression, predicted clinical response to pembrolizumab in head and neck squamous cell carcinoma, suggesting that clinical responses to PD-1 blockade may be related partly to blockade of PD-1/PD-L2 interactions (29). It also remains to be determined how much of PD-L2′s tumor cell–intrinsic function is relevant clinically in the context of MAPKi resistance. Further mechanistic understanding and biomarker discovery of PD-L2′s tumor cell–intrinsic function would facilitate this evaluation.

The recognition of tumor cell–intrinsic transcriptomic reprogramming and its potential immune-suppressive, immune-editing or prometastatic impacts early during therapy calls for preemptive, up-front combinatorial therapies. Although the intrinsic requirement of tumor cell mesenchymal transition for metastatic capability has been called into question (30–32), its importance in early adaptive resistance to targeted therapy is strongly suggested by our study. We also observed, in On-Tx clinical tumors, MAPKi inducing tumor cell–intrinsic (in human melanoma cell lines and the CD45-negative fraction of YUMM1.7 tumors) upregulation of multiple factors that could promote metastasis. One example is the fatty-acid receptor CD36, which has been nominated recently (33) as a critical driver of metastasis-initiating cells. Tumor neoantigens can result from somatic mutations or reexpression of cancer testis antigens. MAPKi-induced epigenomic reprogramming raises the possibility of cancer testis antigens as endogenous immune or exogenous immunotherapeutic targets. Furthermore, immune-suppressive candidates include tumor cell–intrinsic elaboration of TGFβ2, VEGFA/C, and IFN upregulation. Interestingly, chronic IFN stimulation of tumor cells has been shown to induce T-cell inhibitory receptors and reduce sensitivity to aPD-1 blockade (34). Also, EGFR, VEGF, and other triggers of STAT3 signaling could promote immunosuppression (35). Loss of PTEN or induction of AKT has been associated with decreased T-cell infiltration into tumors and T cell–mediated tumor killing through expression of immunosuppressive cytokines (36). Importantly, MAPK-redundant resistance that results from transcriptomic reprogramming invariably upregulates AKT signaling. Thus, further analysis of this data set and exploitation of these experimental models should generate additional combinatorial targets with the potential to tighten the bottleneck of melanoma evolution in response to MAPKi therapy.

All melanoma tumors (pretreatment, on-treatment) were obtained with the approval of institutional review boards (IRB) and patients' consent. mRNAs from 24 snap-frozen specimens were extracted by the mirVana RNA Isolation Kit (Life Technologies) and paired-end sequenced with read length of 2 × 100 bps (Illumina HiSeq2000, unstranded TruSeq RNA Sample Preparation Kit v2), except for patient #9 (2 × 150 bps, Illumina HiSeq3000, Kapa Stranded mRNA-seq kit) and patient #5 (for which microarray was used). Also, RNA-seq data of 22 independent specimens (15) were analyzed in parallel. Selected tumor samples were also subjected to whole-exome sequencing (n = 8 pairs) to verify tumor content and to methylome profiling (n = 4 pairs) to compare with methylomic profiles of cell lines.

All cell lines were routinely tested for Mycoplasma, and cell line and subline identities have been ensured by RNA-seq and the GenePrint 10 system (Promega) at routine intervals during the course of this study for banking and experimental studies. All M series cell lines were established from patient-derived tumors at the University of California, Los Angeles. SKMEL28 was obtained from Dr. Alan Houghton (between 2008 and 2010). YUMM1.7 was obtained from Dr. Marcus Bosenberg (between 2014 and 2016). All cell lines were maintained in DMEM high glucose with 10% heat-inactivated FBS (Omega Scientific) and 2mmol/L glutamine in a humidified, 5% CO₂ incubator. To derive DTPP clones, parental melanoma cells seeded at low density were treated with drugs as indicated every 2 to 3 days for 3 to 6 weeks, and highly proliferative colonies were ring-isolated and expanded. SDR and DDR sublines were derived previously as stated in Supplementary Methods. ShPD-L2 and vector control (Thermo Fisher) were packaged into lentiviral particles for infection. PD-1–GFP and its empty lentiviral vector were purchased from Origene. See Supplementary Methods for details on all antibody, small-molecule inhibitor, and growth factor treatments.

Cells were lysed in RIPA buffer (Sigma) with protease inhibitor cocktail (Roche) and phosphatase inhibitor cocktails (Santa Cruz Biotechnology) for Western blotting. Paraffin-embedded formalin-fixed sections were subjected to heat-induced antigen retrieval. The biotin–streptavidin system (Vector Laboratories) and DAB were used to visualize human Ki-67, and Alexa Fluor conjugates (Life Technologies) and the TSA plus fluorescein system (Perkin Elmer) were used to detect human PDGFRβ, EGFR, and c-MET, and mouse Ki-67, PD-L2, CD8, and CD11c. For immunocytochemistry, cells were fixed with 4% paraformaldehyde and stained with Alexa Fluor 555 (for Ki-67) and Alexa Fluor 488 (for CD271) conjugates (Life Technologies). Nuclei were counterstained by hematoxylin or DAPI. Fluorescent images were digitized using the Applied Imaging Ariol SL-50 scanner at 20× magnification. Quantification of CD8 and Ki-67 signals was performed using the Definiens Tissue Studio 4.3 (Definiens, Inc.) according to the manufacturer's instructions. Please refer to Supplementary Methods for primary antibody information. Cell-surface PD-L1 and PD-L2 were detected by flow cytometry using APC-conjugated anti–PD-L1 and PE-conjugated anti–PD-L2 antibodies (BioLegend).

Senescence was assessed by senescence-β-galactosidase staining kit (Cell Signaling Technologies). CSFE (Molecular Probes) dilution was monitored by flow cytometry to follow cell division. Clonogenic and cell-cycle assays are described in Supplementary Methods. For vital imaging, cells were plated onto gridded-dishes (Sigma) and imaged every 2 days at five predesignated areas.

YUMM1.7 cells (1 million per flank, two flanks per mouse) were subcutaneously implanted into C57BL/6 mice. Tumors were treated with BRAFi when the volumes reached to 150–200 mm³. aPD-L2 antibodies and isotype control were administered intraperitoneally. Tumors were dissociated into single cells and stained with a αCD45 antibody. The CD45-negative population was sorted for RNA extraction. See Supplementary Methods for additional details.

FASTQ reads were mapped to the UCSC hg19 or mm9 (for the YUMM1.7 mouse melanoma study) reference genome using Tophat2. Normalized expression levels of genes (and transcripts in the case of mH2A1.1 and mH2A1.2) were expressed in FPKM values as generated by cuffquant and cuffnorm. Because we observed overexpression of the last two exons in the Pdcd1lg2 gene without the accompanying expression of the other exons in the YUMM1.7 tumors treated with BRAFi for 30 days, we excluded the reads mapping to the last two exons when we quantified full-length Pdcd1lg2 mRNA expression in the YUMM1.7 RNA-seq dataset. A gene was defined as differentially expressed when its expression increased or decreased by at least 2-fold for cell lines or 1.5-fold (|log₂ FC| ≥ 0.585) for tumors (including YUMM1.7 tumors). To overcome noise in differential expression values caused by extremely low FPKM levels, we added a pseudo-FPKM value of 0.1 to all expression values. Recurrence analysis of up-expression or down-expression of genes was performed for cell lines/YUMM1.7 tumors/patient tumors with respect to their parental/vehicle-treated/pretreatment lines or tumors, respectively. We used the ESTIMATE algorithm (37) to calculate the extent of immune infiltration of the baseline and On-Tx tumors based on gene expression data. Tumor purities were called using WES analysis by Sequenza (38). The purities of samples without WES data were estimated by the fraction of Sanger-derived BRAF^MUT peaks compared with the wild-type peaks, assuming heterozygous BRAF mutation. Paired gene set enrichment was performed as described (2); single-sample gene set enrichment was quantified using GSVA (39). GO term enrichment was computed based on the online functional annotation tool ENRICHR (40). RNA-seq data have been made available through the Gene Expression Omnibus (GEO) at the accession number GSE75313.

Differentially methylated CpG sites in tumor/cell line samples (profiled by Illumina 450K Methylation arrays) were computed as described (2). We intersected the proliferative and invasive H3K27ac peak regions (24) with the Illumina 450K Methylation probes. For each peak region with more than one CpG probe, we calculated the median of the methylation changes across the CpG probes. Using the independently derived H3K27ac signals and methylation changes from our studies, we computed the Pearson correlation between the fold differences of H3K27ac signals (invasive vs. proliferative cultures) and % methylation changes (in Rr lines or On-Tx tumors vs. their respective P lines or pretreatment tumors) of the CpG sites within the H3K27ac peaks. To analyze differential mRNA expression (between Rr and P lines) within regions of differential H3K27ac and methylation, we selected those H3K27ac peaks overlapping with one or more differentially methylated CpG sites (q-value ≤ 0.05, |Δβ| ≥ 10%) and whose nearest genes were differentially expressed (≥1.5-fold up-expression or down-expression) in more than half of the Rr lines. These genes were grouped into four different categories (hypomethylation or hypermethylation with mRNA up-expression or down-expression).

P.O. Scumpia is a consultant/advisory board member for EMD Serono and Pfizer. D.B. Johnson is a consultant/advisory board member for BMS, Merck, Incyte, and Genoptix. R.S. Lo served on the advisory boards of Shire, Novartis, and Amgen. No potential conflicts of interest were disclosed by the other authors.

Conception and design: C. Song, A. Hong, W. Hugo, R.S. Lo

Development of methodology: C. Song, M. Piva, L. Sun, A. Hong, S. Lomeli, J. Qian, C.C. Yu, W. Hugo, R.S. Lo

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): C. Song, M. Piva, L. Sun, G. Moriceau, S. Lomeli, J. Qian, C.C. Yu, R. Damoiseaux, M.C. Kelley, K.B. Dahlman, J.A. Sosman, D.B. Johnson, A. Ribas, W. Hugo, R.S. Lo

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): C. Song, M. Piva, L. Sun, A. Hong, X. Kong, S. Lomeli, J. Qian, C.C. Yu, P.O. Scumpia, A. Ribas, W. Hugo, R.S. Lo

Writing, review, and/or revision of the manuscript: C. Song, M. Piva, L. Sun, J.A. Sosman, D.B. Johnson, W. Hugo, R.S. Lo

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): L. Sun, H. Zhang, R. Damoiseaux, M.C. Kelley, W. Hugo, R.S. Lo

Study supervision: W. Hugo, R.S. Lo

We thank G. Bollag (Plexxikon Inc.) for providing PLX4032/4720, M.W. Bosenberg for YUMM1.7, and H. Shi for technical assistance. We thank all patients who donated tissues for this research. Patient-informed consent was obtained for the research performed in this study.

This work has been funded by the Burroughs Wellcome Fund (to R.S. Lo), the NIH (1R01CA176111 to R.S. Lo; P01CA168585 to A. Ribas and R.S. Lo; K12CA0906525 to D.B. Johnson), the Ressler Family Foundation (to R.S. Lo and A. Ribas), the Melanoma Research Alliance (to R.S. Lo, D.B. Johnson, and W. Hugo), the Ian Copeland Melanoma Fund (to R.S. Lo), the SWOG/Hope Foundation (to R.S. Lo and A. Ribas), the Steven C. Gordon Family Foundation (to R.S. Lo), the American Skin Association (W. Hugo), the American Association for Cancer Research-Amgen, Inc. Fellowship in Clinical/Translational Cancer Research (16-40-11-HUGO to W. Hugo), the Department of Defense Horizon Award (to A. Hong), the Dermatology Foundation (to G. Moriceau), a National Cancer Center Aggressive Cancer Research Postdoctoral Fellowship (to J. Qian), the ASCO Conquer Cancer Career Development Award (to D.B. Johnson), and the American Cancer Society Research Professorship (to J.A. Sosman).