. 2016 Apr 8;352(6282):189-96.

doi: 10.1126/science.aad0501.

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

Itay Tirosh¹, Benjamin Izar², Sanjay M Prakadan³, Marc H Wadsworth 2nd³, Daniel Treacy¹, John J Trombetta¹, Asaf Rotem⁴, Christopher Rodman¹, Christine Lian⁵, George Murphy⁵, Mohammad Fallahi-Sichani⁶, Ken Dutton-Regester⁷, Jia-Ren Lin⁸, Ofir Cohen¹, Parin Shah⁹, Diana Lu¹, Alex S Genshaft³, Travis K Hughes¹⁰, Carly G K Ziegler¹⁰, Samuel W Kazer³, Aleth Gaillard³, Kellie E Kolb³, Alexandra-Chloé Villani¹, Cory M Johannessen¹, Aleksandr Y Andreev¹, Eliezer M Van Allen⁴, Monica Bertagnolli¹¹, Peter K Sorger¹², Ryan J Sullivan¹³, Keith T Flaherty¹³, Dennie T Frederick¹³, Judit Jané-Valbuena¹, Charles H Yoon¹¹, Orit Rozenblatt-Rosen¹, Alex K Shalek¹⁴, Aviv Regev¹⁵, Levi A Garraway¹⁶

Affiliations

¹ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA.
² Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. Department of Medical Oncology, Dana-Farber Cancer Institute, Harvard Medical School, Boston, MA 02215, USA. Center for Cancer Precision Medicine, Dana-Farber Cancer Institute, Boston, MA 02215, USA. bizar@partners.org aregev@broadinstitute.org levi_garraway@dfci.harvard.edu.
³ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. Institute for Medical Engineering and Science, Massachusetts Institute of Technology (MIT), Cambridge, MA 02139, USA. Department of Chemistry, MIT, Cambridge, MA 02142, USA. Ragon Institute of Massachusetts General Hospital, MIT and Harvard University, Cambridge, MA 02139, USA.
⁴ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. Department of Medical Oncology, Dana-Farber Cancer Institute, Harvard Medical School, Boston, MA 02215, USA. Center for Cancer Precision Medicine, Dana-Farber Cancer Institute, Boston, MA 02215, USA.
⁵ Department of Pathology, Brigham and Women's Hospital, Harvard Medical School, Boston, MA 02115, USA.
⁶ Program in Therapeutic Sciences, Department of Systems Biology, Harvard Medical School, Boston, MA 02115, USA.
⁷ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. Department of Medical Oncology, Dana-Farber Cancer Institute, Harvard Medical School, Boston, MA 02215, USA. Department of Genetics and Computational Biology, QIMR Berghofer Medical Research Institute, Brisbane, Queensland, Australia.
⁸ HMS LINCS Center and Laboratory of Systems Pharmacology, Harvard Medical School, Boston, MA 02115, USA.
⁹ Department of Medical Oncology, Dana-Farber Cancer Institute, Harvard Medical School, Boston, MA 02215, USA.
¹⁰ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. Institute for Medical Engineering and Science, Massachusetts Institute of Technology (MIT), Cambridge, MA 02139, USA. Ragon Institute of Massachusetts General Hospital, MIT and Harvard University, Cambridge, MA 02139, USA. Division of Health Sciences and Technology, Harvard Medical School, Boston, MA 02115, USA.
¹¹ Department of Surgical Oncology, Dana-Farber Cancer Institute, Harvard Medical School, Boston, MA 02215, USA. Department of Surgical Oncology, Brigham and Women's Hospital, Harvard Medical School, Boston, MA 02115, USA.
¹² Program in Therapeutic Sciences, Department of Systems Biology, Harvard Medical School, Boston, MA 02115, USA. HMS LINCS Center and Laboratory of Systems Pharmacology, Harvard Medical School, Boston, MA 02115, USA. Ludwig Center at Harvard, Boston, MA 02215, USA.
¹³ Division of Medical Oncology, Massachusetts General Hospital Cancer Center, Boston, MA 02114, USA.
¹⁴ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. Institute for Medical Engineering and Science, Massachusetts Institute of Technology (MIT), Cambridge, MA 02139, USA. Department of Chemistry, MIT, Cambridge, MA 02142, USA. Ragon Institute of Massachusetts General Hospital, MIT and Harvard University, Cambridge, MA 02139, USA. Division of Health Sciences and Technology, Harvard Medical School, Boston, MA 02115, USA. Department of Immunology, Massachusetts General Hospital, Boston, MA 02114, USA.
¹⁵ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. Department of Biology and Koch Institute, MIT, Boston, MA 02142, USA. Howard Hughes Medical Institute, Chevy Chase, MD 20815, USA. bizar@partners.org aregev@broadinstitute.org levi_garraway@dfci.harvard.edu.
¹⁶ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. bizar@partners.org aregev@broadinstitute.org levi_garraway@dfci.harvard.edu.

PMID: 27124452
PMCID: PMC4944528
DOI: 10.1126/science.aad0501

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

Itay Tirosh et al. Science. 2016.

. 2016 Apr 8;352(6282):189-96.

doi: 10.1126/science.aad0501.

Authors

Affiliations

¹ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA.
² Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. Department of Medical Oncology, Dana-Farber Cancer Institute, Harvard Medical School, Boston, MA 02215, USA. Center for Cancer Precision Medicine, Dana-Farber Cancer Institute, Boston, MA 02215, USA. bizar@partners.org aregev@broadinstitute.org levi_garraway@dfci.harvard.edu.
³ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. Institute for Medical Engineering and Science, Massachusetts Institute of Technology (MIT), Cambridge, MA 02139, USA. Department of Chemistry, MIT, Cambridge, MA 02142, USA. Ragon Institute of Massachusetts General Hospital, MIT and Harvard University, Cambridge, MA 02139, USA.
⁴ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. Department of Medical Oncology, Dana-Farber Cancer Institute, Harvard Medical School, Boston, MA 02215, USA. Center for Cancer Precision Medicine, Dana-Farber Cancer Institute, Boston, MA 02215, USA.
⁵ Department of Pathology, Brigham and Women's Hospital, Harvard Medical School, Boston, MA 02115, USA.
⁶ Program in Therapeutic Sciences, Department of Systems Biology, Harvard Medical School, Boston, MA 02115, USA.
⁷ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. Department of Medical Oncology, Dana-Farber Cancer Institute, Harvard Medical School, Boston, MA 02215, USA. Department of Genetics and Computational Biology, QIMR Berghofer Medical Research Institute, Brisbane, Queensland, Australia.
⁸ HMS LINCS Center and Laboratory of Systems Pharmacology, Harvard Medical School, Boston, MA 02115, USA.
⁹ Department of Medical Oncology, Dana-Farber Cancer Institute, Harvard Medical School, Boston, MA 02215, USA.
¹⁰ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. Institute for Medical Engineering and Science, Massachusetts Institute of Technology (MIT), Cambridge, MA 02139, USA. Ragon Institute of Massachusetts General Hospital, MIT and Harvard University, Cambridge, MA 02139, USA. Division of Health Sciences and Technology, Harvard Medical School, Boston, MA 02115, USA.
¹¹ Department of Surgical Oncology, Dana-Farber Cancer Institute, Harvard Medical School, Boston, MA 02215, USA. Department of Surgical Oncology, Brigham and Women's Hospital, Harvard Medical School, Boston, MA 02115, USA.
¹² Program in Therapeutic Sciences, Department of Systems Biology, Harvard Medical School, Boston, MA 02115, USA. HMS LINCS Center and Laboratory of Systems Pharmacology, Harvard Medical School, Boston, MA 02115, USA. Ludwig Center at Harvard, Boston, MA 02215, USA.
¹³ Division of Medical Oncology, Massachusetts General Hospital Cancer Center, Boston, MA 02114, USA.
¹⁴ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. Institute for Medical Engineering and Science, Massachusetts Institute of Technology (MIT), Cambridge, MA 02139, USA. Department of Chemistry, MIT, Cambridge, MA 02142, USA. Ragon Institute of Massachusetts General Hospital, MIT and Harvard University, Cambridge, MA 02139, USA. Division of Health Sciences and Technology, Harvard Medical School, Boston, MA 02115, USA. Department of Immunology, Massachusetts General Hospital, Boston, MA 02114, USA.
¹⁵ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. Department of Biology and Koch Institute, MIT, Boston, MA 02142, USA. Howard Hughes Medical Institute, Chevy Chase, MD 20815, USA. bizar@partners.org aregev@broadinstitute.org levi_garraway@dfci.harvard.edu.
¹⁶ Broad Institute of MIT and Harvard, Cambridge, MA 02142, USA. bizar@partners.org aregev@broadinstitute.org levi_garraway@dfci.harvard.edu.

PMID: 27124452
PMCID: PMC4944528
DOI: 10.1126/science.aad0501

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.

PubMed Disclaimer

Figures

**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).

See this image and copyright information in PMC

Cited by

Clustering single-cell RNA sequencing data via iterative smoothing and self-supervised discriminative embedding.
Xie J, Ruan S, Tu M, Yuan Z, Hu J, Li H, Li S. Xie J, et al. Oncogene. 2024 Jul;43(29):2279-2292. doi: 10.1038/s41388-024-03074-5. Epub 2024 Jun 4. Oncogene. 2024. PMID: 38834657
Egr2 drives the differentiation of Ly6C^hi monocytes into fibrosis-promoting macrophages in metabolic dysfunction-associated steatohepatitis in mice.
Iwata A, Maruyama J, Natsuki S, Nishiyama A, Tamura T, Tanaka M, Shichino S, Seki T, Komai T, Okamura T, Fujio K, Tanaka M, Asano K. Iwata A, et al. Commun Biol. 2024 Jun 3;7(1):681. doi: 10.1038/s42003-024-06357-5. Commun Biol. 2024. PMID: 38831027 Free PMC article.
Single cell RNA-seq: a novel tool to unravel virus-host interplay.
Jogi HR, Smaraki N, Nayak SS, Rajawat D, Kamothi DJ, Panigrahi M. Jogi HR, et al. Virusdisease. 2024 Mar;35(1):41-54. doi: 10.1007/s13337-024-00859-w. Epub 2024 Mar 9. Virusdisease. 2024. PMID: 38817399 Review.
N-terminal tagging of RNA Polymerase II shapes transcriptomes more than C-terminal alterations.
Callan-Sidat A, Zewdu E, Cavallaro M, Liu J, Hebenstreit D. Callan-Sidat A, et al. iScience. 2024 May 7;27(6):109914. doi: 10.1016/j.isci.2024.109914. eCollection 2024 Jun 21. iScience. 2024. PMID: 38799575 Free PMC article.
Temporal Single Cell Analysis of Leukemia Microenvironment Identifies Taurine-Taurine Transporter Axis as a Key Regulator of Myeloid Leukemia.
Rodems BJ, Sharma S, Baker CD, Kaszuba CM, Ito T, Liesveld JL, Calvi LM, Becker MW, Jordan CT, Ashton JM, Bajaj J. Rodems BJ, et al. bioRxiv [Preprint]. 2024 May 14:2024.05.11.593633. doi: 10.1101/2024.05.11.593633. bioRxiv. 2024. PMID: 38798540 Free PMC article. Preprint.

See all "Cited by" articles

References

1. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 2011;144:646–674. - PubMed
1. Meacham CE, Morrison SJ. Tumour heterogeneity and cancer cell plasticity. Nature. 2013;501:328–337. - PMC - PubMed
1. Hodi FS, et al. Improved Survival with Ipilimumab in Patients with Metastatic Melanoma. N Engl J Med. 2010;363:711–723. - PMC - PubMed
1. Brahmer JR, et al. Phase I study of single-agent anti-programmed death-1 (MDX-1106) in refractory solid tumors: safety, clinical activity, pharmacodynamics, and immunologic correlates. J Clin Oncol Off J Am Soc Clin Oncol. 2010;28:3167–3175. - PMC - PubMed
1. Brahmer JR, et al. Safety and Activity of Anti–PD-L1 Antibody in Patients with Advanced Cancer. N Engl J Med. 2012;366:2455–2465. - PMC - PubMed

Publication types

Actions
Actions
Actions

MeSH terms

Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions
Actions

Substances

Actions
Actions

Grants and funding

LinkOut - more resources

Full Text Sources
Other Literature Sources
Medical
- MedlinePlus Health Information
Research Materials
- NCI CPTC Antibody Characterization Program
Miscellaneous
- NCI CPTAC Assay Portal

Save citation to file

Email citation

Add to Collections

Add to My Bibliography

Your saved search

Create a file for external citation management software

Your RSS Feed

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

Affiliations

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

Authors

Affiliations

Abstract

Figures

Similar articles

Cited by

References

Publication types

MeSH terms

Substances

Grants and funding

LinkOut - more resources

Full Text Sources

Other Literature Sources

Medical

Research Materials

Miscellaneous

Abstract

Figures

Similar articles

Cited by

References

Publication types

MeSH terms

Substances

Related information

Grants and funding

LinkOut - more resources

Full Text Sources

Other Literature Sources

Medical

Research Materials

Miscellaneous