Linking Transcriptomic and Imaging Data Defines Features of a Favorable Tumor Immune Microenvironment and Identifies a Combination Biomarker for Primary Melanoma

Melanoma is the most aggressive and lethal form of skin cancer (1). Even when completely surgically resected, tumors less than 5 millimeters in thickness are commonly lethal due to tumor seeding into the lymphatics and systemic circulation (2, 3). Once metastatic disease develops, survival rates range from 20% to 40% at 5 years with state-of-the-art immunotherapy (4, 5). While these rates are disappointing, they represent a vast improvement over survival rates at 5 years prior to immunotherapy, which were under 5% (6). Today, immunotherapy is given as the standard of care to virtually all patients diagnosed with metastatic melanoma (7–9).

Early-stage melanoma presents a more difficult clinical management challenge because the side effects of treatment are less justifiable in the absence of high risk of death. Under 5% of patients with stage I melanoma will die of disease, while 10-year survival rates range from 88% to 75% for stage II disease and from 88% to 24% for stage III disease (2). Prediction of recurrence would thus greatly facilitate management of early-stage disease.

Melanoma is an immunogenic tumor, and interactions between tumor cells and the tumor immune microenvironment (TIME) lead to tumor-driven modulation of immune response (10). It has been hypothesized that the immune system is a key determinant of melanoma recurrence and metastasis. Tumor-infiltrating lymphocytes (TIL) have been studied extensively and shown to have positive prognostic value in melanoma (11–13). However, TILs are not currently included in the American Joint Committee on Cancer staging system for melanoma, and their utility is limited by a lack of phenotypic, functional, and precise locational information (10, 14). Furthermore, tumor-associated macrophages (TAM) within the TIME have been shown to correlate with poor prognosis and survival in melanoma and many other cancers (15, 16). Thus, understanding the melanoma TIME is crucial for both the development of new therapeutic interventions, as well as for stratification of patients for adjuvant therapy.

In the previous work, we defined and validated in two independent cohorts a 53-immune gene melanoma immune profile (MIP) predictive of nonprogression using NanoString transcriptomic analysis (17, 18). MIP allows for stratification of patients with stage II-III melanoma into recurrence risk groups. Adding to this work, we utilized quantitative multiplex immunofluorescence (qmIF) to analyze the tumor microenvironment in patients with stage II-III melanoma and found that a high density of cytotoxic lymphocytes (CTL) and a low density of macrophages in the stroma correlates with improved survival. Furthermore, we found that the distance between CTLs and macrophages is associated with poor survival (19), highlighting the importance of both the composition and spatial localization of cells within the melanoma TIME.

We now expand our analysis of the TIME of primary stage II-III melanoma by combining NanoString and qmIF data. Using advanced statistical techniques, we analyze the clustering, cellular composition, and marker proteins of both primary and metastatic melanomas. In particular, to quantify the degree of spatial proximity among immune and tumor cells, we assess clustering using the pair correlation function (PCF). This technique was originally used in statistical physics to model particle interactions (20) and we recently adapted it to study differences in spatial structures in immunofluorescence imaging of glioblastomas (21). The PCF computes the average number of pairs of two given cell types separated by a given distance (normalized by the number expected if all cells were randomly positioned), giving a robust, quantitative metric of spatial interactions. This more detailed characterization may serve to guide the development of targeted and immune therapies for melanoma.

This study was approved by Columbia University Irving Medical Center's (CUIMC) Institutional Review Board. A database of patients with stage II-III melanoma at CUIMC (New York, NY) was retrospectively created and analyzed using qmIF on mantra and RNA expression using NanoString for MIP (Table 1, mantra: n = 104 and MIP: n = 78; refs. 18, 19). To correlate the results of qmIF with RNA expression, we analyzed the patients who overlap between the mantra and the MIP cohorts (Table 1, Combination: n = 53). Furthermore, a subcohort of 13 patients from the mantra cohort, named the macrophage survey subcohort (Supplementary Table ST1), was selected on the basis of high density of CD68⁺HLA-DR⁻ macrophages for additional phenotypic characterization. Finally, a database of patients with metastatic melanoma at CUIMC was created by screening the records of 86 patients with metastatic melanoma, of whom 57 were excluded because of lack of availability of solid tumor samples prior to immunotherapy treatment or insufficient clinical information. Tissue was requested for the 29 patients remaining, then slides were cut and reviewed with a dermatopathologist. A total of 9 patients were excluded at this stage due to lack of sufficient tumor or scattered tumor in the lymph node. The remaining 20 patients were stained for inForm analysis. Six of these patients were excluded because of poor staining quality or tissue destruction, leaving a final total of 14 patients analyzable on inForm (Supplementary Table ST2).

QmIF imaging was obtained from 104 stage II-III primary melanomas from CUIMC in the mantra cohort (Table 1; ref. 19). Five-micron thick full-section slides were stained using Opal qmIF for DAPI, CD3 (T lymphocytes; clone LN10; Leica; 1:200 dilution), CD8 (cytotoxic T lymphocytes, or CTLs; clone 4B11; Leica; ready to use: RTU), CD68 (macrophages; clone KP1; Biogenex; RTU), SOX10 (tumor; clone BC34; Biocare; 1:300), HLA-DR (antigen presentation and activation; clone LN-3; Abcam; 1:200 dilution), and Ki67 (proliferation; clone MIB1; Abcam; RTU). H&E slides were reviewed by a dermatopathologist to determine representative areas for multispectral image capture at ×20 magnification using mantra. Regions of interest (ROI) were approximately 700 × 520 μm. Multispectral images were then analyzed using inForm software and inForm data from each patient was processed in separate software designed in RStudio (version 0.99.896; https://github.com/thmshrt/transform_essential). In this software, images were combined and analyzed to concatenate variables and determine density and distance of distinct phenotypes. Statistical analyses were then completed on GraphPad Prism, R Version 3.3.1, and Python 3.5. For visualization, intensity values of each of the 5 macrophage markers were first normalized to have a SD of 1. These continuous values were processed with T-distributed stochastic neighbor embedding (tSNE; perplexity = 30), a dimensionality reduction technique that places cells with similar marker profiles close to each other (22). Binarized (±) expression (derived with thresholds using inForm) for each marker was then plotted onto the tSNE figures. Coexpression of markers with HLA-DR was assessed with the χ² test. In addition, the macrophage survey subcohort (n = 13) was stained for DAPI, CD68, HLA-DR, CSF1R (inflammation; clone SP211; LifeSpan Biosciences; 1:750 dilution), CD163 (M2 macrophages; clone 10D6; BioCare Medical; 1:300 dilution), CD33 (myeloid-derived suppressor cells, or MDSCs; clone PWS44; Leica; RTU) and MHC-I (antigen presentation; clone MEM-E/07; Thermo Fisher Scientific; 1:400 dilution) and analyzed in the same way as described above.

To assess clustering in qmIF images, PCFs were computed using the spatstat R package (23). This method compares the observed probability of two cell types being separated by a given distance to that expected by chance. A higher value for PCF indicates a greater degree of clustering at that radius, and the AUC of the function can be used as an overall statistic of clustering. Inhomogeneous cross-PCFs were calculated up to a radius of 50 microns, provided that there was a minimum of 10 immune cells in the image. One sample was excluded for lack of tumor cells. The isotropic edge correction and a normalization power of 2 were applied. The AUC for each PCF was computed as a summary statistic for quantifying the degree of clustering, and was compared across recurrence status with the Mann–Whitney U test. The PCFs are graphically displayed by the point-wise median across samples, with shaded 95% confidence intervals (CI) obtained via bootstrapping.

For each of the 78 patients in the MIP cohort (Table 1, MIP), formalin-fixed paraffin-embedded (FFPE) blocks were cut in 10-μm sections to provide 250 mm² of tissue. RNA extraction was performed and RNA was processed using a NanoString assay, which measures expression of 64 target genes [MIP (53 genes) plus 11 IFN-related and 17 housekeeping genes; ref. 24]. Following hybridization and purification of target–probe complexes, digital counts for each gene-specific target RNA were acquired on the nCounter detection analyzer and normalized in nSolver software (NanoString). MIP was computed for each patient in this cohort using the original signature's gene coefficients, then was correlated with disease-specific survival (DSS; refs. 17, 18).

We assessed the similarity of each sample's gene expression profile to the cellular subtype signatures from LM22, the reference standard for CIBERSORT (25). A total of 20 genes overlapped between the NanoString panel and LM22, which was a not sufficient number of genes to deconvolute with the CIBERSORT algorithm (25). Instead, we calculated the cosine similarity of the samples with each of the 22 cell types from LM22, thereby assessing the extent (from −1 to 1) that the matching genes from each sample resemble each cell type. For the creation of the heatmap, raw gene expression values were log-transformed and scaled to have a mean of 0 and a SD of 1. Hierarchical complete-linkage clustering was performed on both rows and columns with a cosine metric. We assessed enrichment of recurrence among this clustering by dividing the patients at the first branch point into the two largest clusters, followed by a Fisher exact test. Differential expression was evaluated for all 63 genes with t test with Benjamini–Hochberg (BH) correction, and visualized in a volcano plot. Top genes were noted if they had greater than +1.0 log₂-fold change.

To assess the ability of the CD8/CD68 ratio to predict recurrence, we performed a Mann–Whitney U test and visualized the prediction using a ROC curve. We tested the correlation of the top immune genes (>1 log₂-fold change) with the CD8/CD68 ratio using Spearman correlation. We combined the CD8/CD68 ratio with MIP using the log-rank (Mantel–Cox) test to compare low- to medium-risk, low- to high-risk, and medium- to high-risk patients. HRs and 95% CIs were calculated using the Mantel–Haenszel test.

A total of 8 patients from the 104-patient mantra cohort were selected on the basis of CD68⁺HLA-DR⁻ phenotype and CD8 expression (Supplementary Table ST3). Of this subcohort, 4 patients had a high density of stromal HLA-DR⁻ macrophages and close proximity of HLA-DR⁻ macrophages to CTLs, and 4 patients had a high CTL/macrophage ratio with low HLA-DR⁻ macrophage density, as determined by previous qmIF. For digital spatial profiling (DSP), FFPE slides were stained with antibodies conjugated to UV-photocleavable DNA barcodes and specific to 34 proteins. Twelve 600-μm ROIs with high macrophage density were selected for each patient. ROIs were selected in tumor/stroma areas containing less than 50% stroma. ROIs were analyzed using UV excitation, which releases DNA barcodes for downstream quantitation on the NanoString nCounter platform. Protein expression was compared between patients, and statistical analyses were performed using a Mann–Whitney U test.

The mantra cohort consists of 104 patients, 64 of whom had available DSS data, defined as known cause and date of death and/or documented clinical follow-up (19). The MIP cohort consisted of 81 patients, all of whom had available DSS (18). The combination cohort consisted of 53 patients, all of whom had available DSS. In the mantra cohort, 15.4% of patients had known recurrence, 34.6% had no recurrence, and 50% had unknown recurrence status. In the MIP cohort, 29.6% of patients recurred, 61.7% did not recur, and 8.6% had unknown recurrence status. Finally, in the combination cohort of 53 patients, 30.2% of patients recurred, 60.4% of patients did not recur, and 9.4% of patients had unknown recurrence status. The combination cohort consisted of 75.5% males and 24.5% females, and the median patient age was 68. 56.6% of patients had a tumor on their trunk or head, and 72.9% of patients were stage II. Median tumor depth in this cohort was 2.7 mm, and ulceration was present in 60.4% of the cases. Median follow-up was 62 months, and 25 of the patients died of any cause, with 14 having died from melanoma. The 13-patient macrophage survey subcohort stained for macrophage phenotypes consisted of 46.2% of patients with a high CD8/CD68 ratio, 38.5% of patients with a low HLA-DR⁻ macrophage density, and 61.5% of patients with a large distance between CTLs and HLA-DR⁻ macrophages (Supplementary Table ST1). The metastatic melanoma cohort consisted of 14 patients, of whom 1 had stage IIIC disease, 3 had stage IV M1a or M1b disease, 1 had stage IV M1c disease, 8 had stage IV M1d disease, and 1 had mucosal melanoma. The tissue origins of the biopsies of metastases used for staining were divided into 50% soft tissue or skin, 35.7% brain, 7.1% lung, and 7.1% lymph node (Supplementary Table ST2). Finally, the cohort used for DSP analysis consisted of 50% of patients with a high CD8/CD68 ratio, 50% of patients with a low HLA-DR⁻ macrophage density, and 62.5% of patients with a large distance between CTLs and HLA-DR⁻ macrophages (Supplementary Table ST3).

First, we evaluated the mantra cohort for distant metastatic recurrence using the 52 patients with known recurrence status and found that, in addition to correlating with DSS and overall survival (OS), as previously published, the CD8/CD68 ratio is also predictive of recurrence (n = 52, AUC = 0.684, P = 0.0350 by Mann–Whitney U test, Supplementary Fig. S1). Next, we applied qmIF to compare density of immune phenotypes in samples of primary (P) melanoma (stage II-III, mantra cohort) to those with metastatic (M) lesions. We determined the proportions of tumor, macrophage, CD8⁺ T cells, and CD8⁻ T cell as a fraction of the total number of cells as identified by DAPI in each sample, and concatenated these proportions with HLA-DR and Ki67 expression (Fig. 1A). We further quantified and compared the densities of CD8⁺ T cells, macrophages, total immune cells (sum of CD3⁺ T cells plus CD68⁺ macrophages), and CD8/CD68 ratios in primary (P) and metastatic (M) patients (Fig. 1B). As a whole, there was a greater proportion of immune cells in P samples compared with M samples (P = 0.0017, Mann–Whitney U test, Fig. 1B), with specifically more CTLs in primary samples (P = 1.8e-6). Next, we separated the primary melanoma cohort into nonrecurrent (NR, n = 36) and recurrent (R, n = 16) patients, finding that NR samples had a higher CD8/CD68 ratio relative to R patients (P = 0.024). Interestingly, there was no significant difference between the R and M samples (P = 0.52; Fig. 1B).

In further analysis, we evaluated nuclear size by comparing nuclear areas and found that tumor cells in M samples had larger nuclear areas compared with their P counterparts (P < 0.001, Mann–Whitney U test, Fig. 1C). This held true both for Ki67⁻ and Ki67⁺ subpopulations. However, there were no significant differences in the sizes of the nuclei of CTLs between the cohorts (Fig. 1C).

As a second element of our composite analysis of the TIME, we aimed to evaluate an IFN signature in our tumor samples from the 74 patients with known recurrence status in the MIP cohort while including 11 additional IFN-related genes for comparison. This was done for the purpose of correlating the gene expression data with qmIF data. We profiled RNA expression across 64 immune-related genes (53 genes from MIP and 11 additional IFN-related genes) and found that expression was associated with recurrence status using hierarchical clustering (P = 0.049, Fisher exact test, Fig. 2A). To identify the cellular source of this association, we assessed the degree of similarity between the transcriptomic signatures in our samples and the CIBERSORT cell reference signatures (25). In evaluation of the signatures for macrophages of M1 and M2 polarity, we found that R patients had a decreased M1 signature (P = 0.042, Fig. 2B top left) but an increased M2 signature (P = 0.048, Fig. 2B, top right) compared with NR patients. We then compared the CIBERSORT signatures with findings from qmIF. We found that the CD8 T cell signature correlated with the proportion of CTLs in qmIF (P = 0.010 Fig. 2B, bottom left) while the M1 signature correlated with density of HLA-DR⁺ macrophages (P = 0.022, bottom right).

Next, we evaluated differential expression of the 64 genes and found that nearly all genes were increased in patients who did not die of melanoma (Fig. 2C). After filtering for genes that were both significant following BH correction and also had greater than +1.0 log₂-fold change in differential expression, we identified 9 genes, all of which are part of MIP (Table 2). We also found that the top 9 genes do not overlap with the 11 additional IFN-related genes and further that these 11 genes did not improve prediction accuracy. Finally, we evaluated the correlation of these 9 genes with stromal density of CD3, CD3⁺CD8⁺, and CD68⁺, as well as the ratio of CD8/CD68 (Table 2). We found that CD3⁺ cells in the stroma correlated with CXCL9, CD8A, CXCR3, CXCL11, CCL5, and CD2 while CD3⁺CD8⁺ cells in the stroma correlated with the same genes minus CD2. Stromal CD68 cells and the CD8/CD68 ratio had no correlation with the top 9 genes from MIP.

Given that genes from MIP did not directly correlate with the CD8/CD68 ratio, we chose to combine these two biomarkers to substratify the patients. We evaluated the ability of the combination of MIP and the CD8/CD68 ratio to predict DSS using the 53-patient combination cohort (Table 1). We stratified patients into high-, medium-, and low-risk groups using their categorization as high or low risk in each of the datasets based on our previously published criteria (18, 19). High-risk patients had poor MIP and CD8/CD68 profiles; medium-risk patients had a poor MIP profile and a favorable CD8/CD68 profile; and low-risk patients had favorable MIP and CD8/CD68 profiles. We found that combining MIP and the CD8/CD68 ratio effectively distinguished low- and medium-risk patients from high-risk patients (P = 0.003; HR: 7.53; 95% CI: 1.98–28.74 and P = 0.026; HR: 3.51; 95% CI: 1.16–10.61, respectively, Fig. 3A and Supplementary Table ST4). Representative images of each of these risk groups are also shown in Fig. 3B–D. In addition, by creating a medium-risk group we found that the combination biomarker improves on the previous biomarkers, MIP and CD8/CD68 ratio, individually (P = 0.030; HR:3.84; 95% CI: 1.14–12.96 and P = 0.044; HR: 3.04; 95% CI: 1.03–8.99, respectively; Supplementary Table ST4). Finally, we found using univariable Cox analysis in medium- and high-risk patients (n = 42) that the high-risk combination biomarker specifies poor prognosis (P = 0.035; HR: 3.29; 95% CI: 1.09–9.96), whereas standard predictors of prognosis were not found to be significant (Supplementary Table ST5).

Next, we extended the spatial analysis of qmIF by using PCF. We found that CD8⁺ T cells clustered more strongly with tumor cells in the nonrecurrent samples compared with recurrent samples (P = 0.0086, Mann–Whitney U test, Fig. 4A, left). This effect was similar when comparing the clustering of tumor cells with HLA-DR⁺ macrophages (P = 0.016, Fig. 4A, right), but not with HLA-DR⁻ macrophages (data not shown). We provide a representative depiction of the lowest and highest PCF values for each of the tumor–CTL (Fig. 4B, left) and tumor–Macrophage (Fig. 4B, right) interactions.

We next sought to compare the AUCs of these PCFs against gene expression, finding significant positive correlations between the tumor–CTL clustering with CXCL11 (P = 0.016, Spearman correlation, Bonferroni correction) and MFGE8 (P = 0.036), and between the tumor–macrophage clustering with TAP2 (P = 0.035; Supplementary Fig. S2).

Having previously found that high density of HLA-DR⁻ macrophages predicts poor DSS (19), we next investigated the properties of these macrophage subpopulations. We performed additional staining for CD68, HLA-DR, CD163, CSF1R, MHC1, and CD33 over a total of 5519 CD68⁺ macrophages identified across the 13 patients in the macrophage survey subcohort (Fig. 5A). We visualized the concatenation of each of these markers using pie charts comparing HLA-DR⁻ and HLA-DR⁺ macrophages (Fig. 5B). We next visualized this data using tSNE to compare the expression of these markers on HLA-DR⁻ and HLA-DR⁺ macrophages. HLA-DR expression was largely independent of CD163 (P = 0.797, Fig. 5C). However, HLA-DR expression was correlated with MHC1 expression (P = 1.4e-102; OR, 3.4; χ² test), whereas lack of HLA-DR expression was correlated with CD33 expression (P = 1.7e-53; OR, 0.22; Fig. 5C). To further study the cells expressing CSF1R and CD33, we evaluated the concatenation of CSF1R and CD33 with CD163, finding that HLA-DR⁻ cells were more likely to also be CD163⁺CSF1R⁺ compared with HLA-DR⁺ cells (P < 0.001, Fig. 5C). Furthermore, in comparison with HLA-DR⁺ cells, HLA-DR⁻ cells were also more frequently CD163⁺CD33⁺ (P < 0.001, Fig. 5C).

In addition to qmIF analysis, we further evaluated 8 primary melanoma cases using DSP. We found that the ROI from patients with higher stromal density of HLA-DR⁻ macrophages had a lower immune infiltration overall as assessed by quantitation of CD45 per ROI (P < 0.0001, Fig. 5D). In addition, CD3 per ROI was higher in patients with lower stromal density of HLA-DR⁻ macrophages (P < 0.0001, Fig. 5D). Interestingly, patients with higher density of HLA-DR⁻ macrophages also had a significantly higher ratio of CD4 to CD45 (P < 0.0001), but a similar CD8A to CD45 ratio (Supplementary Fig. S3). PD-L1 and PD1 levels per CD45 were significantly higher in patients with higher HLA-DR⁻ macrophage density (P < 0.0001 and P = 0.0002, respectively, Fig. 5D).

Innovations in digital pathology methods allow for quantitation of cell phenotypes and reproducible characterization of morphologic and genomic features, and provide new tools for clinical pathology (18, 19, 26, 27). These analyses enhance personalized medicine and allow for rapid development and implementation of new digital pathology-based biomarkers. Here, we show that nuclear area and cellular clustering correlate with clinical features and propose a biomarker combining transcriptomic and qmIF data.

Current pathology methods such as TIL evaluation and tumor grading correlate with clinical outcome in melanoma, but are generally qualitative and therefore difficult to standardize (2, 13, 14). In previous work, we quantitatively evaluated density and spatial relationships of immune infiltrates in primary melanoma specimens using qmIF, finding that the ratio of CD8 to CD68 can be used as a biomarker for DSS and OS in stage II-III primary melanoma (19). In this work, we find that the CD8 to CD68 ratio also predicts recurrence in primary early-stage melanoma. Furthermore, we evaluate the same qmIF panel in metastatic cases and find that the CD8 to CD68 ratio in metastatic lesions (M) is similar to that in primary melanomas that distantly recur (R). Both M and R have a lower CD8 to CD68 ratio than nonrecurrent primary melanomas (NR). This suggests that the TIME in primary melanomas with high metastatic potential is similar to that in metastatic lesions. These findings build on prior work showing that the ratio of CD8 to CD68 has potential as a favorable prognostic feature in primary melanoma tumors (19). Phenotyping both lymphocytes and myeloid cells within the TIME can improve upon quantitation of lymphocyte populations alone (28).

More precise phenotyping of immune cells and correlation of these findings with clinical features will likely lead to the development of more accurate biomarkers. Myeloid cells that are CD68⁺ are predominately macrophages and have been shown to be an unfavorable feature in several tumor types (29). Macrophage populations, however, are notoriously diverse and can exert both protumorigenic and proimmune functions (30, 31). Previously, we distinguished HLA-DR⁻ macrophages from HLA-DR⁺ macrophages as indicative of a poor prognosis (19). HLA-DR expression indicates an activated state and is broadly associated with an M1, proimmune, phenotype (32, 33). Here, we further characterize the HLA-DR⁻ population as containing suppressive populations expressing CD163 and CSF1R. These findings suggest that CD68⁺HLA-DR⁻ cells are suppressive, consistent with our hypothesis that these cells exert a protumorigenic effect. As a future research direction, we hope to further explore the three-cell interaction between these macrophages, T cells, and tumor cells.

QmIF allows not only for calculation of the density of different immune cell phenotypes but also analysis of cell size and distances between cells (34–36). First, we compare nuclear size in NR, R, and M lesions. Interestingly, we find that nuclei in metastatic lesions are larger than in primary lesions regardless of expression of the proliferation marker, Ki67. However, the size of the nuclei in primary lesions is similar regardless of recurrence status. From a clinical pathology perspective, findings suggest that in the absence of an epidermal component, increased nuclear size may favor a melanoma metastasis. However, larger studies are needed.

Distance between cells can be used to define spatial relationships of immune cells with tumor cells (37). Previously, we found that simple nearest neighbor distance of CD8⁺ T cells to CD68⁺HLA-DR⁻ macrophages predicted poor prognosis (19). Here, we further evaluate T cells and macrophages using clustering algorithms and find that even after accounting for variations in cell density there are patterns of positioning among immune and tumor cells that are related to prognosis. Specifically, clustering of tumor cells with CD8⁺ T cells correlates with nonrecurrence. In addition, although our prior study found that density of HLA-DR⁺ macrophages alone did not predict prognosis, we find that clustering of HLA-DR⁺ macrophages with tumor cells correlates with nonrecurrence. The importance of spatial proximity of immune cells to each other and to tumor cells implies a mechanism by which their infiltration affects disease progression. These findings are consistent with the hypothesis that immune surveillance is operative in primary melanoma tumors and that proximity of tumor cells to immune cells of defined phenotypes impacts melanoma progression. Spatial location of cells may be impacted by local hypoxia or other features and future studies will be useful to further explore metabolic mechanisms underlying spatial distributions of immune cells within melanoma and other cancers (38, 39).

In addition to our qmIF analysis, we also previously defined and validated a transcriptomic 53 immune gene panel, MIP, predictive of recurrence in the same cohort of patients with primary stage II-III melanoma (17, 18). Here, we evaluate transcriptomic profiling in addition to qmIF immune phenotyping and find consistency between the defined cell density and immune gene signatures. Interestingly, we find that density of CD8⁺ T cells and HLA-DR⁺ macrophages correlates with defined CIBERSORT signatures (25) for CD8 T cells and M1 macrophages, respectively.

Furthermore, combining the CD8 to CD68 ratio with MIP allows for development of a composite biomarker that quantitatively assesses the entire specimen. This composite biomarker identifies three separate risk groups, most importantly a high-risk population that may benefit from adjuvant studies of combination immunotherapies. Patient stratification is urgently needed in early-stage melanoma because current protocols support treatment of all stage III patients despite the fact that 5-year survival is 83% in stage IIIB patients and 93% in stage IIIA patients. Meanwhile, survival rates in stage II populations are also generally favorable, ranging from 82% to 94%. Nonetheless melanoma, when it recurs is devastating. Therefore, identifying high-risk groups would facilitate application of combination immunotherapy in the adjuvant setting by allowing stratification for clinical trials. In addition, there are data to suggest that a favorable tumor immune microenvironment predisposes to response to PD-1 blockade (40). Therefore, patients at intermediate risk might be more likely to respond to therapy both in the adjuvant setting and at recurrence while low-risk patients could avoid immunotherapy entirely. Thus, our composite biomarker findings, while exploratory and based on small populations, suggest that further work to quantitatively characterize and correlate phenotypic and genomic features of the TIME using FFPE tissues is likely to improve personalized care of patients with early-stage melanoma.

R.D. Gartrell-Corrado was a paid consultant for Northwest Biotherapeutics and reports receiving a commercial research grant from NanoString. E.M. Rizk reports receiving travel support from NanoString. K.M. Komatsubara is an employee/paid consultant for Roche/Genentech. R.D. Carvajal is an employee/paid consultant for Immunocore, I-Mab, InxMed, Merck, Pierre Fabre, PureTech, Sanofi, and Sorrento. R. Rabadan is an unpaid consultant for AmiedBio and CNIO. Y.M. Saenger reports receiving a commercial research grant from Amgen and other commercial research support from Regeneron and has an immediate family member who has ownership interest (including patents) with Wasaba. No potential conflicts of interest were disclosed by the other authors.

The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH. The funding sources had no role in preparation of the manuscript or the decision to submit for publication.

Conception and design: R.D. Gartrell-Corrado, A.X. Chen, Y.M. Saenger

Development of methodology: R.D. Gartrell-Corrado, A.X. Chen, A. Monod, R. Rabadan, Y.M. Saenger

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): R.D. Gartrell-Corrado, M.H. Bogardus, T.D. Hart, A.M. Silverman, C.-A.Y. Bayan, L.W. Barker, K.M. Komatsubara, R.D. Carvajal, Y.M. Saenger

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): R.D. Gartrell-Corrado, A.X. Chen, E.M. Rizk, D.K. Marks, T.D. Hart, A.M. Silverman, K.M. Komatsubara, R. Chang, Y.M. Saenger

Writing, review, and/or revision of the manuscript: R.D. Gartrell-Corrado, A.X. Chen, E.M. Rizk, D.K. Marks, M.H. Bogardus, A.M. Silverman, C.-A.Y. Bayan, G.G. Finkel, L.W. Barker, R.D. Carvajal, B.A. Horst, Y.M. Saenger

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): R.D. Gartrell-Corrado, E.M. Rizk, D.K. Marks, M.H. Bogardus, C.-A.Y. Bayan, B.A. Horst, Y.M. Saenger

Study supervision: A. Monod, R. Rabadan, Y.M. Saenger

This study was supported by the NIH through grant numbers UH2CA21814901A1 (to Y.M. Saenger) and KL2TR001874 (to R.D. Gartrell-Corrado). A.X. Chen is funded by the Medical Scientist Training Program (T32GM007367). Yvonne Saenger is also supported by funding from the Melanoma Research Alliance and by an Irving Assistant Professorship at Columbia University's NIH/NCATS CTSA Program hub: UL1TR001873. Robyn Gartrell-Corrado is also supported by Swim Across America. R. Rabadán, A. Monod, and A. Chen have been supported by the NCI Center for Topology of Cancer Evolution and Heterogeneity (U54CA193313). A. Monod wishes to acknowledge the support of the New Frontiers in Research Fund–Fonds Nouvelles Frontières en Recherche (SSHRC-NFRF-FNFR, Government of Canada) NFRE-2018-00431.

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.