Abstract
Purpose: Understanding tumor heterogeneity is an important challenge in current cancer research. Transcription and epigenetic profiling of cultured melanoma cells have defined at least two distinct cell phenotypes characterized by distinctive gene expression signatures associated with high or low/absent expression of microphthalmia-associated transcription factor (MITF). Nevertheless, heterogeneity of cell populations and gene expression in primary human tumors is much less well characterized.
Experimental Design: We performed single-cell gene expression analyses on 472 cells isolated from needle biopsies of 5 primary human melanomas, 4 superficial spreading, and one acral melanoma. The expression of MITF-high and MITF-low signature genes was assessed and compared to investigate intra- and intertumoral heterogeneity and correlated gene expression profiles.
Results: Single-cell gene expression analyses revealed varying degrees of intra- and intertumor heterogeneity conferred by the variable expression of distinct sets of genes in different tumors. Expression of MITF partially correlated with that of its known target genes, while SOX10 expression correlated best with PAX3 and ZEB2. Nevertheless, cells simultaneously expressing MITF-high and MITF-low signature genes were observed both by single-cell analyses and RNAscope.
Conclusions: Single-cell analyses can be performed on limiting numbers of cells from primary human melanomas revealing their heterogeneity. Although tumors comprised variable proportions of cells with the MITF-high and MITF-low gene expression signatures characteristic of melanoma cultures, primary tumors also comprised cells expressing markers of both signatures defining a novel cell state in tumors in vivo. Clin Cancer Res; 23(22); 7097–107. ©2017 AACR.
In melanoma, diagnosis and characterization of the primary tumor rely heavily on histology analysis and immunostaining with a limited number of markers. Additional molecular techniques are therefore required to refine and complement these approaches. Little is known concerning cell heterogeneity in primary lesions that cannot be assessed by currently used techniques or whether the presence and prevalence of cells with MITF-high or MITF-low gene signatures can be used to predict future tumor evolution. Here, we show that single cell gene expression profiling on needle biopsies from primary melanoma lesions coupled with RNAscope hybridization provides unprecedented characterization of their heterogeneity and gene expression profiles. Our approach opens up new possibilities to correlate the presence of cells expressing genes of the MITF-high or MITF-low signatures in the primary lesion with the future evolution of the disease and patient outcome.
Introduction
Understanding tumor heterogeneity is an important challenge in current cancer research. In melanoma, cells with differing invasive, proliferative, and tumor-initiating potential have been defined on the basis of the characteristics of established lines and short-term primary cultures. Meta-analysis of gene expression and/or epigenetic profiling of hundreds of melanoma cell lines or primary cultures identified gene expression signatures and profiles of open and active chromatin regions that characterize two distinctive cell states, proliferative or invasive (1–3), also designated as MITF-high/AXL-NF-κB-low or AXL-NF-κB-high/MITF-low (4, 5) based on the expression level of microphthalmia-associated transcription factor (MITF). Although both cell types proliferate in vitro, MITF-low cells display in addition elevated motile and invasive capacity and higher tumor-forming capacity when reinjected as xenografts. Furthermore, MITF-low cells with tumor-initiating properties arise spontaneously in cultures of MITF-high cells (6).
Decoding the epigenetic landscapes of multiple primary cultures of metastatic melanoma showed that MITF and SOX10 are major drivers of the MITF-high state. MITF and SOX10 interact physically and functionally with the PBAF and NURF remodeling complexes to establish the epigenetic landscape of MITF-high cells and coregulate genes involved in cell cycle, metabolism, and invasion (7–12). In contrast, the TEAD and JUN transcription factors are important drivers of the MITF-low state (3). Immunostaining of human melanoma sections showed heterogeneous MITF expression with MITF-high cells intermixed with MITF-low cells that express high levels of the transcription factors POU3F2/BRN2 (13) and/or GLI2 (14). Meta-analyses of melanoma transcriptomes also define a class of tumors with high expression of MITF and its target genes (15).
The abovementioned evidence is consistent with at least two broad groups of MITF-high or MITF-low/negative cell types with distinct phenotypes and gene expression signatures. Melanoma development seems to involve dynamic switching between these two phenotypic states (16) driven by the microenvironment and an integrated stress response involving inhibition of the translation initiation factor eIF2B and repression of MITF expression by ATF4 (17).
We previously used single-cell gene expression analyses to analyze heterogeneity in melanoma cells cultured under different conditions in vitro or as xenograft tumors in mice (18). Here, we extended this approach to primary human melanoma whose molecular analysis is limited by the necessity to preserve intact biopsies for histopathology analysis. As single-cell gene expression analyses require only small numbers of cells, these can be isolated as microbiopsies without compromising subsequent histologic evaluation.
Materials and Methods
Patients and tumors
This study was conducted following approval by the University of Strasbourg Medical faculty ethics board and good clinical practice. Biopsies from epidermal melanomas were obtained from patients with written informed consent. Following biopsies, the pathologist's reports confirmed each sample as malignant melanoma.
Cell lines
501 Mel cells were obtained from Dr. Colin Goding and MM099 and MM047 from Dr. Jean-Christophe Marine. All cell lines were mycoplasma tested using an EZ-PCR Kit (BI Industries), and RNA was prepared after less than 5 passages following defreezing.
Single-cell qRT-PCR from primary melanoma
Needle biopsies were made from patient lesions immediately before surgical resection and immediately dissociated as previously described for mouse xenograft tumors by incubation in HBSS (Sigma-Aldrich) supplemented with collagenase IV (10 mg/mL, Eurobio), dispase II (1 mg/mL, Sigma), DNase I (200 IU/mL, Roche), 75 μmol/L CaCl2, and 125 μmol/L MgCl2, for 45 minutes at 37°C. Cells were filtered through a 100-μm pore size filter (Dutscher). After centrifugation at 100 × g, for 7 minutes at 4°C, pellets were dissolved with HBSS buffer containing 200 IU/mL of DNase I and 125 μmol/L of MgCl2. Cell suspensions were centrifuged at 100 × g for 7 minutes and resuspended in red blood cell lysis buffer (0.15 mol/L NH4Cl, 10 mmol/L KHCO3, and 100 μmol/L EDTA). After centrifugation at 100 × g for 3 minutes, the cells were finally resuspended in HBSS buffer. Cells were captured using the C1 Single-Cell Auto Prep System using the 10-17 μm array (Fluidigm), followed by reverse transcription and preamplification according to the Fluidigm's instructions. Single-cell gene expression experiments were performed using Fluidigm's M96 quantitative PCR (qPCR) DynamicArray microfluidic chips. A 2.25 μL aliquot of amplified cDNA was mixed with 2.5 μL of 2× SsoFast EvaGreen Supermix with Low ROX (Bio-Rad, PN 172-5211) and 0.25 μL of “sample loading agent” (Fluidigm, PN 100-3738), and then inserted into one of the chip “sample” inlets. A total of 100 μmol/L of mixed forward and reverse primers were diluted at 1:10 ratios with TE. Then, 2.5 μL of diluted primers was mixed with 2.5 μL of Fluidigm “Assay Loading Reagent” and individually inserted into the chip “assay” inlets. Samples and probes were loaded into M96 chips using a HX IFC Controller (Fluidigm) and then transferred to a Biomark real-time PCR reader (Fluidigm) following the manufacturer's instructions.
Single-cell qRT-PCR data analysis
Basic analyses and heatmap generation were performed as described previously (18). Initial data analysis of the cycle threshold (Ct) values was done with the “Fluidigm Real-time PCR analysis” software and further data analysis and graphics were performed using R software. Complement of Ct values was defined as expression threshold et = Cmax – Ct = 30 – Ct (19). Absent values were replaced by 0. Then, et values were displayed on heatmap images following clustering of both genes and cells calculated with the unsupervised unweighted pair group method with arithmetic mean (UPGMA) and Euclidean distance. To compare data from different experiments, we converted et values by their ranks using the following method: in each cell, the values xi were sorted in ascending order x1 < x2 <…< xi <…< xN (with N the number of genes). Then, each value xi > 0 was replaced by its rank ri: rN = N, rN−1 = N−1, … ri = i. All values xi = 0 (absent values) were kept as 0. These ranks were used to perform clustering, whereas original et values were displayed on the heatmap images. These ranks were also used in the boxplot figures and to compute t-distributed stochastic neighbor embedding (tSNE). tSNE was performed using the sincell Bioconductor package and a perplexity value of 5 or 50. For correlation of gene expression, Pearson correlation coefficients were calculated on each pair of genes and displayed on heatmap images. Clustering of genes based on correlation coefficients was performed with the unsupervised UPGMA and distance = 1−|Pearson correlation coefficient|. To measure the variability of gene expression within a tumor, we calculated for each gene the median absolute deviation (MAD) of et values in all cells. These values were displayed on a heatmap following clustering of genes and tumors calculated with UPGMA and Euclidean distance. Bulk RNA sequencing (RNA-seq) from 501Mel, MM099, and MM047 cells was performed and analyzed as described previously (9).
RNAscope
mRNAs for MITF, BIRC3, SOX10, and SOX9 in sections from human melanomas and cultured cells were detected with RNAscope assay (Advanced Cell Diagnostics, ACD) according to the manufacturer's protocols. Briefly, patient sections were deparaffinized, incubated with hydrogen peroxide at room temperature for 10 minutes, boiled with target retrieval reagent for 15 minutes, and then treated with protease plus reagent at 40°C for 30 minutes. The sections were hybridized with Hs-MITF probe (ACD, cat. no. 310951), Hs-BIRC3-C2 probe (ACD, cat. no. 417591-C2), Hs-SOX10 probe (ACD, cat. no. 484121), and Hs-SOX9-C2 probe (ACD, cat. no. 404221-C2) at 40°C for 2 hours. Hybridization signals were amplified and visualized with RNAscope Multiplex Fluorescent Reagent Kit v2 (ACD, cat. no. 323100). Images were captured with a fluorescent (Leica DM4000B) or a confocal (Leica DMI6000) microscope.
Results
Analyses of gene expression in single cells from primary melanoma
To investigate and quantify heterogeneity in primary human melanoma tumors by single-cell gene expression, we analyzed 472 cells from needle microbiopsies of primary lesions of 5 patients. Four patients displayed superficial spreading melanoma on the back, chest or leg and the fifth displayed an acral melanoma on the sole of the foot. In each case, needle biopsies were removed immediately before their surgical resection. The biopsy was rapidly dissociated into a single-cell suspension and cells were captured as described previously (Fig. 1A; ref. 18). Gene expression was analyzed by qRT-PCR with primers for genes of the MITF-high or MITF-low expression signatures (2), genes involved in EMT, stem cells markers (20), cell-cycle markers, signaling molecules, and MITF cofactors (Supplementary Dataset S1; refs. 9–11).
Comparing expression of the MITF-high and MITF-low signature genes in bulk RNA-seq data from MITF-high 501Mel cells and MITF-low MM099 and MM047 cells (3) confirmed their distinctive profiles. MITF and its target genes were strongly expressed in 501Mel cells, while MITF-low signature genes were highly expressed in MM099 and MM047 cells (Supplementary Fig. S1A). We note however that the levels of expression of MITF-low signature genes differed between the MM099 and MM047 cells.
Tumor 1 came from an 82-year-old male patient with a large prominent melanoma on the chest (Table 1) that originated from a pigmented area and had grown persistently over a 6-year period. Despite the large size of the lesion, MRI scan showed no metastases at the time of surgery and in the following 2 years. Ninety-six cells were captured and analyzed by unsupervised clustering of both genes and cells. Three major classes of cells with different MITF expression levels were observed (Fig. 1B). The first displayed high MITF expression along with many of its target genes, such as CEACAM1, TYR, KAZN, and MLANA. Within this class, a subset of cells showed much lower expression of a group of genes comprising SOX10, PAX3, ZEB2, SLUG and CDK2. The second displayed generally lower MITF expression along with MITF-negative cells accompanied by lowered expression of MITF target genes and genes coregulated with SOX10 and CDK2. The third showed mixed levels of MITF expression, but were distinguished from the second class by the stronger expression of the SOX10-CDK2 gene cluster. The low or absent expression of CDK2 and CCND1 and the strongly reduced expression of signaling molecules like BRAF and PTEN, the chromatin remodelers SMARCA4 (BRG1) and BPTF suggested that class II cells corresponded to non- or slow cycling cells. In contrast, expression of MITF-low signature genes, such as GLI2, ZEB1, EGF, MET, MYOF, and BIRC3 was generally low in all classes. Thus, no strong upregulation of “invasive”-type markers was seen in the MITF-low cells.
Tumor . | Mutation . | Location . | Size (mm) . | No. of cells . |
---|---|---|---|---|
1 | BRAF, NRAS, KIT WT | Upper chest | 60 × 70 | 96 |
2 | BRAFV600E | Lower back | 20 × 14 | 96 |
3 | BRAF, NRAS, KIT WT | Lower leg | 60 × 67 | 94 |
4 | NRASQ61L | Upper back | 72 × 42 | 96 |
5(Acral) | BRAFV600E | Sole of foot | 50 × 50 | 90 |
6 | ND | Right thigh | 13 × 11 | N/A |
7 | ND | Upper back | 16 × 12 | N/A |
Tumor . | Mutation . | Location . | Size (mm) . | No. of cells . |
---|---|---|---|---|
1 | BRAF, NRAS, KIT WT | Upper chest | 60 × 70 | 96 |
2 | BRAFV600E | Lower back | 20 × 14 | 96 |
3 | BRAF, NRAS, KIT WT | Lower leg | 60 × 67 | 94 |
4 | NRASQ61L | Upper back | 72 × 42 | 96 |
5(Acral) | BRAFV600E | Sole of foot | 50 × 50 | 90 |
6 | ND | Right thigh | 13 × 11 | N/A |
7 | ND | Upper back | 16 × 12 | N/A |
Abbreviations: NA, nonapplicable; ND, not determined.
Tumor 2 came from a 41-year-old male patient with a primary lesion on the lower back (Table 1). At the time of biopsy, this patient showed internal metastases and died 14 months later. A majority of MITF-low/negative cells was observed forming two clusters with differential expression of a set of genes exemplified by RELB or OCA2 (Fig. 2). Only a small population of MITF-high cells was observed where a subset of its target genes was also upregulated. In contrast, other genes coregulated with MITF in tumor 1, such as, KAZN, MYO1D, and CEACAM1, although more strongly expressed in the MITF-high cells, were also well expressed in MITF-low/negative cells. Only a subset of invasion-associated genes, such as BIRC3, THBS1, COL13A1, and EGFR, was expressed in a majority of cells, whereas ZEB1, CYR61, and ITGA2 were only weakly expressed.
Comparison of the ranked expression profiles of tumors 1 and 2 highlighted their important differences (Supplementary Fig. S2). The expression of a group of genes that clustered with MITF was clearly enriched in tumor 1, whereas multiple MITF-low genes were enriched in tumor 2.
Tumor 3 came from an 85-year-old female patient with a melanoma on the lower leg just below the knee (Table 1) and was characterized by a lack of prominent MITF-high cells and little observable heterogeneity (Supplementary Fig. S3A). Despite the low MITF expression, high expression of MLANA and TYR was observed, while other target genes such as CEACAM1 and BIRC7 were weakly or not expressed.
We analyzed cells from two additional primary lesions. Tumor 4 from a 70-year-old patient with a lesion on the upper back showed low or negative MITF expression, but high expression of MLANA and TYR (Supplementary Fig. S3B). The acral melanoma tumor from the sole of the foot of a 64-year-old female patient displayed a majority of cells with elevated MITF expression and a small population of MITF-low/negative cells (Supplementary Fig. S3C).
In these experiments, almost all of the cells analyzed could be defined as melanoma based on their gene expression profiles and not keratinocytes or infiltrating immune cells. Melanoma cells expressed one or a combination of melanocyte specific genes, such as MITF, SOX10, PMEL, or MLANA. Even cells from the MITF-low/negative tumors 3 and 4 expressed either PMEL, MLANA, or TYR or combinations of these markers confirming their identity as melanoma. This may be explained by the fact that tumors were biopsied at the thickest area and that cells were size selected by the microfluidics to capture 10 to 17 μm cells, thus eliminating any smaller immune infiltrate cells.
Variable expression of genes of the MITF-high and MITF-low signatures contributes to inter- and intratumoral heterogeneity
Nonlinear dimensionality reduction using tSNE using two values of perplexity showed that each tumor segregated separately highlighting their distinctive expression profiles (Fig. 3A and B). Both values of perplexity indicated that tumors 1 and 2 were the most closely related, with a small number of cells from tumor 1 segregating with tumor 2. In addition to intertumoral heterogeneity, this analysis also revealed intratumoral heterogeneity that was particularly evident in tumors 1 and 2 using a perplexity value of 5 (Fig. 3A). The other tumors showed lesser heterogeneity in accordance with what was observed in the heatmaps.
We next asked whether intratumoral heterogeneity was due to the variable expression of a common or distinct set of genes in all tumors. To measure the variability of expression, we calculated the MAD of the expression threshold values in all cells for each gene. Genes with the most variable expression differed from tumor to tumor (Fig. 3C). For example, in tumor 2, this approach defined MITF and its target genes among those that showed highest variability, whereas in the acral tumor, ZEB1, TNC, and DKK3 were among the most variable. In agreement with the heatmap and tSNE analysis showing tumor 3 as the most homogenous, only a few genes showed strong variability in this tumor.
Comparing the ranked expression values between tumors identified genes whose expression showed different degrees of variability. VIM and NANOG were highly expressed in each tumor with little intratumor variation, while S100A and SOX10 expression showed more variability (Fig. 3D). In contrast, MITF, BIRC3, and ZEB1 showed much higher intratumor variability and differing expression between tumors (Fig. 3D). In contrast to NANOG, POU5F1 expression was much more variable showing differential regulation of these pluripotency markers in melanoma. Inter and intratumour heterogeneity therefore resulted from the variable expression distinct sets of MITF-high and MITF-low signature genes.
Cells in primary melanoma simultaneously expressing genes of the MITF-high and MITF-low signatures
We next generated heatmaps and clustering for representative MITF-high and MITF-low signature genes in each tumor and calculated the Pearson correlation coefficients on each pair of genes. In tumor 1, expression of FST, CEACAM1, KAZN, OCA2, and MLANA correlated with MITF, similar to what we previously reported in MITF-high 501Mel cells (Fig. 4A; ref. 18). Their expression also correlated in the subgroup of MITF-high cells in tumor 2 (Fig. 4B). In tumor 3, MLANA expression remained somewhat correlated with MITF, whereas that of KAZN and OCA2 was independent of MITF (Fig. 4C). Moreover, the MLANA/MITF ratio was much higher in MITF-low tumor 3 than in tumors 1 and 2. A similar situation was seen in MITF-low tumor 4 (Fig. 4D). In the acral tumor, expression of MLANA, OCA2, and KAZN correlated with MITF, but correlation between MITF and CEACAM1 was low.
SOX9 expression correlated with that of MITF in all tumors except tumor 4. Similarly, POU3F2 correlated with MITF in all tumors except tumor 2 where it was not expressed in the majority of cells. In contrast, SOX10 expression did not always correlate with MITF. For example, in tumor 1, MITF-high/SOX10-low cells and SOX10-high/MITF-low cells were seen. SOX10 was, on the other hand, most highly correlated with PAX3 and ZEB2 in a majority of tumors.
In tumor 2, ZEB1 was expressed in the MITF-high and low/negative subpopulations contrary to their anticorrelation in cultured cells (21–23). Although ZEB1 was most highly expressed in MITF-low tumor 3, its expression was much lower in MITF-low tumor 4. Similar observations were seen for GLI2 and BIRC3. In tumor 3, GLI2 expression was highest in cells expressing most MITF.
Cells expressing both SOX9 and SOX10 were observed in almost all tumors. Even in tumor 4 that showed higher and homogenous SOX9 expression, SOX10 was coexpressed in many cells. The primers used properly discriminated these two closely related genes in qRT-PCR on RNA from the MITF-high and MITF-low cell lines (Supplementary Fig. S1A; refs. 3, 9). Using the SOX10 primers, high expression in 501Mel cells, but no significant expression in MM047 and MM099 cells was seen (Supplementary Fig. S1B). In contrast, high SOX9 expression was detected in the MITF-low lines, but not 501Mel cells (Supplementary Fig. S1B). Note that as seen in the RNA-seq data (Supplementary Fig. S1A), SOX9 expression was higher in MM047 than in MM099, but much higher in both than in 501Mel. Thus, these primers properly discriminated SOX9 and SOX10.
Cells coexpressing genes of the MITF-high and MITF-low signatures in primary melanoma and cutaneous metastases
As the above results showing coexpression of MITF-high and MITF-low signature genes in the same cells contradicted longstanding observations from cultured melanoma cells, we sought confirmation by an independent technique. Immunostaining is limited by the availability of reliable and specific antibodies, whereas RNAscope is a sensitive and specific technique to investigate gene expression in tumor sections. We employed probes for MITF and SOX10 as markers for the MITF-high state and BIRC3 and SOX9 as markers for the MITF-low state. MITF and SOX10 expression could readily be detected in 501Mel cells, whereas no SOX9 or BIRC3 was seen (Supplementary Fig. S4). Similarly, SOX9 and BIRC3 were detected in MM099 cells, whereas no signal for SOX10 or MITF was seen. These results confirmed the specificity of these probes.
RNAscope is sensitive to RNA degradation, and only the more recent tumors 4 and 5 gave signals. However, as cells coexpressing these markers were seen in all 5 analyzed tumors, we reasoned that it should be a general feature, and we hence performed additional experiments on sections from two more recently isolated primary superficial spreading melanomas (tumors 6 and 7; Table 1).
In tumor 5, regions with cells expressing only MITF or BIRC3 were detected along with regions comprising numerous cells expressing both MITF and BIRC3 (Fig. 5A). In tumor 6, adjacent regions expressing either MITF or BIRC3 were observed, whereas in other regions, numerous coexpressing cells were detected (Fig. 5B). Similarly, in tumors 4 and 6, regions with unique or coexpression of MITF and SOX9 were observed (Fig. 5C). With the combinations of SOX10-BIRC3 and SOX10-SOX9, regions of tumor 6 where cells expressing either one or other were observed, whereas other regions of the tumor comprised cells expressing combinations of SOX10 and BIRC3 or SOX9 and SOX10 (Fig. 5D and E). In other regions, intermixing of adjacent cells expressing either SOX9 or SOX10 was seen (Fig. 5F). These data highlighted the regional heterogeneity of tumors where cells with different profiles were segregated or intermixed. Additional confocal microscopy confirmed that cells expressing both SOX10 and BIRC3 or SOX9 could be clearly identified (Fig. 5G and H). These data showed that primary tumors comprised cells expressing either MITF-high and MITF-low markers and cells that simultaneously expressed markers of both states.
In sections from cutaneous metastases, cells expressing MITF together with either BIRC3 or SOX9 were readily identified (Supplementary Fig. S5A). Cells expressing SOX9 and SOX10 were observed together with other regions with predominantly SOX9-expressing cells or SOX10-expressing cells (Supplementary Fig. S5B and S5C). Interestingly, we also detected gland-like structures comprising many cells coexpressing SOX9 and SOX10, or SOX10 and BIRC3 (Supplementary Fig. S6A–S6C). Histology analyses of sections from the same tumors showed the presence of numerous sweat glands corresponding to those labeled by the SOX10-SOX9/BIRC3-labeled cells (Supplementary Fig. S6D). This was further evidenced by acquisition of brightfield and fluorescent images of the same sections, where labeling of the glandular structures can be clearly seen (Supplementary Fig. S6E). As these structures did not label with MITF, they likely did not correspond to melanoma cells, but rather sweat gland cells previously shown to express SOX9 and SOX10 (24, 25). These cells exist in clearly defined structures that distinguish them from disperse cells of the melanoma.
Together, the above data indicate that both primary melanoma and cutaneous metastases comprise mixtures of MITF-high and MITF-low cells and cells expressing markers of both cell states.
Discussion
Single-cell expression analysis reveals inter- and intratumor heterogeneity in primary melanoma
We show that single-cell gene expression can be used to investigate heterogeneity in primary human melanoma using small numbers of cells without compromising histologic evaluation. Our results revealed an important intertumor heterogeneity where each tumor showed distinctive profiles of gene expression. This heterogeneity resulted from the variable expression of subsets of genes, with each tumor displaying a discrete set of most variable genes. Other genes showed differing degrees of variability, some showing very homogenous expression and others intermediate levels of variability. The signature of each tumor resulted from the combinatorial expression profiles of the most variable genes rather than overall changes in the expression of a majority of genes. Differential expression of MITF and its target genes in addition to that of markers like ZEB1, THBS1, and BIRC3 strongly contributes to the specific expression signature of each tumor.
Comparisons between tumors did not use absolute expression values, but ranked data. For example, VIM was one of the highest ranked genes in a large majority of cells from all tumors with little variation from cell to cell within a given tumor and among tumors. In contrast, MITF was highly ranked in cells where its expression was high, but was ranked low in cells where its expression was low. In this way, the expression of each gene is corrected relative to the others in a given cell and tumor. This analysis minimizes potential batch effects, as it makes no assumption about absolute expression values.
Different tumors also displayed varying degrees of intratumor heterogeneity. The MITF-low tumors 3–4 were more homogenous than tumors 1 and 2 whose heterogeneity was augmented by the presence of MITF-high and MITF-low cells. Tumor 5 showed many MITF-high cells, but heterogeneity was mainly the consequence of a subpopulation of cells expressing higher levels of MITF-low genes, such as ZEB1, MYOF, and BIRC3.
Tirosh and colleagues analyzed gene expression by RNA-seq in single cells from a series of metastatic melanomas and one primary melanoma (26). They identified a set of genes exemplified by PMEL, TYR, and MLANA that correlated with MITF similar to what we describe here. Their analyses showed that metastatic melanomas comprised a mixture of cells with MITF-high or MITF-low signatures and that the proportion of each cell type varied from tumor to tumor. We show that primary melanomas also comprise varying ratios of MITF-high and MITF-low cells with in several cases a majority of MITF-low cells. Tirosh and colleagues however did not note the presence of cells coexpressing MITF-high and MITF-low markers.
While this study was in progress, Fluidigm reported that their arrays were prone to capture of cells as doublets. In control experiments using a mix of 501Mel and 1205Lu cells, we detected around 20% of doublet cells. These experiments were done with high cell densities, whereas experiments from primary tumors were performed with limiting numbers of cells such that we often had to load the cell suspensions twice to fill the arrays. Although we cannot exclude the possibility of doublet cells in the analyses described here, it is important to note that our findings are not based on observations from a small number of cells, but from a large number of cells often in multiple tumor samples. Moreover, if a large number of doublets were present, it would be impossible to cluster them, for example, by MITF expression or to see the correlated gene expression profiles that we describe. Finally, the presence of cells simultaneously expressing the MITF-high and MITF-low markers in primary tumors that could have been due to the artifactual capture of cells of both types was confirmed by RNAscope, an independent technique.
Primary melanomas comprise subpopulations of cells coexpressing genes of the MITF-high and MITF-low signatures
A second important finding of this study is that primary tumors comprise subpopulations of cells simultaneously expressing genes of the MITF-high and MITF-low signatures. Cells with an MITF-high signature were present in tumors 1 and 2, where MITF expression correlated with that of its target genes arguing that functional MITF protein was present. Nevertheless, unlike established MITF-high cells (14), MITF-high primary tumor cells expressed GLI2 in tumors 2 and 3, ZEB1 in tumors 2 and 5, and more significantly BIRC3 or SOX9 in tumors 1, 2, and 5. Indeed, SOX9 expression was highest in the MITF-high subpopulation of tumor 2, and MITF and SOX9 were coexpressed in tumors 1, 2, 3, and 5.
The simultaneous expression of markers of both cell states was confirmed by RNAscope, where we identified cells expressing combinations of MITF/SOX10 and SOX9/BIRC3 both in primary tumors and in cutaneous metastases. RNAscope also identified cells expressing only markers of one state or the other. These cells were in distinct and sometimes adjacent regions or intermixed highlighting the spatial heterogeneity within tumors. Given these considerations, it is evident that the single-cell biopsies represent only a fraction of the entire population and their composition depends on the tumor zone from which they were taken. Nevertheless, many of the profiles we found by single-cell analyses were observed in cells from multiple tumors and must therefore represent more general phenomenon and not just chance observations.
SOX9 and MITF/SOX10 expression is generally mutually exclusive in most cultured melanoma cells (3, 27, 28), although some lines, such as 1205Lu or A375M, showed expression of both SOX9 and SOX10 mRNAs (28). Epigenetic profiling of the SOX9, SOX10, and MITF gene loci in primary MITF-high and MITF-low cultures revealed reciprocal open and closed expression states and H3K27ac levels (3). Immunostaining of melanomas in situ also showed exclusive SOX9 and SOX10 protein expression (27). In contrast, our single-cell analyses and RNAscope identified primary tumor cells expressing MITF and SOX9 or SOX10 and SOX9. Although a translational control of the SOX9 protein would explain the discrepancy with immunostaining, the coexpression of these genes in primary tumors indicates a cell state not typically seen cultured melanoma cells.
Similarly, POU3F2 has been reported to repress MITF expression in cultured melanoma lines (13, 29) yet MITF and POU3F2 were often coexpressed in cells from several primary tumors. On the other hand, Wellbrock and colleagues reported that oncogenic BRAF activates MITF expression via POU3F2 binding to its promoter (30). At present, we cannot distinguish between a positive regulation of MITF by POU3F2 in response to activation of MAP kinase signaling in primary tumors and activation of POU3F2 by MITF. Further experiments will be required to determine the mechanisms that account for this observed correlation.
In primary tumors, MITF and SOX10 were poorly correlated, with SOX10 being highly expressed even in the MITF-low tumors consistent with its role in promoting and maintaining melanoma growth (31). SOX10 expression best correlated with PAX3 and ZEB2 consistent with a transcription-regulatory network active during melanocyte development where they bind and activate the MITF M-promoter (23). This hierarchical relationship is not always conserved in melanoma, as evidenced by populations of PAX3-SOX10-ZEB2-high/MITF-low cells in tumors 2 and 3 and MITF-high/PAX3-SOX10-ZEB2-low cells in tumor 1.
The above observations support the idea that primary and metastatic melanomas comprise not only MITF-high and MITF-low cells, but also subpopulations expressing markers of both signatures. Combinations of the three cell populations may be adjacent or intermixed contributing to the spatial heterogeneity of the tumors. Whether the expression of markers of both signatures represents cells in transition between the two phenotypes or a more stable state specific to tumors in vivo remains to be determined.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Authors' Contributions
Conception and design: M. Ennen, C. Keime, I. Davidson
Development of methodology: M. Ennen, C. Keime, C. Thibault-Carpentier
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M. Ennen, A. Kieny, C. Thibault-Carpentier, F. Margerin-Schaller, D. Lipsker
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M. Ennen, C. Keime, G. Davidson, C. Vagne, I. Davidson
Writing, review, and/or revision of the manuscript: M. Ennen, C. Keime, C. Vagne, D. Lipsker, I. Davidson
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M. Ennen, S. Coassolo, F. Margerin-Schaller
Study supervision: M. Ennen, I. Davidson
Other (design and execution of RNAscope experiments and image analysis): G. Gambi
Acknowledgments
We thank, J-C. Marine for cell lines MM047 and MM099, all the staff of the IGBMC common facilities and the staff of the Strasbourg Hospital Dermatology Clinic, the IGBMC high-throughput sequencing facility, a member of the “France Génomique” consortium (ANR10-INBS-09-08). I. Davidson is an “équipe labellisée” of the Ligue Nationale contre le Cancer.
Grant Support
This work was supported by Ligue National Contre le Cancer (Irwin Davidson), Association pour la Recherché Contre le Cancer, Institut National du Cancer (PAIR 13-002), and Agence National de la Recherche (ANR10-Labex-0030-INRT). M. Ennen was supported by a fellowship from the Ligue Nationale Contre le Cancer.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.