Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq

Abstract

To explore the distinct genotypic and phenotypic states of melanoma tumors, we applied single-cell RNA sequencing (RNA-seq) to 4645 single cells isolated from 19 patients, profiling malignant, immune, stromal, and endothelial cells. Malignant cells within the same tumor displayed transcriptional heterogeneity associated with the cell cycle, spatial context, and a drug-resistance program. In particular, all tumors harbored malignant cells from two distinct transcriptional cell states, such that tumors characterized by high levels of the MITF transcription factor also contained cells with low MITF and elevated levels of the AXL kinase. Single-cell analyses suggested distinct tumor microenvironmental patterns, including cell-to-cell interactions. Analysis of tumor-infiltrating T cells revealed exhaustion programs, their connection to T cell activation and clonal expansion, and their variability across patients. Overall, we begin to unravel the cellular ecosystem of tumors and how single-cell genomics offers insights with implications for both targeted and immune therapies.

Figure 1. Dissection of melanoma with single-cell RNA-seq

(A) Overview of workflow. (B) Chromosomal landscape of inferred large-scale copy number variations (CNVs) distinguishes malignant from non-malignant cells. The Mel80 tumor is shown with individual cells (y-axis) and chromosomal regions (x-axis). Amplifications (red) or deletions (blue) were inferred by averaging expression over 100-gene stretches on the respective chromosomes. Inferred CNVs are concordant with calls from whole-exome sequencing (WES, bottom). (C,D) Single cell expression profiles distinguish malignant and non-malignant cell types. Shown are t-SNE plots of malignant (C, shown are the six tumors each with >50 malignant cells) and non-malignant (D) cells (as called from inferred CNVs as in B) from 11 tumors with >100 cells per tumor (color code). Clusters of non-malignant cells (called by DBScan, (17)) are marked by dashed ellipses and were annotated as T cells, B cells, macrophages, CAFs and endothelial cells, from preferentially expressed genes (Fig. S2, Table S2–3).

Figure 2. Single-cell RNA-seq distinguishes cell cycle and other states among malignant cells

(A) Estimation of the cell cycle state of individual malignant cells (circles) on the basis of relative expression of G1/S (x-axis) and G2/M (y-axis) gene-sets in a low-cycling (Mel79, top) and a high-cycling (Mel78, bottom) tumor. Cells are colored by their inferred cell cycle states, with cycling cells (red), intermediate (light red) and non-cycling cells (black); cells with high expression of KDM5B (Z-score>2) are marked in cyan filling. (B) IHC staining (40x magnification) for Ki67+ cells shows concordance with the signature-based frequency of cycling cells for Mel79 and Mel78 (as for other tumors; Fig S4C). (C) KDM5B/Ki67 staining (40x magnification) in corresponding tissue showing small clusters of KDM5B-high expressing cells negative for Ki67 (Fig. S4). (D) An expression program specific to Region 1 of Mel79, identified on the basis of multifocal sampling. The relative expression of genes (rows) is shown for cells (columns) ordered by the average expression of the entire gene-set. The region-of-origin of each cell is indicated in the top panel (Fig. S5).

Figure 3. MITF- and AXL-associated expression programs vary between and within tumors, and following treatment

(A) Average expression signatures for the AXL program (y-axis) or the MITF program (x-axis) stratify tumors into ‘MITF-high’ (black) or ‘AXL-high’ (red). (B) Single-cell profiles show a negative correlation between the AXL program (y-axis) and MITF program (x-axis) across individual malignant cells within the same tumor; cells are colored by the relative expression of the MITF (black) and AXL (red) programs. Cells in both states are found in all examined tumors, including three tumors (Mel79, Mel80 and Mel81) without prior systemic treatment, indicating that dormant resistant (AXL-high) cells may be present in treatment naïve patients. (C) Mel81 and Mel80 immunofluorescence staining of MITF (green nuclei) and AXL (red), validating the mutual exclusivity among individual cells within the same tumor (Fig. S8). (D) Relative expression (centered) of the AXL-program (top) and MITF-program (bottom) genes in six matched pre-treatment (white boxes) and post-relapse (gray boxes) samples from patients who progressed through RAF/MEK inhibition therapy; numbers at the top indicate patient index. Samples are sorted by the average relative expression of the AXL vs. MITF gene-sets. In all cases, the relapsed samples had increased ratio of AXL/MITF expression compared to their pre-treatment counterpart. This consistent shift of all six patients is statistically significant (P<0.05, binomial test), as are the individual increases in AXL/MITF for four of the six sample pairs (P<0.05, t-test; black and gray arrows denote increases that are individually significant or non-significant, respectively). (E) Quantitative, multiplexed single-cell immunofluorescence for AXL expression (y-axis top) and MAP-kinase pathway inhibition (pERK levels, y-axis) in the example cell lines WM88 and MELHO treated with increasing concentrations (x-axis) of either RAF inhibitor alone (dark gray bars) or a combination of RAF/MEK-inhibitors (light gray bars). We observed an increasing fraction of AXL-high cells (top panels) as well as a dose-dependent decrease of p-ERK (bottom panels) (Figs. S11–12 show additional cell lines).

Figure 4. Deconvolution of bulk melanoma profiles reveals cell-cell interactions

(A) Bulk tumors segregate to distinct clusters on the basis of their inferred cell type composition. Top panel: heat map showing the relative expression of gene sets defined from single-cell RNA-seq as specific to each of five cell types from the tumor microenvironment (y-axis) across 471 melanoma TCGA bulk-RNA signatures (x-axis). Each column is one tumor and tumors are partitioned into 10 distinct patterns identified by K-means clustering (vertical lines and cluster numbers at the top). Lower panels show, from top to bottom, tumor purity estimated by ABSOLUTE (DNA) and RNA-seq analysis (RNA), specimen location (from TCGA), and AXL/MITF scores. Tumors with high abundance of CAFs are correlated with an increased ratio of AXL-program/MITF-program expression (bottom). (B) Inferred cell-to-cell interactions between CAFs and T cells. Scatter plot compares for each gene (circle) the correlation of its expression with inferred T cell abundance across bulk tumors (y-axis, from TCGA transcriptomes) to how specific its expression is to CAFs vs. T cells (x-axis, based on single-cell transcriptomes). Genes that are highly specific to CAFs in a single cell analysis of tumors (red), but also associated with high T cell abundance in bulk tumors (black border) are candidates for CAF/T cell interactions. (C) 90 samples with 80 tumor specimens (black dots) showing a correlation (R=0.86) between C3/CD8 signal by quantitative immunofluorescence and 10 normal control specimens (grey dots) (Fig. S18A–F shows normalization and additional specimens). (D) Correlation coefficient (y-axis) between the average expression of CAF-derived complement factors shown in (B) and that of T cell markers (CD3/D/E/G, CD8A/B) across 26 TCGA cancer types with >100 samples (x-axis, left panel) and across 36 GTEx tissue types with >100 samples (x axis, right panel). Bars are colored on the basis of correlation ranges as indicated at the bottom.

Figure 5. Activation-dependent and independent variation in T-cell exhaustion markers

(A) Single T cell stratification into CD4+ and CD8+ cells (upper panel), CD25+FOXP3+ and other CD4 cells (middle panel) and their inferred activation state (lower panel, from average expression of the cytotoxic and naïve gene-sets in (B)). (B) Average expression of markers of cytotoxicity, exhaustion and naïve cell states (rows) in (left to right) Tregs, CD4+ T cells, and CD8+ T cells; CD4+ and CD8+ T cells are each further divided into five bins by their cytotoxic score (ratio of cytotoxic to naïve marker expression levels), showing an activation-dependent co-expression of exhaustion markers. Bottom: proportion of cycling cells (calculated as in Fig. 2B). Asterisks denote significant enrichment or depletion of cycling cells in a specific subset compared to the corresponding set of CD4+ or CD8+ T cells (P<0.05, hypergeometric test). (C) Immunofluorescence of PD-1 (upper panel, green), TIM-3 (middle panel, red) and their overlay (lower panel) validates their co-expression. (D) Activation-independent variation in exhaustion states within highly cytotoxic T cells. Scatter plot shows the cytotoxic score (x-axis) and exhaustion score (y-axis, average expression of the Mel75 exhaustion program as in Fig. S21) of each CD8+ T cell from Mel75. In addition to the overall correlation between cytotoxicity and exhaustion, the cytotoxic cells can be sub-divided into highly exhausted (red) and lowly exhausted cells (green) based on comparison to a LOWESS regression (black line). (E–F) Relative expression (log₂ fold-change) in high vs. low exhaustion cytotoxic CD8+ T cells from five tumors (x-axis), including 28 genes that were significantly upregulated (P<0.05, permutation test) in high-exhaustion cells across most tumors (E) and 272 genes that were variably associated with high-exhaustion cells across tumors (F). Three independently derived exhaustion gene-sets were used to define high and low exhaustion cells (Mel75, (17, 46, 48)), and the corresponding results are represented as distinct columns for each tumor. (G) Expanded TCR clones. Cells were assigned to clusters of TCR segment usage (black bars; Fig. S23), and cluster size (x-axis) was evaluated for significance by control analysis in which TCR segments were shuffled across cells (grey bars). The percentage of Mel75 cells (y-axis) is shown for clusters of small size (1–4 cells) that likely represent non-expanded cells, medium size (5–6 cells) that may reflect expanded clones (FDR=0.12), and large size that most likely reflect expanded clones (FDR=0.005). (H) Expanded clones are depleted of non-exhausted cells and enriched for exhausted cells. Mel75 cells were divided by exhaustion score into low exhaustion (green, bottom 25% of cells) and medium-to-high exhaustion (red, top 75%). Shown is the relative frequency of these exhaustion subsets (y-axis) in each TCR-cluster group (x-axis, as defined in G), defined as log2-ratio of the frequency in that group compared to the frequency across all Mel75 cells. All values were significant (P<10⁻⁵, binomial test).