Research Article

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

See allHide authors and affiliations

Science  08 Apr 2016:
Vol. 352, Issue 6282, pp. 189-196
DOI: 10.1126/science.aad0501

Single-cell expression profiles of melanoma

Tumors harbor multiple cell types that are thought to play a role in the development of resistance to drug treatments. Tirosh et al. used single-cell sequencing to investigate the distribution of these differing genetic profiles within melanomas. Many cells harbored heterogeneous genetic programs that reflected two different states of genetic expression, one of which was linked to resistance development. Following drug treatment, the resistance-linked expression state was found at a much higher level. Furthermore, the environment of the melanoma cells affected their gene expression programs.

Science, this issue p. 189


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.

Tumors are complex ecosystems defined by spatiotemporal interactions between heterogeneous cell types, including malignant, immune, and stromal cells (1). Each tumor’s cellular composition, as well as the interplay between these components, may exert critical roles in cancer development (2). However, the specific components, their salient biological functions, and the means by which they collectively define tumor behavior remain incompletely characterized.

Tumor cellular diversity poses both challenges and opportunities for cancer therapy. This is exemplified by the varied clinical efficacy achieved in malignant melanoma with targeted therapies and immunotherapies. Immune checkpoint inhibitors can produce clinical responses in many patients with metastatic melanomas (37); however, the genomic and molecular determinants of response to these agents remain incompletely understood. Although tumor neoantigens and PD-L1 expression clearly correlate with this response (810), it is likely that other factors from subsets of malignant cells, the microenvironment, and tumor-infiltrating lymphocytes (TILs) also play essential roles (11).

Melanomas that harbor the BRAFV600E (V600E: Val600→Glu600) mutation are commonly treated with inhibitors of rapidly accelerated fibrosarcoma kinase (RAF) and mitogen-activated protein kinase (MEK), before or after immune checkpoint inhibition. Although this regimen improves survival, virtually all tumors eventually develop resistance to these drugs (12, 13). Unfortunately, no targeted therapy currently exists for patients whose tumors lack BRAF mutations—including NRAS mutant tumors, those with inactivating NF1 mutations, or rarer events (such as RAF fusions). Collectively, these factors highlight the need for a deeper understanding of melanoma composition and its effect on the clinical course.

The next wave of therapeutic advances in cancer will probably be accelerated by technologies that assess the malignant, microenvironmental, and immunologic states most likely to inform treatment response and resistance. Ideally, we would be able to assess salient cellular heterogeneity by quantifying variation in oncogenic signaling pathways; drug-resistant tumor cell subsets; and the spectrum of immune, stromal, and other cell states that may inform immunotherapy response. Toward this end, single-cell genomic approaches enable detailed evaluation of genetic and transcriptional features present in hundreds to thousands of individual cells per tumor (1416). In principle, this approach may allow us to identify all major cellular components simultaneously, determine their individual genomic and molecular states (15), and ascertain which of these features may predict or explain clinical responses to anticancer agents. To explore this question, we used single-cell RNA sequencing (RNA-seq) to examine heterogeneities in malignant and nonmalignant cell types and states and to infer their possible drivers and interrelationships in the complex tumor cellular ecosystem.

Profiles of individual cells from patient-derived melanoma tumors

We measured single-cell RNA-seq profiles from 4645 malignant, immune, and stromal cells isolated from 19 freshly procured human melanoma tumors that span a range of clinical and therapeutic backgrounds (table S1). These included 10 metastases to lymphoid tissues (9 to lymph nodes and 1 to the spleen), 8 to distant sites (5 to subcutaneous or intramuscular tissue and 3 to the gastrointestinal tract), and one primary acral melanoma. Genotypic information was available for 17 of the 19 tumors, of which 4 had activating mutations in BRAF and 5 in NRAS oncogenes; eight patients had BRAF/NRAS wild-type melanomas (table S1).

To isolate viable single cells that are suitable for high-quality single-cell RNA-seq, we developed and implemented a rapid translational workflow (Fig. 1A) (15). We processed tumor tissues immediately after surgical procurement and generated single-cell suspensions within ~45 min, using an experimental protocol optimized to reduce artifactual transcriptional changes introduced by disaggregation, temperature, or time (17). Once in suspension, we recovered individual viable immune (CD45+) and nonimmune (CD45) cells (including malignant and stromal cells) by flow cytometry (fluorescence-activated cell sorting). Next, we prepared cDNA from the individual cells, followed by library construction and massively parallel sequencing. The average number of mapped reads per cell was ~150,000 (17), with a median library complexity of 4659 genes for malignant cells and 3438 genes for immune cells, comparable to previous studies of only malignant cells from fresh glioblastoma tumors (15).

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

(A) Overview of workflow. WES, whole-exome sequencing; RBC, red blood cell; FACS, fluorescence-activated cell sorting. (B) Chromosomal landscape of inferred large-scale CNVs allows us to distinguish malignant from nonmalignant 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 WES (bottom). (C and D) Single-cell expression profiles allow us to distinguish malignant and nonmalignant cell types. Shown are t-SNE plots of malignant [(C), shown are the six tumors, each with >50 malignant cells] and nonmalignant (D) cells [as called from inferred CNVs as in (B)] from 11 tumors with >100 cells per tumor (see color code below the panels). Clusters of nonmalignant cells [called by DBScan (17, 19)] are marked by dashed ellipses and were annotated as T cells, B cells, macrophages, CAFs, and endothelial (Endo.) cells from preferentially expressed genes (fig. S2 and tables S2 and S3). NK, natural killer cells.

Single-cell transcriptome profiles distinguish cell states in malignant and nonmalignant cells

We used a multistep approach to distinguish the different cell types within melanoma tumors on the basis of both genetic and transcriptional states (Fig. 1, B to D). First, we inferred large-scale copy number variations (CNVs) from expression profiles by averaging expression over stretches of 100 genes on their respective chromosomes (15) (Fig. 1B). For each tumor, this approach revealed a common pattern of aneuploidy, which we validated in two tumors by bulk whole-exome sequencing (WES) (Fig. 1B and fig. S1A). Cells in which aneuploidy was inferred were classified as malignant cells (Fig. 1B and fig. S1).

Second, we grouped the cells according to their expression profiles (Fig. 1, C and D, and fig. S2). To do this, we used nonlinear dimensionality reduction [t-distributed stochastic neighbor embedding (t-SNE)] (18), followed by density clustering (19). Generally, cells designated as malignant by CNV analysis formed a separate cluster for each tumor (Fig. 1C), suggesting a high degree of intertumor heterogeneity. In contrast, the nonmalignant cells clustered by cell type (Fig. 1D and fig. S2), independent of their tumor of origin and metastatic site (fig. S3). Clusters of nonmalignant cells were annotated as T cells, B cells, macrophages, endothelial cells, cancer-associated fibroblasts (CAFs), and natural killer cells on the basis of their preferentially or distinctively expressed marker genes (Fig. 1D, fig. S2, and tables S2 and S3).

Analysis of malignant cells reveals heterogeneity in cell cycle and spatial organization

We next used unbiased analyses of the individual malignant cells to identify biologically relevant melanoma cell states. After controlling for intertumor differences (17), we examined the six top components from a principal component analysis (PCA) (table S4). The first component correlated highly with the number of genes detected per cell and probably reflects technical aspects, whereas the other five significant principal components highlighted biological variability.

The second component (PC2) was associated with the expression of cell cycle genes (Gene Ontology project: “cell cycle” P < 10−16; hypergeometric test). To characterize cycling cells more precisely, we used gene signatures that have previously been shown to denote G1/S or G2/M phases in both synchronization (20) and single-cell (16) experiments in cell lines. Cell cycle phase–specific signatures were highly expressed in a subset of malignant cells, distinguishing cycling cells from noncycling cells (Fig. 2A and fig. S4A). These signatures revealed variability in the fraction of cycling cells across tumors (13.5% on average, ±13 SD) (fig. S4B), allowing us to designate both low-cycling (1 to 3%, e.g., Mel79) and high-cycling tumors (20 to 30%, e.g., Mel78), consistent with Ki67+ staining results (Fig. 2B and fig. S4C).

Fig. 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 tumor (Mel79, top) and a high-cycling tumor (Mel78, bottom). Cells are colored by their inferred cell cycle states: cycling cells, red; intermediate, pink; and noncycling cells, gray. Cells with high expression of KDM5B (z score > 2) are shown in green. N denotes number of cells. (B) Immunohistochemistry staining (40× 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 and Ki67 staining (40× magnification) in corresponding tissue showing small clusters of KDM5B-high expressing cells negative for Ki67 (fig. S4). DAPI, 4′,6-diamidino-2-phenylindole. (D) An expression program specific to region one 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).

A core set of cell cycle genes was induced (fig. S4D, red dots; and table S5) in both low-cycling and high-cycling tumors, with one notable exception: cyclin D3 (CCND3), which was induced in cycling cells only in high-cycling tumors (fig. S4D). In contrast, KDM5B (JARID1B) showed the strongest association with noncycling cells (Fig. 2A, green dots), mirroring Patel et al.’s findings in glioblastoma (15). KDM5B encodes a H3K4 histone demethylase associated with a subpopulation of slow-cycling and drug-resistant melanoma stemlike cells (21, 22) in mouse models. Immunofluorescence (IF) staining validated the presence and mutually exclusive expression of KDM5B and Ki67. KDM5B-expressing cells were grouped in small clusters, consistent with observations in mouse and in vitro models (21) (Fig. 2C and fig. S4E).

Two principal components (PC3 and PC6) primarily segregated different malignant cells from one treatment-naïve tumor (Mel79). In this tumor, we analyzed 468 malignant cells from four distinct regions after surgical resection (fig. S5A). We identified 229 genes with higher expression in the malignant cells of region one compared with those of other tumor regions [Fig. 2D, false discovery rate (FDR) < 0.05; and table S6]. A similar expression program was found in T cells from region one (fig. S6 and table S6), suggesting a spatial effect that influences multiple cell types. The genes with the highest preferential expression in region one are also generally coexpressed across melanoma tumors profiled in bulk in The Cancer Genome Atlas (TCGA) (23) (fig. S6). Many of these genes encode immediate early-activation transcription factors linked to inflammation, stress responses, and a melanoma oncogenic program (e.g., ATF3, FOS, FOSB, JUN, JUNB). Several of these transcription factors (e.g., FOS, JUN, NR4A1/2) are regulated by cyclic adenosine monophosphate (cAMP) and cAMP response element–binding protein signaling, which has been implicated as a mitogen-activated protein kinase (MAPK)–independent resistance module in BRAF-mutant melanomas treated with RAF and MEK inhibition (24). Other top genes differentially up-regulated in region one included those involved in survival (MCL1), stress responses (EGR1/2/3, NDRG, HSPA1B), and NF-κB signaling (NFKBIZ), which has also been associated with resistance to RAF and MEK inhibition (25). Immunohistochemistry analysis confirmed the elevated NF-κB and JunB levels in cells of region one compared with cells in the other regions of this tumor (fig. S5B).

Heterogeneity in the abundance of a dormant, drug-resistant melanoma subpopulation

Collectively, the above observations imply that pretreatment melanoma tumors may harbor subsets of malignant cells that are less likely to respond to targeted therapy. The transcriptional programs associated with principal components PC4 and PC5 were highly correlated with expression of the MITF gene (microphthalmia-associated transcription factor), which encodes the master melanocyte transcriptional regulator and a melanoma lineage-survival oncogene (26). Scoring genes by their correlation to MITF across single cells, we identified a “MITF-high” program consisting of MITF itself and several MITF target genes, including TYR, PMEL, and MLANA (table S7). A second transcriptional program, negatively correlated with the MITF program and with PC4 and PC5 (Pearson correlation P < 10−24), included AXL and NGFR (p75NTR), a marker of resistance to various targeted therapies (27, 28) and a putative melanoma cancer stem cell marker (29), respectively (table S8). Thus, these transcriptional programs resemble reported (25, 3032) MITF-high, as well as MITF-low and AXL-high (“AXL-high”), transcriptional profiles that can distinguish melanoma tumors, cell lines, and mouse models. Notably, the AXL-high program has been linked to intrinsic resistance to RAF and MEK inhibition (25, 30, 31).

Although at the bulk tumor level each melanoma could be classified as MITF-high or AXL-high (Fig. 3A), at the single-cell level every tumor contained malignant cells corresponding to both transcriptional states. Using single-cell RNA-seq to examine each cell’s expression of the MITF and AXL gene sets, we observed that MITF-high tumors, including treatment-naïve melanomas, harbored a subpopulation of AXL-high melanoma cells that was undetectable through bulk analysis, and vice versa (Fig. 3B). The malignant cells thus spanned the continuum between AXL-high and MITF-high states in every investigated tumor (Fig. 3B and fig. S7). We performed IF staining to further validate the mutually exclusive expression of the MITF-high and AXL-high programs in cells from the same bulk tumors (Fig. 3C and fig. S8).

Fig. 3 MITF- and AXL-associated expression programs vary between and within tumors, as well as after 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) categories. (B) Single-cell profiles show a negative correlation between the AXL program (y axis) and the 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 IF 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 genes (top) and MITF program genes (bottom) in six matched pretreatment (white boxes) and postrelapse (gray boxes) samples from patients who progressed through therapeutic RAF and MEK inhibition. Numbers at the top indicate patient index. Samples are sorted by the average relative expression of the AXL versus MITF gene sets. In all cases, the relapsed samples had an increased ratio of AXL-to-MITF expression compared with their pretreatment counterparts. This consistent shift of all six patients is statistically significant (P < 0.05, binomial test), as are the individual increases in AXL and MITF for four of the six sample pairs (P < 0.05, t test; black and gray arrows denote increases that are individually significant or nonsignificant, respectively). (E) Quantitative, multiplexed single-cell IF for AXL expression (top y axes) and MAP kinase pathway inhibition (p-ERK levels, bottom y axes) in the example cell lines WM88 and MELHO treated with increasing concentrations (x axis) of either a RAF inhibitor alone (dark gray bars) or a combination of RAF and 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 and S12 show results for additional cell lines).

Because malignant cells with AXL-high and MITF-high transcriptional states coexist in melanoma, we hypothesized that treatment with RAF and MEK inhibitors would increase the prevalence of AXL-high cells after the development of drug resistance. To test this, we analyzed RNA-seq data from a cohort (13) of six paired BRAFV600E melanoma biopsies taken before treatment and after resistance to single-agent RAF inhibition (vemurafenib; 1 patient) or combined RAF and MEK inhibition (dabrafenib and trametinib; 5 patients), respectively (tables S9 and S10). We ranked the 12 transcriptomes on the basis of the relative expression of all genes in the AXL-high program compared with those in the MITF-high program. In each pair, we observed a shift toward the AXL-high program in the drug-resistant sample [Fig. 3D; P < 0.05 for same effect in six of six paired samples, binomial test; P < 0.05 for four of six individual paired-sample comparisons shown by black arrows (17)]. RNA-seq data from an independent cohort (33) also showed that a subset of drug-resistant samples exhibited increased expression of the AXL program (fig. S9). Other genes previously implicated in resistance to RAF and MEK inhibition were also increased in a subset of the drug-resistant samples. PDGFRB (platelet-derived growth factor receptor β) (34) was up-regulated in a similar subset as the AXL program, whereas MET (33) was up-regulated in a mutually exclusive subset (fig. S9), suggesting that AXL and MET may reflect distinct drug-resistant states.

To further assess the connection between the AXL program and resistance to RAF and MEK inhibition, we studied single-cell AXL expression in 18 melanoma cell lines from the Cancer Cell Line Encyclopedia (35) (table S11). Flow cytometry analysis revealed a wide distribution of the proportion of AXL-positive cells, from <1 to 99% per cell line, which correlated with bulk mRNA levels and was inversely associated with sensitivity to small-molecule RAF inhibition (table S11).

We treated 10 cell lines (17) with increasing doses of a combination of RAF and MEK inhibitors (dabrafenib and trametinib) and found an increase in the proportion of AXL-positive cells in 6 cell lines initially composed of a small (<3%) pretreatment AXL-positive population (fig. S10A). In contrast, cell lines with an intrinsically high proportion of AXL expression showed modest or no changes (fig. S10B). We obtained similar results by multiplexed quantitative single-cell IF, which also demonstrated that the increased fraction of AXL-positive cells after inhibition of RAF and MEK is associated with rapid decreases in extracellular signal–regulated kinase (ERK) phosphorylation (reflecting MAP kinase signaling inhibition) (Fig. 3E and figs. S11 and S12). In summary, both melanoma tumors and cell lines demonstrate drug-resistant tumor cell subpopulations that precede treatment and become enriched after MAP kinase–targeted treatment.

Nonmalignant cells and their interactions within the melanoma microenvironment

Various nonmalignant cells make up the tumor microenvironment. The composition of the microenvironment has an important effect on tumorigenesis and also in the modulation of treatment responses (1). Tumor infiltration with T cells, for example, is predictive for the response to immune checkpoint inhibitors in various cancer types (36).

To resolve the composition of the melanoma microenvironment, we used our single-cell RNA-seq profiles to define distinct expression signatures of each of five distinct nonmalignant cell types: T cells, B cells, macrophages, endothelial cells, and CAFs. Because our signatures were derived from single-cell profiles, we could avoid confounders and ensure that each signature is determined by cell type–specific profiles (17). Next, we used these signatures to infer the relative abundance of those cell types in a larger compendium of tumors (17) (Fig. 4A and fig. S13). We found a strong correlation (correlation coefficient R ~ 0.8) between our estimated tumor purity and that predicted from DNA analysis (37) (Fig. 4A, first lane below the heat map).

Fig. 4 Deconvolution of bulk melanoma profiles reveals cell-to-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 represents one tumor, and tumors are partitioned into 10 distinct patterns identified by k-means clustering (vertical lines and cluster numbers at the top). Endo, endothelial cells; Macro., macrophages. (Lower panels, 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 a high abundance of CAFs are correlated with an increased ratio of AXL-to-MITF expression (bottom). LN, lymph node. (B) Inferred cell-to-cell interactions between CAFs and T cells. The 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 the specificity of its expression in CAFs (black) versus T cells (gray) (x axis, based on single-cell transcriptomes). Genes that are highly specific to CAFs in a single-cell analysis of tumors but are also associated with high T cell abundance in bulk tumors (red) are candidates for interaction between CAF cells and T cells. (C) Of the 90 samples, 80 tumor specimens (black dots) show a correlation (R = 0.86) between C3 and CD8 signals, as analyzed by quantitative IF. Ten normal control specimens (gray dots) are also shown (fig. S18, A to 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 (Genotype-Tissue Expression Project) tissue types with >100 samples (x axis, right panel). Bars are colored on the basis of correlation ranges, as indicated at the bottom.

We partitioned 471 tumors from TCGA into 10 distinct microenvironment clusters on the basis of their inferred cell type composition (Fig. 4A). Clusters were mostly independent of the site of metastasis (Fig. 4A, second lane), with some exceptions (e.g., clusters 8 and 9). Next we examined how these different microenvironments may relate to the phenotype of the malignant cells. In particular, CAF abundance is predictive of the AXL-MITF distinction, with CAF-rich tumors expressing the AXL-high signature (Fig. 4A, bottom lane). Interestingly, an AXL-high program was expressed by both melanoma cells and CAFs. However, we distinguished AXL-high genes that are preferentially expressed by melanoma cells (“melanoma-derived AXL program”) from those that are preferentially expressed by CAFs (“CAF-derived AXL program”). Both sets of genes were correlated with the inferred CAF abundance in tumors from TCGA (fig. S14) (38). Furthermore, the MITF-high program, which is specific to melanoma cells, was negatively correlated with inferred CAF abundance. Taken together, these results suggest that CAF abundance may be linked to preferential expression of the AXL-high over the MITF-high program by melanoma cells. Thus it is possible that specific tumor-CAF interactions may shape the melanoma cell transcriptome.

Interactions between cells play crucial roles in the tumor microenvironment (1). To assess how cell-to-cell interactions may influence tumor composition, we searched for genes expressed by cells of one type that may influence or reflect the proportion of cells of a different type in the tumor (fig. S15). For example, we searched for genes expressed primarily by CAFs (but not T cells) in single-cell data that correlated with T cell abundance (as inferred by T cell–specific genes) in bulk tumor tissue from the TCGA data set (23). We identified a set of CAF-expressed genes that correlated strongly with T cell infiltration (Fig. 4B, red circles). These included known chemotactic (CXCL12 and CCL19) and immune-modulating (PD-L2) genes, which are expressed by both CAFs and macrophages (fig. S16). A separate set of genes, exclusively expressed by CAFs, that correlated with T cell infiltration (fig. S16) included multiple complement factors [C1S, C1R, C3, C4A, CFB, and C1NH (SERPING1)]. Notably, these complement genes were specifically expressed by freshly isolated CAFs but not by cultured CAFs (fig. S17) or macrophages (fig. S16). These findings are intriguing, as studies have implicated complement activity in the recruitment and modulation of T cell–mediated antitumor immune responses [in addition to their role in augmenting innate immunity (39)].

We validated a high correlation (R > 0.8) between complement factor 3 (C3) levels (one of the CAF-expressed complement genes) and infiltration of CD8+ T cells. We performed dual IF staining and quantitative slide analysis of two tissue microarrays with a total of 308 core biopsies, including primary tumors, metastatic lesions, normal skin with adjacent tumor, and healthy skin controls (Fig. 4C and fig. S18) (17). To test the generalizability of the association between CAF-derived complement factors with T cell infiltration, we expanded our analysis to bulk RNA-seq data sets across all TCGA cancer types (Fig. 4D). Consistent with the results in melanoma, complement factors correlated with inferred T cell abundance in many cancer types and more highly than in normal tissues (e.g., R > 0.4 for 65% of cancer types but only for 14% of normal tissue types). Although correlation analysis cannot determine causality, this indicates a potential in vivo role for cell-to-cell interactions.

Diversity of tumor-infiltrating T lymphocytes and their functional states

The activity of TILs, particularly CD8+ T cells, is a major determinant of successful immune surveillance. Under normal circumstances, effector CD8+ T cells exposed to antigens and costimulatory factors may mediate lysis of malignant cells and control tumor growth. However, this function can be hampered by tumor-mediated T cell exhaustion, such that T cells fail to activate cytotoxic effector functions (40). Exhaustion is promoted through the stimulation of coinhibitory checkpoint molecules on the T cell surface (PD-1, TIM-3, CTLA-4, TIGIT, LAG3, and others) (41); blockade of checkpoint mechanisms has shown clinical benefit in subsets of melanoma and other malignancies (3, 10, 42, 43). Although checkpoint ligand expression (e.g., PD-L1) and neoantigen load clearly contribute (9, 44, 45), no biomarker has emerged that reliably predicts the clinical response to immune checkpoint blockade. We reasoned that single-cell analyses might yield features to elucidate response determinants and possibly identify new immunotherapy targets.

Thus, we analyzed the single-cell expression patterns of 2068 T cells from 15 melanomas. We identified T cells and their main subsets [CD4+, regulatory T cells (Tregs), and CD8+] on the basis of the expression levels of their respective defining surface markers (Fig. 5A, top, and table S12). Within both the CD4+ and CD8+ populations, a PCA distinguished cell subsets and heterogeneity of activation states on the basis of the expression of naïve and cytotoxic T cell genes (Fig. 5, A and B, and fig. S19).

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

(A) Single T cell stratification into CD4+ and CD8+ cells (top panel), CD25+FOXP3+ and other CD4 cells (middle panel), and their inferred activation state [lower panels, from average expression of the cytotoxic and naïve gene sets in (B)]. Th, T helper cells; Tregs; regulatory T cells. (B) (Top) Average expression of markers of cytotoxicity (Cyto.), exhaustion (Exhau.), and naïve cell states (rows) in (from left to right) Tregs, CD4+ T helper 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 activation-dependent coexpression 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, as compared with the corresponding set of CD4+ or CD8+ T cells (P < 0.05, hypergeometric test). (C) IF analysis of PD-1 (top, green), TIM-3 (middle, red), and their overlay (bottom) validates their coexpression. (D) Activation-independent variation in exhaustion states within highly cytotoxic T cells. The 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 subdivided into cells with high (red) and low exhaustion (green), based on comparison to a LOWESS (locally weighted scatter plot smoothing) regression (black line). (E and F) Relative expression (log2 fold-change) in high- versus low-exhaustion cytotoxic CD8+ T cells from five tumors (x axis), including 28 genes that were significantly up-regulated (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 (dark gray bars) (fig. S23), and cluster size (x axis) was evaluated for significance by control analysis in which TCR segments were shuffled across cells (light gray bars). The percentage of Mel75 cells (y axis) is shown for clusters of small size (one to four cells) that probably represent nonexpanded cells, medium size (five or six cells) that may reflect expanded clones (FDR = 0.12), and large size (more than six cells) that most likely reflect expanded clones (FDR = 0.005). (H) Expanded clones are depleted of nonexhausted cells and enriched for exhausted cells. Mel75 cells were divided according to exhaustion score into categories of 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 the log2 ratio of the frequency in that group compared with the frequency across all Mel75 cells. All values were significant (P < 10−5, binomial test).

Next we sought to determine the exhaustion status of each cell from the expression of key coinhibitory receptors (PD1, TIGIT, TIM3, LAG3, and CTLA4). In several cases, these coinhibitory receptors were coexpressed across individual cells; we validated this phenomenon for PD1 and TIM3 by IF staining (Fig. 5C). However, exhaustion gene expression was also highly correlated with the expression of both cytotoxicity markers and overall T cell activation states (Fig. 5B). This observation resembles an activation-dependent exhaustion expression program, such as those reported previously (46, 47). Accordingly, expression of coinhibitory receptors (alone or in combinations) may not be sufficient by itself to characterize the salient functional state of tumor-associated T lymphocytes in situ or to distinguish exhaustion from activation.

To define an activation-independent exhaustion program, we leveraged single-cell data from CD8+ T cells sequenced in a single tumor (Mel75, 314 cells). These data allowed cytotoxic and exhaustion programs to be deconvolved. Specifically, PCA of Mel75 T cell transcriptomes identified a robust expression module that consisted of all five coinhibitory receptors and other exhaustion-related genes, but not cytotoxicity genes (fig. S21 and table S13).

We used the Mel75 exhaustion program, along with previously published exhaustion programs (46, 48), to estimate the exhaustion state of each cell. An exhaustion state was defined as high or low expression of the exhaustion program relative to that of the cytotoxicity genes (Fig. 5D) (17). Accordingly, we defined exhaustion states in Mel75 and in four additional tumors with the highest number of CD8+ T cells (68 to 214 cells per tumor). We identified the top preferentially expressed genes in high-exhaustion cells compared with low-exhaustion cells (both defined relative to the expression of cytotoxicity genes). This allowed us to define a core exhaustion signature across cells from various tumors.

Our core exhaustion signature yielded 28 genes that were consistently up-regulated in high-exhaustion cells of most tumors, including coinhibitory (TIGIT) and costimulatory (TNFRSF9/4-1BB and CD27) receptors (Fig. 5E and table S14). In addition, most genes that were significantly up-regulated in high-exhaustion cells of at least one tumor had distinct associations with exhaustion across the different tumors (Fig. 5F, 272 of 300 genes with P < 0.001 by permutation test; fig. S22, A and B; and table S14). These tumor-specific signatures included variable expression of known exhaustion markers (table S14) and could be linked to immunotherapeutic response or might reflect the effects of previous treatments. For example, CTLA-4 was highly up-regulated in exhausted cells of Mel75 and weakly up-regulated in three other tumors but was completely decoupled from exhaustion in Mel58. Interestingly, Mel58 was derived from a patient with an initial response and subsequent development of resistance to CTLA-4 blockade with ipilimumab (Fig. 5F and fig. S22C). Another variable gene of interest was the transcription factor NFATC1, previously implicated in T cell exhaustion (49). NFATC1 and its target genes were preferentially associated with the activation-independent exhaustion phenotype in Mel75 (fig. S22, D and E), suggesting a potential role of NFATC1 in the underlying variability of exhaustion programs among patients.

Finally, we explored the relationship between T cell states and clonal expansion. T cells that recognize tumor antigens may proliferate to generate discernible clonal subpopulations defined by an identical T cell receptor (TCR) sequence (50). To identify potentially expanded T cell clones, we used RNA-seq reads that map to the TCR to classify single T cells by their isoforms of the V and J segments of the α and β TCR chains, and we searched for enriched combinations of TCR segments. Most observed combinations were found in few cells and were not enriched. However, approximately half of the CD8+ T cells in Mel75 had one of the seven enriched combinations identified (FDR = 0.005) and thus may represent expanded T cell clones (Fig. 5G and fig. S23). This putative T cell expansion was also linked to exhaustion (Fig. 5H), such that low-exhaustion T cells were depleted in expanded T cells (TCR clusters with more than six cells) and enriched in nonexpanded T cells (TCR clusters with one to four cells). In particular, the nonexhausted cytotoxic cells are almost all nonexpanded cells (Fig. 5H). Overall, this analysis suggests that single-cell RNA-seq may allow for the inference of functionally variable T cell populations that are not detectable with other profiling approaches (fig. S24). This knowledge may empower studies of tumor response and resistance to immune checkpoint inhibitors.


Our analysis has uncovered intra- and interindividual, spatial, functional, and genomic heterogeneity in melanoma cells and associated tumor components that shape the microenvironment, including immune cells, CAFs, and endothelial cells. We identified a cell state in a subpopulation of all melanomas studied that is linked to resistance to targeted therapies, and we used a variety of approaches to validate the presence of a dormant drug-resistant population in a number of melanoma cell lines.

By leveraging single-cell profiles from a few tumors to deconvolve a large collection of bulk profiles from TCGA, we discovered different microenvironments associated with distinct malignant cell profiles. We also detected a subset of genes expressed by one cell type (e.g., CAFs) that may influence the proportion of other cell types (e.g., T cells); this indicates the importance of intercellular communication for tumor phenotype. Putative interactions between stromal-derived factors and immune cell abundance in melanoma core biopsies suggest that future diagnostic and therapeutic strategies should account for tumor cell composition rather than bulk expression. Furthermore, our data suggest potential biomarkers for distinguishing exhausted and cytotoxic T cells that may aid in selecting patients for immune checkpoint blockade.

Although future work is necessary to clarify the interplay between these cell types and functional states in space and time, the ability to carry out a number of highly multiplexed single-cell observations within a tumor allows us to identify meaningful cell subpopulations and gene expression programs that may inform both the analysis of bulk transcriptional data and precision treatment strategies. Conceivably, single-cell genomic profiling may soon enable a deeper understanding of the complex interplay among cells within the tumor ecosystem and its evolution in response to treatment, thereby providing a versatile new tool for future translational applications.

Supplementary Materials

Materials and Methods

Figs. S1 to S24

Table S1 to S16

References (5167)

References and Notes

  1. Materials and methods are available as supplementary materials on Science Online.
Acknowledgments: We thank M. Singer, A. Anderson, V. K. Kuchroo, and A. Aguirre for fruitful discussions; S. Barthel and T. Schatton for providing the Tim-3 antibody; Q. Zhan for preparing IF staining; and J. Thurman for guidance on complement biology. L.A.G., B.I., S.M.P., A.Re., O.R.-R., A.K.S., I.T., M.H.W.II, The Broad Institute, Brigham and Women’s Hospital, Dana-Farber Cancer Institute, MIT, and the president and fellows of Harvard College have filed a patent application (BI-2015/077) that relates to tumor and microenvironment gene expression, compositions of matter, and methods of use thereof. L.A.G., B.I., A.Ro., and the Dana-Farber Cancer Institute have filed a patent application (DFCI 2105.001) that relates to cancer-patient–derived tumor dissociation for biological analysis. L.A.G was supported by National Cancer Institute (NCI) grants P01CA163222 and R35CA197737, the Dr. Miriam and Sheldon Adelson Medical Research Foundation, the Melanoma Research Alliance, and the Ludwig Center at Harvard Medical School. L.A.G. is a member of the scientific advisory board for Warp Drive; a consultant for Novartis, Bayer Oncology, and Foundation Medicine, and an equity holder in Foundation Medicine. A.Re. was supported by funds from the Howard Hughes Medicine Institute, the Klarman Cell Observatory, STARR Cancer Consortium, NCI grant 1U24CA180922, Koch Institute support (core) grant P30-CA14051 from NCI, and the Broad Institute. A.Re. is a scientific advisory board member for ThermoFisher Scientific and Syros Pharmaceuticals and a consultant for Driver Group. E.M.V.A. is a consultant for Roche Ventana, Takeda, and Third Rock Ventures. A.K.S. was supported by the Searle Scholars Program, the Beckman Young Investigator Program, and the NIH New Innovator Award (DP2 OD020839). B.I. was supported by the Wong Family Award for Translational Oncology of the Dana-Farber Cancer Institute. I.T. was supported by a Human Frontier Science Program long-term fellowship, a Rothschild fellowship, STARR Cancer Consortium, and an Integrative Cancer Biology Program grant (U54CA112962). I.T., A.K.S., and O.R.-R. were supported by the Broad Institute. O.R.-R. was supported by a grant from the Next Generation Fund at the Broad Institute of MIT and Harvard. S.W.K. was supported by a fellowship from the NSF Graduate Research Fellowships Program. M.F.-S. was supported by the NCI (grant K99CA194163), and P.K.S. was supported by the NIH (grant P50GM107618) and the Ludwig Center at Harvard. Processed single-cell and bulk RNA-seq data are available through the Gene Expression Omnibus (accession numbers GSE72056 and GSE77940). Raw RNA-seq and WES data will be available through dbGAP (the database of Genotypes and Phenotypes). The raw human data are available at the Data Use Oversight System (DUOS) ( with accession number DUOS-000002.

Correction (1 February 2019): The raw human data have been deposited in the Broad Institute's Data Use Oversight System (DUOS) repository. The accession number has been added to the Acknowledgments.

Stay Connected to Science

Navigate This Article