Single-Cell Characterization of the Immune Microenvironment of Melanoma Brain and Leptomeningeal Metastases

Melanoma metastasizes to many organs, with common sites being the brain, skin, lungs, liver, and bones (1, 2). CNS involvement is clinically evident in over 40%–60% of metastatic melanoma patients, with rates being as high as 75% at autopsy (3). In about 5% of cases, melanoma cells metastasize to the leptomeninges, the subarachnoid space, and cerebrospinal fluid (CSF). Patients with leptomeningeal melanoma metastases (LMM) have the worst prognosis and are characterized by rapid disease progression (mean survival 8–10 weeks) and a death from neurologic causes (4–6). Recent advances in the development of targeted therapies and immunotherapies have led to improved survival for most melanoma patients with advanced disease (7–10). Despite this, responses are often heterogeneous within individual patients, with some lesions responding well while others do not (11). Patients on BRAF-MEK inhibitor therapy and immunotherapy frequently show progression in the brain even when the extracranial disease is well controlled (12–15). There is also evidence that responses to targeted therapy and immunotherapy are shorter in duration in melanoma brain metastases (MBM) than those at extracranial sites (16, 17). In patients with LMM, outcomes after either targeted therapy or immunotherapy are particularly poor (6, 18).

The reasons underlying organ-specific patterns of tumor progression and therapeutic response are not well understood. Although it is possible that organ-specific metastasis may result from the acquisition of distinct oncogenic mutations, it is clear that each metastatic site differs in its immune microenvironment (19–21). This likely reflects the need to keep inflammatory responses localized (such as in the skin, lungs, and gut) and to respond to distinct organ-specific pathogenic insults (22). At this time, very little is known about the immune microenvironments of MBM or LMM. Whereas patients with MBM can respond favorably to both targeted therapy and immunotherapy, patients with LMM respond very poorly, suggestive of significant differences between the immune microenvironment of these two CNS compartments. Although the majority of patients with LMM progress rapidly, a small proportion of individuals do well on therapy and can live for many years after their diagnosis (23). The biology underlying these extreme responses has never been explored, but could give critical new insights that allow improved therapies to be developed. In order to better understand these differences, we have undertaken a single-cell RNA-sequencing (scRNA-seq) analysis of 43 clinical melanoma specimens from 26 individuals taken from three distinct metastatic sites: the skin, brain, and the leptomeninges/CSF. Many of these samples were serial and were obtained pre- and post-immune or targeted therapy. We developed a novel computational platform, ISCVA (Interactive Single-Cell Visual Analytics), described below, that allowed visualization and analysis of single-cell transcriptomic data. Our analyses revealed the tissue microenvironments (TME) of brain and leptomeningeal metastatic sites to be very distinct and show differential responses to systemic therapy. Correlation analyses identified a hitherto uncharacterized population of dendritic cells (DC3) that were associated with a more favorable immune environment and improved overall survival (OS) in melanoma patients.

This study was conducted in accordance with recognized ethical guidelines (e.g., Declaration of Helsinki, CIOMS, Belmont Report, U.S. Common Rule). Forty-one human melanoma specimens from 24 patients, and 2 CSF samples from 2 non-LMM patients were procured under protocols approved by the Institutional Review Board (MCC50102, MCC50103, MCC50172, MCC19332, MCC18597, MCC19441/ICAF Pro000024988). Upon collection, samples were immediately placed on ice and transferred for processing. For CSF specimens, cells were quantified and directly loaded for library preparation. Samples from brain and skin metastases were digested using the Miltenyi Tissue Dissociation Kit (Bergisch Gladbach) and enriched for live cells using FACS sorting prior to quantification and loading. Another RNA-seq data set from a Moffitt cohort of 135 melanoma patients was analyzed under protocol MCC#19147 (24) and further specimens used for validation studies were obtained and analyzed under protocols MCC21044 and MCC20779, all approved by the Institutional Review Board.

A single-cell suspension from each tissue was quantified and analyzed for viability using the Nexcelom Cellometer K2 and then loaded onto the 10X Genomics Chromium Single-Cell Controller for scRNA-seq library preparation (10X Genomics). Briefly, the single cells, reagents, and 10X Genomics gel beads were encapsulated into individual nanoliter-sized Gelbeads in Emulsion (GEMs) and then reverse transcription of poly-adenylated mRNA was performed inside each droplet. The cDNA was amplified, purified, and cDNA libraries were then prepared in bulk reactions using the 10X Chromium Single-Cell 3′ Library Prep Kit. Approximately 50,000 to 1,000,000 mean sequencing reads per cell were generated on the Illumina NextSeq 500 instrument using v2.5 flow cells. Demultiplexing, barcode processing, alignment, and gene counting were performed using the 10X Genomics CellRanger software.

To facilitate the rapid analysis of single-cell data sets, we developed a new computational tool consisting of two major components. The first component comprised of a collection of Bash and R scripts [utilizing many of the widely used algorithms in the single-cell community, including Seurat for general processing (25), SingleR for cell type recognition (26), and single-cell signature explorer for gene set signature scoring (27)] that processed the scRNA-seq data offline, and a second web-based component (built with state-of-the-art web front-end technologies, including react.js from Facebook, tensorflow.js from Google, and Plotly.js) that allows convenient real-time interactive exploration and ad hoc analysis. The heterogeneity analyses implemented in SinCHet (28) were also performed as part of the analytical modules. A node.js backend was also created to serve the on-demand queries of the web application, allowing for real-time interactive investigation of genes expressed in selected samples or subsets of cells. The webtool (along with all of the data analyzed herein) can be accessed at http://iscva.moffitt.org. The full data set is also available at GEO (GSE174401). Two publicly available data sets, one from melanoma clinical samples and another from peripheral blood of healthy donors, are included as example data sets (29, 30).

We used the R package Seurat to process the aggregated 31,075-cell transcript count matrix generated from the 10X Cell Ranger pipeline. A total of 2,367 cells with >20% reads in mitochondria were filtered from further analyses. A two-stage clustering was performed. The first stage was for broad cell type identification while the second stage was for cell subpopulation identification. Regularized negative binomial regression, along with a second linear regression against mitochondria read percentage as implemented in SCTransform, was used to adjust for cell-to-cell technical variations. We used the unsupervised clustering algorithm Louvain in the principal component analysis (PCA) space to identify cell groups among all cells that passed QC. In addition, the R package SingleR was used to generate cell type candidates for each cell using multiple reference panels. Each cell group is then assigned a broad cell category based on the combination of SingleR predictions—mostly using the BlueprintEncode panel and literature review. Ambient contamination of immunoglobulin transcripts was detected in four of the samples. To remove the effect of these contaminants, the entire analysis, except SingleR, was repeated with the expression levels of these genes zeroed out for all cells that were not categorized as pDC, B, or plasma cells. We applied a second stage of clustering and curation to respectively map out the substructure of the myeloid, lymphoid, and neuronal cell populations. The myeloid cell analysis included all cells that fell under the broad categories of pDC and macrophage/monocyte/DC, while the lymphoid analysis included all T/NK cells. The second-stage analyses generated unsupervised clustering of the subpopulations of cells and then annotated these clusters according to their distinguishing gene expression, using published markers on myeloid and lymphoid cell subpopulations as a guide. In the myeloid analysis, Infomap was used instead of Louvain as an alternative clustering method to separate the cells at a finer detail, with attention to known DC subpopulations. The template for generating an R markdown report with summary statistics of identified cell subpopulations, lists of differentially expressed genes, and heat map was also included. All backend scripts, R templates, and the javascripts are available for download: https://github.com/zhihua-chen/iscva/. For validation purposes, scRNA-seq data of 4,645 cells from 19 melanoma patients from a data set previously published by Tirosh and colleagues (30) were loaded into ISCVA. Briefly, the QC and two-stage clustering described above was performed to identify broad cell types in general, and detailed subpopulations within myeloid and lymphoid cells as described above.

Human FFPE tissue samples were stained using the PerkinElmer OPAL TM Automation IHC kit on the BOND RX autostainer (Leica Biosystems) following the automated OPAL IHC procedure (PerkinElmer) with DAPI stain (Akoya Biosciences) and antibodies against LAMP3 (Abcam), IDO1 (Origene), and CCL19 (Invitrogen). Slides were imaged using the Vectra3 Automated Quantitative Pathology Imaging System. Multilayer TIFF images were exported from InForm (PerkinElmer) and loaded into HALO (Indica Labs) for quantitative image analysis. The tissue was segmented into individual cells using the DAPI and a positivity threshold for each DC3 marker was determined based on published staining patterns and intensity for that specific antibody. LAMP3⁺ IDO1⁺ CCL19⁺ cells were identified as DC3s. Immune phenotyping of melanoma brain tissues was performed following the same protocol, using antibodies against CD4 (Cell Marque), CD8 (DAKO/Agilent), Foxp3 (Abcam), and CD68 (CST).

C57BL/6J mice were injected with 50,000 SM1 melanoma cells (passage #5) obtained from Dr. Eric Lau (Moffitt Cancer Center) via stereotactic injection into the brain and then treated with isotype control or anti–PD-1 antibody (200 μg/100 μL, IP; clone RMP1-14; Bio X Cell, NH) every 5 days. At endpoint, tumors were harvested and digested using the Tissue Dissociation Kit. Single-cell suspensions were stained for immunophenotyping using two panels of antibodies for T cells and myeloid cells (31) and analyzed by flow cytometry, as described previously (31). Statistical analysis of data was carried out using a one-sided Mann–Whitney test. Cells were kept for a maximum of 10 passages (total). Cell lines were routinely tested for Mycoplasma (every three months) and were authenticated (6-month intervals) by STR authentication (last date of stock testing: December 18, 2019). All the protocols were reviewed and approved by Institutional Animal Care and Use Committee at University of South Florida (approval R IS00004882).

Fresh human melanoma tumors were harvested from the leptomeninges at autopsy and from the brain during surgery. Tumors were digested using the Tissue Dissociation Kit. Single-cell suspension of tumor harvested from the leptomeninges was stained using antibodies against CD45, CD3, CD4, CD8, CD163, CD25, CD127, and Live/Dead stain (Supplementary Table S1). Single-cell suspensions of MBM tumors were stained using antibodies against CD45, CD19, CD3, CD4, HLA-DR, CD163, CD25, CD127, CD14, and Live/Dead stain (Supplementary Table S1). Samples were analyzed by flow cytometry.

Proportions of cells were calculated by the number of cells in each identified population divided by the total number of cells from a patient sample, if not specified otherwise. Kruskal–Wallis analyses were performed to compare the size of cell populations among the three sites. Correlation among each cell proportions of all samples was assessed using Spearman correlation. To evaluate potential associations between the proportion of each cell population and patient OS outcomes, Cox regression was performed with site of metastasis as a covariate in the model to account for the known survival differences among different metastatic sites. When multiple samples were available, mean cell population proportion was used in the analyses. The optimal cut-point for the proportion of cell population of interest in the survival analyses was chosen using the maximally selected rank statistic in the “coin” R package and visualized using Kaplan–Meier curves. Log-rank tests were performed to compare survival difference between two groups.

The method SingleCellSignalR (32) was used to investigate the cell–cell interactions mediated by ligand–receptor complexes between DC3 cells and CD4 or CD8 subpopulations using scRNA-seq data. Briefly, SingleCellSignalR infers the ligand (L) and receptor (R) interactions between two cell types using a regularized LRscore of the 3,251 ligand–receptor interactions from its compiled and curated LR databases. LRscore is a scaled product score of average ligand expression in one cell type and average receptor expression in other cell type. The heat map visualized the interactions with LRscore higher than 0.5 between DC3 and CD4 or CD8 subpopulations, further filtered by DC3-specific markers. The interaction score is a product of the number of all interactions with LRscore higher than 0.5 and the sum of these interactions' LR scores.

The level-3 RNA-seq expression data and clinical covariates from TCGA melanoma patient samples were used (N = 459 patients). LAMP3, IDO1, IDO2, and CCL19 were identified as robust markers of DC3 cells with high specificity both in our own data set and the Tirosh melanoma patient data set. To investigate the association between the DC3 population and patient outcome in the TCGA data set, we performed Cox regression model between the expression level (in log₂ scale) of each of the DC3-specific genes and patients' OS outcome. PCA was performed to generate a signature to summarize gene expression of these four markers. The results of survival analyses using the median as the cut-point for each gene and the first principal component (PC1) were visualized using Kaplan–Meier curves. Log-rank tests were performed to compare survival difference between two groups. OS is defined from the time of sample collection to death and censored at the last follow-up. The same analysis was also performed on an additional RNA-seq data set from a Moffitt melanoma patient cohort (N = 135; ref. 24).

The mouse experiment was carried out as described previously (31). SW1 (passage #10) melanoma cells were obtained from Dr. Eric Lau (Moffitt Cancer Center) in 2017. STR validation and Mycoplasma testing were performed as for the SM1 cells (see above). Briefly, 1.5 × 10⁶ SW1 melanoma cells were subcutaneously injected into the flanks of C3H/HeNCrl mice. Mice were administered anti–PD-1 antibody (200 μg/100 μL, IP) or IgG isotype control every 5 days. On day 10, the immunotherapy was stopped and the mice received either vehicle or ceritinib (25 mg/kg) + trametinib (1 mg/kg) in combination via oral gavage daily. At endpoint or by day 41, tumors were measured and collected for analysis by scRNA-seq for the presence of DC3 cells. Statistical analysis of tumor size differences was carried out using a parametric unpaired t test.

To interrogate the tumor microenvironment (TME) at a single-cell resolution, we performed droplet-based scRNA-seq on 43 patient specimens. These were derived from 8 punch biopsies of skin metastases from melanoma patients enrolled on a phase I clinical trial of vemurafenib–cobimetinib–XL888 (including 5 serial specimens from two patients, NCT02721459), 14 surgical specimens from melanoma patients undergoing resection of brain metastases (including 2 serial specimens from one patient), and 19 serial specimens of CSF from 6 patients with LMM (Fig. 1A). CSF from 2 non-melanoma patients without LMM was included as a control. Among the patients with LMM, the majority of samples (n = 15) were from patients with a poor outcome, but one set of four samples was from an extreme outlier who responded to systemic therapy and survived for over 38 months after diagnosis of LMM. All the skin and LMM samples harbored a BRAF V600E mutation. Among the brain metastasis samples, 8/14 harbored a BRAF V600E mutation, with one BRAF D594N mutant, one BRAF V600K mutant, one BRAF K601E/GNAS R201P mutant, one NRAS Q61K mutant, and two without any identifiable driver mutation (Supplementary Table S2).

To analyze the scRNA-seq data, we developed our own analytics tool, ISCVA, that used a two-stage architecture. The first stage processed the scRNA-seq data. The second stage allowed for the interactive selection and visualization of subpopulations of cells and genes in both UMAP and t-distributed stochastic neighbor embedding (t-SNE) projections (Fig. 1A), which revealed the melanoma samples to consist of a diverse landscape of malignant and host cell types (Fig. 1B). We used the unsupervised clustering algorithm Louvain in the PCA space to find cell groups among all cells that passed QC, which were assigned a broad cell category based on the majority vote using the SingleR predictions (BlueprintEncode panel) followed by manual curation from the literature. A total of 10 broad cell types (melanoma cells, macrophage/monocyte/DCs, T/NKs, pDCs, fibroblasts, B cells, plasma cells, neurons, endothelial and other cells) were identified at this level. Analysis of all 43 samples revealed the melanoma cells to be the most diverse population of cells and that individual patients' cells formed their own distinct clusters, whereas the stromal and immune cell subtypes from different patients tended to cluster together (Fig. 1B). Analysis of the cell type composition for each site also revealed evident differences in the TME composition (Fig. 1C). Cells from the LMM samples were more abundant in T/NK cells and had the highest proportion of myeloid/monocytic cells. Brain metastasis samples tended to have the highest proportion of melanoma cells per sample, lower T/NK cell numbers, and the presence of plasma cells, B cells, and neurons (as well as other neuronal cells). Samples from skin metastases frequently contained fibroblasts and the highest numbers of B cells (Fig. 1C). At the individual sample level, considerable variability was seen in the cell composition. Where serial samples were available from individual patients, the TME composition changed in response to therapy and disease progression (Fig. 1D).

Melanoma is uniquely sensitive to immunotherapy approaches. The nature and extent of immune infiltration (aka a “hot” or “cold” immune microenvironment) is thought to be predictive of outcome and response to immunotherapy (33–35). Important differences were noted in the T and NK cell profiles from each of the three metastatic environments (Fig. 2A). Notably, lymphocytes from brain metastases tended to have more overlap with those from the skin compared with the CSF. A high-level overview of the lymphocyte population indicated distinct patterns of gene expression between the three sites of metastasis (Fig. 2B). To understand the landscape of lymphocytes in the different TMEs in finer detail, we applied a second stage of clustering and curation to map the substructure of the cell populations, including all of the T/NK cells (Fig. 2C and D). A breakdown of the lymphocyte subsets revealed there to be a total of 17 subgroups, including 8 subclusters of CD4 cells [with one cluster of CD4 regulatory T cells (Treg)], 6 clusters of CD8 T cells, one cluster of gamma-delta (γδ) T cells, and two subgroups of NK cells (Fig. 2C and Supplementary Table S3 show how gene expression differentiates each cluster). To further determine the likely activation status of each identified T/NK cell subset, we next interrogated the data for known activation/exhaustion markers (Fig. 2E).

Among the major populations of T cells identified, clusters #1–4 and 7 expressed activation markers (cluster #3 being follicular CD4 T helper cells). Clusters #5–6 expressed low levels of activation markers and apoptosis genes and were judged to be inactive/exhausted (Fig. 2C and E). Cluster #8 was characterized as CD4 Tregs. Cluster #9 represented γδ T cells. Among the CD8 T cells, cluster #10 was dysfunctional and clusters #11–15 expressed activation markers. Clusters #16 and #17 were NK cells. It is worth noting that gene-expression profiles associated with classic T-cell activation markers varied among the three metastatic sites for several T and NK cell subclusters including clusters #1, 3, 6, 8, and 13–17 (Supplementary Fig. S1A).

A high-level view of the immune landscape demonstrated the LMM samples to have the highest proportion of CD4 T cells and the smallest numbers of CD8 cells. The MBM specimens had the lowest proportion of CD4 and NK cells, but a greater proportion of CD8 T cells (Fig. 2F). A more detailed breakdown revealed the greatest heterogeneity to be within the CD4 T cells, with a distinct organ-specificity observed (Fig. 2G). The LMM specimens were characterized by the highest numbers of exhausted and apoptotic CD4 cells. The brain samples contained activated CD4 T cells, with the skin also showing enrichment for the greatest proportion of proliferating CD4 T cells (Fig. 2F and G). Overall, there was greater overlap in the CD4 T-cell clusters between the skin and brain compared with brain and CSF (Fig. 2F and G). A similar heterogeneity was also observed in the CD8 T-cell proportions, with the LMM samples being highly enriched for Tem CD8 cells, whereas the brain and skin showed the highest proportion of cells in the highly activated cluster #13 (Fig. 2E–G). Analysis of the B cells showed the greatest cell numbers to be present in the skin samples, with the lowest numbers found in the LMM samples (Fig. 1C). A more in-depth analysis of the B-cell and plasma subsets showed the presence of activated B cells at all three sites (Supplementary Fig. S1B and S1C). The skin B-cell populations were enriched for markers of follicular B cells that were less prominent in the LMM and brain samples. Meanwhile, the B cells in the brain tumors expressed markers associated with memory B cells (Supplementary Fig. S1B and S1C). Overall, the lymphocyte landscape of the MBM was similar to that of the skin metastases, whereas the LMM was unique and characterized by more dysfunctional populations of CD4 and CD8 T cells.

The expansion of myeloid cell populations in tumors, including DCs, macrophages, monocytes, and myeloid-derived suppressor cells (MDSC), is a dynamic process that can lead to either antitumor or protumorigenic effects (36). A total of 32 clusters were found in the myeloid analysis, each annotated to one of seven broad subtypes: plasmacytoid DCs (pDC: #1), conventional DC types 1 and 2 (cDC1, cDC2: #2, #3), DC cluster 3 (DC3: #4), macrophages (#5), monocyte-derived DCs (moDC: #6), and myeloid-derived suppressor (MDSC: #7)-like cells (Fig. 3A–D; Supplementary Table S4). The expression of macrophage-specific markers was heterogeneous across the three metastatic sites (Fig. 3E; Supplementary Fig. S2A). It was noted that the CSF specimens contained the most diverse macrophage subpopulations and that these mostly expressed genes associated with alternative macrophage activation (Supplementary Fig. S2B–S2D). Expression of markers associated with specific DC subtypes was more homogeneous among the three metastatic sites (Supplementary Fig. S2A). The cDCs were of particular interest, as scRNA-seq has proven instrumental in defining their landscape and has identified several new subclasses (37). DC3s remain inadequately characterized and were identified through multistep clustering, which identified them as being distinct from the other DC cell subgroups, with a gene-expression profile that matched those previously identified as DC3s in human and mouse lung cancers (Fig. 3F and G; Supplementary Fig. S2A; refs. 37, 38).

Again, we noted diversity in the myeloid cell composition depending upon the metastatic site. The high-level overview identified a larger proportion of monocytic MDSC-like cells in the LMM and brain samples, higher percentages of macrophages in the brain and skin metastases, and a higher proportion of DCs in the LMM and skin (Fig. 3H). A more detailed breakdown showed the LMM samples to be high in cDC2s, and to have the highest proportion of moDCs along with the smallest proportion of macrophages of the three organ sites (Fig. 3I). The MBM specimens were high in immunosuppressive myeloid cell populations and low on possible antigen-presenting DCs. This site was very low in all types of DCs (cDC1, cDC2, and DC3) and had the highest proportion of MDSC-like cells, along with large numbers of macrophages (Fig. 3H and I). Among the myeloid cell populations, there was closer general alignment between the LMM and brain than the brain and the skin (Fig. 3H and I).

Among our cohort of 6 LMM patients, we had 5 individuals with the usual dismal outcome and 1 patient (patient F) who lived for over 38 months following their diagnosis of LMM, giving us a unique opportunity to define the differences between responders and nonresponders. Multiple CSF samples were available from each patient (often after therapy), along with two control CSF specimens from non-melanoma patients without any cancer in their leptomeninges (Fig. 4A; Supplementary Fig. S3A; Supplementary Table S2). It was noted that the exceptional responder had a lymphocyte profile closer to that of the non-tumor CSF controls and was characterized by a high proportion of CD4 T cells, CD8 Tem cells, and a higher proportion of DCs (pDCs, cDC2, and DC3) and moDCs, B cells, and plasma cells compared with the non-cancer controls (Fig. 4B–D; Supplementary Fig. S3A). By contrast, the poorly responding patients had a CSF environment that had high levels of inactivated CD4 T cells, MDSC-like cells, fewer pDCs, and more cDC2s. Neuronal cells and endothelial cells were also detected in their CSF (which may have been indicative of tissue damage associated with metastasis), which was lacking in the non-cancer controls and the exceptional responder.

An analysis of these serial CSF samples revealed new insights into the immune responses associated with LMM (Fig. 4D–F; Supplementary Fig. S3). Samples from the exceptional responder (patient F) revealed slow disease progression on dabrafenib/trametinib therapy, characterized by decreasing numbers of cDC2s and NK cells and increasing numbers of CSF-associated macrophages (samples 29, 23, and 30 in Fig. 4D), suggesting their potential role in regulating responses to therapy. At day 714, asymptomatic progression of disease was suspected, and the patient was started on nivolumab monotherapy and binimetinib–encorafenib. This switch in therapy was associated with increased infiltration of CD8 Tem cells, NK cells, and a decrease in the proportion of macrophages (sample 26; Fig. 4D) and the patient remained stable for >1 year following this intervention. In a different (poorly responding) individual (patient D), increasing tumor burden (between samples 15 and 16) was associated with decreased total immune cell numbers and higher numbers of moMDSCs (Fig. 4E). Treatment with radiation decreased melanoma cell numbers in the CSF and altered the myeloid composition to mostly moDCs; continued disease progression was associated with increased numbers of inactivated T cells (Fig. 4E). In another patient with a poor outcome (patient B; Fig. 4F), initial analysis of the CSF following the LMM diagnosis revealed the presence of inactivated T cells and a high percentage of macrophages. Following nivolumab treatment (sample 8), the myeloid compartment underwent a shift to increased DC numbers (pDC, cDC2, moDC) and MDSCs along with a continued inactivation of the T-cell compartment (Fig. 4F). At autopsy (sample 9), the patients' CSF showed high numbers of melanoma cells and macrophages. Similar patterns were seen in other poor responders (Supplementary Fig. S3B), with nivolumab treatment improving the myeloid compartment (increased cDC2s) but not positively impacting the T cells (e.g., the lymphocyte compartment was again characterized by inactivated T cells; Supplementary Fig. S3B). There are several potential limitations to sampling CSF, including low cell counts and the possibility that these cells do not accurately represent the adherent tumor at the leptomeninges. There is also the likelihood of a brain metastasis shedding cells into the subarachnoid space or ventricular CSF that may complicate analyses. Although we were not able to compare the transcriptional profiles of melanoma cells floating in the CSF to those adhering to the leptomeninges, we did perform an orthogonal validation of the immune cell type proportions between matched samples of CSF and specimens collected from the leptomeninges at autopsy (patient B). These analyses confirmed concordant cell numbers between the scRNA-seq of CSF and multiplexed immunofluorescence analyses of the thoracic spinal cord, demonstrating the potential of CSF sampling in the interrogation of the TME from LMM patients (Fig. 4G). Together, these results demonstrated that good responders with LMM had a more favorable immune environment that could be shaped following treatment with immune-checkpoint inhibitors, whereas poor responders were characterized by inactivated T cells and a less favorable myeloid compartment with more macrophages and MDSCs.

We next characterized how interpatient variability and therapeutic intervention shaped the immune environment of MBM (Fig. 5A–C). The samples were divided into those from individuals who had craniotomies (i) prior to any systemic therapy, (ii) after targeted therapy alone (BRAF-MEK inhibitor combination), (iii) after immune-checkpoint inhibitor therapy (anti–PD-1 and/or anti–CTLA-4), or (iv) after immune-checkpoint inhibitor therapy and targeted therapy (Fig. 5A–C; Supplementary Fig. S4; Supplementary Table S2). Although the patient immune landscapes were very heterogeneous, some trends were noted. Many pretreatment MBM samples had very few total immune cells (e.g., virtually no lymphocytes or myeloid cells) and were mostly tumor cells (e.g., MB-04, MB-05, MB-08, and MB-15; Fig. 5C). In most cases, systemic targeted therapy or immunotherapy was associated with increased numbers of infiltrating immune cells. Treatment with targeted therapy was associated with an increased ratio of CD8 T cells to CD4 T cells, which was less frequently seen in pretreatment specimens or following immunotherapy (Fig. 5C; Supplementary Fig. S4). A history of immunotherapy was associated with a more diverse lymphocyte landscape than seen in pretreatment samples, or in those from patients who had received only targeted therapy (Fig. 5C). The myeloid cell compartments of the immunotherapy and targeted therapy-treated brain metastases were very similar. Other striking differences included the targeted therapy–treated tumors being high in neuronal cells, whereas the immunotherapy-treated tumors had a high proportion of B cells and plasma cells, suggesting a sustained humoral antitumor immune response (Fig. 5C; Supplementary Fig. S4B). We then looked at the history of prior radiotherapy in these patients and noted that radiation-treated patients had fewer myeloid cells and a much less diverse myeloid compartment (Fig. 5C). To determine the fidelity of our scRNA-seq, we performed flow cytometry on eight craniotomy specimens (for which sample remained) and found good agreement in the numbers of multiple cell types between the two approaches, including NK cells, CD4 T cells, CD8 T cells, DCs, monocytes, macrophages, and B cells (Fig. 5D).

Our analyses suggested that many drug-naïve MBMs had little immune infiltrate and that this increased in response to systemic therapy. To explore this in more detail, we analyzed a rare case in which two consecutive craniotomies were obtained from an MBM patient pre-systemic therapy and then after a short course of systemic BRAF-MEK inhibitor therapy (2 months of therapy 1 year prior to the second craniotomy) followed by anti–PD-1 therapy (nivolumab 3 weeks prior to the second craniotomy). These matched samples clearly demonstrated very few NK cells, CD4 and CD8 cells, DCs, or macrophages (but high levels of MDSCs) in the pretherapy sample and a dramatic increase in immune infiltrate (T cells, NK cells, DCs, macrophages) after systemic therapy (Fig. 6A and B). Fibroblasts were also detected (possibly from the scarring following the first surgical procedure), as were limited numbers of B cells. These analyses suggest that systemic treatment can increase the immune infiltration of otherwise immunologically “cold” MBMs. To confirm this in a more controlled setting, we performed preclinical experiments in which mice with MBM (BRAF-mutant mouse SM1 cells injected intracranially) were treated with anti–PD-1 therapy (Fig. 6C), and the extent of tumor immune infiltrate was measured. It was noted that systemic anti–PD-1 therapy increased multiple immune cell types, including activated CD4 and CD8 T cells, MDSCs, Tregs, and DCs, consistent with the scRNA-seq findings in patient specimens (Fig. 6D). The observation that some MBMs had little immune infiltrate prior to systemic therapy suggested that post-therapy immune cells may have trafficked from the periphery to the brain. To determine the congruence between extracranial and intracranial immune responses we next compared the immune infiltrate from MBM to that of melanoma skin metastases. It was noted that although lymphocyte and myeloid cell numbers were significantly lower in the brain compared with the skin, the phenotypic make-up of the lymphocyte compartment was quite similar (despite the samples being from different individuals; Fig. 6E and F). The proportions of myeloid cells differ slightly, with the brain having a greater proportion of macrophages and a much-reduced percentage of DCs (Fig. 6E and F).

Having demonstrated that systemic therapy markedly altered the immune environment of LMM, and MBM samples, we next defined how each of the cell subtypes affected the other cell types in the tumor landscape (Fig. 7A). Among the DCs, cDC1, cDC2, and DC3 correlated with reduced numbers of melanoma cells, macrophages, endothelial cells, and fibroblasts. DC3s showed the strongest correlation with a diverse array of CD4 and CD8 T cells, with cDC2s showing a more restricted pattern of T cell correlation. It was noted that the presence of DC3 cells was associated with higher expression of MHC class I and II molecules in melanoma cells (Supplementary Fig. S5). As DC3s remain poorly characterized, we next performed an in silico analysis of ligand–receptor interactions between DC3s and multiple T-cell subsets, and identified high interaction scores between DC3s and two subclasses of activated tissue-resident memory CD8 T cells and CD4 follicular helper cells (Supplementary Fig. S6), suggesting that DC3s may be driving an active T-cell response. DC3s were detected in 10/43 samples, from 8/24 melanoma patients analyzed (33%), but were not observed in the CSF from non-LMM individuals (Supplementary Fig. S3A; Supplementary Table S2). Multiplexed IF analysis was undertaken to confirm the presence of this rare cell type across multiple MBM specimens that had matched scRNA-seq data (Fig. 7B; Supplementary Fig. S7). There was good agreement between the scRNA-seq and multiplex IF data, particularly as these cells were of very low frequency (typically 1–68 cells per sample; Fig. 7B). Interestingly, our analyses of our limited patient cohort indicated that the presence of DC3s might be associated with better patient outcomes (based on best clinical response observed to systemic therapy; Fig. 7B). We performed an association analysis to investigate whether the presence of each cell population was associated with survival outcomes in melanoma patients (N = 24 patients), using site of metastasis as a covariate in a Cox proportional hazard model. We identified the DC3 cells, along with cluster #1 of CD4 T cells and γδ T cells, to be associated with longer OS times (Supplementary Table S5). Because the DC3s were rare, we calculated the proportion of these cells in the total myeloid cell population and used this proportion in the survival analyses to obtain a reliable estimate of hazard ratio. The survival analyses showed the higher mean DC3 proportion to be associated with better survival outcome (HR = 0.009, P = 0.019). We then interrogated two additional independently published scRNA-seq data sets: one from human PBMCs and another from 19 melanoma patients (29, 30) using the same multistep clustering approach. In the Tirosh melanoma validation data set, we identified 10 T/NK and 25 myeloid cell subpopulations (Supplementary Figs. S8 and S9). We identified DC3s in 21% of the 19 melanoma patients, but not in the normal PBMC population. A robust gene-expression signature was associated with the DC3 cell population consisting of IDO1, IDO2, CCL19, and LAMP3 (Fig. 3D; Supplementary Fig. S10). These 4 genes were positively correlated with one another in the melanoma TCGA data set (n = 495) and PCA was performed to generate a signature to summarize the gene expression of these 4 markers (Supplementary Fig. S10A). Given the known survival differences between LMM patients and other melanoma patients, survival differences were visualized and the optimal cut-points determined separately for LMM and other patients of our scRNA-seq cohort (Fig. 7C; Supplementary Fig. S11). We evaluated the expression of these 4 genes and showed them to predict for better OS in the melanoma TCGA data set both individually and through use of PC1 signature (Fig. 7D; Supplementary Fig. S12). Analysis of an additional bulk RNA-seq data set from a Moffitt cohort containing 135 melanoma patients further confirmed the positive correlation among the 4 genes and the presence of DC3 cell markers to be associated with better OS in melanoma patients (Fig. 7E; Supplementary Fig. S10C).

To address if the presence of DC3s correlated with favorable therapeutic response in preclinical models, we analyzed the scRNA-seq data from a recently published data set from our lab in which a syngeneic mouse melanoma model was treated with (i) targeted therapy, (ii) immunotherapy, or (iii) immunotherapy and targeted therapy in sequence (31). It was noted that the best antitumor responses were seen to the sequence of immunotherapy followed by targeted therapy and that these were the only tumor specimens to harbor any DC3 cells (Fig. 7F). Taken together, these data highlight a role for DC3s in positively regulating the melanoma immune microenvironment that are associated with improved survival across multiple independent patient cohorts.

Melanoma is a heterogeneous disease, with differences in the speed of tumor progression and therapeutic responses frequently observed between different metastatic sites within the same patient (11, 39). There is evidence that baseline and therapy-induced levels of immune infiltration contribute to this heterogeneity through mechanisms that remain poorly defined (33, 34, 40).

Until now, the immune environment of the LMM has never been explored at a single-cell resolution. Our in-depth analysis of the immune landscape of 2 sites of melanoma CNS metastasis revealed some striking new findings. A first key observation was that the phenotypic make-up of the CD4 and CD8 T-cell populations in the brain showed a high degree of overlap with those of skin metastases (albeit with much fewer total T-cell numbers in the brain vs. skin), whereas the gene-expression profiles of CD4 and CD8 T cells from LMM patients were quite distinct. The immune environment of LMM was highly immune suppressed and enriched for populations of T cells that were predicted to be exhausted or inactivated. Although CSF was long considered to be an immunologic desert, normally containing less than 5 white blood cells per mm³ (41), there are a few reports of immune responses in the leptomeninges, such as the presence of clonally expanded CD8 T cells in the CSF of patients with Alzheimer's disease (42).

The CSF of one individual with LMM who survived long term (>38 months following diagnosis) differed from that of patients with poor survival and had an immune landscape closer to that of CSF from non-LMM patients, although with a higher overall level of immune infiltrate. Evidence of a more active immune environment in this individual was supported by analysis of serial CSF specimens, in which responses to anti–PD-1 and BRAF/MEK inhibitor therapy were accompanied by an increase in T-cell diversity, NK cells, and a reduction in suppressive macrophages. The likely role of macrophages in contributing to an immunosuppressive environment of the CSF space was suggested by our correlation analyses, which showed their presence to be negatively correlated with cDC1s, CD8 Tem cells, CD8 cycling Trm cells, γδ T cells, and multiple clusters of CD4 T cells. Recent preclinical work has supported a role for PD-1 inhibition in shaping myeloid lineages away from immunosuppressive phenotypes into differentiation programs that favor systemic tumor immunity mediated through DCs and effector monocytes (43). In contrast, systemic therapy in LMM patients with a poor outcome led to changes in the myeloid cell compartment but did little to improve the lymphocyte response. In patients with rapidly progressing LMM, the balance of the CD4 T cells shifted toward phenotypes that were exhausted and nonfunctional. Our previous studies on the LMM microenvironment demonstrated the CSF to contain high levels of complement, other innate immune components, and immunosuppressive growth factors such as transforming growth factor (TGF)-β (44). It is therefore possible that circulating factors and suppressive myeloid cells work together to reprogram the T cells in the leptomeningeal microenvironment into a dysfunctional state.

There is already evidence from studies on autoimmunity and from high-dimensional analyses of T-cell trafficking that individual organs harbor unique populations of T cells with distinct cytokine secretion profiles (45, 46). The similarity of the T-cell populations between the skin and brain metastases was intriguing and suggested that the T-cell population infiltrating MBM may have been educated extracranially, before migrating to the brain. This idea was supported by preclinical studies in mouse models demonstrating that systemic PD-1 therapy increased the infiltration of multiple immune subsets into the brain. The observation that many MBMs contained virtually no DCs or other antigen-presenting cells prior to systemic therapy, likely contributed to the lack of immune response at the intracranial site. Studies on the immune environment of matched cranial and extracranial melanoma metastases have shown that brain metastasis samples have less TCR diversity than the extracranial samples (20). Despite this, some patients did have significant immune infiltration in their brain metastases, despite no history of prior systemic therapy. These differences could be a reflection of differences in blood–brain barrier integrity between patients or possibly even the anatomic location of the tumor within the brain, with some sites offering the potential for better immune priming.

Our single-cell data set gave us a unique opportunity to better understand how the presence of individual immune cell types shaped the entire tumor immune environment. Recent studies have implicated increased B-cell numbers, associated with tertiary lymph node–like structures, and the expansion of large clones of effector memory CD8 T cells as potential biomarkers of immunotherapy response in melanoma (33, 35, 47). Association analyses allowed us to identify a minor subset of DC3s that positively shaped the tumor microenvironment and predicted for better OS in three independent patient data sets. We showed that the interactive tool, ISCVA, along with the multistage clustering approach, was able to identify this small subpopulation of DC3 cells, at a similar frequency, in a second independent melanoma single-cell data set (30). Of note, our study and the prior study of Tirosh and colleagues used different analysis platforms (10X Genomics vs. SMART-seq), validating the broad utility of our ISCVA analysis pipeline in analyzing single-cell data sets and robustness of the identified signatures across platforms.

DCs are the most potent subtype of antigen-presenting cells and are highly effective at activating naïve T cells. Our scRNA-seq analysis has identified four distinct subsets of DCs, including pDCs, cDC1, cDC2, and DC3 (37). Among these, DC3s remain relatively poorly characterized but are known to exhibit an “activated” DC phenotype (37, 48). In prostate cancer, the activation of both cDC1s and DC3s correlates with increased patient survival following treatment with ipilimumab and GVAX (49). In terms of gene-expression profile, DC3s show some overlap with cDC1-expressed genes such as BATF3 and IRF8, but do not express others, including XCR1 and CLEC9A (37). We found that the DC3 cells in melanoma patient specimens expressed high levels of the costimulatory molecules CD80 (B7-1) and CD86 (B7-2) and their presence was associated with reduced numbers of melanoma cells and immunosuppressive stromal and myeloid cells, including endothelial cells, fibroblasts, and macrophages. DC3s were also correlated with increased numbers of B cells, plasma cells, follicular helper CD4 T cells, NK cells, pDCs, and activated clusters of CD8 T cells. To better understand the role of DC3s in our melanoma specimens, we inferred potential cell–cell interactions using the scRNA-seq data and uncovered ligand–receptor pairs between DC3s and the T-cell subsets. These analyses suggested that DC3s were likely to interact with some of the most active subsets of CD8 T cells, and that the presence of DC3s in the melanoma predicted for high levels of MHC expression in the tumor cells. These data support the idea that the presence of DC3s help to drive a sustained immune response that correlated with better survival outcomes. Our observations in the clinical samples were supported by preclinical experiments in syngeneic mouse models that showed better therapeutic responses to be associated with the accumulation of DC3s in the tumors.

Analysis and visualization of complex single-cell RNA-seq data encompassing 20–30 thousand genes within tens of thousands of cells at a glance is a challenge, particularly for non-experts (such as biologists and clinicians). The most commonly used platform, Loupe, developed and provided by 10X Genomics, is powerful, but does not provide cell type annotation. Here, we developed our own ISCVA with a two-stage architecture. The first stage utilizes recently developed algorithms, which have been peer reviewed and published for single-cell data processing, QC, cluster analyses, and cell typing and commonly used by the single-cell research community [e.g., Seurat for general processing (25), SingleR for cell type recognition (26), and single-cell signature explorer for gene set signature scoring (27, 50)]. The second stage contains web-based components, built with state-of-the-art web front-end technologies, including react.js from Facebook, tensorflow.js from Google, and Plotly.js that allows convenient real-time interactive exploration and ad hoc analysis. Our platform, which is publicly available online (on a purpose-built website), allows investigators with no prior training to visualize multiple parameters including specific genes, sample types, and cell types from their single-cell RNA-seq data in different projections (e.g., t-SNE and UMAP). Figures can be easily exported for presentation and publication and indeed this platform was used to generate many of the figures in this manuscript.

Our study does have a number of limitations. Our data are focused on a cohort of patients who are mostly BRAF mutant (21/24) and it is possible that patients with other genetic drivers of melanoma may have different cellular landscapes. We were also limited in analyzing skin metastases as our only site of extracranial metastasis, and it is likely that other extracranial organ sites may have different immune/cellular landscapes. Furthermore, because the patients who typically undergo surgical resection of brain metastases have tumors that are not responding to therapy, have hemorrhaged, or have significantly increased in size, it is possible that our MBM cohort is skewed toward patients with worse clinical outcomes. It is worth noting, however, that the OS for our patient cohort ranges from 7.9 months to 75.7+ months, with the median being around 27.6 months, which is similar to that of other MBM patients undergoing treatment with systemic targeted and immune-checkpoint therapies (51). Furthermore, our focus on the floating melanoma cells in the CSF of LMM patients does not address the phenotypic or transcriptional properties of the larger proportion of melanoma cells that adhere to the surfaces of the leptomeninges. Despite these caveats, our studies have provided unique new insights into the immune environment of melanoma metastases in these two CNS compartments and in particular how systemic therapy helps shape these. We further identified a novel population of DC3s that correlated with improved survival and therapy response in our patient cohort and in two additional validation data sets. This work represents the first atlas of two distinct CNS sites of melanoma metastasis, offers important new insights into the biology of this devastating disease, and provides a novel data analysis and visualization platform for interrogating scRNA-seq data sets.

I. Smalley reports grants from NIH during the conduct of the study. A. Sarnaik reports grants, personal fees, and nonfinancial support from Iovance Biotherapeutics Inc; grants from Provectus Inc; and personal fees from Guidepoint Inc, Defined Health Inc, Gerson Lehrman Group Inc, Physicians' Education Resource LLC, Medscape Inc, and Medstar Health outside the submitted work; in addition, A. Sarnaik has a patent for culture of tumor-infiltrating lymphocytes from tumor digest, filed March 24, 2021, US Patent Application No. 17/279,327 pending. V.K. Sondak reports other support from Neogene Therapeutics, as well as personal fees from Merck, BMS, Eisai, Replimune, Regeneron, and Novartis outside the submitted work. Z. Eroglu reports grants and personal fees from Novartis; personal fees from Genentech, Array, Natera, OncoSec, Regeneron and SunPharma; and grants from Pfizer outside the submitted work. P.A. Forsyth reports other support from AbbVie, BTG, Inovio, Novocure, Boehringer Ingelheim, NCI Neuro-Oncology Branch Peer Review, NCRI, NIH, Novellus, Physician Sciences Oncology Network, and Ziopharm, as well as grants from CDMRP research grant, Department of Defence, Moffitt Center of Excellence Celgene project, NIH/NCI, Pfizer, State of FL Bankhead Coley, NIH/NCI 1R21, Moffitt CBMM, Moffitt Foundation Research Acceleration Fund/Schulze, and FL Academic Cancer Center Alliance outside the submitted work. K.S.M. Smalley reports personal fees from Elsevier outside the submitted work. No disclosures were reported by the other authors.

I. Smalley: Conceptualization, supervision, writing–original draft, project administration, formal analysis, funding acquisition, investigation, writing–review and editing. Z. Chen: Data curation, software, formal analysis, validation, investigation. M. Phadke: Investigation. J. Li: Data curation, software, formal analysis. X. Yu: Data curation, software, formal analysis. C. Wyatt: Investigation. B. Evernden: Resources. J.L. Messina: Resources, validation. A. Sarnaik: Resources. V.K. Sondak: Resources, formal analysis, writing–review and editing. C. Zhang: Formal analysis, methodology. V. Law: Resources, formal analysis. N. Tran: Resources. A. Etame: Resources, writing–review and editing. R.J.B. Macaulay: Resources, formal analysis, validation. Z. Eroglu: Resources. P.A. Forsyth: Conceptualization, resources, investigation, writing–review and editing. P.C. Rodriguez: Conceptualization, formal analysis, investigation, writing–review and editing. Y.A. Chen: Conceptualization, resources, data curation, software, formal analysis, supervision, methodology, writing–original draft, writing–review and editing. K.S.M. Smalley: Conceptualization, resources, formal analysis, supervision, funding acquisition, investigation, writing–original draft, project administration, writing–review and editing.

We would like to thank the patients and their families for their valuable contributions to this study. We would like to acknowledge Paige Carbon for technical assistance with the sample collection. We thank Carlos Moran Segura, Jonathan Nguyen, and the Advanced Analytical and Digital Pathology team at Moffitt Cancer Center for technical assistance with tissue staining and analysis. We would like to thank Jodi Kroeger and the Flow Cytometry Core at Moffitt Cancer Center for technical assistance with immunophenotyping analyses. This work was supported by NIH grants P50 CA168536, R21 CA198550, R21 CA216756 (to K.S.M. Smalley), K99 CA226679 (to I. Smalley), the Department of Defense W81XWH1810268 (to K.S.M. Smalley), the Melanoma Research Foundation (to K.S.M. Smalley), and Moffitt DRP Innovation Core Project Award (to Y.A. Chen and X. Yu). The Molecular Genomics, Flow Cytometry, Advanced Analytical and Digital Pathology and Bioinformatics and Biostatistics Core at Moffitt are supported in part by the NCI through a Cancer Center Support Grant (P30-CA076292) and the Moffitt Foundation. We would like to thank Lesa Kennedy and Bill Christy for their generous support of the Moffitt Melanoma and Skin Cancer Center of Excellence that significantly contributed to these studies.

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.