Mapping phenotypic heterogeneity in melanoma onto the epithelial-hybrid-mesenchymal axis

Epithelial to mesenchymal transition (EMT) is a well-studied hallmark of epithelial-like cancers that is characterized by loss of epithelial markers and gain of mesenchymal markers. Melanoma, which is derived from melanocytes of the skin, also undergo phenotypic plasticity toward mesenchymal-like phenotypes under the influence of various micro-environmental cues. Our study connects EMT to the phenomenon of de-differentiation (i.e., transition from proliferative to more invasive phenotypes) observed in melanoma cells during drug treatment. By analyzing 78 publicly available transcriptomic melanoma datasets, we found that de-differentiation in melanoma is accompanied by upregulation of mesenchymal genes, but not necessarily a concomitant loss of an epithelial program, suggesting a more “one-dimensional” EMT that leads to a hybrid epithelial/mesenchymal phenotype. Samples lying in the hybrid epithelial/mesenchymal phenotype also correspond to the intermediate phenotypes in melanoma along the proliferative-invasive axis - neural crest and transitory ones. As melanoma cells progress along the invasive axis, the mesenchymal signature does not increase monotonically. Instead, we observe a peak in mesenchymal scores followed by a decline, as cells further de-differentiate. This biphasic response recapitulates the dynamics of melanocyte development, suggesting close interactions among genes controlling differentiation and mesenchymal programs in melanocytes. Similar trends were noted for metabolic changes often associated with EMT in carcinomas in which progression along mesenchymal axis correlates with the downregulation of oxidative phosphorylation, while largely maintaining glycolytic capacity. Overall, these results provide an explanation for how EMT and de-differentiation axes overlap with respect to their transcriptional and metabolic programs in melanoma.

Introduction

Epithelial to mesenchymal transition (EMT) is a well-characterized phenomenon involved in multiple axes of cancer progression, such as metastasis and drug resistance. EMT is commonly associated with morphological changes, functional changes (increased migration, invasion, and immune invasion) (1–3) and molecular changes, including upregulation of EMT markers and transcription factors (TFs), such as VIM, ZEB1, SNAI1 and TWIST1. While the phenomenon of EMT has largely been characterized for epithelial cancers (such as breast cancer and lung adenocarcinoma), similar molecular, functional and morphological changes have also been observed in non-epithelial cancers, such as sarcomas (4, 5), glioblastoma (6), myeloma (7), lymphoma (8, 9), leukemia (10, 11) and melanoma (12) in preclinical and clinical settings.

Treatment of melanoma tumors harboring BRAFV600E mutation often involves targeted therapy strategies that inhibit BRAF or MEK signaling. While these targeted agents provide clinical benefit to melanoma patients, resistance to these therapies is common. Therapy-resistant melanomas often undergo de-differentiation, which is characterized by loss of melanocytic markers such as MLANA, TRPM1 and TYR and gain of invasive molecular markers such as c-JUN, NGFR and ZEB1 (13–16). The de-differentiation trajectory of melanoma cells is characterized by a transition along the proliferation-invasion axis, from a melanocytic phenotype to an undifferentiated phenotype while passing through the intermediate transitory and neural crest stem cell-like (NCSC) phenotypes (Figure 1A). This trajectory is the reverse of the differentiation that occurs during melanocyte development, where undifferentiated tissue in the embryonic neural plate give rise to highly migratory and mesenchymal neural crest cells, some of which differentiate into melanocytes upon reaching the epidermis (17). Therapy resistant melanoma is also commonly associated with a mesenchymal-like phenotype with more invasive and aggressive features (13, 16, 18–20). These relationships between de-differentiation, invasion, and EMT pathways in response to therapy suggest EMT and de-differentiation programs in melanoma may be linked.

FIGURE 1

Figure 1 Mapping phenotypic heterogeneity in melanoma onto the EMT axis. (A) A schematic representation. Volcano plots depicting Spearman’s correlation coefficients and -log₁₀(p-value) of 78 datasets for the Verfaillie proliferative and invasive gene set with (B) 76GS EMT scoring metric, and with (C) KS EMT scoring metric (D) Boxplot depicting range of correlation coefficients for KS and 76GS with Verfaillie invasive and proliferative gene sets. Volcano plots depicting the Spearman’s correlation coefficient and -log₁₀(p-value) of 78 datasets for Verfaillie proliferative and invasive gene set with (E) Epithelial gene set (E scores) and (F) Mesenchymal gene set (M scores). (G) Boxplot depicting range of correlation coefficients for E and M scores with Verfaillie invasive and proliferative gene sets. Inset labelled “Significant” is calculated as the fraction of datasets (out of 78) which show a significant correlation trend (r < - 0.36 or r > 0.36, p < 0.05). The absolute number of significant points (datasets) for the specified cut-off is indicated in brackets. “Proliferative” and “Invasive” labels represent the percentage of significant correlations that are between the EMT score and proliferative score or invasive score, respectively.

The similarity between EMT and de-differentiation programs extends beyond cell-intrinsic alterations and impacts cell-extrinsic changes as well. EMT often leads to varied extracellular matrix (ECM) stiffness and density (21–23) and altered cell-matrix and cell-cell interactions (24, 25). In melanoma, acquisition of de-differentiated and invasive phenotypes is often accompanied with changes in composition and physical properties of ECM, and modified cell-matrix interactions and cell morphology (26–28). Increased expression of matrix metalloproteases (MMPs), immune evasion (characterized by both signaling-mediated immune suppression (e.g. by TGF-ß release) and prevention of immune cell entry into tumors by dense collagen matrix/low α-SMA expression), increased inflammatory markers (such as TNF-α, NF-kB and AP-1) and cytoskeleton remodeling have been closely linked to the acquisition of an invasive phenotype and loss of melanocytic differentiation regulator MITF (29–34). All of these changes are reported with EMT progression as well in multiple epithelial cancers (35–37). Such extensive similarity between EMT and de-differentiation programs in cancer-microenvironment cross-talk and niche construction underscore the potential of common regulatory pathways involved in both EMT and de-differentiation.

Another common feature that links EMT in epithelial cancers to de-differentiation in melanoma is the presence of intermediate or hybrid phenotypes. Hybrid epithelial/mesenchymal (E/M) cells express molecular and functional characteristics of both epithelial (high proliferation and cell-cell adhesion, low invasion) and mesenchymal (low proliferation and cell-cell adhesion, high invasion) cells (38). On the other hand, melanoma intermediate phenotypes, which comprise transitory and neural crest-like stem cell-like (NCSC) phenotypes, exhibit combined features of proliferative and invasive phenotypes (39, 40) (Figure 1A). Gene regulatory networks for EMT and melanoma provide a mechanistic basis for explaining the existence of these hybrid/intermediate states (41, 42). An overlap in key regulators and stabilizers for hybrid E/M phenotypes and melanoma phenotypes (such as ZEB1, NFATC2, CDH1, SNAI2, NRF2) suggest common regulatory links (13, 43–49). For instance, SNAI2, a stabilizer of the hybrid E/M phenotype, is a key regulator of the NCSC phenotype and metastasis in melanoma, suggesting its involvement in regulating the intermediate phenotypes in melanoma as well (45, 49). However, certain regulators show opposite trends in melanoma and EMT. For instance, ZEB2 is considered an inducer of EMT in epithelial cancers, but in the context of melanoma, it inhibits the mesenchymal phenotype (19, 50). Other molecules that show opposite effects include KLF4 (51, 52) and TFAP2A (53, 54). Thus, understanding the mechanistic underpinning of how the de-differentiation and EMT programs are linked can help decipher reasons for the similarities and differences between these pathways across cancers.

In this study, we map the de-differentiation axis in melanoma (also called proliferative-invasive/P-I axis) to the EMT axis using previously defined scoring metrics (3, 55–57). We compare the extent to which a gain in a mesenchymal signature corresponds to a loss in the epithelial signature during de-differentiation of melanoma. By deciphering the interdependencies between de-differentiation and mesenchymal programs, the differences in molecular regulation between EMT and de-differentiation can be better understood. We have identified that the mesenchymal program, but not the epithelial program, is closely linked to de-differentiation. Although the mesenchymal signature enrichment shows a strong negative correlation with a differentiated/melanocytic transcriptional program, it does not increase monotonically during de-differentiation. This non-monotonic trend is also captured by metabolic programs associated with EMT, such as glycolysis and HIF1α, but not with metabolic programs associated with differentiation/melanocytic genes, such as the MITF-regulated OXPHOS pathway. Our results indicate that phenotypic heterogeneity in melanoma occurs along a proliferative-invasive axis that correlates with a “one-dimensional EMT” in which cells transition along a mesenchymal axis without an alteration in epithelial phenotype. Deciphering such inter-connections among multiple axes of plasticity in a cancer cell population may guide potent combinatorial therapeutic strategies aimed at controlling transitions to a more hybrid cell type with combined features of both proliferation and invasion.

Materials and methods

Software and datasets

Publicly available datasets from Gene Expression Omnibus (GEO), The Cancer Genome Atlas (TCGA), Cancer Cell Line Encyclopedia (CCLE- Broad Institute) (58), and National Cancer Institute-60 (NCI-60) databases were analyzed. Microarray data were downloaded from GEO using GEOquery Bioconductor R package. All analyses done on R version 4.1.0. ggplot2, and ggpubr R packages were used to create and customize plots.

Pre-processing of datasets

Microarray datasets, with un-mapped probe IDs, were pre-processed by mapping the probe IDs onto their gene symbols using the relevant platform annotation table. In the case of multiple probes mapping to the same gene, the mean expression of all the probes was considered for that gene. For non-normalized RNA-Seq datasets TPM normalization followed by log2 transformation with an offset value of 1 was used.

ssGSEA

Single-sample Gene Set Enrichment Analysis, an extension of Gene Set Enrichment Analysis (GSEA) (56, 59), calculates separate enrichment scores for each sample and a gene set. Each score represents the degree to which genes in a gene set are up or down-regulated in a sample. We calculated ssGSEA scores for the Verfaillie proliferative and Verfaillie invasive gene sets (60), Hoek proliferative and Hoek invasive gene sets (39), the epithelial (E) and mesenchymal (M) gene sets of the EM tumor gene signature genes and cell lines gene signatures in the KS scoring metric (57), and the Tsoi melanocytic, transitory, NCSC, and undifferentiated gene set (40) (Table 1).

TABLE 1

Table 1 List of scores used for quantifying various axes of heterogeneity.

Calculation of EMT scores

We calculated EMT scores of datasets using four metrics- 76 Gene Signature (76GS), Kolmogorov -Smirnov test (KS), E score and M score (Table 1). 76GS and KS were calculated as defined earlier (1, 55, 57). 76GS score is a metric for how epithelial a sample is, i.e., higher scores reflect greater association with an epithelial phenotype. The KS score is a metric for how mesenchymal a sample is. The higher the KS score of a sample, the greater is its association with a mesenchymal phenotype. While 76GS scores do not have a pre-defined range of scores, KS scores lie within a +1 to -1 range. The E and M scores are ssGSEA scores for epithelial and mesenchymal gene lists, respectively, for the KS scoring metric (3). For calculating KS, E and M scores, datasets were classified based on whether the samples were derived from cell-lines or tumors and the appropriate gene sets were used.

Correlations

All correlation values were calculated using Spearman’s correlation coefficient, unless mentioned otherwise. Spearman’s correlation coefficient method generates a coefficient ranging between –1 to +1, where +1 indicates a strong positive correlation, and –1 indicates a strong negative correlation between two variables. It determines the correlation between any monotonically related variables- linear or non-linear. Correlations with R >0.36 and p<0.05 are considered significant.

Moving window average

A moving window average is used to quantify the gradient for a variable along a given axis. A window covering 60% of the entire range of the axis is created and the average value of the variable for all samples in the window is calculated. Then the window is then shifted by 1% and the average is re-calculated. This is iteratively repeated to cover the entire range. The slope of the averages determines the magnitude and direction of the gradient.

Conditional probabilities

Once the cell lines were sorted into their respective phenotypes and the conditional probabilities were obtained, the statistical significance and p-values for the conditional probabilities were calculated using the one-proportion Z test.

The z-score was calculated using the equation

$$z=\frac{\widehat{p}-p_0}{\sqrt{\frac{p_0\left(1-p_0\right)}{n}}}$$

where $\widehat{p}$ is the observed proportion, p₀ is the null probability, and n is the sample size. The obtained value of z was then converted into the corresponding p-value using the standard normal distribution. If the obtained p-value < 0.05, it was considered significant.

Assigning phenotypes to samples

In order to identify samples belonging to the 4 phenotypes (melanocytic, transitory, NCSC and undifferentiated), we calculated ssGSEA scores based on gene sets for each of these phenotypes (40). Samples lying in the top 10% scores were assigned that particular phenotype. Taking a cut-off value of less than 10% would can enable only one point being selected for each phenotypes in datasets having less than 20 samples (e.g. int in Figure 4D, GSE101434) while increasing this threshold might lead to non-specific phenotype cells being selected in larger datasets. Thus, 10% was chosen as an optimal cut-off.

Metabolic scores

The oxidative phosphorylation (OXPHOS) and glycolysis (Glyco) scores in our study were calculated using ssGSEA carried out with the corresponding hallmark gene sets for these pathways [obtained from Molecular Signature Database (MSigDB) (61)]. The HIF-1 signature - which is a surrogate for glycolysis - was quantified based on a method previously reported (64). This method uses expression levels of their downstream target genes to capture the respective enzyme activities. A total of 59 downstream genes for HIF-1 were used and the scores were obtained using the Singscore method performed on these gene sets (62, 65). The fatty acid oxidation (FAO) scores were calculated based on the equation previously reported (63) which uses expression levels of 14 FAO enzyme genes.

Results

Enrichment of mesenchymal genes can capture the extent of de-differentiation in melanoma

To test whether EMT and de-differentiation in melanoma programs are correlated with one another, we used previously-defined EMT scores – KS and 76GS (55, 57) – and ssGSEA scores for Verfaillie proliferative and invasive (60) and Hoek proliferative and invasive melanoma gene sets (39) and investigated their correlation coefficients across 78 datasets. Additionally, to dissect the contributions of epithelial and mesenchymal gene set separately, we calculated the ssGSEA scores (56, 59) for corresponding gene sets individually too (57), referred here as E and M scores, respectively (3). A sample with a higher 76GS or E score is more epithelial while a higher KS or M score refers to more mesenchymal phenotype. Thus, given the overlap between mesenchymal and invasive programs, we expected invasive scores to correlate positively with KS and M scores and negatively with 76GS and E scores. We also expected opposite trends for proliferation scores: negative correlations with KS and M scores and positive correlation with 76GS and E scores. We visualized the relationships between these pathways as volcano plots in which each dot corresponds to a dataset analysed. For positively-correlated metrics, we expect the majority of data sets to lie in the top right rectangle, while those displaying a significant negative correlation are expected to lie in the top left rectangle.

A total of 34 out of 78 datasets (43.59%) showed a significant negative correlation (r < - 0.36, p < 0.05) between 76GS and one of the two Verfaillie (proliferative, invasive) scores (66). In 30 out of those 34 datasets (88.23%), 76GS scores correlated negatively with invasive scores, while in remaining 4 datasets (11.76%), 76GS scores correlated negatively with proliferative scores (Figure 1B, left). Similarly, among 45 datasets that showed a positive correlation (r > 0.36, p < 0.05) between 76GS scores and one of Verfaillie scores, 38 (84.4%) cases had a positive correlation between 76GS and proliferative scores, and in the remaining seven datasets, 76GS scores correlated positively with invasive scores (Figure 1B, right). Overall, both the scoring metrics (76GS and KS) displayed correlations with Verfaillie and Hoek proliferative and invasive scores across the 78 datasets to support a relationship between E/M plasticity and the proliferative/invasive axis (Figures 1B–D, S1A–C).

Because gain of mesenchymal features is reported more commonly in melanoma as compared to loss of epithelial features, we decoupled the epithelial and mesenchymal components of the scoring metrics (E and M scores, respectively). The KS method provides separate information on genes that are associated with an epithelial phenotype and those with a mesenchymal state. Using the genes from the KS scoring method we segregated the genes and calculated individual ssGSEA scores for epithelial and mesenchymal gene lists and re-evaluated their correlation with proliferative and invasive scores in melanoma. While epithelial genes continued to show random distributions of samples throughout the plot, mesenchymal genes showed clear segregation of proliferative and invasive scores based on Spearman’s correlation coefficients, i.e., invasive scores were positively correlated with M score while proliferative scores were negatively correlated with the M scores (Figures 1E-G, S1D–F). This observation suggests that mesenchymal genes, but not epithelial genes, can capture the phenotypic heterogeneity displayed by melanoma along the proliferative-invasive axis.

To provide further support for these observations, we focused only on Verfaillie gene sets, since they have levels of overlap with gene sets for the intermediate phenotypes that were previously identified (40) (Figure S1G). Thus, a continuous scoring metric defined for the Verfaillie gene set is expected to be more sensitive for capturing intermediate phenotypes as compared to the Hoek gene set.

Because correlation coefficients only provide an overall trend in data, we wished to determine how proliferative and invasive scores vary along the E and the M axis. For this purpose, we generated two dimensional EMT plots of the data sets in which E and M scores are represented along each of the two axes. These plots display the relative position of a sample along an epithelial or mesenchymal axis (3, 56, 59). We then overlay information on the proliferative and invasive scores for each sample. As expected, across various datasets, proliferative and invasive scores for samples had a stronger visible gradient along the M axis as compared to the E axis (Figures 2A–B). To quantify this gradient, we used a rolling window to estimate the increase of average proliferative and invasive scores across the E and M axis. For this, we start with a rolling window covering 60% of the entire range along a given axis and calculate the average proliferative (P) or invasive (I) score within that window. Then the window is shifted by 1% and the average is re-calculated. This process is repeated until the entire range is covered, and the change in averages is plotted. For an axis that strongly correlates with the change in scores, we expect a steeper slope. The nature of a slope (positive or negative) is determined by the correlation between the axis and the score. Both axes trend in the expected direction, with a positive slope for invasive scores and negative slope for proliferative scores along the M axis and vice versa for the E axis (Figure 2C). This analysis also reveals that the M axis has a steeper curve than the E axis for both P and I scores. These results suggest that proliferative-invasive heterogeneity in melanoma can be considered as a “one-dimensional form” of EMT where the mesenchymal program enrichment increases as cells become more invasive, but the epithelial program need not be suppressed concomitantly (Figures 2, S2), as often tacitly assumed for the case of EMT. Other non-epithelial cancers, such as glioblastoma (GBM) and sarcoma, also display larger variation along the M-score axis as compared to the E-score axis, suggesting that “one-dimensional EMT” might not be specific to melanoma alone (Figures S3A–B). Moreover, we also observe that more de-differentiated phenotypes in sarcoma display higher M scores, while in GBM a switch from proneural to mesenchymal phenotypes is clearly visualised along the M-score axis. Thus, phenotypic plasticity along a mesenchymal axis in non-epithelial cancers can take various trajectories.

FIGURE 2

Figure 2 Scoring metrics based on mesenchymal genes capture de-differentiation better than metrics based on epithelial genes. Two dimensional EMT plots of different types of datasets- (i) GSE7127 (63 melanoma cell lines - microarray), ii. CCLE (59 cell lines - microarray), iii.GSE4843 (45 tumor samples - microarray), iv.GSE65904 (214 tumor samples - microarray),v. GSE72056 (1257 single-cell tumor samples), vi.GSE81383 (307 single-cell tumor sample) depicting the variation of (A) Proliferative scores along the E and M score axes. (B) Invasive scores along the E and M score axes. (C) Quantifying the proliferative and invasive score gradient along the E-M axes using a rolling window.

The mesenchymal axis follows a non-monotonic relationship with de-differentiation

Because the M score axis was able to capture the phenomenon of de-differentiation quantified by continuous scoring metrics, such as the proliferative and invasive scores, we next tested if the discretized phenotypes also arrange themselves in order of appearance along the two dimensional EMT plane. The classification of samples into four categories - melanocytic, transitory, NCSC and undifferentiated (in order of increasing de-differentiation) - for GSE80829, GSE10916, GSE4843, GSE7127 and GSE116237 was done as previously defined (15, 42). Along the proliferative-invasive plane, samples displayed a strong negative relationship between the two scores, i.e., proliferative scores of samples decreased as their invasive score increased. The four phenotypes also appeared in the expected order (18, 40), with the melanocytic samples having the highest proliferative scores and lowest invasive scores, and the undifferentiated samples displaying the lowest invasive scores and highest proliferative scores (Figure 3A). However, the two dimensional EMT plane failed to resolve the four phenotypes in terms of these four phenotypes showing non-overlapping scores. Since the E score axis performed poorly previously (Figures 1E, G) in these metrics, we quantified the ability of M score axis alone to resolve the four phenotypes by quantifying the conditional probability of a sample to belong to the intermediate phenotypes (transitory and NCSC), given that they lie in an intermediate M score range. Interestingly, samples with intermediate M scores were significantly likely to belong to the transitory phenotype (Figures 3B, S3C, Table 2). However, the probability of these samples to belong to the NCSC phenotype was negligible. In some datasets (GSE7127, GSE116237), the melanocytic phenotype was also significantly enriched in the intermediate M score populations. However, the melanocytic phenotype cells were enriched in the bottom M score population as well, and were not uniquely present in the intermediate score range like the transitory phenotype cells (Figures S3D–F).

FIGURE 3

Figure 3 Variation of the four molecular phenotype scores along the epithelial, mesenchymal, proliferative, and invasive axes. (A) Plotting samples classified into four phenotypes onto the E-M, proliferative-invasive score axes. (B) Venn diagram depicting the intersection of the four phenotype scores of samples and intermediate M scores. p represents p-value for the conditional probability that a sample belongs to the phenotype given that they lie in the intermediate M score range.

TABLE 2

Table 2 Conditional probabilities for a sample belonging to a particular phenotype given that it lies in the intermediate M score range.

To further dissect the relationship between the four phenotypes and the M score axis, we quantified the change in M score with respect to the invasive scores for the four phenotypes. To identify the four phenotypes, we used ssGSEA scores for gene sets defined for each of the four phenotypes (40). The top 10% of samples that had the highest scores for a particular gene set, were assigned the label of that particular phenotype. We observed that in these samples there was a non-monotonic increase in M scores as invasive score/de-differentiation increased. As samples progressed from NCSC to undifferentiated, M scores either decreased (Figures 4C–E) or remained the same (Figures 4A–B, F). In the context of melanocyte development, neural crest cells are precursors for melanocytes with high migratory potential and high levels of EMT markers (17, 67, 68). Thus, the non-monotonic increase in the mesenchymal program seen here is reminiscent of the differentiation of melanocytes.

FIGURE 4

Figure 4 The mesenchymal axis follows a non-monotonic relationship with de-differentiation. Plotting M scores against invasive scores for different phenotypes along the P-I axis in many datasets: (A) GSE7127 (B) GSE158607 (C) GSE80829 (D) GSE101434 (E) GSE65904 (F) GSE19234.

Metabolic reprogramming along the proliferative-invasive axis in melanoma

The EMT status of epithelial cancer cells is often associated with distinct metabolic programs. Generally speaking, EMT is negatively correlated with the enrichment of oxidative phosphorylation (OXPHOS) and fatty acid oxidation (FAO), but positively correlated with glycolysis (62). In melanoma, the proliferative state is associated with high levels of OXPHOS and the invasive phenotype is associated with high levels of glycolysis (69–72), reinforcing the commonalities between these two different instances of phenotypic plasticity. Computational analysis has suggested the existence of four metabolic sub-populations (63): 1) OXPHOS-high/glycolysis-low, 2) OXPHOS-low/glycolysis-high, 3) OXPHOS- low/glycolysis-low, and 4) OXPHOS high/glycolysis-high.

To assess whether the OXPHOS-glycolysis metabolism axis can be mapped onto the proliferation-invasion axis, we calculated Spearman’s correlation coefficients between the metabolic scores (OXPHOS and glycolysis) and the de-differentiation scores (proliferative and invasive scores) (Figures 5A–C) across the 78 datasets. In 38 out of 78 datasets where the OXPHOS scores correlate significantly with proliferative scores, 34 datasets show a positive correlation. Similarly, among 43 datasets showing a significant correlation of OXPHOS scores with invasive scores, all of them showed negative correlation. Thus, overall, OXHOS scores corelated positively with proliferative scores and negatively with invasive scores (Figure 5A). Glycolysis scores, on the other hand, did not show a clear relationship with EMT status, with a subset of datasets showing trends in both the directions (positive and negative correlation) both for proliferative and invasive scores (Figure 5B). This difference is reminiscent of prior observations for the association of EMT with OXPHOS and glycolysis in which glycolysis is only moderately correlated with EMT status, but OXPHOS is consistently negatively correlated with EMT (62). This trend is substantiated by observations that in cases where OXPHOS is positively correlated with proliferative scores or negatively correlated with invasive scores, glycolysis scores do not show any particular direction of enrichment with either proliferative or invasive axes (Figure 5C).

FIGURE 5

Figure 5 Mapping metabolic programs associated with EMT onto the de-differentiation program axes. Volcano plots depicting Spearman’s correlation coefficient and -log₁₀(p-value) of 78 datasets for (A) Hallmark OXPHOS and Verfaillie gene set. (B) Hallmark glycolysis and Verfaillie gene set. (C) Spearman’s correlation coefficient between OXPHOS and Glycolysis and Verfaillie scores. (D) Hallmark OXPHOS and Tsoi gene set. (E) Hallmark glycolysis and Tsoi gene set. (F). Spearman’s correlation coefficient between OXPHOS and Glycolysis and Tsoi scores. N represents number of samples present in a given quadrant.

We next sought to dissect whether intermediate melanoma phenotypes might be enriched for a specific metabolic profile. To investigate this trend, we calculated the Spearman’s correlation coefficients for metabolic scores and ssGSEA scores for gene signatures corresponding to each of the four molecular phenotypes of melanoma (Figure 5D–F). OXPHOS showed a clear shift from datasets with a significant positive correlation with a melanocytic phenotype to a significant negative correlation for the undifferentiated phenotype (Figure 5D). On the contrary, glycolysis scores do not show a clear shift from negative to positive correlations with de-differentiation (Figure 5E). Similar to the non-monotonic trend observed for M-scores, the glycolysis scores show the strongest positive correlation trends for the NCSC phenotype. Undifferentiated phenotype scores have a mixture of positively correlated and negatively correlated datasets with respect to glycolysis scores. Put together, these observations suggest that the regulatory modules controlling the switch to glycolysis are likely linked to the mesenchymal program rather than the de-differentiation one. On the other hand, regulatory modules for OXPHOS are likely to be closely linked to the melanocytic differentiation program. This trend is in accordance with experimental evidence that suggests that OXPHOS in melanoma cells is regulated by PGC1α, a downstream target of MITF, a key regulator of melanocyte differentiation (73, 74). Interestingly, fatty acid oxidation, which is also directly controlled by MITF via SCD (75), also displays trends similar to OXPHOS (Figure S4A) while a HIF1α signature, that is commonly associated with the invasive phenotype follows a non-linear trend similar to glycolysis (Figure S4B), suggesting that it is linked to the mesenchymal program rather than the de-differentiation program.

Discussion

De-differentiation in melanoma occurs in response to targeted therapy. This process may be mediated by transitions across a spectrum of phenotypes in which melanocytic cells treated with BRAF/MEK inhibitors pass through a transitory phenotype, followed by the NCSC phenotype, before becoming completely un-differentiated (15, 18, 40). This trajectory is accompanied by loss of a proliferative signature and gain of invasive characteristics. Here, we decipher the relationship between de-differentiation and EMT in melanoma. These processes are often considered to co-occur during drug treatment (14, 16, 34); however, comparison of EMT and de-differentiation scores reveal that the two processes may be more closely related to the mesenchymal program rather than the loss of an epithelial-like state or an EMT program per se. This observation is reminiscent of previous results in breast cancer and melanoma in which regulatory genes for the mesenchymal and de-differentiated phenotypes overlapped while those corresponding to epithelial and differentiated (melanocytic) phenotypes did not overlap and were tissue-specific (76). Previous pan-cancer studies have also highlighted that downregulation of epithelial components and upregulation of mesenchymal features need not always be as strongly coupled as often assumed (77, 78). Moreover, differences along these two axes need not be necessarily reflected at a transcriptional level (79). Together, these observations highlight the need to analyze epithelial and mesenchymal axes independently, rather than as a conventional single metric for EMT.

Our results also indicate that metabolic programs can be linked either with the de-differentiation program or the mesenchymal program. OXPHOS and fatty acid oxidation are both indirectly regulated by MITF. In the case of OXPHOS, MITF regulates PGC-1α (74); in the fatty acid oxidation pathway, MITF regulates SCD (75). MITF, which controls both metabolic pathways, decreases with increasing de-differentiation. This trend is explained by the decline in MITF associated with de-differentiation, in accordance with the MITF rheostat model (80). On the other hand, glycolysis and HIF-1α signatures seem to be co-regulated with the mesenchymal program. Previous studies in epithelial cancers have shown how well-established EMT transcription factors (EMT-TFs) regulate the metabolic profile of a cell and cause a switch to glycolysis (also called Warburg effect) (81). Consistently, neural crest cells also display decay of glycolytic capabilities as they differentiate into melanocytes (82). Our analysis suggests that the metabolic state of a cell is closely linked to the transcriptional program governing it at a given time point. Thus, de-differentiation captures the transcriptional and metabolic states observed during melanocyte development.

Although our study focuses on melanoma, EMT-like phenotypic switching is also characteristic of other non-epithelial cancers and de-differentiation of melanocytes independent of malignant transformation. De-differentiation of melanocytes into pluripotent stem cells demonstrated a reduction in expression levels of E-Cadherin, an epithelial marker, and similarities to mesenchymal stem cells (83). Molecular subtypes of glioblastoma multiforme (GBM), a non-epithelial cancer, include the pro-neural, classical, and mesenchymal phenotypes, which exist along a spectrum of worsening prognosis (84). Single-cell analysis reveals that these molecular subtypes recapitulate neurodevelopmental trajectories, with proneural cells forming a major composition of proliferative glial progenitor-like cells (85, 86). A proneural-to-mesenchymal transition (PMT) is characterized by an increase in mesenchymal markers, such as SNAI1 and ZEB1. Similarly, glioma stem cells (GSCs) exist as proneural GSCs and mesenchymal GSCs, which can give rise to the complete spectrum of intra-tumor heterogeneity, including the classical phenotype (87), reminiscent of epithelial and mesenchymal CSCs reported in breast cancer (88). Moreover, samples belonging to the classical subtype are depleted of pro-neural GSCs and enriched for mesenchymal GSCs, possibly suggesting that mesenchymal GSCs are more likely to give rise to the classical subtype. This trend strengthens the semi-independent nature of EMT and stemness as seen in epithelial cancers (78). Another study in GBM cell lines reports that loss of N-cadherin (a well-established mesenchymal marker) increases invasiveness (89), reinforcing the trends that increased migration and invasion is not always an inevitable consequence of carcinoma-associated EMT (90). These scenarios of non-overlapping behaviors in terms of invasiveness, stemness and EMT, seen both for epithelial and non-epithelial cancers, advocate for improving existing therapeutic strategies by targeting multiple axes of cellular plasticity simultaneously rather than individually.

Our study focuses on the overlap between the de-differentiation and the EMT axis during drug treatment in melanoma samples. However, de-differentiation is not the only trajectory altered by drug treatment. Cells can follow multiple paths to therapy resistance, one of which is by attaining a hyper-pigmented phenotype (15, 91, 92). The mapping of these trajectories and states to the E-M axis remains to be studied. In addition, another axis of cellular plasticity commonly associated with EMT is immune suppression and immune evasion. Previous studies have shown that the expression levels of programmed death-ligand 1 transmembrane protein (PD-L1) – a driver of immune evasion - does not increase monotonically with EMT (3). Consistently, in melanoma, the expected trend of worse response to anti-PD-1 therapy with increasing de-differentiation is not observed; rather, results from the CheckMate 038 clinical trial indicate that the NCSC phenotype is associated with a better outcome to immune checkpoint blockade therapy as compared to the melanocytic phenotype (93). The extent of overlap between the axes of EMT, immune evasion, and de-differentiation require further study to design temporally-sequenced effective combination therapies that can shift the differentiation and EMT status of melanoma toward a less invasive and more immune activated state. Recent in vitro investigations in melanoma have shown proof-of-principle evidence of phenotypic plasticity driven drug resensitization as a mechanism underlying the beneficial impact of intermittent therapy (94).

Despite providing the abovementioned insights, our study suffers from various limitations. First, no mechanism-based models have been developed to gain insights into the emergent dynamics for the observed trends. A better understanding of the dynamics can help identify more effective therapeutic strategies by fine-tuning the interval, sequence, and dosage for combinatorial and/or sequential therapeutic strategies (95). Second, our analysis only characterizes phenotypes at a transcriptomic level, although preliminary investigation supports consistent trends at a proteomic level too (Figure S5). Third, due to limited availability of longitudinal transcriptomic data for varying treatment durations, our analysis is not restricted to time-resolved data exclusively. Preclinical data shows that short duration of drug treatment can induce a NCSC phenotype that is highly mesenchymal (14, 16), while prolonged treatment (8-12 weeks) can drive an undifferentiated phenotype. Our study indicates that a prolonged treatment can induce further de-differentiation but not always a concomitant increase in mesenchymal status, a prediction that needs detailed experimental validation. However, this observation of the NCSC phenotype being the most mesenchymal is in accordance with melanocyte development. Neural crest cells are progenitors of melanocytes that undergo EMT during development to delaminate and migrate from the neural tube to the epidermis, where they lose their EMT signature and differentiate into melanocytes (17, 67, 68). Thus, the non-monotonic variation in EMT during development (the initial increase during migration followed by decrease during differentiation) can be possibly recapitulated during treatment-induced de-differentiation. We propose that the often-presumed overlap between the mesenchymal and invasive axes may arise from the lack of information for longer time scales (since most in vitro drug treatment studies are performed in under three weeks), and often held assumptions about linearly increasing trends. However, increasing evidence suggests that maximum stemness is associated with hybrid E/M phenotypes rather than ‘extreme’ mesenchymal or epithelial phenotypes, suggesting that many such associations among axes of plasticity can be non-monotonic in nature (96–98).

Overall, our transcriptomic data-based analysis highlights the partially overlapping nature of EMT with molecular attributes of de-differentiation and metabolism during drug treatment in melanoma. We provide a framework for studying multiple intertwined axes of plasticity and heterogeneity (EMT, metabolic reprogramming, proliferative-invasive status) and identifying the degree to which these axes overlap.

Code Availability

All codes and scores generated for this paper can be found at:

https://github.com/csbBSSE/EMT_Melanoma

https://github.com/priyanka8993/EMT_score_calculation

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

MP and GR performed research, analyzed data and wrote the first draft of the manuscript. PT, NA, SM, AR and DB performed research. JS and MKJ conceived research and worked on manuscript revisions. MKJ supervised research. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by Ramanujan Fellowship awarded to MKJ by Science and Engineering Research Board (SERB), Department of Science and Technology (DST), Government of India (SB/S2/RJN-049/2018) and by Infosys Young Investigator award to MJ supported by Infosys Foundation, Bangalore. MP is supported by KVPY fellowship (DST). JS is supported by NIH 1R01CA233585-03, Department of Defense W81XWH-18-1-0189, and the Triangle Center for Evolutionary Medicine.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2022.913803/full#supplementary-material

References

1. Chakraborty P, George JT, Tripathi S, Levine H, Jolly MK. Comparative study of transcriptomics-based scoring metrics for the epithelial-Hybrid-Mesenchymal spectrum. Front Bioengineer Biotechnol (2020) 8:220. doi: 10.3389/fbioe.2020.00220/bibtex

CrossRef Full Text | Google Scholar

2. Hanahan D, Weinberg RA. Hallmarks of cancer: The next generation. Cell (2011) 144(5):646–74.

PubMed Abstract | Google Scholar

3. Sahoo S, Nayak SP, Hari K, Purkait P, Mandal S, Kishore A, et al. Immunosuppressive traits of the hybrid Epithelial/Mesenchymal phenotype. Front Immunol (2021) 12:797261/bibtex. doi: 10.3389/fimmu.2021.797261/bibtex

CrossRef Full Text | Google Scholar

4. Jolly MK, Ware KE, Xu S, Gilja S, Shetler S, Yang Y, et al. E-Cadherin Represses Anchorage-Independent Growth in Sarcomas through Both Signaling and Mechanical Mechanisms E-Cadherin Represses Anchorage-Independent Growth. Molecular Cancer Research (2019) 17(6):1391–402.

PubMed Abstract | Google Scholar

5. Somarelli JA, Shetler S, Jolly MK, Wang X, Bartholf Dewitt S, Hish AJ, et al. Mesenchymal-epithelial transition in sarcomas is controlled by the combinatorial expression of microRNA 200s and GRHL2. Molecular and cellular biology (2016) 36(19):2503–13.

PubMed Abstract | Google Scholar

6. Siebzehnrubl FA, Silver DJ, Tugertimur B, Deleyrolle LP, Siebzehnrubl D, Sarkisian MR, et al. The ZEB1 pathway links glioblastoma initiation, invasion and chemoresistance. EMBO Mol Med (2013) 5(8):1196–212. doi: 10.1002/emmm.201302827

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Roccaro AM, Mishima Y, Sacco A, Moschetta M, Tai YT, Shi J, et al. CXCR4 regulates extra-medullary myeloma through epithelial-Mesenchymal-Transition-like transcriptional activation. Cell Rep (2015) 12(4):622–35. doi: 10.1016/j.celrep.2015.06.059

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Lemma S, Karihtala P, Haapasaari KM, Jantunen E, Soini Y, Bloigu R, et al. Biological roles and prognostic values of the epithelial–mesenchymal transition-mediating transcription factors twist, ZEB1 and slug in diffuse large b-cell lymphoma. Histopathology (2013) 62(2):326–33. doi: 10.1111/his.12000

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Sánchez-Tilló E, Fanlo L, Siles L, Montes-Moreno S, Moros A, Chiva-Blanch G, et al. The EMT activator ZEB1 promotes tumor growth and determines differential response to chemotherapy in mantle cell lymphoma. Cell Death Different (2013) 21(2):247–57. doi: 10.1038/cdd.2013.123

CrossRef Full Text | Google Scholar

10. Stavropoulou V, Kaspar S, Brault L, Sanders MA, Juge S, Morettini S, et al. MLL-AF9 expression in hematopoietic stem cells drives a highly invasive AML expressing EMT-related genes linked to poor outcome. Cancer Cell (2016) 30(1):43–58. doi: 10.1016/j.ccell.2016.05.011

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Wu S, Du Y, Beckford J, Alachkar H. Upregulation of the EMT marker vimentin is associated with poor clinical outcome in acute myeloid leukemia. J Trans Med (2018) 16(1):1–9. doi: 10.1186/s12967-018-1539-y/figures/7

CrossRef Full Text | Google Scholar

12. Kahlert UD, Joseph J, Kruyt FAE. EMT- and MET-related processes in nonepithelial tumors: importance for disease progression, prognosis, and therapeutic opportunities. Mol Oncol (2017) 11(7):860–77. doi: 10.1002/1878-0261.12085

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Denecker G, Vandamme N, Akay Ö., Koludrovic D, Taminau J, Lemeire K, et al. Identification of a ZEB2-MITF-ZEB1 transcriptional network that controls melanogenesis and melanoma progression. Cell Death Different (2014) 21(8):1250–61. doi: 10.1038/cdd.2014.44

CrossRef Full Text | Google Scholar

14. Fallahi-Sichani M, Becker V, Izar B, Baker GJ, Lin J-R, Boswell SA, et al. Adaptive resistance of melanoma cells to RAF inhibition via reversible induction of a slowly dividing de-differentiated state. Mol Syst Biol (2017) 13:905. doi: 10.15252/msb

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Rambow F, Rogiers A, Marin-Bejar O, Aibar S, Femel J, Dewaele M, et al. Toward minimal residual disease-directed therapy in melanoma. Cell (2018) 174(4):843–55. doi: 10.1016/j.cell.2018.06.025

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Ramsdale R, Jorissen RN, Li FZ, Al-Obaidi S, Ward T, Sheppard KE, et al. The transcription cofactor c-JUN mediates phenotype switching and BRAF inhibitor resistance in melanoma. Sci Signaling (2015) 8(390):ra82. doi: 10.1126/scisignal.aab1111/suppl_file/8_ra82_sm.pdf

CrossRef Full Text | Google Scholar

17. Dupin E, le Douarin NM. Development of melanocyte precursors from the vertebrate neural crest. Oncogene 2003 22:20 (2003) 22(20):3016–23. doi: 10.1038/sj.onc.1206460

CrossRef Full Text | Google Scholar

18. Su Y, Wei W, Robert L, Xue M, Tsoi J, Garcia-Diaz A, et al. Single-cell analysis resolves the cell state transition and signaling dynamics associated with melanoma drug-induced resistance. Proc Natl Acad Sci USA (2017) 114(52):13679–84. doi: 10.1073/pnas.1712064115

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Vandamme N, Denecker G, Bruneel K, Blancke G, Akay O, Taminau J, et al. The EMT transcription factor ZEB2 promotes proliferation of primary and metastatic melanoma while suppressing an invasive, mesenchymal-like phenotype. Cancer Res (2020) 80(14):2983–95. doi: 10.1158/0008-5472.can-19-2373

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Wouters J, Kalender-Atak Z, Minnoye L, Spanier KI, de Waegeneer M, Bravo González-Blas C, et al. Robust gene expression programs underlie recurrent cell states and phenotype switching in melanoma. Nat Cell Biol (2020) 22(8):986–98. doi: 10.1038/s41556-020-0547-3

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Deng Y, Chakraborty P, Jolly MK, Levine H. A theoretical approach to coupling the epithelial-mesenchymal transition (EMT) to extracellular matrix (ECM) stiffness via LOXL2. Cancers 2021 (2021) 13(7):1609. doi: 10.3390/cancers13071609

CrossRef Full Text | Google Scholar

22. Fattet L, Jung HY, Matsumoto MW, Aubol BE, Kumar A, Adams JA, et al. Matrix rigidity controls epithelial-mesenchymal plasticity and tumor metastasis via a mechanoresponsive EPHA2/LYN complex. Dev Cell (2020) 54(3):302–16.e7. doi: 10.1016/j.devcel.2020.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Kumar S, Das A, Sen S. Extracellular matrix density promotes EMT by weakening cell-cell adhesions. Mol Biosyst (2014) 10(4):838–50. doi: 10.1039/c3mb70431a

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Bianchi A, Gervasi ME, Bakin A. Role of β5-integrin in epithelial-mesenchymal transition in response to TGF-β. Cell cycle (Georgetown, Tex.) (2010) 9(8):1647–59. doi: 10.4161/cc.9.8.11517

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Kilinc AN, Han S, Barrett LA, Anandasivam N, Nelson CM. Integrin-linked kinase tunes cell-cell and cell-matrix adhesions to regulate the switch between apoptosis and EMT downstream of TGFβ1. Mol Biol Cell (2021) 32(5):402–12. doi: 10.1091/mbc.E20-02-0092

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Kaur A, Ecker BL, Douglass SM, Kugel CH, Webster MR, Almeida F, et al. Remodeling of the collagen matrix in aging skin promotes melanoma metastasis and affects immune cell motility. Cancer Discov (2019) 9(1):64–81. doi: 10.1158/2159-8290.CD-18-0193

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Long JE, Wongchenko MJ, Nickles D, Chung W-J, Wang B, Riegler J, et al. Therapeutic resistance and susceptibility is shaped by cooperative multi-compartment tumor adaptation. Cell Death Different (2019) 26(11):2416–29. doi: 10.1038/s41418-019-0310-0

CrossRef Full Text | Google Scholar

28. Spoerri L, Tonnessen-Murray CA, Gunasingh G, Hill DS, Beaumont KA, Jurek RJ, et al. Phenotypic melanoma heterogeneity is regulated through cell-matrix interaction-dependent changes in tumor microarchitecture. BioRxiv (2021) 06(09)2020:141747. doi: 10.1101/2020.06.09.141747

CrossRef Full Text | Google Scholar

29. Dilshat R, Fock V, Kenny C, Gerritsen I, Lasseur RMJ, Travnickova J, et al. Mitf reprograms the extracellular matrix and focal adhesion in melanoma. ELife (2021) 10:e63093. doi: 10.7554/elife.63093

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Jensen C, Madsen DH, Hansen M, Schmidt H, Svane IM, Karsdal MA, et al. Non-invasive biomarkers derived from the extracellular matrix associate with response to immune checkpoint blockade (anti-CTLA-4) in metastatic melanoma patients. J Immunother Cancer (2018) 6(1):1–10. doi: 10.1186/s40425-018-0474-z/figures/4

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Kim MH, Kim J, Hong H, Lee S-H, Lee J-K, Jung E, et al. Actin remodeling confers BRAF inhibitor resistance to melanoma cells through YAP/TAZ activation. EMBO J (2016) 35(5):462–78. doi: 10.15252/embj.201592081

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Lal G, Contreras PG, Kulak M, Woodfield G, Bair T, Domann FE, et al. Human melanoma cells over-express extracellular matrix 1 (ECM1) which is regulated by TFAP2C. PLoS One (2013) 8(9):e73953. doi: 10.1371/journal.pone.0073953

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Miskolczi Z, Smith MP, Rowling EJ, Ferguson J, Barriuso J, Wellbrock C. Collagen abundance controls melanoma phenotypes through lineage-specific microenvironment sensing. Oncogene (2018) 37(23):3166. doi: 10.1038/S41388-018-0209-0

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Riesenberg S, Groetchen A, Siddaway R, Bald T, Reinhardt J, Smorra D, et al. MITF and c-jun antagonism interconnects melanoma dedifferentiation with pro-inflammatory cytokine responsiveness and myeloid cell recruitment. Nat Commun (2015) 6:8755. doi: 10.1007/s10911-010-9177-x

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Radisky ES, Radisky DC. Matrix metalloproteinase-induced epithelial-mesenchymal transition in breast cancer. J Mammary Gland Biol Neoplas (2010) 15(2):201–12. doi: 10.1007/s10911-010-9177-x/figures/3

CrossRef Full Text | Google Scholar

36. Suarez-Carmona M, Lesage J, Cataldo D, Gilles C. EMT and inflammation: inseparable actors of cancer progression. Mol Oncol (2017) 11(7):805–23. doi: 10.1002/1878-0261.12095

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Tripathi SC, Peters HL, Taguchi A, Katayama H, Wang H, Momin A, et al. Immunoproteasome deficiency is a feature of non-small cell lung cancer with a mesenchymal phenotype and is associated with a poor outcome. Proc Natl Acad Sci USA (2016) 113(11):E1555–64. doi: 10.1073/pnas.1521812113/suppl_file/pnas.1521812113.sd07.xlsx

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Jolly MK, Murphy RJ, Bhatia S, Whitfield HJ, Redfern A, Davis MJ, et al. Measuring and modelling the epithelial- mesenchymal hybrid state in cancer: Clinical implications. Cells Tiss Organ (2022) 211(2):110–33. doi: 10.1159/000515289

CrossRef Full Text | Google Scholar

39. Hoek KS, Eichhoff OM, Schlegel NC, Döbbeling U, Kobert N, Schaerer L, et al. In vivo switching of human melanoma cells between proliferative and invasive states. Cancer Res (2008) 68(3):650–6. doi: 10.1158/0008-5472.can-07-2491

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Tsoi J, Robert L, Paraiso K, Galvan C, Sheu KM, Lay J, et al. Multi-stage differentiation defines melanoma subtypes with differential vulnerability to drug-induced iron-dependent oxidative stress. Cancer Cell (2018) 33(5):890–904. doi: 10.1016/j.ccell.2018.03.017

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Jolly MK, Tripathi SC, Jia D, Mooney SM, Celiktas M, Hanash SM, et al. Stability of the hybrid epithelial/mesenchymal phenotype. Oncotarget (2016) 7(19):27067–84. doi: 10.18632/oncotarget.8166

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Pillai M, Jolly MK. Systems-level network modeling deciphers the master regulators of phenotypic plasticity and heterogeneity in melanoma. IScience (2021) 24(10):103111. doi: 10.1016/j.isci.2021.103111

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Bocci F, Tripathi SC, Vilchez Mercedes SA, George JT, Casabar JP, Wong PK, et al. NRF2 activates a partial epithelial-mesenchymal transition and is maximally present in a hybrid epithelial/mesenchymal phenotype. Integr Biol (2019) 11(6):251–63. doi: 10.1093/intbio/zyz021

CrossRef Full Text | Google Scholar

44. Caramel J, Papadogeorgakis E, Hill L, Browne GJ, Richard G, Wierinckx A, et al. A switch in the expression of embryonic EMT-inducers drives the development of malignant melanoma. Cancer Cell (2013) 24(4):466–80. doi: 10.1016/j.ccr.2013.08.018

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Gupta PB, Kuperwasser C, Brunet JP, Ramaswamy S, Kuo WL, Gray JW, et al. The melanocyte differentiation program predisposes to metastasis after neoplastic transformation. Nat Genet (2005) 37(10):1047–54. doi: 10.1038/ng1634

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Jessen C, Kreß JKC, Baluapuri A, Hufnagel A, Schmitz W, Kneitz S, et al. The transcription factor NRF2 enhances melanoma malignancy by blocking differentiation and inducing COX2 expression. Oncogene (2020) 39(44):6841. doi: 10.1038/S41388-020-01477-8

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Perotti V, Baldassari P, Molla A, Nicolini G, Bersani I, Grazia G, et al. An actionable axis linking NFATc2 to EZH2 controls the EMT-like program of melanoma cells. Oncogene 2019 38:22 (2019) 38(22):4384–96. doi: 10.1038/s41388-019-0729-2

CrossRef Full Text | Google Scholar

48. Subbalakshmi AR, Kundnani D, Biswas K, Ghosh A, Hanash SM, Tripathi SC, et al. NFATc acts as a non-canonical phenotypic stability factor for a hybrid Epithelial/Mesenchymal phenotype. Front Oncol (2020) 10:553342. doi: 10.3389/fonc.2020.553342

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Subbalakshmi AR, Sahoo S, Biswas K, Jolly MK. A computational systems biology approach identifies SLUG as a mediator of partial epithelial-mesenchymal transition (EMT). Cells Tiss Organ (2022) 211(6):105–18. doi: 10.1159/000512520

CrossRef Full Text | Google Scholar

50. Vandewalle C, Comijn J, de Craene B, Vermassen P, Bruyneel E, Andersen H, et al. SIP1/ZEB2 induces EMT by repressing genes of different epithelial cell–cell junctions. Nucleic Acids Res (2005) 33(20):6566–78. doi: 10.1093/nar/gki965

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Subbalakshmi AR, Sahoo S, McMullen I, Saxena AN, Venugopal SK, Somarelli JA, et al. KLF4 induces mesenchymal–epithelial transition (MET) by suppressing multiple EMT-inducing transcription factors. Cancers (2021) 13(20):5135. doi: 10.3390/cancers13205135

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Zhang D, Lin J, Chao Y, Zhang L, Jin L, Li N, et al. Regulation of the adaptation to ER stress by KLF4 facilitates melanoma cell metastasis via upregulating NUCB2 expression. J Exp Clin Cancer Res (2018) 37(1):1–14. doi: 10.1186/s13046-018-0842-z

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Campbell NR, Rao A, Hunter M, Sznurkowska MK, Briker L, Zhang M, et al. Cooperation between melanoma cell states promotes metastasis through heterotypic cluster formation. Dev Cell (2021) 56(20):2808–25.e10. doi: 10.1016/j.devcel.2021.08.018

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Dimitrova Y, Gruber AJ, Mittal N, Ghosh S, Dimitriades B, Mathow D, et al. TFAP2A is a component of the ZEB1/2 network that regulates TGFB1-induced epithelial to mesenchymal transition. Biol Direct (2017) 12(1):8. doi: 10.1186/s13062-017-0180-7/figures/5

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Byers LA, Diao L, Wang J, Saintigny P, Girard L, Peyton M, et al. An epithelial-mesenchymal transition gene signature predicts resistance to EGFR and PI3K inhibitors and identifies axl as a therapeutic target for overcoming EGFR inhibitor resistance. Clin Cancer Res (2013) 19(1):279–90. doi: 10.1158/1078-0432.ccr-12-1558

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci USA (2005) 102(43):15545–550. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Tan TZ, Miow QH, Miki Y, Noda T, Mori S, Huang RY, et al. Epithelial-mesenchymal transition spectrum quantification and its efficacy in deciphering survival and drug responses of cancer patients. EMBO Mol Med (2014) 6(10):1279–93. doi: 10.15252/emmm.201404208

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Barretina J, Caponigro G, Stransky N, Venkatesan K, Margolin A, Kim S, et al. The cancer cell line encyclopedia enables predictive modelling of anticancer drug sensitivity. Nature (2012) 483(7391):603–7. doi: 10.1038/nature11003

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Barbie DA, Tamayo P, Boehm JS, Kim SY, Moody SE, Dunn IF, et al. Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nat (2009) 462(7269):108–12. doi: 10.1038/nature08460

CrossRef Full Text | Google Scholar

60. Verfaillie A, Imrichova H, Atak ZK, Dewaele M, Rambow F, Hulselmans G, et al. Decoding the regulatory landscape of melanoma reveals TEADS as regulators of the invasive cell state. Nat Commun (2015) 6:6683. doi: 10.1038/ncomms7683

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database hallmark gene set collection. Cell Syst (2015) 1(6):417–25. doi: 10.1016/j.cels.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Muralidharan S, Sahoo S, Saha A, Chandran S, Majumdar SS, Mandal S, et al. Quantifying the patterns of metabolic plasticity and heterogeneity along the epithelial–Hybrid–Mesenchymal spectrum in cancer. Biomolecules (2022) 12(2):297. doi: 10.3390/biom12020297

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Jia D, Paudel BB, Hayford CE, Hardeman KN, Levine H, Onuchic JN, et al. Drug-tolerant idling melanoma cells exhibit theory-predicted metabolic low-low phenotype. Front Oncol (2020) 10:1426. doi: 10.3389/fonc.2020.01426

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Yu L, Lu M, Jia D, Ma J, Ben-Jacob E, Levine H, et al. Modeling the genetic regulation of cancer metabolism: Interplay between glycolysis and oxidative phosphorylation. Cancer Res (2017) 77(7):1564–74. doi: 10.1158/0008-5472.can-16-2074/652665/am/modeling-the-genetic-regulation-of-cancer

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Foroutan M, Bhuva DD, Lyu R, Horan K, Cursons J, Davis MJ. Single sample scoring of molecular phenotypes. BMC Bioinf (2018) 19(1):1–10. doi: 10.1186/s12859-018-2435-4/figures/2

CrossRef Full Text | Google Scholar

66. Taylor R. Interpretation of the correlation coefficient: A basic review. J Diagn Med Sonogr (1990) 6(1):35–9. doi: 10.1177/875647939000600106

CrossRef Full Text | Google Scholar

67. Tang Y, Durand S, Dalle S, Caramel J. EMT-inducing transcription factors, drivers of melanoma phenotype switching, and resistance to treatment. Cancers (2020) 12(8):2154. doi: 10.3390/cancers12082154

CrossRef Full Text | Google Scholar

68. Wessely A, Steeb T, Berking C, Heppt MV. How neural crest transcription factors contribute to melanoma heterogeneity, cellular plasticity, and treatment resistance. Int J Mol Sci (2021) 22(11):5761. doi: 10.3390/ijms22115761

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Abildgaard C, Guldberg P. Molecular drivers of cellular metabolic reprogramming in melanoma. Trends Mol Med (2015) 21(3):164–71. doi: 10.1016/j.molmed.2014.12.007

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Bettum IJ, Gorad SS, Barkovskaya A, Pettersen S, Moestue SA, Vasiliauskaite K, et al. Metabolic reprogramming supports the invasive phenotype in malignant melanoma. Cancer Lett (2015) 366(1):71–83. doi: 10.1016/j.canlet.2015.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Gelato KA, Schöckel L, Klingbeil O, Rückert T, Lesche R, Toedling J, et al. Super-enhancers define a proliferative PGC-1α-expressing melanoma subgroup sensitive to BET inhibition. Oncogene 2018 37:4 (2017) 37(4):512–21. doi: 10.1038/onc.2017.325

CrossRef Full Text | Google Scholar

72. Laurenzana A, Chillà A, Luciani C, Peppicelli S, Biagioni A, Bianchini F, et al. uPA/uPAR system activation drives a glycolytic phenotype in melanoma cells. Int J Cancer (2017) 141(6):1190–200. doi: 10.1002/ijc.30817

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Haq R, Shoag J, Andreu-Perez P, Yokoyama S, Edelman H, Rowe GC, et al. Oncogenic BRAF regulates oxidative metabolism via PGC1α and MITF. Cancer Cell (2013) 23(3):302–15. doi: 10.1016/j.ccr.2013.02.003

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Vazquez F, Lim JH, Chim H, Bhalla K, Girnun G, Pierce K, et al. PGC1α expression defines a subset of human melanoma tumors with increased mitochondrial capacity and resistance to oxidative stress. Cancer Cell (2013) 23(3):287–301. doi: 10.1016/j.ccr.2012.11.020

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Vivas-García Y, Falletta P, Liebing J, Louphrasitthiphol P, Feng Y, Chauhan J, et al. Lineage-restricted regulation of SCD and fatty acid saturation by MITF controls melanoma phenotypic plasticity. Mol Cell (2020) 77(1):120–137.e9. doi: 10.1016/j.molcel.2019.10.014

PubMed Abstract | CrossRef Full Text | Google Scholar

76. Klinke DJ, Torang A. An unsupervised strategy for identifying epithelial-mesenchymal transition state metrics in breast cancer and melanoma. IScience (2020) 23(5):101080. doi: 10.1016/j.isci.2020.101080

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Cook DP, Vanderhyden BC. Context specificity of the EMT transcriptional response. Nat Commun (2020) 11(1):2142. doi: 10.1038/s41467-020-16066-2

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Sahoo S, Ashraf B, Duddu AS, Biddle A, Jolly MK. Interconnected high-dimensional landscapes of epithelial–mesenchymal plasticity and stemness in cancer. Clin Exp Metast (2022) 1:1–12. doi: 10.1007/S10585-021-10139-2

CrossRef Full Text | Google Scholar

79. Norgard RJ, Pitarresi JR, Maddipati R, Aiello-Couzo NM, Balli D, Li J, et al. Calcium signaling induces a partial EMT. EMBO Rep (2021) 22(9):e51872. doi: 10.15252/embr.202051872

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Rambow F, Marine JC, Goding CR. Melanoma plasticity and phenotypic diversity: Therapeutic barriers and opportunities. Genes Dev (2019) 33(19–20):1295–1318. doi: 10.1101/gad.329771.119

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Youssef KK, Nieto MA. Glucose metabolism takes center stage in epithelial-mesenchymal plasticity. Dev Cell (2020) 53(2):133–5. doi: 10.1016/J.Devcel.2020.03.021

PubMed Abstract | CrossRef Full Text | Google Scholar

82. Zheng X, Boyer L, Jin M, Mertens J, Kim Y, Ma L, et al. Metabolic reprogramming during neuronal differentiation from aerobic glycolysis to neuronal oxidative phosphorylation. ELife (2016) 5(JUN2016):e13374. doi: 10.7554/elife.13374

PubMed Abstract | CrossRef Full Text | Google Scholar

83. Vidács DL, Veréb Z, Bozó R, Flink LB, Polyánka H, Németh IB, et al. Phenotypic plasticity of melanocytes derived from human adult skin. Pigment Cell & Melanoma Res (2022) 35(1):38–51. doi: 10.1111/pcmr.13012

PubMed Abstract | CrossRef Full Text | Google Scholar

84. Fedele M, Cerchia L, Pegoraro S, Sgarra R, Manfioletti G. Proneural-mesenchymal transition: Phenotypic plasticity to acquire multitherapy resistance in glioblastoma. Int J Mol Sci (2019) 20(11):2746. doi: 10.3390/ijms20112746

CrossRef Full Text | Google Scholar

85. Couturier CP, Ayyadhury S, Le PU, Nadaf J, Monlong J, Riva G, et al. Single-cell RNA-seq reveals that glioblastoma recapitulates a normal neurodevelopmental hierarchy. Nat Commun (2020) 11(1):3406. doi: 10.1038/s41467-020-17186-5

PubMed Abstract | CrossRef Full Text | Google Scholar

86. Phillips HS, Kharbanda S, Chen R, Forrest WF, Soriano RH, Wu TD, et al. Molecular subclasses of high-grade glioma predict prognosis, delineate a pattern of disease progression, and resemble stages in neurogenesis. Cancer Cell (2006) 9(3):157–73. doi: 10.1016/j.ccr.2006.02.019

PubMed Abstract | CrossRef Full Text | Google Scholar

87. Wang L, Babikir H, Müller S, Yagnik G, Shamardani K, Catalan F, et al. The phenotypes of proliferating glioblastoma cells reside on a single axis of variation. Cancer Discov (2019) 9(12):1708–19. doi: 10.1158/2159-8290.CD-19-0329

PubMed Abstract | CrossRef Full Text | Google Scholar

88. Liu S, Cong Y, Wang D, Sun Y, Deng L, Liu Y, et al. Breast cancer stem cells transition between epithelial and mesenchymal states reflective of their normal counterparts. Stem Cell Rep (2013) 2(1):78–91. doi: 10.1016/j.stemcr.2013.11.009

CrossRef Full Text | Google Scholar

89. Camand E, Peglion F, Osmani N, Sanson M, Etienne-Manneville S. N-cadherin expression level modulates integrin-mediated polarity and strongly impacts on the speed and directionality of glial cell migration. J Cell Sci (2012) 125(4):844–57. doi: 10.1242/jcs.087668

PubMed Abstract | CrossRef Full Text | Google Scholar

90. Schaeffer D, Somarelli JA, Hanna G, Palmer GM, Garcia-Blanco MA. Cellular migration and invasion uncoupled: Increased migration is not an inexorable consequence of epithelial-to-Mesenchymal transition. Mol Cell Biol (2014) 34(18):3486–99. doi: 10.1128/mcb.00694-14/suppl_file/zmb999100574so2.pdf

PubMed Abstract | CrossRef Full Text | Google Scholar

91. Goyal Y, Dardani IP, Busch GT, Emert B, Fingerman D, Kaur A, et al. Pre-determined diversity in resistant fates emerges from homogenous cells after anti-cancer drug treatment. BioRxiv (2021) 2021:471833. doi: 10.1101/2021.12.08.471833

CrossRef Full Text | Google Scholar

92. Su Y, Ko ME, Cheng H, Zhu R, Xue M, Wang J, et al. Multi-omic single-cell snapshots reveal multiple independent trajectories to drug tolerance in a melanoma cell line. Nat Commun (2020) 11(1):2345. doi: 10.1038/s41467-020-15956-9

PubMed Abstract | CrossRef Full Text | Google Scholar

93. Kim YJ, Sheu KM, Tsoi J, Abril-Rodriguez G, Medina E, Grasso CS, et al. Melanoma dedifferentiation induced by IFN-γ epigenetic remodeling in response to anti–PD-1 therapy. J Clin Invest (2021) 131(12). doi: 10.1172/jci145859

PubMed Abstract | CrossRef Full Text | Google Scholar

94. Kavran AJ, Stuart SA, Hayashi KR, Basken JM, Brandhuber BJ, Ahn NG. Intermittent treatment of BRAF V600E melanoma cells delays resistance by adaptive resensitization to drug rechallenge. Proc Natl Acad Sci (2022) 119(12):e2113535119. doi: 10.1073/PNAS.2113535119

PubMed Abstract | CrossRef Full Text | Google Scholar

95. Goldman A, Majumder B, Dhawan A, Ravi S, Goldman D, Kohandel M, et al. Temporally sequenced anticancer drugs overcome adaptive resistance by targeting a vulnerable chemotherapy-induced phenotypic transition. Nat Commun (2015) 6(1):2015. doi: 10.1038/ncomms7139

CrossRef Full Text | Google Scholar

96. Grosse-Wilde A, D’Hérouël AF, McIntosh E, Ertaylan G, Skupin A, Kuestner RE, et al. Stemness of the hybrid Epithelial/Mesenchymal state in breast cancer and its association with poor survival. PLoS One (2015) 10(5):e0126522. doi: 10.1371/journal.pone.0126522

PubMed Abstract | CrossRef Full Text | Google Scholar

97. Kröger C, Afeyan A, Mraz J, Eaton EN, Reinhardt F, Khodor YL, et al. Acquisition of a hybrid E/M state is essential for tumorigenicity of basal breast cancer cells. Proc Natl Acad Sci USA (2019) 116(15):7353–62. doi: 10.1073/pnas.1812876116

PubMed Abstract | CrossRef Full Text | Google Scholar

98. Pasani S, Sahoo S, Jolly MK. Hybrid E/M phenotype(s) and stemness: A mechanistic connection embedded in network topology. J Clin Med (2021) 10(1):60. doi: 10.3390/jcm10010060

CrossRef Full Text | Google Scholar