Mapping heterogeneity in patient-derived melanoma cultures by single-cell RNA-seq

Abstract

Recent technological advances in single-cell genomics make it possible to analyze cellular heterogeneity of tumor samples. Here, we applied single-cell RNA-seq to measure the transcriptomes of 307 single cells cultured from three biopsies of three different patients with a BRAF/NRAS wild type, BRAF mutant/NRAS wild type and BRAF wild type/NRAS mutant melanoma metastasis, respectively. Analysis based on self-organizing maps identified sub-populations defined by multiple gene expression modules involved in proliferation, oxidative phosphorylation, pigmentation and cellular stroma. Gene expression modules had prognostic relevance when compared with gene expression data from published melanoma samples and patient survival data. We surveyed kinome expression patterns across sub-populations of the BRAF/NRAS wild type sample and found that CDK4 and CDK2 were consistently highly expressed in the majority of cells, suggesting that these kinases might be involved in melanoma progression. Treatment of cells with the CDK4 inhibitor palbociclib restricted cell proliferation to a similar, and in some cases greater, extent than MAPK inhibitors. Finally, we identified a low abundant sub-population in this sample that highly expressed a module containing ABC transporter ABCB5, surface markers CD271 and CD133, and multiple aldehyde dehydrogenases (ALDHs). Patient-derived cultures of the BRAF mutant/NRAS wild type and BRAF wild type/NRAS mutant metastases showed more homogeneous single-cell gene expression patterns with gene expression modules for proliferation and ABC transporters. Taken together, our results describe an intertumor and intratumor heterogeneity in melanoma short-term cultures which might be relevant for patient survival, and suggest promising targets for new treatment approaches in melanoma therapy.

Keywords: melanoma; single cell transcriptome sequencing; stem cells.

Figure 1. Dissecting heterogeneity in a BRAF/NRAS wild type patient-derived melanoma culture revealed by single-cell RNA-seq

(A) A similarity network illustrates the diversity of single-cell expression landscapes. Three major groups of cells color-coded in green, red and purple are identified. Expression landscapes are visualized as mean SOM expression portraits of each group. Group-specific modules (A–D) of co-regulated genes are identified. (B) Pairwise correlation heatmap quantifies the mutual similarity of single-cell expression landscapes as Pearson's correlation coefficient (r). An anti-correlation for group 1 cells compared to cells of group 2 and 3 is discovered. (C) Mean SOM expression portrait of all cells illustrates the distribution of main expression spots across the SOM landscape. Spot-specific genes which could be assigned to different functional signatures (Proliferation, Oxphos, Pigmentation and Stroma) are identified. Top 5 genes of each spot defining the assignment to a functional group are shown for each spot. A refined analysis discovered subgroups (a–d) which are associated with different functions, e.g., different phases of cell cycle and ABC-membrane transporters. (D) Number of over-expression spots per cell increases from group 1 to group 3 revealed by the spot number distribution. Quality of cluster assignment is demonstrated by a silhouette plot. Genbank database accession number: GSE81383.

Figure 2. Patient-derived BRAF/NRAS wild type melanoma culture resembles partially primary and metastatic melanoma gene expression signatures

(A) Heatmap showing correlation of independent studies [–20] and our gene expression modules detected by single-cell RNAseq. Four of the tested gene expression data sets (Highgrade Up, Relapse Up, Metastasis Up and Pigmentation Type) show upregulation in cells of either group 1 or 2. Each column represents a single cell. (B) GSZ-profiles of identified gene sets illustrate the proportion of their similarity with single cells of our study. Each bar represents a single cell. (C) Selected gene sets are mapped on the SOM portrait to identify additional substructures. A distinct expression of genes between the Highgrade Up and the Relapse Up type can be recognized within spot A (cell cycle signature) indicated by red arrows. (D) Comparison of single-cell expression data from the patient-derived melanoma culture with melanoma expression data of an independent study by mapping top 100 genes of each expression spot (A–D) against melanoma data by Raskin and co-workers [21]. Cell cycle genes (spot A) are upregulated in metastatic melanoma and in high-grade melanomas.

Figure 3. Kinome analysis of the BRAF/NRAS wild type culture reveals potential targets for clinical treatments

(A) Gene expression levels of 4 cyclin-dependent kinases (CDK1, CDK2, CDK4 and CDK6) are analyzed for each of the 3 identified main groups. Data are given as standard box plots. CDK2 and CDK4 are consistently highly expressed in all groups. Positions of each CDK gene are shown on the mean SOM expression landscape indicated by arrows. (B) Growth inhibition curves are shown for treatments of a patient-derived melanoma short-term culture using inhibitors directed against CDK2 (purvanolol), CDK4 palbociclib, activated BRAF (vemurafenib), ERK1/2 (SCH772984), and MEK1/2 (cobimetinib).

Figure 4. Expression of marker genes of melanoma stem or initiating cells (MMIC)

(A) Heatmap showing scaled expression of MMIC marker genes [23] (e.g., ABCB5, KDM5B and ABCC4), various ABC-transporters and ALDH genes. Additional ABC-transporters and ALDH genes shown are all co-expressed in spot E (Supplementary Figure S5). Rows represent single cells. (B) Expression profile of spot E comprising 16 ABC-transporter and 8 ALDH genes. Selected expression portraits of cells expressing these genes are shown. (C) Expression profile of spot F where KDM5B (JARID1B) is heterogeneously expressed across single cells analyzed. Selected expression portraits of cells expressing KDM5B (JARID1B) are shown.

Figure 5. Comparative analysis reveals intertumor heterogeneity

(A) A similarity network was calculated and reveals separate clusters of BRAF/NRAS wild type (wt/wt), BRAF mutant/NRAS wild type (mut/wt) and BRAF wild type/NRAS mutant (wt/mut) cells, respectively. A replicate analysis of wt/wt cells was performed and fits into the respective wt/wt cluster. Cells in a highly proliferative state were highlighted in each cluster. (B) Pairwise correlation heatmap quantifies the mutual similarity of single-cell expression landscapes as Pearson's correlation coefficient (r). (C) Comparative SOM analysis between the different cultures generated five main expression signatures. Three expression signatures characterize each individual cell culture which represents intertumor heterogeneity. (D) Integrated proliferation signature which is shared between all cell cultures analyzed and group-specific expression signature for group 3 cells of the wild type short-term culture. Genbank database accession number: GSE81383.