Proximity of immune and tumor cells underlies response to BRAF/MEK-targeted therapies in metastatic melanoma patients

Acquired resistance to BRAF/MEK-targeted therapy occurs in the majority of melanoma patients that harbor BRAF mutated tumors, leading to relapse or progression and the underlying mechanism is unclear in many cases. Using multiplex immunohistochemistry and spatial imaging analysis of paired tumor sections obtained from 11 melanoma patients prior to BRAF/MEK-targeted therapy and when the disease progressed on therapy, we observed a significant increase of tumor cellularity in the progressed tumors and the close association of SOX10+ melanoma cells with CD8+ T cells negatively correlated with patient’s progression-free survival (PFS). In the TCGA-melanoma dataset (n = 445), tumor cellularity exhibited additive prognostic value in the immune score signature to predict overall survival in patients with early-stage melanoma. Moreover, tumor cellularity prognoses OS independent of immune score in patients with late-stage melanoma.

for immune cell activation 5,6 ), or association with clinical parameters pre-and post-treatment (Fig. 1). We observed a significant increase in SOX10 + melanoma cells in the progressed tumors compared to the tumors prior to BRAF/MEK-inhibition from the same patient (Fig. 2a). While the frequency of preexisting SOX10 + melanoma cells prior to treatment was not associated with PFS, a greater increase in SOX10 + melanoma cells from baseline to progression was associated with shorter PFS under BRAF/MEK inhibition treatment (Spearman r = −0.842, p = 0.004) (Fig. 2b). We speculated that melanoma cells may suppress immune cell responses for immune escape and disease progression. However, there was no alteration in the number of CD8 + T cells or CD11c + DCs in the progressed tumors compared to the tumors prior to therapy ( Supplementary Fig. 1). Neither the preexisting frequency nor change of frequency of CD8 + T cells or CD11c + DCs was associated with PFS ( Supplementary Fig. 1). Intercellular sensing and communication (e.g., soluble cytokines/ chemokines) requires close proximity (~40 μm) of interacting cells 7 . We next counted the SOX10 + melanoma cells that paired with CD8 + T cells and normalized the pairs to the total CD8 + T-cell number since CD8 + T-cell number was not altered due to therapy. The data show that those patients with a high number of SOX10 + melanoma cells at progression and a higher number of SOX10 & CD8 + T-cell pairs had a significantly shorter PFS (Spearman r = −0.75, p = 0.025). This effect was most pronounced for patients where melanoma cells were paired within 45 μm of CD8 + T cells (Spearman r = −0.783, p = 0.017) (Fig. 2c). No significant differences were observed in CD40 and CD80 content. There is an extensive data overlap between patients with and without prior immunotherapy or immunotherapy-chemotherapy combination before the BRAF/MEK-therapy, suggesting our observation is likely independent of prior chemo-/immuno-therapy status (Supplementary Fig. 2). Given the uneven distribution of SOX10 + cells in the tumor sections, visual estimates of percentage tumor nuclei in the region-of-interests (ROIs) were determined by pathologists who were blinded to treatment ( Supplementary Fig. 3). These ROIs (1-3 per tumor) which exhibited viable tumor cells and reliable SOX10 expression were selected by the pathologists. Consistent with whole-slide SOX10 quantitation, the change of %SOX10 + cells in ROIs confirmed a negative correlation with patient's PFS (Spearman r = −0.733, p = 0.0311). We also observed a moderate positive correlation between the visually estimated percentage tumor nuclei (generally >60%) and the computer-assisted SOX10 + cell density quantitation (generally 10-60%). Nevertheless, there was no alteration of visually estimated percentage tumor nuclei in the progressed tumors compared to the tumors prior to BRAF/ MEK-inhibition. As such, percentage tumor nuclei does not correlate to the patient's PFS. Visual estimation alone is known to be variable and inaccurate 8 . Our whole-slide MxIHC-based SOX10% results determined by computer-assisted cell density quantitation are in line with the observations of AstroPath multispectral imaging platform and may highlight the benefit of computer-assisted cell density quantitation to capture treatmentinduced responses and antibody cocktail staining of melanoma cell markers to reduce false-negative signals 9 .
Immune score based on the ESTIMATE method was previously generated from 11 different tumor types 10 and confirmed in melanoma 11 to infer tumor purity, stromal and immune cell admixture from expression data in The Cancer Genome Atlas (TCGA). We further showed that OS was associated with immune score stratified by disease stage in melanoma (n = 445) (Fig. 3a). In contrast, increased SOX10 z-scores were associated with decreased OS of patients with late-stage (n = 192; hazard ratio [HR] = 1.186; 95% confidence interval [CI] = 1.006-1.397; p = 0.03), but not early-stage melanoma (n = 217; HR = 1.203; 95% CI = 1.008-1.436; p = 0.192) (Fig. 3b). Notably, compared to the main effect of immune score only, the multivariable analysis presented an additive interaction of SOX10 with an immune score to predict OS in patients with early-stage (likelihood ratio [LR] test, p = 0.005), but not late-stage melanoma (LR test, p = 0.554) (Fig.  3c, d). Together, our data revealed significant interaction and an additive effect of the prognostic value of SOX10, confirmed by another melanoma marker MLANA ( Supplementary Fig. 4), with an immune score in early-but not late-stage melanoma. Thus, for early-stage patients with a low immune score, the OS of melanoma patients with tumors of high SOX10 (Wald test, p = 0.0012) or MLANA (Wald test, p = 0.0051) was worse than that of those with low levels of SOX10 or MLANA. For late-stage patients, increased SOX10 (log-rank test, p = 0.03) or MLANA (log-rank test, p = 0.018) prognoses worse OS independent of immune score.
In this retrospective study, we performed MxIHC and subsequent whole-tumor imaging spatial analyses of tumors obtained from 11 melanoma patients prior to BRAF/MEK-targeted therapy and after the disease progressed. The data revealed that at the time of progression, the close association of SOX10 + melanoma cells with CD8 + T cells negatively correlated with PFS of melanoma patients and may indicate a potential mechanism for acquired resistance to BRAF/MEK-targeted therapy. The observed increase of tumor cellularity in progressed melanoma tumors could be either a result or a cause of the escape from BRAF/MEKinhibition. Additional on-treatment biopsies, as well as assessment of the genetic profile of the tumor cells that expanded posttreatment, would be essential to further explore whether there was an expansion of a treatment-resistant clone that led to poor PFS or was it the density of the tumor cells in relation to the immune effector cells per se that was associated with poorer PFS. Since SOX10 was associated with OS of late-stage melanoma patients in the TCGA dataset, it is plausible that the SOX10 expression level is a prognostic marker in melanoma, and is not linked to any specific treatment response. Using a preclinical mouse model, we showed the detrimental role of PD-L1:PD-1 axis between melanoma cell:Tc-cell in melanoma. Besides cell-cell interaction mechanisms, tumor cells may secrete immunosuppressive cytokines/chemokines and/or compete with immune cells in the microenvironment for the components required for their own metabolism, further inhibiting immune cell functions 16 . Indeed, CD8 + T cells in primary melanoma were spatially distant from proliferating (Ki67 + ) tumor cells compared to nondividing (Ki67 − ) tumor cells, suggesting that rapidly growing primary tumors may suppress and/or exclude CD8 + T cells or fail to produce factors that recruit these cells into the tumor 17 .
We acknowledge that there are several limitations to our study. First, this study was conducted at a single center with a small  sample size, which may result in institution-specific biases. Second, the clinical data were assessed retrospectively and not in a controlled, prospective fashion. Finally, the melanoma patients received BRAF inhibitor with or without MEK inhibitor. However, we observed a consensus association between outcomes and tumor cellularity irrespective of therapy type. The prognostic value of tumor cellularity was confirmed in the TCGAmelanoma dataset (n = 445). Further research is needed to understand the crosstalk between rapidly proliferating melanoma cells and T cells to decipher the biological mechanisms underlying the acquired resistance of BRAF/MEK-targeted therapy.

Patient material
Institutional IRB approval and written informed consent from all patients were obtained before study initiation. All patient donors signed informed consent before providing tissue samples. Patient samples were collected on a tissue-collection protocol approved by the Vanderbilt University IRB. Paired advanced melanoma samples from 11 patients (pre-and posttreatment of BRAF/MEK-targeted therapy) were collected as part of NCT01205815 clinical trial 5 .
Multiplex immunohistochemistry (MxIHC) assessment and spatial analysis  18 . Computational single-cell segmentation/identification was performed on the Ilastik probability maps using CellProfiler 19 and MatLab. Pairwise cell analysis was performed by calculating the Euclidian distance between identified single-cell centroids.

Immune score estimation
Immune scores were derived from RNA-seq data in TCGA Skin Cutaneous Melanoma based on the ESTIMATE algorithm, which were previously generated from 11 different tumor types (not including melanoma) via an R library of estimate package that is available on https://bioinformatics. mdanderson.org/estimate/rpackage.html 10 .

Flow cytometric analysis
The details of staining and flow cytometry analysis protocols is according to our previously published methodology 5

Statistical analysis
For BRAF/MEK-targeted therapy clinical data, treatment effects in standard two-group paired experiments were compared using the Wilcoxon matched-pairs signed-rank test. A Spearman's rank correlation test was used to evaluate the association between PFS and each variable. A scatter plot was shown with a quadratic regression and its 95% confidence interval for visualization. For univariable analysis of the TCGA Skin Cutaneous Melanoma data, survival curves were estimated using the Kaplan-Meier method and compared between groups including the following variables: gender, stage (I + II vs. III + IV), immune score (low vs. high), SOX10 (low vs. high), MLANA (low vs. high), and CD8 (low vs. high) with the log-rank test. The log-rank test compares the entire survival experience between groups and to see whether the survival curves are identical or not. Please note that median split was used for tuning a continuous variable into a binary variable to visualize the difference between groups. In the survival (Kaplan-Meier) plot, the hazard ratio (HR) of death with 95% CI was reported per interquartile range (IQR) change in continuous predictors based on univariable Cox regression. Taking the survival plot by SOX10 in Fig. 1 as an example, for the unit of one interquartile change of the SOX10 z-score, the risk of death increased by 32.1% (HR = 1.321). If the HR is less than 1, for example, 0.454 (like survival plot by the immune score at the early stage in Fig. 3a), the risk of death falls by 54.6% per IRQ change in immune score. When comparing two groups (such as gender), the HR was used to show which groups are more likely to experience an event (death) first. Wilcoxon rank-sum test was used to test the difference in gene expression profiles among four subgroups including SOX10 hi CD8 hi , SOX10 hi CD8 low , SOX10 low CD8 hi , and SOX10 low CD8 low . P value was adjusted with Bonferroni correction for pairwise comparisons of four subgroups. Redundancy analysis was used to determine how well each variable [immune score, sex, disease stage (I + II vs III + IV), etc.] could be explained by other predicted from the remaining variables. There were no redundant variables detected. Multiple imputation analysis was performed to account for missing data with ten repetitions. It was predetermined that variables with more than 20% missing values were excluded from the analysis (such as BRAF mutations, NRAS mutation, height, and weight). The race was not included in this report because 97% were Caucasian (see in Supplementary Table 2). Variable selection was conducted using regularization methods (elastic-net penalty). (1) Multivariable Cox (proportional hazards) regression was performed for investigating the association between overall survival and predictors and was used to estimate the adjusted hazard ratios. With a continuous predictor, the HR indicates the change in the risk of death per IQR change. Diagnostic tests revealed significant non-proportionality in the disease stage. To ameliorate this violation of modeling assumptions, Cox regression models were stratified by early stage and late stage. (2) The likelihood ratio (LR) test determines the goodness-of-fit between the full model which included interactions between immune score and each marker and the reduced model which contained main effects. A statistically significant interaction implies that the effect of score changes with the level of marker (and the effect of marker differs by the level of score). (3) The concordance index (c index) was used to measure of predictive accuracy and to assess the added value of a new marker (SOX10 and MLANA) where c index is the probability of concordance between the predicted and the observed survival. The c index with bootstrap approach was used to adjust for optimism/overfitting in measures of predictive ability for internal validation. For evaluating the association between OS and an interaction effect of SOX10 and CD8, multivariable cox proportional hazards regression stratified by the disease stage was performed where the SOX10 z-score and log2 of CD8 were analyzed to reduce skewness from the predictors. Moreover, multivariable linear regression was used to evaluate the interaction effect of SOX10 z-score and log2 of CD8 on each log2 of gene expressions.
For mouse studies, the progression of tumor volume (mm 3 ) over time among groups of mice with or without therapies were compared with a linear mixed-effects regression model to take into account the correlation structure with the repeated measures data within a mouse. A square root or a natural log transformation was implemented to better meet the normality assumptions. The likelihood ratio test was performed to identify statistically significant time by treatment effect. A statistically significant interaction implies that the magnitude of treatment differences depends on the actual day of measurement. Using model-based (least-square) means, the average tumor growth between treatment groups was estimated and compared with the Wald test. R version 4.0.4 was used for statistical analysis.

Reporting Summary
Further information on research design is available in the Nature Research Reporting Summary linked to this article.

DATA AVAILABILITY
The TCGA Skin Cutaneous Melanoma (Firehose Legacy) data that support the findings of this study are available from the website [cBioPortal for Cancer Genomics] (https://www.cbioportal.org/study/summary?id=skcm_tcga). De-identified BRAF/ MEK-targeted therapy clinical data are available on request due to privacy restrictions. Data of patient characteristics and response to BRAF-targeted therapy have been presented in aggregate form in Supplementary Table 1. All the data that support the findings of this study are available on request from the corresponding author A.R. (ann.richmond@vanderbilt.edu). The data were not publicly available due to potential compromise of research participant privacy. No data use agreement is required for aggregate data, though access to individual patient-level data requires a