A core GRN explains the patterns of non-genetic heterogeneity along the proliferative-invasive axis in melanoma cells

Phenotypic plasticity in cancer cells is driven by a myriad of complex molecular interactions that drive coordinated changes in molecular and functional aspects of various biological pathways. Identifying a minimal regulatory network that can represent many such experimentally reported interactions and whose emergent dynamics can recapitulate the hallmark patterns of phenotypic plasticity and heterogeneity is, therefore, a first important step to elucidate systems-level behavior, while effectively approximating the biological system at hand. To examine the plasticity of melanoma cells along the proliferative-invasive axis, we identified a regulatory network that involves transcription factors implicated in a PIT—MITF, JUN, ZEB1, SOX10, and SOX9. MITF can drive the proliferative phenotype at least partly through LEF1, and transactivate its own promoter.24 25 It can also suppress JUN, a regulator of the invasive cell state, through binding to its enhancer.26 Furthermore, JUN can suppress the activity of MITF directly and/or indirectly,26 and self-activate transcriptionally through forming Jun/ATF-227 heterodimers. Thus, MITF and JUN are engaged in mutual antagonism. Similarly, MITF and ZEB1 can form a ‘toggle switch’,28 such that MITF directly represses ZEB1,29 and ZEB1 can bind to the MITF promoter region and inhibit MITF.30 Similar to MITF and JUN, ZEB1 can autoactivate via multiple feedback loops.31 SOX10 has been reported in maintaining a proliferative state,32 but it has also been shown to promote invasive features in melanoma.33 34 SOX10 inactivation leads to SOX9 upregulation, and SOX9 can downregulate SOX10 expression by binding to SOX10 promoter, thus indicating mutual antagonism between them.32 SOX9 can also be involved in upstream enhancer mediated positive autoregulation.35 Finally, SOX10 has been shown to activate the expression of MITF by directly binding to its promoter.36 37 Together, these interactions constitute a GRN incorporating the complex crosstalk among diverse players of PIT or its reverse IPT (figure 1A).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1

Non-genetic heterogeneity along the proliferative-invasive spectrum enabled by an underlying gene regulatory network in melanoma. (A) Gene regulatory network to model the proliferative-invasive heterogeneity in cell states in melanoma. The red hammerheads represent inhibitory links, and the black arrows represent activating connections. (B) Hierarchically clustered heatmap of simulated steady states allowed by the gene regulatory network and qualitative classification of the four emergent cell states. (C) Scatter plot showing all the steady states projected onto the proliferative-invasive plane defined by sum of expression of MITF, SOX10 (proliferative score) and those of SOX9, ZEB1 and JUN (invasive score). The steady states are colored by cluster labels obtained from hierarchical clustering. (D) Scatter plot of the association between the proliferative (MITF+SOX10) and invasive (SOX9+ZEB1+JUN) scores for the CCLE group of skin cancer cell lines (top). Scatter plot showing association between the proliferative (MITF+SOX10) and invasive (SOX9+ZEB1+JUN) scores for clinical samples from TCGA cohort of SKCM patients (bottom). (E) Volcano plot for the results of meta-analysis of melanoma datasets, accounting for the associations between the proliferative and invasive scores; n_negative and n_positive denote the number of datasets (out of 32) that are correlated negatively (Spearman correlation coefficient <−0.3; p<0.05) and positively (Spearman correlation coefficient >0.3; p<0.05) (left). Scatter plot for the comparisons of the proliferative-invasive correlations based on the five gene signature and the gold standard Verfaillie signatures. Points colored in red are datasets that show significant correlations in both metrics (right). (F) (i) Scatterplot of single cell RNA-Seq data showing each cell of each cell line projected on the proliferative-invasive plane based on the imputed expression of five modeled master regulators. For each cell line, the dominant phenotype is binarized and the corresponding abundances have been reported. (ii) Boxplots of cells categorized by the dominant binary phenotypes based on the five gene signature and their invasive (left) and proliferative (right) ssGSEA scores based on Verfaillie signatures. CCLE, cancer cell line encyclopedia; TCGA, The Cancer Genome Atlas.

We simulated the dynamics of this GRN using the random circuit perturbation (RACIPE) framework to identify the steady state expression patterns that represent the ‘possibility space’ of phenotypes enabled by this GRN. The RACIPE framework simulates the dynamics of a GRN through a set of coupled ordinary differential equations (ODEs), with each ODE for a node in the GRN, and kinetic parameters corresponding to production, degradation and regulation being sampled from biologically relevant ranges.38 The ODEs are then solved across diverse parameter sets and initial conditions to identify the ensemble of steady states. This ensemble is indicative of possible phenotypes allowed by the corresponding network topology. For the GRN considered here, we first plotted the ensemble of steady state solutions as a heatmap (figure 1B). On hierarchical clustering, we observed that the invasive genes JUN, ZEB1 and SOX9 were largely coexpressed and so were the proliferative genes MITF and SOX10. Qualitatively speaking, this heatmap revealed four phenotypes—proliferative (MITF^hi/SOX10^hi), undifferentiated (SOX9^hi with JUN and/or ZEB1 high), neural crest-like (SOX10^hi with JUN and/or ZEB1 high) and intermediate (MITF^hi/SOX9^hi) (online supplemental figure S1A). Coexpression of SOX10 along with invasive markers, such as ZEB1 and/or JUN, has been shown to be necessary to maintain the neural crest-like phenotype with invasive characteristics in development.39 40 To classify the ensemble of steady states into cellular phenotypes more systematically, we defined two scores: proliferative score (sum of z-normalized gene expression values of proliferative master regulators—MITF and SOX10) and invasive score (sum of z-normalized gene expression values of invasive marker genes—SOX9, JUN and ZEB1). We then projected all the steady state solutions onto a scatter plot of corresponding proliferative and invasive scores (figure 1C). The four phenotypes observed in the heatmap can also be recapitulated as distinct clusters in this proliferative-invasive plane, with the proliferative and invasive scores showing a strong negative correlation (ρ=−0.837; p<10⁻⁵). Similarly, on coloring the steady states by a ‘neural crest-like score’ (z-normalized expression values of SOX10+JUN+ZEB1), we observed the neural crest-like phenotype present amidst the proliferative-invasive spectrum (online supplemental figure S1B). These phenotypes identified by our dynamical model are in accordance with previously reported multistage differentiation and dedifferentiation states, which define melanoma subtypes with the neural crest-like phenotype as an intermediate cell state between the terminal melanocytic proliferative state and the undifferentiated invasive cell state.41

To cross-validate the observed trend from model simulations, we analyzed a cohort of skin cancer cell lines included in the cancer cell line encyclopedia (CCLE)42 (figure 1D; top). In support of the model simulations, we found that proliferative and invasive scores were strongly negatively correlated in the CCLE skin cancer cell line cohort (ρ=−0.558; p<10⁻⁵). This pattern is also reiterated in the primary tumor cohort from The Cancer Genome Atlas (TCGA) (ρ=−0.427; p<10⁻⁵) (figure 1D; bottom). Furthermore, in both CCLE and TCGA samples, we observed that the proliferative score defined based on MITF and SOX10 levels correlated strongly positively with the enrichment of a much larger gene set defined for proliferative status in melanoma.43 Similarly, the invasive score derived from the model (defined as the sum of expression of SOX9, JUN and ZEB1) correlated strongly with the enrichment of the Verfaillie gene set used to define invasive status in melanoma (online supplemental file 1). These observations strongly suggest that the five genes involved in the GRN model seem sufficient to capture the major patterns in the proliferative and invasive nature of melanoma samples, both for cell lines and patient tumors. As expected, the proliferative and invasive scores defined based on Verfaillie gene expression signatures negatively correlated with each other in both CCLE and TCGA (online supplemental figure S1E), with the extent of antagonism similar to that seen using our five gene signature (figure 1D).

To ensure the reliability of the chosen five genes for our regulatory network, we compared the correlation of the five identified master regulators (MITF, SOX10 for proliferative; JUN, ZEB1, SOX9 for invasive) against a distribution of correlation coefficients of random combinations of transcription factors (negative control) with the Verfaillie proliferative and invasive signatures in both CCLE skin cancer cell line group and TCGA SKCM patient cohort (online supplemental figure S2). We observed that the MITF-SOX10 pair and the JUN-ZEB1-SOX9 triplet are generally more strongly correlated to the Verfaillie proliferative and invasive signatures, respectively, compared with any random combination of two or three transcription factors respectively in the CCLE and TCGA datasets (online supplemental figure S2A,B). This suggests that these sets of transcription factors are key drivers of proliferative-invasive phenotypic regulation in melanoma. Furthermore, we found that the trend remains consistent for the MITF-SOX10 pair even when the background reference distribution is created by sampling pairs of transcription factors exclusively from the Verfaillie proliferative gene list (online supplemental figure S2C,D). However, a similar analysis for the triplet of invasive genes did not show very strong trends and specificity in CCLE and TCGA datasets (online supplemental figure S2C,D). To improve the specificity of our analysis, we added TCF4 and IRF1 to the signature, as they were among the top correlated transcription factors, which significantly improved the observed trends (online supplemental figure S2E).

To provide further validation of our observations in the CCLE and TCGA cohorts, we analyzed 32 publicly available melanoma bulk transcriptomic datasets (online supplemental table S1) and observed a significant negative correlation (ρ<−0.3, p<0.05) between the defined proliferative and invasive score as the predominant trend (figure 1E; left). In 27 out of the 32 datasets, our proliferative score defined based on our minimalistic GRN (MITF+SOX10) correlated positively with the Verfaillie proliferative score (online supplemental figure S3A). Similarly, in 16 out of 32 datasets, our invasive score (=SOX9+ JUN+ ZEB1) correlated positively with Verfaillie invasive score (online supplemental figure S3B). These results indicate that our minimal GRN is sufficient to capture the proliferative-invasive antagonism, observed across many datasets when compared against the well-established Verfaillie proliferative and invasive signatures (figure 1E, online supplemental figure S3C). Finally, we estimated the strength of antagonism between proliferative and invasive scores vis-à-vis that between the Verfaillie proliferative and invasive gene signatures.43 We found that the anticorrelation between the two scores defined based on our GRN is largely consistent (in 16 out of 32 datasets) with the Verfaillie scores defined proliferative-invasive spectrum (figure 1E; right). The antagonism between the proliferative and invasive scores, as well as the positive correlation between the invasive score and Verfaillie invasive signature in both TCGA and CCLE datasets, were also seen when using the updated invasive score (including IRF1, TCF4, SOX9, JUN, ZEB1) (online supplemental figure S3D,E). The updated invasive score also showed similar trends in meta-analysis of melanoma datasets as noted earlier (compare online supplemental figure S3F, i with figure 1E). Additionally, 25 out of 32 datasets showed a positive correlation (ρ>0.3, p<0.05) between the updated five gene invasive score and the Verfaillie invasive score (online supplemental figure S3F ii). Thus, our model observations can recapitulate the patterns of phenotypic heterogeneity across many publicly available melanoma datasets, including in vitro, in vivo and patient sample data sets.

To analyze phenotypic heterogeneity at a single-cell level in melanoma with respect to the five genes in the core GRN, we analyzed single-cell transcriptomic profiles of patient-derived melanoma cultures (GSE134432) along the proliferative-invasive axis.9 We projected 39,263 cells belonging to 10 genetically homogeneous patient-derived MM lines on the proliferative-invasive plane (figure 1F,i). A density histogram for the expression of the five core genes, fitted with a kernel density estimate, revealed a largely bimodal distribution (online supplemental figure S4A) just as in our simulations (online supplemental figure S1A). On binarization of gene expression levels in each melanoma cell (see the Methods section), we quantified the gene state as high (1) or low (0). We noticed that a large majority (>92%) of the cells for specific cell lines can be categorized into either of the four observed phenotypes that is, the phenotypes of the cells were found to be consistent within their cell line. The cells belonging to the cell lines—MM001, MM011, MM031, MM074 and MM087—can be assigned to a proliferative subtype with high levels of MITF and SOX10 and low levels of SOX9, JUN and ZEB1 ((MITF, SOX10, SOX9, JUN, ZEB1)=11000). In contrast, MM047, MM029 and MM099 cell lines are typically invasive, with low levels of MITF and SOX10 and high levels of SOX9, JUN and ZEB1 high ((MITF, SOX10, SOX9, JUN, ZEB1)=00111). Interestingly, A375 can be classified as a neural crest like cell type with high levels of SOX10 and ZEB1 and low levels of MITF, SOX9 and JUN ((MITF, SOX10, SOX9, JUN, ZEB1)=01 001). Moreover, the cell line MM057 belonged to an intermediate phenotype with a subset of both proliferative and invasive factors being relatively high. This cell state can also be found as a subset of the proliferative phenotype in our simulation results (figure 1B), further underscoring the applicability of a parameter ensemble modeling approach to capture the ‘possible cell states’ allowed by a GRN. Similar patterns were revealed when cell lines were projected in a two-dimensional proliferative and invasive score plane (online supplemental figure S4B). After the classification of these cell-lines into various phenotypes along the proliferative-invasive spectrum based on the five GRN model, we tested our classification against the activity-based score of Verfaillie proliferative and invasive signatures (figure 1F, ii). The phenotype corresponding to a high expression of MITF and SOX10 and low expression of SOX9, JUN and ZEB1 (11000) had the highest Verfaillie proliferative score and the least Verfaillie invasive score. Conversely, the 00111 phenotype where the invasive players are more highly expressed compared with the proliferative ones had low Verfaillie proliferative scores and high Verfaillie invasive scores. The other phenotypes where one or more of both proliferative and invasive genes are high were ranked as intermediate with respect to the extreme 00111 and 11000 phenotypes, marking the ends of the Verfaillie proliferative-invasive axis. Similarly, the cells we categorized as a neural crest-like phenotype based on MITF, SOX10, SOX9, JUN, ZEB1 ((MITF, SOX10, SOX9, JUN, ZEB1)=01 001) had enrichment of activity of the neural crest associated gene signature (online supplemental figure S4C).

Overall, we demonstrate that emergent dynamics of a minimal five-gene-based GRN can explain the hallmark patterns underpinning phenotypic heterogeneity in melanoma cell lines and tumors, both at the bulk and single-cell transcriptomic levels, along the proliferative-invasive spectrum.

Perturbations to the GRN can cause non-genetic phenotypic changes along the proliferative-invasive spectrum during melanoma disease progression

To understand how perturbing our core GRN can impact proliferative-invasive characteristics, we simulated two scenarios: (A) downregulation (DE) of SOX10 (figure 2A) and (B) OE of SOX9 (figure 2E). DE of SOX10 led to an increase in the frequency of the invasive phenotype and a decrease in the proliferative phenotype, when compared with the unperturbed wild-type network. We also found a decreased frequency of the neural-crest phenotype and an increased proportion of the intermediate phenotype (figure 2A). The model simulations are substantiated by in vitro bulk transcriptomic analysis of the effects of SOX10 shRNA-mediated knockdown in M010817 human melanoma cells (GSE37059).32 SOX10 knockdown led to a decrease in levels of the MITF and SOX10-based proliferative score, indicative of an enrichment of the invasive state (figure 2B, left). Simultaneously, these cells have a higher invasive score, as defined by expression levels of SOX9, JUN and ZEB1 (figure 2B, right) as well as with an invasive score defined by the expression of SOX9, JUN, ZEB1, IRF1 and TCF4 (online supplemental figure S5A). Furthermore, cells with SOX10 knockdown displayed perturbation of the regulon-based proliferative and invasive scores as compared with control cells (online supplemental figure S5A). We next interrogated the impacts of SOX10 knockdown in the A375 neural crest-like melanoma cell line (GSE180568).44 Similar to M010817 melanoma cells (figure 2B), SOX10 knockdown in A375 cells led to an increase in the invasive score and a slight decrease in the proliferative scores (figure 2C, online supplemental figure S5B) (GSE180568). We also found that the increased invasiveness was accompanied by a loss of the neural crest-like signature, indicative of a transition from a neural crest-like phenotype to an undifferentiated phenotype, as predicted by the in silico perturbation simulations (figure 2A). To assess the dynamic trajectory of transition from a more proliferative to a more invasive phenotype, we next analyzed single cell RNA-seq data of SOX10 knock down time-course experiments of three available cell lines (MM057, MM074 and MM087), all of which exhibit a proliferative phenotype (figure 1F,i). On projecting the MM057 control and SOX10-siRNA treated cells (GSE134432) on the proliferative/invasive plane, we observed a gradual shift of cells from proliferative to invasive. By 72 hours of SOX10 knockdown MM057 cells were as invasive as the canonical invasive cell lines (MM047, MM099 or MM029), with low levels of MITF and SOX10 and high levels of SOX9, JUN and ZEB1 high (00111) (compare figure 2D with figure 1F,i). Similar trajectories and shifts in phenotypic distributions to varying degrees were observed for MM074 and MM087 cell lines on knock down of SOX10 (online supplemental figure S5C).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2

In silico and experimental perturbations of melanoma cell lines reveal shifts along the proliferative-invasive axis. (A) Bar plots showing proportions of steady state solutions in control cases compared with SOX10 downregulation across the four phenotypes. Error bars represent standrd deviation (SD) of n=3 technical RACIPE replicates. (B) Bar plots of the experimentally observed significant changes (demarcated by *) in the proliferative (defined as the sum of z-normalized expression of MITF and SOX10) (left) and invasive scores (defined as the sum of z-normalized expression of JUN, ZEB1 and SOX9) (right) on SOX10 knockdown (GSE37059). (C) Bar plots showing experimentally observed significant changes in the Invasive score (defined as the sum of z-normalized expression of JUN, ZEB1 and SOX9) (left) and neural crest-like ssGSEA score (Tsoi et al) (right) (GSE180568). (D) Scatter plot of single cell RNA-Seq data for the transition of SOX10 knockdown cells in MM057 cells at 24 hours, 48 hours and 72 hours in comparison to control data (GSE134432). (E) Bar plots showing proportions of steady state solutions in control cases compared with SXO9 over expression across the four phenotypes. Error bars represent SD of n=3 technical RACIPE replicates. (F) Bar plots of the experimentally observed changes in the proliferative (defined as the sum of z-normalized ssGSEA scores of MITF and SOX10 regulons) (left) and invasive scores (defined as the sum of z-normalized ssGSEA scores of JUN, ZEB1 and SOX9 regulons) on SOX9 over expression (right) (GSE57463). *Represents a statistically significant difference in the levels based on Student’s t-test. RACIPE, random circuit perturbation.

Next, we simulated SOX9 OE, which led to an increased frequency of the invasive phenotype (figure 2E). Concurrently, both the proliferative and neural crest-like phenotypes experienced a decrease in their percentage distribution, with a mild increase in that of the intermediate phenotype. To validate this model prediction, we analyzed publicly-available bulk transcriptomic profiles of M010817 human melanoma cells in which SOX9 has been overexpressed (GSE57463).45 Although we did not observe a marked change in the proliferative score based on the expression of MITF and SOX10 alone (online supplemental figure S5D), we found a decrease in the proliferative nature (figure 2D, left), defined by the activity score of the MITF and SOX10 regulon and an increase in invasiveness (figure 2D, right), as defined by the activity score of SOX9, JUN, and ZEB1. The invasive score (JUN, SOX9 and ZEB1 expression) showed a significant increase as reflected in the activity scores of the regulons of the corresponding genes (online supplemental figure S5D), supporting our model prediction. A similar trend of increasing in invasive score on SOX9 OE was also recapitulated using the updated invasive score (online supplemental figure S5D). Together, these analyses reveal that perturbations in levels of various nodes in our minimal GRN is capable of recapitulating phenotypic switching along the proliferative-invasive axis in melanoma.

Phenotypic plasticity along the proliferative-invasive axis can alter PD-L1 levels in melanoma cells

After validating the ability of our GRN to predict phenotypic plasticity in melanoma cells, we sought to investigate interconnection between phenotypic plasticity along the proliferative-invasive spectrum and the immunosuppressive traits of melanoma cells. To understand immune evasion in melanoma, we coupled immune-effector molecules: PD-L1 (suppressor on adaptive immune system16) and IFNγ signaling (figure 3A). We extended our GRN based on following experimental observations: (1) IFNγ is a key regulator of PD-L1 levels,21 (2) IFNγ regulates differentiation state of melanoma cells,21 46 47 (3) IFNγ can transcriptionally activate PD-L1 via NF-κB and/or JAK-STAT-IRF1 pathway,48 49 and (4) IFNγ can suppress MITF expression by inducing CBR and STAT1 association, thereby inhibiting CREB binding to MITF promoter.50 Likewise, PD-L1 is connected to the GRN through the following experimentally observed connections: (1) MITF can degrade PD-L1 via SA-49,51 (2) JUN activates the expression of PD-L1 through STAT3,52 and (3) and ZEB1 can upregulate PD-L1 expression53 and promote immune escape in melanoma.54 Furthermore, we created a variant version of the regulatory network including IRF1 (online supplemental figure S6A). IRF1 is known to be a potent activator of PD-L155 and is also known to be upregulated by IFNγ signaling.56 SOX10 can also down regulate IRF1 in the context of melanoma cells.57 Thus, PD-L1 and IFNγ signaling crosstalk extensively with multiple drivers of proliferative and invasive phenotypes in melanoma.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3

Steady state dynamics of emergent cellular phenotypes on addition of the PD-L1-IFNγ signaling motif to the gene regulatory network. (A) Gene regulatory network coupling the proliferative-invasive regulators with the PD-L1-IFNγ signaling motif in melanoma. The red hammerheads represent inhibitory links, and the black arrows represent activations. (B) Scatter plot showing all the steady states projected onto the proliferative-invasive plane defined by expression of MITF, SOX10 (proliferative score) and SOX9, ZEB1 and JUN (invasive score). The steady states are colored by cluster labels obtained from hierarchical clustering; (C) (i) Scatter plot showing all the steady states projected onto the proliferative-invasive plane colored based on PD-L1 levels. (ii) Strip plot showing the PD-L1 steady state levels for the four phenotypes. The horizontal lines mark the stratification of PD-L1 levels into low, mid and high regions. (D) (i) Scatterplot of single cell RNA-Seq data showing each cell of each cell line projected on the invasive-proliferative score and the imputed PD-L1 expression plane. (ii) Boxplots of cells categorized by the dominant binary phenotypes based on the five gene signature and their imputed PD-L1 levels. (E) (i) Scatter plot showing increase in PD-L1 levels of SOX10 knockdown cells in MM074 cell line as they transition from a proliferative phenotype to an invasive phenotype along 24 hours, 48 hours and 72 hours time course single cell RNA-seq data in comparison to control data (GSE134432). (ii) Boxplots of cells categorized by different time points quantifying their imputed PD-L1 levels. (F) (i) Scatterplot of single cell RNA-seq data projecting cells of two cell lines—A375 (red and orange corresponding to the resistant and sensitive clones, respectively) and 451Lu (green and blue corresponding to the resistant and sensitive clones, respectively) on the proliferative-invasive plane and (ii) Scatterplot of single cell RNA-seq data projecting Vemurafenib sensitive and resistant A375 clones with the difference of invasive and proliferative scores on the x-axis and the PD-L1 levels on the y-axis. (G) Barplot showing the proliferative scores (left), invasive scores (middle) and PD-L1 levels (right) across control cells, relapsed cells, and cells retransplanted into mice after relapse (GSE40213). *Represents a statistically significant difference in the levels based on Student’s t-test.

As before, we used RACIPE to generate the ensemble of steady-state solutions, enabled by this extended GRN. We plotted all the resultant steady state values on the proliferative-invasive plane after performing hierarchical clustering (figure 3B). We observed similar phenotypic clusters as observed for the earlier GRN (figure 1C), demonstrating a strong proliferative-invasive antagonism (ρ=−0.817; p<10⁻⁵). On overlaying PD-L1 levels on the proliferative-invasive plane, we noted that undifferentiated/invasive and neural crest subtypes tended to have high PD-L1 expression whereas intermediate and proliferative cell types had relatively low levels of PD-L1 (figure 3C,i). To quantify the association of PD-L1 for these phenotypes, we first evaluated the histogram of PD-L1 distribution, which showed distinct trimodally (online supplemental figure S6B). Further analysis revealed that neural crest and undifferentiated states express high to intermediate levels of PD-L1, in contrast to proliferative and intermediate phenotypes, for which PD-L1 expression is primarily in the low range (figure 3C,ii). The variant network including IRF1 (online supplemental figure S6A) also showed very similar phenotypic clusters (online supplemental figure S6C) as observed for the network shown in figure 3A.

To validate these model predictions for association of proliferative-invasive status with PD-L1 levels, we investigated the CCLE group of skin cancer cell lines. We used two different metrics to quantify the proliferative-invasive status: (A) the five genes included in the GRN (=SOX9+JUN+ZEB1–SOX10–MITF), and b) the Verfaillie gene lists (= ssGSEA(Verfaillie invasive geneset)–ssGSEA (Verfaillie proliferative geneset). In CCLE, we observed weak positive associations of PD-L1 expression with both of these metrics (online supplemental figure S6D). Similar observations were seen for TCGA cohort of SKCM patients (online supplemental figure S6E). These trends remained consistent even with the new metric that included updated invasive score including IRF1 and TCF4 (online supplemental figure S6F). Further, single cell RNA-seq analysis of melanoma cell lines revealed that more invasive cell lines were more likely to express higher levels of PD-L1 compared with the more proliferative cell lines (figure 3D, online supplemental figure S7A). Further, among invasive cell lines, we observed PD-L1 expression levels to be quite heterogeneous with cells exhibiting low, medium, and high levels (figure 3D,i). This pattern is reminiscence of our model simulations where the neural-crest like and invasive phenotypes expressed low, intermediate as well as high levels of PD-L1, while the proliferative phenotype had predominantly low PD-L1 levels (figure 3C,ii).

Next, we examined if induction of an invasive state via perturbations to the GRN caused concomitant changes in PD-L1 expression levels. First, we analyzed single cell RNA-seq datasets for siRNA-mediated knockdown of SOX10 in MM057, MM074 and MM087 melanoma lines (figure 3E). PD-L1 levels significantly correlated to the proliferative-invasive status of MM074 cells (figure 3E,i), where PD-L1 expression levels were markedly elevated 72 hours after SOX10 knockdown as compared with the non-transfected controls (figure 3E,ii). Similar trends, although weaker, were noted in MM057 and MM087 cells (online supplemental figure S7B), suggesting a ‘semi-independent’ behavior between PD-L1 expression and the proliferative-invasive status. This pattern denotes that while PIT alters PD-L1 expression, these levels can be modified by other signaling pathways.

PIT in melanoma cells is also often associated with resistance to targeted therapies. More invasive cells are more resistant to targeted therapies such as BRAF/MEK inhibitors.13 58 Based on this association, we hypothesized that targeted therapy may lead to changes in the GRN that would impact levels of PD-L1. To provide support for this hypothesis, we analyzed single cell RNA-seq data from paired parental and vemurafenib-resistant clones of A375 and 451Lu melanoma cells (GSE108383) 59(figure 3F,i–ii). On projecting the two cell lines onto a proliferative-invasive plane based on the five genes in the GRN, we observed that the A375 cell line was comparatively more invasive in nature compared with the 451Lu cell line. The parental 451Lu cell line was predominantly proliferative based on the five gene signature ((MITF, SOX10, SOX9, JUN, ZEB1)=11000) while the A375 cell line was predominantly neural crest-like ((MITF, SOX10, SOX9, JUN, ZEB1)=01 001). The vemurafenib-resistant A375 clone was distinctively more invasive compared with its parental counterpart, with decreased SOX10 and increased ZEB1 levels, suggesting that emergence of adaptive resistance to vemurafenib may involve a transition from a neural crest to invasive phenotype in melanoma (figure 3F,i, online supplemental figure S7C). Conversely, vemurafenib-resistant 451Lu cells did not show such a distinctive phenotypic switch, as the parental and the resistant clones were both proliferative based on the five gene signature. The phenotypic switch from neural crest-like to invasive in A375 cells was accompanied by an increase in PD-L1 levels (figure 3F,ii). No such change was observed in 451Lu cells that did not undergo a switch in phenotypes during the evolution of resistance to vemurafenib (online supplemental figure S7D). These trends reinforce our model predictions that a switch along the proliferative-invasive spectrum can contribute to increased PD-L1 expression.

We next analyzed transcriptomic profiles for control, relapsed and retransplanted relapsed tumors for the HCmel3 mouse melanoma model (GSE40213)8 (figure 3G). We found that the control and the retransplanted tumor samples were more proliferative and less invasive compared with the relapsed melanoma tumors. The levels of PD-L1 levels mirrored the invasive nature of the tumor samples. More importantly, it shows that the increase in PD-L1 levels was reversible as melanoma cells underwent IPT (figure 3G). Inflammation-induced reversible dedifferentiation has been reported to cause melanoma cells to resist T-cell therapy.8 Overall, our results show that PD-L1 expression in melanoma cells is reversibly associated, at least in part, with the gain of a more invasive phenotype, thereby providing further support for the role of phenotypic plasticity along the proliferative-invasive axis in mediating upregulation of PD-L1.

Role of IFNγ signaling in induction of PD-L1 levels and dedifferentiation of melanoma cells

A positive association of PD-L1 and IFNγ has been extensively reported experimentally.60–62 We tested whether our model could capture this association. Our simulation results suggest that phenotypes with high or intermediate IFNγ signaling are likely to have high PD-L1 expression, thus highlighting a strong positive influence of IFNγ on PD-L1 (ρ=0.75; p<10⁻⁵) (figure 4A). To get further supporting evidence for this trend, we analyzed bulk transcriptomic data of melanoma cell lines treated with IFNγ. We found that on treatment, both the Hallmark IFNγ response pathway activity and PD-L1 expression was enriched, indicating an association between IFNγ signaling and PD-L1 levels (GSE154496)63 (figure 4B). As expected, the PD-L1 expression levels were correlated positively with IFNγ activity levels both in CCLE and TCGA samples (figure 4C, online supplemental figure S8A).

• Download figure
• Open in new tab
• Download powerpoint

Figure 4

Control of PD-L1 levels and the differentiation state of melanoma cells by IFNγ signaling. (A) Scatterplot showing the association between PD-L1 levels and IFNγ signaling at steady state levels. (B) Paired plot showing changes in levels of (i) Hallmark IFNγ signaling and (ii) PD-L1 levels on treatment of melanoma cells with IFNγ. *Represents a statistically significant difference in the levels based on a paired Student’s t-test. (C) Scatter plots showing association between Hallmark IFNγ response (x-axis) and PD-L1 expression (y-axis) in the TCGA melanoma cohort. (D) (top) Scatterplot showing the spread of steady state solutions with invasive—proliferative score on the x-axis and IFNγ signaling on the y-axis. The steady states have been colored by the corresponding four discretized phenotypes. The horizontal line represents the stratification between low and high IFNγ signaling regions. (Bottom) Bar plots showing the fraction of steady state solutions that belong to each of the four phenotypes segregated based on the level of IFNγ signaling. The error bars represent the SD based on n=3 RACIPE replicates. *Represents a statistically-significant difference in the levels based on a Student’s t-test. (E) Scatterplot showing associations between the (i) Verfaillie Invasive—Proliferative scores and (ii) five-gene signature based invasive-proliferative scores with Hallmark IFNγ signaling. (F) Scatterplot showing associations between the Verfaillie Invasive—Proliferative scores with Hallmark IFNγ signaling in the TCGA melanoma cohort. (G) Scatterplot of single cell RNA-Seq data showing individual cells from each cell line projected on the invasive-proliferative and Hallmark IFNγ signaling planes. Each color represents a different cell line. CCLE, cancer cell line encyclopedia; RACIPE, random circuit perturbation; TCGA, The Cancer Genome Atlas.

Next, we investigated how PIT status associates with IFNγ signaling. Our simulations suggest that, similar to PD-L1, IFNγ signaling is higher in neural crest-like and undifferentiated phenotypes (online supplemental figure S8B,C). Furthermore, while all four phenotypes are observed at lower IFNγ levels, neural crest-like and undifferentiated phenotypes predominate at higher levels of IFNγ (figure 4D). This trend is reminiscent of experimental data suggesting that IFNγ signaling can induce dedifferentiation in melanoma.47 In CCLE and TCGA samples, we observed strong positive associations between the proliferative-invasive status and IFNγ signaling for both the reduced five-gene signature or Verfaillie-based metric of proliferative-invasive status (figure 4E,F). This positive association was also recapitulated when considering the updated invasive score and IFNγ signaling (online supplemental figure S8D). Additionally, treatment of melanoma cells with IFNγ signaling led to significant, but variable decreases in proliferative scores and an increase in invasive scores (online supplemental figure S8E). Similar results were observed when we analyzed a group of n=8 melanoma cell lines treated with either TNF or IFNγ (online supplemental figure S8F–H) (GSE152755).47 While treatment of TNF and IFNγ caused expected significant increases in the corresponding pathways, PD-L1 levels were also induced, especially for IFNγ-treated cells (Fig S8F). However, the extent of dedifferentiation in either case was variable (online supplemental figure S8G,H). This is consistence with reports where transcription factors such as JUN have been shown to have a relatively modest role in IFNγ induced expression of PD-L1 in melanoma cells.49 Analysis of single cell RNA-seq data of ten patient-derived MM lines (GSE134432), revealed a strong positive correlation between intrinsic IFNγ signaling levels and invasive/proliferative scores (p=0.609; p<10⁻⁵) (figure 4G). Together, these results show that increased IFNγ signaling in melanoma can affect levels of PD-L1 and/or differentiation status of melanoma cells via the underlying GRN.

Combinatorial control of PD-L1 levels by IFNγ signaling and proliferative-invasive status of cells underlies heterogeneous expression patterns of PD-L1 in melanoma

Our model simulations indicated that proliferative and intermediate states typically have lower PD-L1 levels, but neural crest-like and undifferentiated phenotypes have either intermediate or high levels of PD-L1 (figure 5A, left). The highest levels of PD-L1 observed were also associated with a higher activation of IFNγ signaling (figure 5A, right). These observations suggested that PD-L1 levels could be upregulated through two possible paths: (A) increase in invasive nature and (B) enhanced IFNγ signaling. We further quantified the contribution of these two paths and their interdependency.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5

Analysis of coupled dynamics for PD-L1, proliferative-invasive phenotype, and IFNγ signaling. (A) Associations between invasive—proliferative score and PD-L1 levels colored by the four phenotypes (left) and IFNγ signaling levels (right). (B) Stacked bar plots showing the fraction of steady states belonging to low (gray), mid (orange), or high (red) PD-L1 levels in either low IFNγ signaling levels (left) or high IFNγ signaling levels (right). The main proportion of steady states belonging to each group are indicated. (C–E) Bar plots showing the relative fold changes in the mRNA expression levels of MITF, JUN, NGFR and PD-L1 after 48 hours in non-treated (NT), TNFα/TGFβ treated, IFNγ treated and TNFα/TGFβ /IFNγ treatment conditions in melanoma cells (C) SK-MEL-5 (D) MALME-3M (E) A375. *p<0.05, **p<0.01, ***p<0.001, ****p<0.0001 represents a statistically significant difference in the levels based on paired t-test. (F) Scatterplot of single cells from melanoma patients before immune checkpoint inhibitors therapy projected on the Invasive-Proliferative score and Hallmark IFNγ signaling (top) and box plots of their corresponding PD-L1 levels in a patient specific grouping sorted ascending order according to median expression values (bottom). Each color in scatterplot corresponds to a particular patient sample.

Thus, we segregated the scenarios belonging to high versus low IFNγ signaling and subsequently calculated the conditional probability of each of the four phenotypes to display low, medium, and high levels of PD-L1. For low IFNγ signaling, proliferative and intermediate phenotypes were most likely (>95%) to express low PD-L1 levels, but for neural crest-like and invasive phenotypes, PD-L1 levels were either medium (49%–51%) or high (32%–40%) (figure 5B, left). Importantly, while proliferative and intermediate phenotypes were more homogeneous in terms of PD-L1 levels, the neural crest-like and invasive ones exhibited higher variance. This behavior recapitulates the trends seen in a single-cell RNA-seq cohort of melanoma cell lines spread across the proliferative-invasive spectrum (figure 3D,i). Together, these results suggest that in absence of IFNγ signaling, a PIT can increase PD-L1 levels. For high IFNγ signaling, amplified heterogeneity in PD-L1 is seen for proliferative and intermediate phenotypes. Conversely, invasive and neural crest-like are most likely (>98%) to display high PD-L1 and thus more homogeneous (figure 5B, right). Importantly, we notice that a substantial proportion of cells (>25%), irrespective of their proliferative-invasive status, can exhibit high PD-L1 levels, under scenarios of high IFNγ signaling (compare the proportion of PD-L1 high cells (red) across phenotypes in figure 5B, right vs left). This trend is consistent with experimental reports showing PIT inducing factors such as JUN having a relatively minor influence on IFNγ-induced PD-L1 expression.49 Put together, our simulations highlight that increased IFNγ signaling, and a more invasive nature can both augment PD-L1 levels.

To experimentally test these model predictions, we chose three melanoma cell lines—SK-MEL-5, MALME3M (both representative of a proliferative phenotype) and A375 (representative of the neural-crest phenotype) based on their gene expression profiles in the CCLE group of melanoma cell lines. SK-MEL-5 and MALME3M have higher proliferative scores compared with their invasive scores and lie toward the proliferative end of the spectrum on a two-dimensional proliferative-invasive plane of the Verfaillie gene sets (online supplemental figure S9A). A375 cells, however, are more toward the center of this phenotypic spectrum (online supplemental figure S9A, left) and exhibit higher levels of neural crest markers compared with MALME3M and SK-MEL-5 (online supplemental figure S9A, right). To test our model predictions, we treated each of these cell lines with either a combination of TNFα and TGFβ to induce a PIT or we treated the cell lines with IFNγ as an independent modulator of PD-L1 levels and assessed the mRNA levels of MITF (proliferative marker), JUN (invasive marker), NGFR (neural-crest marker) and PD-L1 levels after a duration of 48 hours. Treatment of MALME3M and SK-MEL-5 cells with both TNFα and TGFβ led to significantly lower levels of MITF and higher levels of JUN, indicative of a less proliferative and more invasive cell state (figure 5C,D). These changes were also accompanied by significant increase in NGFR levels, indicating the transition from a proliferative to a neural crest/invasive phenotype, as well as enhanced PD-L1 levels (figure 5C,D, right). Treatment with IFNγ alone, on the other hand, did not cause such drastic changes for proliferative (MITF), invasive (JUN) and neural crest (NGFR) markers (figure 5C,D). One explanation for this modest impact of IFNγ alone on these markers can be its epigenetic mode of action, which are often discernible over longer time scales.64 Despite minimal impact on PIT, IFNγ treatment alone induced PD-L1 levels to an even higher extent as compared with TNFα/TGFβ treatment (figure 5C,D, right), thus demonstrating that IFNγ can induce PD-L1 levels efficiently with minimal effects on PIT at relatively short timescales. On a combinatorial treatment of TNFα/TGFβ and IFNγ, PD-L1 levels increased synergistically in both proliferative cell lines, potentially showing an emergent relationship between the two modes of induction of PD-L1 levels (figure 5C,D, right). This observation resonates well with our simulation results that under conditions of low IFNγ levels, proliferative cells are predominantly (97.35%; figure 5B, left) low in terms of PD-L1 levels, but under scenarios of high IFNγ levels and a switch (induced by TNFα/TGFβ) to more neural crest like or undifferentiated phenotype, most cells (>98.1%; figure 5B, right) express high PD-L1.

Next, we checked the effects of individual and combinatorial treatments of TNFα/TGFβ and IFNγ on a neural crest cell line, A375. We observed a robust decrease of proliferative (MITF) and upregulation of invasive (JUN) markers on TNFα/TGFβ treatment (figure 5E). IFNγ treatment alone inhibited MITF, but had modest effects on levels of JUN, consistent with earlier reports.49 Earlier observations report that neural-crest phenotype has the highest NGFR expression65 and thus we would expect that a transition to a more undifferentiated phenotype would lead to a decrease in NGFR levels. Validating this hypothesis, we observed robust decrease of NGFR mRNA expression (figure 5E), indicating that the cells that were originally of the neural crest phenotype transitioned to a more undifferentiated invasive phenotype. This was again accompanied by robust upregulation of PD-L1 levels in TNFα/TGFβ and IFNγ treated cases with the combinatorial treatment (figure 5E, right), showing the highest increase as predicted by our model.

Finally, to determine if these trends were observed in melanoma cells from treatment naïve patients (GSE115978),66 we assessed how PD-L1 levels depended on both the extent of PIT and activity of IFNγ signaling. Typically, cells from each patient sample cluster together in this plot of invasive score and IFNγ signaling activity (figure 5F, top), showcasing greater interpatient phenotypic variability than intrapatient variability. Across the patient samples, we observe scenarios of varying degrees of invasive score and IFNγ signaling, enabling us to evaluate how these both axes control PD-L1 levels.

First, we focused on pretreatment samples. Compared with other cells, Mel129pa cells (pink cluster) were less invasive and had lower IFNγ signaling (figure 5F, top). These cells had the lowest PD-L1 levels (figure 5F, bottom), thus validating our model predictions that proliferative phenotype with low IFNγ signaling has minimal PD-L1 expression (figure 5B, left). These analyses support a model in which at least one of the two axes is needed to increase PD-L1 levels: PIT or IFNγ signaling. Conversely, patient cells that were enriched for both IFNγ signaling and an invasive phenotype (Mel97—blue cluster; Mel89—purple cluster) had the highest PD-L1 levels (figure 5F). This observation also recapitulates our model predictions that neural crest-like or invasive states are highly likely to have upregulated PD-L1 levels (figure 5B, right). Next, we focused on samples which had only one of the two axes being active: PIT and IFNγ signaling. Samples with high IFNγ signaling but low invasive scores (Mel103—brown cluster; Mel 71—orange cluster) had higher levels of PD-L1 than the ones with high invasive score but with low IFNγ signaling (Mel81—green cluster) (figure 5F). However, both these scenarios had lower PD-L1 levels compared with the ‘double positive’ high invasive and high IFNγ activity subset (figure 5F). Therefore, our dynamical model can reproduce the origins of PD-L1 heterogeneity in melanoma as a function of contributions from PIT and IFNγ signaling.

Exploring the dual routes of therapeutic resistance in melanoma: pre-existing and induced heterogeneity pathways and varied trajectories

To better decipher the complex interplay between PD-L1 expression, PIT, and IFNγ signaling in melanoma under the influence of both targeted therapies and immune checkpoint inhibitor therapies, we analyzed transcriptomic profiles of WM989 cells treated with vemurafenib (a potent, highly selective inhibitor of mutated BRAF V600E) after 48 hours of treatment, and 7 days postdevelopment of therapy resistance (GSE97681).67 We observed that melanoma cells post vemurafenib treatment at both these time points were more invasive and had enhanced IFNγ signaling activity (figure 6A). This combined increase along both IFNγ signaling axis and proliferative-invasive axis is accompanied by an increase in PD-L1 expression levels. A short-term response (48 hours of treatment with vemurafenib) was found to be sufficient to induce a PIT (figure 6A). This trajectory of drug resistance was also seen in other studies in other datasets for the same cell line (WM989) treated either treated with vemurafenib (GSE161299) (online supplemental figure S9B) alone or in combination with pinometostat, a clinical-stage DOT1L inhibitor (GSE161298) (online supplemental figure S9C).

• Download figure
• Open in new tab
• Download powerpoint

Figure 6

Role of pre-existing and induced heterogeneity in development of drug resistance to both targeted therapy and immunotherapy. (A) Box plot illustrating the changes in proliferative-invasive status (top), IFNγ signaling (middle) and PD-L1 expression levels (bottom), following treatment with vemurafenib for 48 hours and 7 days. (B) (i) Scatterplot of single-cell RNA sequencing (scRNA-seq) data for MEK inhibitor-treated A375 cell line, comparing untreated and vemurafenib-treated conditions, projected on the proliferative-invasive plane based on the imputed expression of Verfaillie signatures. (ii) Scatterplot depicting IFNγ signaling levels in the A375 cell line under untreated and vemurafenib-treated conditions. The red box denotes the transcriptional state which was present in untreated but was absent in Vemurafenib resistant clones. (C) (i) Scatterplots of all A375 cell lines under different MEK inhibitor treatment conditions, projected on proliferative-invasive planes, based on the imputed expression of Verfaillie signatures. (ii) Scatterplot illustrating the differences in IFNγ signaling levels under various MEK inhibitor treatment conditions. (D) Scatterplot showing immune checkpoint inhibitor pretreatment (black) and post-treatment (red) single cells projected on the proliferative-invasive plane. Arrows represent hypothesized directions of drift of melanoma cells during adaptation to therapy. (E) Contour showing immune checkpoint inhibitor pretreatment (black) and post treatment (red) single cells projected on the proliferative-invasive spectrum and Hallmark IFNγ signaling. Arrows represent hypothesized directions of drift of melanoma cells due to adaptation to therapy. *Denotes a statistically significant difference in levels based on Student’s t-test.

Next, we investigated the single cell transcriptomics profiling of 31,000 BRAF-mutant A375 cells in 4 treatment conditions: untreated/sensitive cells, cells resistant to vemurafenib, cells double-resistant to vemurafenib and cobimetinib, cells double-resistant to vemurafenib and trametinib (GSE164897).68 We observed that the untreated cells exhibited two distinct subpopulations—one of which exhibited a transcriptional state present only in the untreated cells (figure 6B, red box), and the other overlapped with the state seen in treated cells, indicating pre-existing heterogeneity in terms of varying invasive nature. In the treated population, cells displayed a predominantly more invasive state (figure 6B,i). This change could indicate two possible scenarios that are not necessarily mutually exclusive: (A) pre-existing heterogeneity in proliferative-invasive status of cells may have provided more invasive cells with a selective advantage during the evolution of drug resistance to vemurafenib, (B) drug treatment induced a subpopulation of cells to alter phenotypes. Furthermore, these cells had lower levels of IFNγ signaling (figure 6B,ii) despite having a higher invasive score. This trend agrees with another instance (GSE108383) of vemurafenib resistant clones of A375 and 451Lu cells characterized by reduced IFNγ signaling compared with their sensitive counterparts (online supplemental figure S9D). We also observed that a slight increase in PD-L1 levels post-treatment in vemurafenib-resistant cells, and a positive correlation with the invasive status (online supplemental figure S9E). These results indicated that in some adaptation strategies, the positive association between PIT and IFNγ signaling can be compromised, while the levels of PD-L1 and the extent of PIT can increase concurrently.

On the other hand, in case of combinatorial drug treatment, that is, treatment with BRAF inhibitor vemurafenib in combination with MEK1 inhibitor trametinib or cobimetinib, we observed the emergence of another phenotype that was even more invasive than the vemurafenib-resistant one (figure 6C,i). This state was earlier not observed in untreated samples, and was accompanied by an increase in IFNγ signaling (figure 6C,ii), suggesting an induction to a distinct new transcriptional state not previously present in the pretreatment setting. Such new attractor states enabled by the induction through combined BRAF/MEK1 inhibition may also contribute to the development of resistance and tumor relapse. Furthermore, the cell-state transition trajectory(ies) during the evolution of resistance to a targeted therapy combination could depend on cell type itself. For example, a combinatorial treatment of MEK and CDK4/6 inhibitors over a prolonged period (~1 month) promoted four different cell populations (GSE230538),69 each with distinct cell state trajectories on the proliferative-invasive plane. We observed that the cell lines which were drug-naïve on treatment over 1, 4 and 33 days may or may not show a reversal to their original proliferative-invasive status during the evolution of resistance. For example, MEL-JUSO and M20 became significantly more proliferative as well as significantly more invasive on becoming resistant while the cell lines IPC-298 and SK-MEL-30 showed an increase in both proliferative and invasive scores but returned to a proliferative-invasive state like their treatment naïve case at the end of 33 days of drug treatments (online supplemental figure S9F).

We next investigated changes in transcriptional state along the proliferative-invasive spectrum and IFNγ signaling on cells becoming resistant to immune checkpoint inhibitors (ICIs). We analyzed single-cell transcriptomics data from patient-derived melanoma cells that are resistant to immune checkpoint inhibitor (GSE115978) therapy.66 Projecting the pretreatment and post-treatment single-cell profiles along the axes of invasive nature, PD-L1 expression and IFNγ signaling activity, we observed two distinct trajectories that can correspond to adaptive response to immune checkpoint inhibitor therapy. In one trajectory, melanoma cells become more invasive, as shown by a distinct increase in invasive score and a concomitant decrease in proliferative scores (figure 6D). These cells also have reduced IFNγ signaling activity (figure 6E). This trajectory was recently reported in experimental data demonstrating that melanoma patients responding to anti-PD-1 therapy exhibit a more proliferative gene expression signature, while the non-responders were enriched in invasive and neural crest-like phenotypes.70 Another trajectory that cells can follow is cells exhibiting a more proliferative nature, without a noticeable change in IFNγ signaling activity (figure 6D,E). A similar trajectory is reported in response to anti-BRAF therapy, where cells can exhibit a hyperdifferentiated phenotype.71 72 Finally, we analyzed post-treatment samples, and observed diverse trajectories of drug adaptation. We plotted the single-cell RNA-seq data from melanoma patients’ post-therapy on two-dimensional plane of proliferative-invasive score and IFNγ signaling (online supplemental figure S9G, top). Similar to treatment naïve samples (figure 5F), interpatient variability dominates over intrapatient variability across both the axes (online supplemental figure S9G). We observed that the least invasive cells had the highest expression of PD-L1 levels (Mel88, Mel78, Mel98 - orange, blue and pink clusters, respectively). However, these cell clusters were also higher in IFNγ signaling compared with the more invasive clusters, such as Mel110 (purple) and Mel106 (brown), which had relatively lower levels of PD-L1 (online supplemental figure S9G, bottom). These results indicate that IFNγ signaling levels may be a dominant determinant of PD-L1 levels as compared with the proliferative-invasive status of cells in immune checkpoint inhibitor-resistant melanoma cells. Such co-occurrence of adaptive drug-tolerant states is reflective of how isogenic cells exposed to identical therapeutic stress can diverge toward different ‘attractors’ in the phenotypic landscape, driving non-genetic heterogeneity. Together, our results highlight the importance of both pre-existing heterogeneity and cell-state transition driven heterogeneity along the different trajectories on the IP and IFNγ signaling axes, as two distinct but not necessarily exclusive routes to emergence of drug resistance in cancer cells. However, the exact underlying mechanisms of resistance, and the interplay between these two routes, require further investigation and more mechanistic modeling, indicating the future scope of our study.

Overall, our dynamical model, coupled with single-cell transcriptomic analysis of melanoma samples, indicate that both IFNγ signaling activity and PIT can govern PD-L1 levels in a combinatorial manner in treatment-naïve conditions and specific cases of resistance to targeted and/or immune checkpoint inhibitor therapies. Moreover, diverse cellular trajectories of adaptive response to anti-BRAF or immune checkpoint inhibitor therapy can alter how PIT and IFNγ-induced dedifferentiation shape the PD-L1 heterogeneity patterns in a tumor.