|Home | About | Journals | Submit | Contact Us | Français|
Transcriptomic phenotypes defined for melanoma have been reported to correlate with sensitivity to various drugs. In this study, we aimed to define a minimal signature that could be used to distinguish melanoma sub-types in vitro, and to determine suitable drugs by which these sub-types can be targeted. By using primary melanoma cell lines, as well as commercially available melanoma cell lines, we find that the evaluation of MLANA and INHBA expression is as capable as one based on a combined analysis performed with genes for stemness, EMT and invasion/proliferation, in identifying melanoma subtypes that differ in their sensitivity to molecularly targeted drugs. Using this approach, we find that 75% of melanoma cell lines can be treated with either the MEK inhibitor AZD6244 or the HSP90 inhibitor 17AAG.
Melanoma is a heterogeneous tumor, as can be defined either by mutations or by transcriptomic analyses. More than 80% of all melanoma cells harbor a mutated gene within the MAPK pathway, all which are mutually exclusive1, 2. Patients with BRAF mutations are candidates for Vemurafenib, a selective BRAF inhibitor3, and the administration of a MEK inhibitor (Trametinib) in these patients has shown to increase survival time4, demonstrating the critical function of the MAPK pathway. In addition, close to 60% of melanoma contain mutations within the CDKN2A pathway related genes. Although these gene mutations are also mutually exclusive within the pathway, they overlap with MAPK related mutations, indicating the presence of at least two distinct pathways that are deregulated in this tumor type. However, the mutational landscape as defined by these and other mutations has been reported not to correlate with phenotypes that are based on the transcriptome of the same tumors1, 2. That gene expression based phenotypes reveal heterogeneity existing among melanoma tumors is a well-established fact1, 2, 5–7. Current data indicates the existence of at least 3 melanoma sub-groups as can be defined by transcriptomic analyses, with a preferential up-regulation of either differentiation, proliferation, or immune related genes1, 2. Personalized treatment, such as CTLA4 and PD-1 blockade, a most promising treatment for advanced melanoma8, is likely to benefit from these classifications as they might help design such therapy on an individual basis9. However, melanoma cells per se show transcriptome based heterogeneity when studied in vitro as well, and this is important as it has been shown to correlate with drug sensitivity10. Thus, complete cure of melanoma might have to be pursued from two arms, one based on in vivo effects, and therefore the microenvironment of the tumor cells, and the other based on in vitro defined mechanisms related to drug sensitivity. Switching between phenotypes has been observed for melanoma in vivo, as well as in vitro, and can result in altered drug sensitivity, irrespective of mutational status for melanoma cells5, 11; suggesting that a combination treatment consisting of two drugs affecting each one of the two phenotypes might be ultimately necessary for efficient tumor eradication. Drug resistance has traditionally been associated with a stem-cell-like phenotype and a more mesenchymal state, as compared to an epithelial one in almost all tumors12, 13. The “stemness” of melanoma cells is thought to exist in a dynamic state, as melanoma cells with stem-cell markers have been reported to arise from cells lacking these markers, stem-cell like melanoma cells might not be required for tumor initiation, and the stem- and non-stem-cell phenotypes can switch dynamically7, 14. Similar to the plasticity reported for stemness in melanoma, epithelial-to-mesenchymal transition (EMT) has also been implied as a factor driving metastasis and causing drug resistance in melanoma15, 16. Thirdly, melanoma specific in vitro analyses aiming primarily to explain differing metastatic potential generated gene lists that could robustly classify melanoma cell lines or tumor cells into either a “proliferative” or an “invasive” character6, 17, 18; and switching among these phenotypes was also shown as well5, 14. Metastasis, drug resistance and invasion have all been attributed both to the stem cell phenotype, as well as to EMT and have been observed in those classes defined by various researchers for melanoma. In this study we asked if there was overlap between melanoma cells classified by these various signatures and if such melanoma subtypes would have unique responses to targeted therapy. We also aimed the identification of a minimum number of biomarkers that could help classify melanoma subtypes with differing drug sensitivity profiles.
Primary melanoma cell lines used in this study were previously described19, 20. All cell lines (primary and commercial) were incubated in humidified incubators supplied with 5% CO2 and cultured in RPMI1640 medium supplemented with 10% FBS, 1% 10K/10K Penicillin-Streptomycin (P/S Solution, Lonza, Basel, Switzerland) and 1% 200mM L-Glutamine (L-Glu, Lonza, Basel, Switzerland). Two exceptions where commercial cell lines Mewo and SK-MEL-28 which were cultured in DMEM with the same supplements. Melanocytic progeny of primary melanoma cell lines (PrMCLs) were investigated by staining for S100, MART-1, gp100 and GD2 expression using polyclonal rabbit anti-S100, monoclonal A-103, HMB-45 (Invitrogen, Carlsbad, CA, USA) and R-24 (a kind gift of Dr. Aliza Adler) antibodies, and flow cytometry based (FACS) detection of HMW-MMA (CSPG4) (R&D Systems Minneapolis, MN, USA) (Supplementary Table 1). Samples were analyzed by FACS on a Becton Dickinson LSR II cytometer. Data analysis was performed using FCS express V.3 (De Novo Software).
Quantitative profiling of INHBA and MLANA by q-RT-PCR was done using a 7500 Real-time PCR System (Applied Bio-systems, Carlsbad, CA, USA) using SYBR® Green PCR Master Mix (Catalog number: 4309155, Applied Biosystems, Carlsbad, CA). PCR reactions were run at cycling conditions according to manufacturer’s instructions. GAPDH was used as the endogenous control in all reactions. PCR primers were designed using NCBI Primer–Blast Tool (http://www.ncbi.nlm.nih.gov/tools/primer-blast/). Primer sequences were: MLANA_F1-ATTAAGGAAGGTGTCCTGTGCC, MLANA_R1-CAGCCGTGGTGTAAGAGTGG, INHBA_F1-GTTTGCCGAGTCAGGAACAGC, INHBA_R1-GAAGAGCCAGACTTCTGCACG, GAPDH_F1-AAAGGGTCATCATCTCTGCCCCC and GAPDH_R1-GGCATTGCTGATGATCTTGAGGCTG. Relative expression values of INHBA and MLANA were calculated by the ddCt method as previously described21.
AZD0530, 17AAG and AZD6244 were purchased from Selleck Chemicals. Drugs were screened on transparent 96-well-plates in which cancer cell lines were seeded at 1000 cells/well 96 hours before cytotoxicity measurements. Drugs were screened using 10 different doses ranging from 0.001µM to 50µM except 10a for which 6 different doses ranging from 0.625µM to 20µM were tested. All stock solutions were prepared in 0.1% DMSO. IC50 values were calculated as previously described22. Cytotoxicity was measured with CellTiter-Glo® Luminescent Cell Viability Assay (Promega, Madison WI, USA) after 72 hours of drug treatment.
Normalized activity area or IC50 values of 24 different commercial drugs were downloaded from CCLE database (https://portals.broadinstitute.org/ccle/home) for in silico drug screening analyses.
EMT and CSC scores were calculated based on an R script23 (http://www.R-project.org/) which runs a principal component analyses (PCA) with the following probeset expression values: EMT scoring was based on probeset expression values of E-CAD (201131_s_at) and VIM (201426_s_at); CSC scoring was using CD271 (205858_at), CD20 (228599_at, 228592_at, 217418_x_at, 210356_x_at), CD133 (204304_s_at) and ALDH1 (212224_at), all which were previously published melanoma stem-cell markers24. CCLE and primary melanoma cell lines were classified as PRO (proliferative), INT (intermediate) and INV (invasive) according to Pro (r) and Inv (r) values generated using the online HOPP tool and the associated gene list18. CCLE cell lines with a Pro (r) value of more than 0.652 and an Inv (r) value of less than 0.171 were designated as “PRO”, those with a Pro (r) value between 0.041 to 0.630 and an Inv (r) value between 0.254 to 0.600 were designated as “INT” and cells with a Pro (r) value of less than −0.200 and an Inv (r) value of more than 0.445 were designated as “INV” cells. For primary cell lines, designations were as follows: PRO: Pro (r) > 0.533 and Inv (r) < 0.162; INT: Pro (r): 0.186 to 0.588 and Inv (r): 0.221 to 0.581; INV: Pro (r) < 0.03 and Inv (r) > 0.554.
E-MTAB-783 and GSE36133 commercial cancer cell line datasets were normalized separately using RMA normalization without background correction. Differentially expressed probesets were determined between melanoma cell lines and other solid tumor cell lines separately in both datasets according to p-value and fold change (FC) cutoffs, 0.0001 (Bonferroni corrected) and 4 (in either direction), respectively. Significant probesets common to both datasets were determined using VENNY “Oliveros, J.C. (2007–15) An interactive tool for comparing lists with Venn Diagrams. (http://bioinfogp.cnb.csic.es/tools/venny/index.html).” These 223 probesets constitute the “gene list #1”. Standardized expression values of these probesets obtained from the C-MAP microarray data of primary melanoma cell lines20 were used for hierarchical clustering analysis using Euclidean distance as the similarity metric and complete linkage as the clustering method in Cluster3.0 software. Results were visualized in JavaTreeView. Two main clusters were generated with this approach. One of the clusters showed upregulation of genes which were also upregulated in commercial melanoma cell lines compared with other solid cancer cell lines. Second cluster of primary cell lines showed upregulation of genes that were downregulated in commercial melanoma cell lines, when compared to other solid cancer cell lines. In the second round of analyses, the most differentially expressed genes between two main clusters were determined using the same significance cutoffs used in generating “gene list#1”. This resulted in “gene list#2”. INHBA and MLANA were selected among the most significantly differentially expressed genes from gene list#2 for further validation of phenotypic features of the two main clusters.
Melanoma cell lines derived from primary tumor tissues were used throughout this study. Cells were initially tested for melanoma markers S100, MART-1 and gp100 to differentiate them from non-melanoma cells20. As a further measure of specificity, gene expression data obtained from all cell lines via C-MAP technology was used to generate two gene lists representing “melanoma” (MEL) and “fibroblasts” (FIB) like qualities using a previously defined algorithm18. Gene expression based correlation measures obtained for melanoma cell lines, as well as for normal skin fibroblasts and cancer-associated fibroblasts show that this approach is an effective means for separating melanoma from fibroblasts. As shown in Figure 1, gene expression based analysis revealed that only 4 out of 36 melanoma cell lines show an inverse correlation with the MEL signature, while only 6 out of 65 fibroblasts have a FIBr value lower than 0.5. The 4 melanoma cell lines that can be thus designated “fibroblast-like” are however, clearly melanoma cell lines as judged by morphology. Other primary cell lines which were negative for S100, MART-1 or gp100, and which lacked a melanoma-like morphology were also identified as non-melanoma and fibroblast-like by this analysis and were excluded from our studies (not shown). Similar to our analyses, “fibroblast-like” cells, as defined by low MELr and high FIBr values, are present in two out of three datasets generated from melanoma primary cell lines (Supplementary Fig. 1a). We consistently observed that such cells grow much slower, compared to those with high MELr scores, explaining why they are extremely rare among commercially available melanoma cell lines (CMCL, Supplementary Fig. 1b). The average daily proliferation rate of cells with low INVMELr scores (<0.4) was 18%, compared to 65% for cells with high INVMELr scores (>0.4); similarly cells with low FIBr scores (<0), had an average proliferation rate of 67%, compared to 20% for those with high FIBr scores (>0).
Drug resistance in melanoma has been reported to associate with qualities determined by phenotypes such as EMT, stemness, and measures pertaining to the cell’s proliferative and invasive capacity. We evaluated each primary melanoma cell line, as well as CMCLs for each of these 3 qualities in silico, to determine any overlap between these definitions, and if one or a combination of these could be used to better explain sensitivity to a molecularly targeted drug. For this purpose, cells were first classified into proliferative (PRO), invasive (INV) or intermediate (INT) categories using an online algorithm based on a previously published gene list18. Each cell was then assigned a score corresponding to its epithelial/mesenchymal state (EMT score), in addition to a score reflecting how stem-like the cell was (CSC score; see methods). Not surprisingly, only one CMCL, and 3 PrMCLs had both epithelial and CSC-like qualities, as these two phenotypes are expected to show mutual exclusivity (Quadrant 1, Figure 2a). Similarly, most INV cells were more mesenchymal than epithelial, while PRO cells were more epithelial. However, while the INT group among CMCLs was mostly of a CSC-like and mesenchymal phenotype (Q3, Figure 2a), the INT group among primary melanoma cell lines had more non-CSC-like and epithelial features (Q4, Figure 2b). We analyzed if cells harboring BRAF (V600) mutations were preferentially of a particular transcriptomic subtype. CSC and EMT scores were not significantly different between mutant and WT cells. However, the PRO scores of mutant cells were significantly larger than that of WT cells, and the inverse was observed for INV scores; p<0.01 and p=0.03, respectively. Suggesting BRAF mutations were not a primary determinant of stemness or the EMT state of cells, but were related to PRO and INV scores.
We then asked if the classification based on the combination of these three transcriptomic groups could predict drug sensitivity better than that which could be predicted by using the “PRO/INV/INT” classification by itself. For this purpose, we correlated INV and PRO values obtained for melanoma cell lines for which there was cytotoxicity data in the CCLE database22. The strongest direct correlation between PRO values and simultaneous negative correlation with INV values was observed for cytotoxicity data from two MEK inhibitors AZD6244 and PD0325901, in line with the correlation of these phenotypes with BRAF mutation status. The exact opposite pattern was observed most prominently for three other drugs, the HSP90 inhibitor 17-AAG, the SRC-inhibitor AZD0530, and a cMET inhibitor PF-2341066 (Supplementary Table 2). As HSP90 and MEK inhibitors have previously been shown to affect opposite phenotypes25, and similar findings exist for SRC and MEK inhibitor combinations10, we decided to focus on 17AAG, AZD6244 and AZD0530. We first determined in silico, how melanoma cells, as defined by the combination of 3 scores, behaved in response to these drugs and whether this classification would help improve over the PRO/INT/INV classification. As shown in Figure 3, PRO cells were primarily sensitive to the MEK inhibitor AZD6244, while INV cells were to 17AAG (HSP90 inhibitor) and AZD0530 (SRK inhibitor). By further classifying cells into quadrants based on the added EMT and stemness score we were unable to identify a subpopulation that was more sensitive to any of the three drugs tested, with the most resistant and sensitive groups identified based on quarters showing similar values as those obtained for PRO, INT and INV groups (Figure 3). This indicates that classifying cells based on a combination of EMT and CSC markers is able to identify distinct cell populations that differ in their chemosensitivity profiles but does not help identify cells that are more sensitive to these drugs, compared to the classification generated by the PRO/INV score.
Next, we wanted to validate the in silico cytotoxicity correlations in vitro using PrMCLs. In contrast to what we observed for CMCLs, PRO, INT and INV groups were not as clearly separated into quadrants in analyses of PrMCLs (Figure 2). Most PRO cells were observed in Q4 (Epithelial and non-CSC-like), together with INT cells, while INV cells were present in quarters 1, 2 and 3, which is a different distribution compared to the classification we observed for CMCLs. Cytotoxicity profiles of PrMCLs however, were quite different between Q-based and PRO/INV based classes (Figure 4). For 17-AAG, although INV cells were the most sensitive, INT cells were generally more resistant than PRO cells. A single cell tested in Q2 was the most sensitive, supporting this observation. For AZD6244, PRO and Q4 cells were relatively sensitive. Although a single cell in Q2 was also sensitive, as this is a sole observation, we do not think it indicates an exception. AZD0530, similar to its in silico data, did not show much variation among groups, although INV cells and that in Q2 were somewhat more sensitive. We then tested a novel drug 10a, that has been shown to be a specific SRC inhibitor26. 10a was effective in particular on PRO cells and those in Q4 (epithelial, non-CSC-like), but also on cells in Q2 (mesenchymal, CSC-like) (Figure 4). There were no statistically significant differences among either quadrant, or PRO/INV based classes for 17AAG and AZD0530 IC50 values. However, there were significant differences for AZD6244 between Q3 vs. Q4 (p=0.03, t-test), and for 10a, in quadrant based (p=0.026) or PRO/INT/INV based (p<0.001) classifications as identified by 1-way ANOVA. This observation reveals that PRO/INV based and quadrant based classifications are distinct, in that the cells categorized by both, differ in their responses to targeted therapy, especially to the SRC inhibitor 10a. Interestingly the sensitivity profile of 10a is the complete inverse of that observed for AZD0530. Our data shows that SRC inhibition by 10a shows a profile almost identical to that of AZD6244 for PrMCLs (Figure 4). As MEK and SRC are shared in most pathways, this suggests the SRC-MAPK pathway is active and responsible for growth primarily in PRO cells and those in quarters 2 and 4. AZD0530 shows minimal activity in both CMCLs and PrMCLs at the concentrations tested and does not happen to have distinctly different effects on subgroups of melanoma. That its activity profile is different from 10a suggests either that AZD0530 is unable to inhibit SRC in these cells or that 10a is inhibiting a non-SRC protein.
We also compared cytotoxicity responses among cells with and without BRAF mutations. We did not observe any significant differences in IC50 values for any of the drugs tested between the two groups, although BRAF mutant cells were generally more sensitive to the MEK inhibitor AZD6244 (data not shown).
All phenotypical categorizations described thus far are based on the evaluation of a relatively large number of genes. However, as the various categories defined by the three different phenotypical markers generate clusters with similar drug sensitivity profiles, we argued that such profiles could be reflected by a much smaller number of gene expression patterns. As a number of genes used in these classifiers are related to differentiation markers of melanoma, as a discovery strategy, we sought to identify genes whose expression could differentiate melanoma from non-melanoma tumors by comparing gene expression profiles of tumors within E-MTAB-783 and GSE36133 datasets. 223 probesets with a Bonferroni corrected p value of 0.0001 or below and a difference of more than 4 folds that were common to both datasets were used to cluster PrMCLs in silico (Supplementary Figure 2). This analysis revealed two distinct clusters (A and B) which were then used to discover genes that were most differentially expressed between these two clusters (Supplementary Figure 2). MLANA and INHBA were selected as representatives of the most upregulated genes in each of the two distinct clusters, especially since they have been reported to identify opposite phenotypes in melanoma: MLANA is critical for the regulation of mammalian pigmentation27, while INHBA is reported to decrease melanin content of B16 murine melanoma cell lines28. Furthermore, INHBA is a member of transforming growth factor superfamily which induces a mesenchymal phenotype by increasing N-Cadherin expression in esophageal carcinoma cell lines29. These two genes are inversely expressed in both commercial as well as in PrMCLs (Supplementary Figure 3). Although most INHBA expressing cells are low expressors of MLANA and vice versa (thus in quadrants 1 and 3, Supplementary Figure 4), there are cells that express both or neither. Both genes have been reported to show a strong correlation between mRNA and protein levels17, 30; we therefore, believe the mRNA based phenotype should be reflected at the protein level as well. As expected, MLANA and INHBA mRNA expression correlated tightly with MITF and AXL, two genes that are known to play critical roles in phenotype switching in melanoma31, 32 (Supplementary Figure 5). In line with BRAF mutations not correlating with previously reported transcriptomic signatures for melanoma, BRAF mutation status of cells does not seem to correlate with MLANA or INHBA either (Supplementary Figure 4).
We observed that MLANA and INHBA expression based classification was successful in revealing drug sensitivity differences for both CMCLs as well as PrMCLs: We combined cells in Q2, 3 and 4 as their sensitivity profiles were similar. We then compared their sensitivity to cells within Q1 (MLANA-high, INHBA-low). This revealed the expected trends for all three drugs. In silico analysis showed that Q1 cells (MLANA-high, INHBA-low) were relatively sensitive to AZD6244 (MEK inhibitor), when they were relatively resistant to 17AAG (HSP90 inhibitor) and AZD0530 (SRC inhibitor), compared to cells in Q2, 3, 4 (Figure 5). We then validated these findings in vitro for PrMCLs, and included the SRC inhibitor 10a to our analysis. As shown in Figure 6, AZD6244 and 10a are both more effective on Q1 (MLANA-high, INHBA-low) cells, while 17AAG and AZD0539 showed the exact opposite trend. The differential effect of AZD6244 was statistically significant. To test if the cells that were resistant to AZD6244 were also sensitive to 17-AAG or 10a, we plotted the drug sensitivity profiles of cells for either of the two drugs. This analysis, however, revealed that while half of either CMCLs or PrMCLs followed this pattern, about a quarter were sensitive to both, and another quarter of the cells were resistant to both drugs analyzed (Supplementary Figure 6). As previously observed, 10a and AZD6244 sensitivity profiles were highly similar (Supplementary Figure 6c). We therefore conclude that phenotyping based on MLANA and INHBA will be able to identify cells that are sensitive to either AZD6244/10a or 17AAG, in about 50% of the cells tested. 25% of cells will be sensitive to both drugs, and a quarter might be resistant to both.
Plasticity of stemness, EMT and INV/PRO phenotypes in melanoma is clearly related to changes in the signaling pathway choices of these cells. We hypothesized that different melanoma cells, as characterized by gene expression profiles reported to relate to distinct phenotypes, would primarily rely on one particular signaling pathway which could be blocked by one of the drugs generated to target specific signaling pathways. We also hypothesized that the various phenotypes defined for melanoma could overlap. Despite the limited number of cell lines used here, we have been able to support our second hypothesis, as we were able to show that cells grouped by utilizing multiple gene lists are not necessarily better characterized in terms of their drug sensitivity profiles, compared to using any given one. If each of these different phenotypes arises due to the utilization of a distinct pathway, then the phenotypic switch should relate to a shift in the signaling pathways the cell uses, resulting in altered drug sensitivity to a molecularly targeted drug. In line with this argument, we were able to generate partial support for our first hypothesis: although we were able to identify distinct drugs what were more effective on cells, as defined by a given phenotype, and that drugs that target different molecules within the same pathway show similar outcomes (shown for 10a and AZD6244), we observed that drugs showing opposite correlations with a given phenotype (AZD6244 or 10a versus 17AAG) did not necessarily have mutually exclusive actions on cells. This is in contrast to the fact that the biomarkers identifying these cells are almost mutually exclusive (MLANA and INHBA). The explanation for this could be that mRNA based biomarkers, however good they may be in identifying cells with homogenous transcriptomic phenotypes, are not necessarily identifying cells that utilize the same survival pathways. This discrepancy could then likely be due to the fact that simultaneous but distinct pathways are active in a given cell and they have to be targeted at the same time for effective control of cell growth. Despite not being statistically significant, the similarity in sensitivity profiles we observe for 17AAG and AZD0530, two drugs that target non-overlapping signal pathways, support this line of thought. It is possible that the ultimate aim of the cell is to alter its transcriptomic profile, generating a given functional phenotype. If various pathways can be utilized for this purpose, then the cell will use one out of the several ones possible. On the other hand, the fact that melanoma cells that develop resistance to B-RAF inhibitors re-activate the MAPK pathway suggests there are situations where the cell relies solely on one pathway33.
Targeted therapy via agents like Dabrafenib (BRAF inhibitor) and Trametinib (MEK inhibitor) significantly increase overall survival of metastatic melanoma patients with BRAF V600E or V600K mutations34, 35. That MEK and RAF inhibitor cytotoxicity profiles are similar according to our in silico analysis is also supportive of this. One reason we could not correlate BRAF mutations to sensitivity to drugs tested here could be the limited number of mutant cells we studied. This is supported by the fact that BRAF V600 mutant cells have much higher PRO scores compared to INV scores, and the fact that PRO/INV scores do correlate with drug sensitivity; especially to that of the MEK inhibitor AZD6244. However, melanoma subclasses as defined by transcriptomic analyses cannot be explained based on their mutational profile1, 2. It is therefore likely that a simple classification based on MLANA and INHBA can help identify tumors that will benefit from one of the drugs analyzed in this study better by itself or in combination with mutation information.
Our results indicate that melanoma cells are possibly utilizing multiple survival pathways and that the best personalized treatment approach will be one that first reveals the signaling pathways critical for the survival in each cell, followed by the utilization of a drug or a combination of drugs that can target all of these pathways at once. Our data shows that a profiling based on gene expression testing of only MLANA and INHBA is as capable of classifying in vitro responses to molecularly targeted drugs. Larger validation experiments in vivo and tests of drug synergy will need to be performed to confirm this observation.