Intratumoral Heterogeneity in Ngfr+ Melanoma Subpopulations Shapes Immune Evasion and Immunotherapy Resistance

Intratumoral Heterogeneity in Ngfr+ Melanoma Subpopulations Shapes Immune Evasion and Immunotherapy Resistance | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Article Intratumoral Heterogeneity in Ngfr + Melanoma Subpopulations Shapes Immune Evasion and Immunotherapy Resistance Daniel Peeper, Sebastiaan Schieven, Joleen Traets, Arno Velds, and 14 more This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-6506453/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Human melanomas dedifferentiate into a neural crest-like cell state when exposed to T cell cytokines or MAPK pathway inhibitors. This transformation is associated with cellular heterogeneity and the emergence of small, therapy-resistant melanoma populations characterized by elevated nerve growth factor receptor (NGFR) expression. However, the extent of this heterogeneity and its impact on immunotherapy response remain unclear. By dissecting intratumor heterogeneity in patient melanomas, we show here that even within NGFR + tumor subpopulations, remarkable phenotypic and functional diversity exists. Combined single-cell RNA sequencing (scRNA-seq) of NGFR + fractions and spatial transcriptomics uncovered pronounced diversity among single-cell clusters, characterized by patient-specific gene regulatory networks (GRNs) and distinct spatial organization. Furthermore, we identify an NGFR + subpopulation marked by co-expression of Platelet-Derived Growth Factor Receptor (PDGFR), which is associated with increased resistance to the T cell cytokines IFNg and TNF. Clinically corroborating these findings, we observed that NGFR + /PDGFR + mesenchymal-like cells are enriched in melanomas infiltrated with active T cells yet failing to respond to immune checkpoint blockade treatment. Our results highlight extreme heterogeneity within human melanoma, which is spatially organized and regulated by patient-specific GRNs, and harboring a distinct subfraction linked to immunotherapy resistance. Biological sciences/Cancer/Skin cancer/Melanoma Biological sciences/Computational biology and bioinformatics/Gene regulatory networks Biological sciences/Immunology/Immune evasion NGFR heterogeneity spatial transcriptomics PDGFR cytokine resistance and immune checkpoint blockade resistance Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 INTRODUCTION Melanoma is a heterogeneous cancer type, in part owing to both its high mutational load and the presence of multiple transcriptional phenotypes 1 – 3 . Originally, melanomas were classified as either proliferative or invasive 4 , 5 , mostly as a function of the expression levels of the key melanocyte transcription factor MITF and its target genes. Later, this model was refined to include a higher phenotypic complexity at the transcriptomic level 6 – 8 . The different melanoma phenotypes can transition from one to the other, for example upon inhibition of the MAPK pathway 8 . We and others have shown that such transitions are associated with cell state-specific gene signatures and marker genes, for example AXL and NGFR 7 , 9 – 11 . NGFR, a cell surface neurotrophic factor receptor marking neural crest stem cells (NCSCs) 7 , 12 , 13 stands out, as it can be co-expressed with either differentiation or dedifferentiation markers in melanoma cells 11 , 12 , 14 , 15 . Its expression can be induced through various mechanisms, including MAPK pathway inhibition 8 , 16 , 17 , TGF-β 18 , in vitro T cell challenge 19 , adoptive cell transfer (ACT) 14 , 15 , immunotherapy 20 and immune cytokines such as TNF 14 , 15 , 21 and IFNg 22 , 23 . This complex regulation of NGFR conceivably contributes to its heterogeneous expression pattern. We previously showed that NGFR + subpopulations commonly pre-exist in patient melanomas already before any treatment, and that they typically manifest in heterogeneous patterns 19 . Furthermore, we demonstrated that NGFR marks melanoma cells that are relatively insensitive to cytotoxic T cells, as well as their cytokines IFNg and TNF, and can be associated with T cell exclusion 19 . In line with these observations, others have revealed NGFR's involvement in rendering melanoma cells resistant to natural killer (NK) cells 24 . Additionally, in a patient case study, the proportion of NGFR-expressing cells increased following immune checkpoint blockade (ICB) 20 . Given that a substantial fraction of melanoma patients exhibits either intrinsic or acquired resistance to immunotherapy, potentially attributed to NGFR expression, we suggested previously that it may be beneficial to target NGFR to overcome resistance 19 . These observations together suggest that the NGFR + melanoma subpopulation can manifest in heterogeneous patterns, express different marker proteins, and is induced by distinct stimuli. However, it is unclear whether it reflects an intrinsically homogeneous cell group or that even this subpopulation is heterogeneous in nature. Clearly, this would have immediate consequences for responses to immune cells and therapy. Therefore, we dissected NGFR heterogeneity from clinical melanoma specimens using different technologies including single-cell RNA-sequencing (scRNA-seq) and spatial transcriptomics, in combination with functional analyses. Specifically, we investigated whether (i) this heterogeneity manifests only between cell groups or that even single cells can co-express different markers; (ii) this pattern of heterogeneity is consistent among patients; (iii) heterogeneous NGFR + cells are spatially organized; (iv) heterogeneous NGFR + subpopulations show differential immune responses; (v) these responses are reflected by differential clinical responses to immunotherapy. MATERIAL AND METHODS Cell lines and cell culture conditions All melanoma cell lines are from the Peeper laboratory cell line stock. Melanoma cell lines were cultured in DMEM (Gibco) with fetal bovine serum (FBS, Sigma), 100 U/ml penicillin and 0.1 mg/ml streptomycin (both Gibco) under standard conditions. All cell lines were authenticated by short tandem repeat (STR) profiling (Promega) and regularly confirmed to be mycoplasma-free by PCR. Western blotting and antibodies Western blotting and antibodies As described previously 25 . Cell pellets were lysed in RIPA buffer (50 mM TRIS pH 8, 150 mM NaCl, 1% Nonidet P40, 0.5% sodium deoxycholate, 0.1% SDS) supplemented with HALT™ protease and phosphatase inhibitor cocktail (100x) (Fisher Scientific, cat # 78444). Lysis was performed for 30 minutes and vortexed every 10 minutes. Samples were centrifuged at maximum speed for 10 minutes. Protein concentration was determined using a Bradford assay (Bio-Rad). Protein concentrations were normalized to each other and mixed with 4x LDS sample buffer (Fisher Scientific, 15484379) containing 10% β-Mercaptoethanol (Merck) (final concentration 2.5%) and incubated for 5 minutes at 95°C. Western blotting was performed with standard techniques using 4–12% Bis-Tris polyacrylamide-SDS gels (NuPAGE, Life Technologies) and nitrocellulose membranes (Whatman, GE Healthcare). Blotting was performed using the iBlot dry blotting system from Invitrogen. Blots were blocked in 4% milk in PBS plus 0.2% Tween 100 and incubated with primary antibody: AXL (1:1,000, C89E7, Cell Signaling), Melan-A (1:10,000, M7196, Dako); PDGFRb (1:1,000, 3169, CST), PDGFRa (1:1,000, D1E1E, Cell Signaling), Tubulin (1:10,000, T9026, Sigma), and NGFR (1:1,000, #8238, Cell Signaling). The following secondary antibodies were used: goat anti-rabbit peroxidase conjugate (1:7,500, G21234) and goat anti-mouse (1:7,500, G21040), both purchased from Invitrogen. Immunoblots were incubated with Clarity™ Western ECL Substrate (cat# 170–5061, Biorad). Luminescence was captured by the Bio-Rad ChemiDoc imaging system. Both 8-bit tiff and 16-bit raw tiff images were used to make the figures. IHC of human melanomas As described previously 25 . The collection and use of human tissue was approved by the Medical Ethical Review Board of the Antoni van Leeuwenhoek. Patients gave informed consent for secondary use of tumor tissue. The study was approved by the Institutional Review Board (IRB). Immunohistochemistry of the FFPE tumor samples was performed on a BenchMark Ultra autostainer (Ventana Medical Systems). Briefly, paraffin sections were cut at 3 µm, heated at 75°C for 28 minutes and deparaffinized in the instrument with EZ prep solution (Ventana Medical Systems). Heat-induced antigen retrieval was carried out using Cell Conditioning 1 (CC1, Ventana Medical Systems) for 32 minutes at 95°C (MITF) and 40 minutes at 95°C (AXL). MITF was detected using clone C5/D5 (1:800 dilution, 32 minutes at 37°C., LSBio), AXL using clone C89E7 (1:100 dilution, 32 minutes at room temperature., Cell Signaling). For AXL, signal amplification was applied using the Optiview Amplification Kit (8 minutes, Ventana Medical Systems). Melan-A was detected using clone A103 (1:40 dilution, 32 minutes at 37°C, Agilent/Dako). NGFR was detected using clone D4B3 (1:400 dilution, 1 hour at room temperature, Cell signaling). For Melan-A, signal amplification was applied using the Optiview Amplification Kit (4 minutes, Ventana Medical Systems). CD3 was stained as described previously (RM-9107-S, Thermo Scientific) 19 . PDGFRb was detected using clone 28E1 (1:50 dilution, 60 minutes at room temperature, Cell Signaling). Stainings were developed using brown (OptiView DAB Detection Kit (Ventana Medical Systems) or red (UltraView Universal Alkaline Phosphatase Red Detection (Roche Diagnostics, Ventana)) visualization. Slides were counterstained with Hematoxylin and Bluing Reagent (Ventana Medical Systems). For a sequential Periodic Acid–Schiff (PAS) staining, slides were removed from the BenchMark Ultra autostainer, rinsed in distilled water, incubated in 0.5% Periodic Acid (7 minutes, VWR) followed by Schiff’s Reagent (30 minutes, CellaVision / RAL Diagnostics) after washing steps in between. Slides were counterstained with Hematoxylin (KliniPath). A PANNORAMIC® 1000 scanner from 3DHISTECH was used to scan the slides at a 40x magnification. Percentage positive tumor cells was quantified. Cytoplasmic and membrane staining were scored for AXL, NGFR, PDGFRβ and Melan-A. Only nuclear MITF staining was scored. Scoring was performed by a certified pathologist. A semi-quantitative scoring was applied to quantify CD3 infiltration. Tumor regions were scored as: not infiltrated (score of 0), lowly infiltrated (a score of 2), moderately infiltrated (a score of 4) and highly infiltrated (a score of 6). The average of all scores was taken if tumors received more than one score. Only intra-tumoral CD3 + T cells were scored. The majority of the tumors studied were untreated/baseline and obtained from surgical resections or using a 14-gauge biopsy needle. One tumor was an anti-PD1 treatment relapse sample (patient tumor 44 in Fig. 1 A and S1 A). Flow cytometry Method for cell surface staining as described previously 25 . Cells were stained with antibodies targeting surface molecules of interest according to manufacturer’s instructions and analyzed on a Fortessa flow cytometer LSR ( BD Bioscience). The following antibodies were used: AXL-PE conjugated antibody (1:200, FAB154P, R&D), PDGFRb-PE conjugated antibody (1:50, FAB1263P, R&D) and NGFR-APC (1:200, 345107, Biolegend) for 20 minutes at 4°C. Data depicted is mean fluorescence intensity (MFI) of sample - MFI of unstained sample, unless otherwise stated. Colony formation assay 100,000 cells per well were seeded in a 12 wells plate. The day after, cells were challenged with cytokines IFNg (300-02, Prepotech) and TNF (11343017, Immunotools). After 5 days, medium was removed and cells were stained with a crystal violet solution containing 0.1% crystal violet (Sigma) and 50% methanol (Honeywell) for 1 hour. For quantification, the crystal violet stain was dissolved in 10% acetic acid (Sigma). Absorbance of this solution was measured on an Infinite 200 Pro spectrophotometer (Tecan) at 595 nm. Single-cell RNA-sequencing sample preparation Collected tumors from patients 35 (lymph node lesion), 37 (lymph node lesion), 44 (thigh lesion within fat tissue) and 45 (skin lesion, Fig. 1 A) were digested in FBS free RPMI medium (Gibco) with 100 U/ml penicillin/0.1 mg/ml streptomycin, collagenase IV (1:50, 17104-019, ThermoFisher Scientific) and pulmozyme (1 mg/ml, Roche) for 30 minutes at 37°C on a rotator. Digest was filtered over a 100 µm filter and frozen down in 90% FBS and 10% DMSO (45-34943, Sigma-Aldrich). Next, digests were thawed in 3 ml cold RPMI medium with 10% FBS and 4 µl benzonase (1:1,000, 70746-3, VWR). Cells were spun for 8 minutes 14,000 rpm at 4°C. Supernatant was removed, and cells were resuspended in 4 ml cold RPMI with 4 µl benzonase (1:1,000). Samples were centrifuged for 8 minutes 14,000 rpm at 4°C followed by resuspension in cold PBS with 1% BSA. The following antibodies were used for staining: NGFR-APC (1:200, 345107, Biolegend) and CD45-Alexa Fluor 488 (1:400, 304019, Biolegend). Cells were incubated in staining solution (i.e., antibodies diluted in PBS + 1% BSA) for 20–30 minutes. Then, cells were washed twice with cold PBS with 1% BSA and centrifuged (8 minutes 14,000 rpm at 4°C). Cells were kept in PBS with 2% of FBS before sorting for NGFR + and CD45 − cells. Cells were sorted by BD FACSAria Fusion. DAPI was used as a life dead marker. Cells were considered NGFR positive if they had a NGFR signal higher than unstained and/or fluorescence minus one (FMO) control. CD45 − /NGFR + cells were used as input for scRNA-seq. Single-cell RNA-sequencing For each sample, the Chromium Controller platform of 10X Genomics was used for single-cell partitioning and barcoding. Each cell’s transcriptome was barcoded during reverse transcription, pooled cDNA was amplified and Single-Cell 3’ Gene Expression were prepared according to the manufacturer’s protocol “Chromium NextGEM Single-Cell 3’ Reagent Kits v3.1” (CG000315, 10X Genomics). Both Single-Cell 3’ Gene Expression libraries were quantified on a 2100 Bioanalyzer Instrument following the manufacturer’s protocol “Agilent DNA 7500 kit” (G2938-90024, Agilent Technologies). These Single-Cell 3’ Gene Expression libraries were combined to create one sequence library pool which was quantified by qPCR, according to manufacturer’s protocol “KAPA Library Quantification Kit Illumina® Platforms” (KR0405, KAPA Biosystems). A NovaSeq 6000 Illumina sequencing system was used for paired end sequencing of the Single-Cell 3’ Gene Expression libraries at a sequencing depth of approximately 40,000-110,000 mean reads per cell. An estimated number of 5,000–10,000 cells per sample were targeted. NovaSeq 6000 paired end sequencing was performed using 28 cycles for Read 1, 10 cycles for Read i7, 10 cycles for Read i5 and 90 cycles for Read 2, using NovaSeq SP Reagent Kit v1.5 (cat# 20028401, Illumina) and NovaSeq S2 Reagent Kit v1.5 (cat# 20028316, Illumina). Single-cell RNA-sequencing data analysis ScRNA-seq libraries were sequenced on an Illumina Novaseq 6000 using paired-end dual index reads. Demultiplexing and FastQ generation was performed using either bcl2fastq version 2.20.0 (Illumina) or BCLconvert version 3.9.3 (Illumina). FastQ data was aligned and quantified using the Cell Ranger package version 6.1.2 (10X Genomics). All analyses on the single-cell data were performed in R version 4.2.3 using Seurat version 4.3.0 26 . After loading the expression matrix DoubletFinder 27 was used to identify the heterotypic doublets. Cells labeled as doublets, having more than 25% mitochondrial reads or less than 1000 detected features were removed from the expression matrix which was subsequently normalized using SCTransform v2, regressing out cell cycle and mitochondrial fraction 28 . After dimensionality reduction and UMAP projection, the cells were clustered using the Louvain algorithm with a resolution of 1.0. ScopeLoomR was used to export and import the data for the SCENIC analysis. To identify malignant cells, each single-cell experiment was subset based on the AUCell score of Jerby-Arnon malignant gene set (> 0.12) 29 . Single-cell regulatory network inference and clustering (SCENIC) SCENIC was run on the raw count matrix after mitochondrial read and doublet removal in addition to malignant cell identification as mentioned above. SCENIC analysis was performed using pySCENIC version 0.12.1. The co-expression modules procedure was run using GRNboost2, regulon prediction was run using the motif ranking databases hg38 500bp up/100bp down and 10kbp up/down ( https://resources.aertslab.org ). To limit potential undesirable effects of the stochastic nature of the gradient-boosting step within SCENIC, the complete pipeline was run ten times on the four patient samples individually (i.e., patients 35, 37, 44 and 45 in Fig. 1 A). The set of predicted regulons were imported into R and filtered based on recurrence of both the detected regulons (6 out of the 10 runs) as well as the predicted target genes (6 out of 10 runs). The filtered regulons were then used as input for AUCell 30 . Regulon activity and differentially activated regulons The activity matrix generated by AUCell was added as an assay to the Seurat object. After scaling, the AUCell matrix was used for dimensionality reduction, UMAP projection and regulon clustering (Louvain). Differentially activated regulons were identified using Seurat’s FindAllMarkers using the Wilcoxon rank sum test requiring an adjusted P value < 0.05, a log fold-change threshold of 0.01 and a minimum positive fraction of 0.95. To calculate the reciprocal overlap of regulons between two clusters across samples, the ratio of the overlapping regulons within each cluster was averaged. Reactome analysis Differentially active regulon intersect of patients 35, 37, 44 and 45 clusters was introduced to Reactome to perform pathway overrepresentation analysis. Negative log10 FDR values were used to plot the results, as they include a corrected over-representation probability. Spatial transcriptomics To perform spatial transcriptomics, we made use of the Visium Spatial Gene Expression for FFPE platform of 10X Genomics. HE staining, tissue adhesion test quality check and all other procedures were performed as explained by Visium Spatial Gene Expression for FFPE from 10X Genomics ( CG000407, CG000408 and CG000409). Patient tumor material in FFPE was used as input material for the Visium. Region of interest was marked and slices were taken for DV200 estimation (Agilent TapeStation). FFPE blocks were chosen for Visium if the DV200 value was equal to or higher than 50% and had NGFR protein expression. Selected tumor quadrants were sliced and placed on the Visium slide. An extra slice of the same quadrant region was taken along for an additional NGFR staining. Visium Spatial Gene Expression libraries were prepared according to the Visium Spatial Gene Expression User Guide for FFPE (CG000407). For each experiment, the Visium Spatial Gene Expression libraries were pooled to create one sequencing pool. The sequencing pool was quantified by qPCR, according to the KAPA Library Quantification Kit Illumina® Platforms protocol (KR0405, KAPA Biosystems). Paired end sequencing was performed on NextSeq 550 and NovaSeq 6000 Systems (Illumina) using Illumina sequencing reagent kits (cat no. 20024906, cat no. 20028401, Illumina), at a sequencing depth of approximately 30,000–85,000 read pairs per tissue covered spot. Visium libraries were sequenced on an Illumina Novaseq 6000 using paired-end dual index reads. Demultiplexing and FastQ generation was performed using Illumina BCLconvert version 3.9.3 (Illumina). FastQ data was quantified and spatially resolved using the Space Ranger package version 1.3.0 (10X Genomics). All analyses on the spatial data were performed in R version 4.2.3 using Seurat version 4.3.0 26 . The expression matrix was normalized using SCTransform v2, clustered using the Louvain algorithm with a resolution of 1.0. Imaging Visium slide The images have been acquired on an ZeissAxiover 200M microscope equipped with a 20x /0.75 NA Plan-Achromat objective, a 0.55 NA condenser and an Axiocam 512 color camera (pixel size 0.246 x 0.246 µm) in the ZEN 2.3 software. Images were stitched and exported as merged RGB color tiff. Gene signatures scoring AUCell was used to score the enrichment of the different melanoma gene signatures 1 , 5 – 8 , 31 , 32 for both the single-cell and the Visium data. For the Pozniak signatures, the top 100 most significant genes were selected for each signature. Of note, we excluded the Pozniak melanocytic signature for the Visium analyses, since it largely consisted of ribosomal genes which are not measured with the Visium FFPE probes set. The AUC value matrix was transformed to a Z-score for each gene-set. Then, the values were averaged for every Louvain cluster and the resulting matrix was shown as a heatmap using an additional scaling step for each row (gene signature) to emphasize the differences between the clusters of each patient tumor. Spatial mapping of single-cell regulon clusters SCENIC derived regulons from patient tumors 35, 37, 44 and 45 were used for scoring regulon activity, Louvain clustering and determining the differentially active regulons (identical to the single-cell method) on the matched Visiums. Next, the differentially active regulons of the clusters from both the scRNA-seq and the Visium 10X platforms were compared to determine the degree of similarity. When a scRNA-seq regulon scored the highest significant overlap with a Visium regulon cluster, the scRNA-seq regulon cluster number was transferred to the Visium regulon cluster. Analysis of patient RNA-sequencing datasets Raw transcript counts from all datasets (Hugo 33 and Riaz 34) were compiled using Kallisto 0.48.0. Samples with abnormal counts, as detected in principal component analysis (PCA) or ≤5% of the median total transcript counts, were excluded. Similarly, transcripts with 0 counts or total counts ≤0.1% of the median were also left out of the analysis. CombatSeq function from the SVA package (3.42.0) was used for dataset-related batch effect removal. Transcript counts were then transformed into TPMs considering the Kallisto produced length matrix. TPMs were then aggregated by gene and scaled for downstream analysis. To limit non-melanoma cell noise, samples were filtered based on purity and tumor markers. The ESTIMATE deconvolution algorithm in the immunedeconv package (2.1.0) was used to predict tumor purity. Samples that had lower than 0.4 estimated tumor purity score and were below 50% median TPM counts for the sum of the expression of S100B , SOX10 , MLANA , MITF genes were filtered out. Next, tumor samples were scored for the following signatures: the IFNg axis was based on the 28 gene signature by Ayers and colleagues 35 . The APM score was calculated using the eight selected genes by Thompson and colleagues 36 . T cell reactivity and T cell signatures were taken from Chow and colleagues 37 . The tertiary lymphoid structure (TLS) chemokines signature was derived from Li and colleagues 38 . For the TNF signature, we used the PID_TNF_PATHWAY 39 . The BATF3 signature was used as described by Hoefsmit and colleagues 40 . Since antigenicity was calculated using different approaches in the considered datasets, a different calculation method was applied. The antigenicity predictor was scaled for each dataset independently, to avoid biases when employing different methods. For the Riaz and Hugo datasets we took their neopeptides that had an affinity ≤500 nM and/or a ranking percentage ≤2%. If expression data was available, only expressed neopeptides were considered. Final Z-scores where then scaled across all patients. The scaled TPM expression of each signature gene is averaged, then the resulting average is scaled across all patients. The Immune Activity Score (IAS) is calculated as the average of Z-scores for all the different signatures, including antigenicity, and aims to give an indication of the overall T cell presence, activity and priming status in the biopsied sample. Based on the median IAS score, we separated patient tumors in high and low immune active respectively (HI_IAS and LO_IAS). Each cohort contained responders and non-responders (R and NR) creating a final classification of four patient groups. All analyses on the combined RNA-seq data were performed in R version 4.1.3 Single-cell analysis of publicly available datasets Single-cell RNA-sequencing data was used from a study by Pozniak and colleagues 31 . Key marker genes for the mesenchymal-like cell state were calculated using Seurat’s (4.4.0) FindMarkers function across all available before treatment (BT) samples and on the SCT data. Those markers with an adjusted P value 1.5 and pct.1 > 0.2 of the cluster’s cells were considered for further enrichment analysis of the combined RNA-seq dataset. Analyses were performed in R version 4.1.3. Statistics To compare two means, a two-tailed Student T test was used. One-way ANOVA test for more than 2 comparisons. For bulk RNA-seq, all signatures, scores and gene expression comparisons were statistically assessed using Wilcoxon rank tests. Statistics were performed by Prism (Graphpad Software Inc., version 9.0) or in R (4.1.3.). A P value of lower than 0.05 was regarded as being statistically significant. Results NGFR-expressing melanomas are highly heterogeneous for differentiation and dedifferentiation markers To determine the degree of heterogeneity within NGFR + human melanomas, we first determined the expression of differentiation and dedifferentiation markers. We stained a panel of 65 clinical human melanoma samples derived from 46 patients for the differentiation markers MITF and Melan-A and the dedifferentiation markers AXL and NGFR 7 – 11 (Figure S1A) . Because we wished to study NGFR + melanoma subpopulations, we focused on 41 (63% of all stained patient tumors) of those, displaying a range of 1-100% NGFR positivity on viable tumor cells ( Fig. 1 A ) . We observed that across these 41 tumors, NGFR correlated inversely with both Melan-A and MITF expression, while positively correlating with AXL, in agreement with previous observations ( Fig. 1 B-D ) 8 , 11 , 15 . MITF and Melan-A mark a differentiated melanocytic phenotype, whereas NGFR marks a neural crest stem cell (NCSC) phenotype and AXL an undifferentiated mesenchymal phenotype ( Fig. 1 E ) 7 , 41 . However, we noted that > 50% of NGFR + melanomas (22 out of 41) also comprised regions positive for differentiation (Melan-A and/or MITF) and dedifferentiation (AXL) markers ( Fig. 1 A-I and S1 B-E). This included patient melanomas #35, 37, 44 and 45, which showed different degrees and patterns of marker heterogeneity and for which sufficient high-quality material was available for several downstream analyses ( Fig. 1 A-D, 1 F-I and S1 B-E ) . Thus, the NGFR + fraction does not represent a single population but instead is heterogeneous, comprising a phenotypically diverse pool of differentiated and undifferentiated melanoma cells, and possibly intermediate cell states. This observation led us to investigate NGFR + heterogeneity in more detail, focusing on single cells, inter-patient variation, spatial organization, functional consequences and clinical responses. Therefore, we set up a pipeline to determine both transcriptional and spatial heterogeneity in NGFR + patient melanomas, combining scRNA-seq on FACS-sorted NGFR + viable cells with spatial transcriptomics on matched formalin-fixed paraffin-embedded (FFPE) samples ( Fig. 1 J ) , which was complemented with functional analyses and clinical responses to immunotherapy. ​​​ ​ ScRNA-seq identifies several transcriptomically distinct NGFR + melanoma phenotypes First, we set out to investigate whether the NGFR + fraction shows heterogeneity in gene expression at the single-cell level. ScRNA-seq was performed on FACS-sorted NGFR + melanoma cells from the above-mentioned four patient melanomas, from whom a total of 24,953 NGFR + /CD45 − cells were isolated, expressing a wide range of cell surface NGFR expression levels ( Fig. 2A-B and S2 A ) . Using Uniform Manifold Approximation and Projection (UMAP) on the gene expression as a dimension reduction approach and Louvain clustering, we observed that NGFR + cells were represented by several different gene expression clusters. Melanoma cells from patients 35, 37 and 44 were grouped into 13 clusters and patient 45 in 17 clusters (Figure S2B) . Expression of several markers, namely SOX10 , S100B , MITF and MLANA , demonstrated that most sorted NGFR + cells were positive for markers from the melanocytic lineage (Figure S2C). Furthermore, almost all sorted cells scored positively for a malignant melanoma signature 29 , confirming that we had selected cancer cells (Figure S2D). Those melanoma cells were then subjected to a second round of UMAP analysis and Louvain clustering, which revealed 12, 12, 10 and 15 clusters for patients 35, 37, 44 and 45, respectively (Fig. 2C) . To determine the phenotypic identities of the NGFR + single cells within each patient tumor, the expression of key melanoma markers 7 – 11 , 32 , 41 – 43 was assessed and a series of melanoma signatures that were previously reported and commonly used were scored 1 , 5 – 8 , 31 , 32 . This was done using AUCell 30 , an algorithm that quantifies the enrichment of input signature genes within the expressed genes of each single cell (Fig. 2D) . All clusters showed absence of PTPRC expression (encoding CD45), confirming the exclusion of immune cells. The signatures were grouped into five main categories: differentiated, undifferentiated, transitory/neural crest-like, stress-like and mitotic. The gene expression clusters from all four patients showed great diversity for these signatures, implying that cells that are positive for the neural crest/stem cell marker NGFR can exist in several different states along the (de)differentiation trajectory. Gene expression analysis of melanoma phenotype switch markers (i.e., MITF , Melan-A , NGFR , AXL , PDGFRA , PDGFRB and EGFR ) confirmed the presence of both differentiated and undifferentiated NGFR + subpopulations (Fig. 2D) . Strikingly, several clusters scored positively for multiple signature categories. For example, patient 35 cluster 9 scored highly positive for undifferentiated, transitory/neural crest-like and stress-like categories, while patient 45 cluster 7, scored positively for differentiated, transitory/neural crest-like and stress-like groups. Together, these results indicate that individual patients harbor an NGFR + fraction comprising subpopulations that exist at different positions along the (de)differentiation trajectory. NGFR + melanoma cells are characterized by patient-specific gene-regulatory networks We next sought to find a mechanistic explanation for the observed great inter-patient diversity in gene expression within the NGFR + cell fractions. We considered the possibility that NGFR marks different melanoma cell states governed by unique master regulators. To investigate this, we employed motif based single-cell regulatory network inference and clustering (SCENIC) on the scRNA-seq data followed by UMAP analysis and Louvain clustering. SCENIC is a computational tool for gene-regulatory network (GRN) inference and cell state discovery 30 , taking a three-step approach: inference of candidate target genes through co-expression inference, followed by motif-based filtering and AUC transcription factor (TF) activity scoring 45 . This allowed us to first identify key underlying GRNs of relative cell states within patients and, second, to interrogate whether these GRNs are shared between patients. SCENIC revealed 9–13 independent cell states or regulon clusters within the different patients ( Fig. 3 A ). Most regulon clusters, including cluster 4 in patient 37 and cluster 3 in patient 44, were characterized by dozens of differentially active TFs ( Fig. 3 B and S3 A). Of note, an overlay of gene expression and regulon clusters revealed that some gene expression clusters were associated with the same regulons. This was evident for regulon cluster 2 in patient 35, which is composed of gene expression clusters 2 and 6; for regulon cluster 4 in patient 37 (gene expression clusters 8, and 9 and part of 1) and cluster 7 in patient 45 (gene expression clusters 12 and part of 0) (Figure S3B) . These results suggest that NFGR marks several unique cell states per patient that are characterized by highly diverse GRNs. Next, the overlap between regulons across the four patients was determined. We found 103–145 regulons per patient, 50 of which were shared among the four patients. Using Reactome enrichment analysis we observed a common enrichment across the four samples for Nerve Growth Factor (NGF)-stimulated transcription, consistent with our NGFR-focused scRNA-seq approach (Figure S3C ). Many other regulons were unique to individual patients, raising the possibility that cell states are established in a patient-dependent fashion. To investigate this further, we determined the average reciprocal overlap of differentially active regulons between clusters across patients. Only significant regulons per cluster were considered as cluster-specific cell state signatures. Little overlap was observed between clusters of the patients investigated ( Fig. 3 D ) . Together, these results indicate that NGFR + melanoma cells are characterized by patient-specific GRNs. Spatial organization of NGFR + subpopulations While highlighting the diversity of melanoma GRNs, the analyses above do not consider how the corresponding cells are spatially organized. We considered this a key parameter to include in our study, given the great intrinsic diversity of NGFR + melanoma fractions as revealed by our immunohistochemical ( Fig. 1 ) and scRNA-seq analyses ( Figs. 2 and 3 ) . Therefore, we used Visium 10X spatial transcriptomics to dissect the spatial and transcriptomic NGFR heterogeneity in the same patient melanomas. For this analysis, selection of tumor quadrant areas was guided by the presence of NGFR protein expression ( Fig. 4 A ). Two slices were taken of the quadrant area of each single FFPE block. One was used for H&E staining on the Visium 10X slide, and the other one was stained for NGFR ( Fig. 4 B and S4 A ) . In agreement with the NGFR IHC staining results ( Fig. 4 B ) , NGFR expression was detected in the Visium 10X quadrants, again showing highly heterogeneous patterns ( Fig. 4 C ) . Furthermore, all samples showed large areas scoring positively for the malignant melanoma signature 29 , indicative of the presence of melanoma cells (Figure S4B) . To investigate the spatial gene expression patterns of the NGFR + tumor regions, gene expression Louvain clustering was performed, followed by melanoma signature scoring and phenotype switch marker gene expression analysis. Among the patients, we identified 10–13 unique gene expression clusters that were spatially organized ( Fig. 4 D ) . Next, we scored these clusters for a series of common signatures representative of different melanoma cell states 1 , 5 – 8 , 31 , 32 . We observed again a high degree of variation in melanoma cell state signatures and phenotype switch marker genes ( Fig. 4 E ) . This result supports our scRNA-seq analyses showing that NGFR + melanoma cells can exist in several different states along the (de)differentiation trajectory (Fig. 2D) . To complement this analysis with spatial annotation, we determined the location of the scRNA-seq regulon clusters. We first generated Visium regulon clusters for each patient. The single cell regulon list was used as input to score the Visium spots, employing AUCell followed by Louvain clustering (Figure S4C) . Next, we highlighted scRNA-seq regulon clusters that most significantly overlapped with Visium regulon clusters in the Visium spots. Whereas most scRNA-seq regulon clusters were spatially organized, there were clear patterns of mutual cell spreading to other compartments, which was apparent for all patients. For example, in patient 44, single cell regulon clusters 2 and 5 formed spatially different compartments, but each formed several satellites in the other cluster ( Fig. 4 F ) . Similarly, satellite formation in other regulon clusters was apparent for clusters 0, 5 and 10 in patient 45, and for clusters 2, 6 and 8 in patient 35. Thus, we observe a pattern, common among melanoma patients, of NGFR + cell states that are mostly spatially organized yet also form satellites in spatial regions dominated by other cell states. PDGFR marks cytokine-resistant and mesenchymal-like NGFR + melanoma cells We previously reported that NGFR + melanoma populations are less sensitive to T cell cytokines 19 . The findings above indicate that even within the melanoma NGFR + cell fraction there is great transcriptional and phenotypic heterogeneity. This raises the possibility that this diversity is associated with functional differences. We thus investigated any functional consequences of NGFR heterogeneity for cytokine sensitivity. A panel of representative NGFR + melanoma cell lines was selected using the established phenotype switch markers Melan-A, AXL, PDGFRa and PDGFRb ( Fig. 5 A and S5 A-C ) . Next, we examined the cytotoxic potential of T cell cytokines IFNg and TNF across this cell line panel. We observed that specifically the cell lines marked by PDGFRa and PDGFRb expression (“PDGFR”, i.e., A875, M080.X1.CL and M019R.X1.CL), showed the highest cytokine resistance ( Fig. 5 B, C and S5 D ) . These results show that while NGFR + melanoma fractions already show reduced sensitivity to cytokines 19 , NGFR + /PDGFR + double-positive tumor cells show an even higher degree of cytokine resistance. To complement these functional and our Visium spatial analyses, we next examined the location of NGFR + /PDGFR + double-positive cells in clinical melanoma samples. We performed NGFR and PDGFRb IHC staining on consecutive FFPE slices combined with Periodic Acid–Schiff (PAS) to distinguish the tumor from the stromal component (i.e., vessels and collagen). PDGFRb expression was detected both on melanoma cells in NGFR + tumor regions and in the vascular area, the latter being indicative of pericytes ( Fig. 5 D ) . To validate these observations, we analyzed a scRNA-seq melanoma cohort 31 . We found that the NGFR + /PDGFR + melanoma cell fraction was highly enriched in the mesenchymal-like cell fraction (“MES”) in baseline melanomas ( Fig. 5 E ) . In the other three patients studied we were unable to convincingly identify NGFR + /PDGFRb + melanoma cells (Figure S5E) , confirming the low frequency of this double-positive population seen in the scRNA-seq analyses (Fig. 2D) . Therefore, we extended the panel of melanoma samples with three additional patients (7b, 36, 43). This analysis validated our finding of NGFR + /PDGFRb + melanoma cells, which were again often located near the stromal component including vessels ( Fig. 5 D ) . Together, PDGFR marks a cytokine-resistant melanoma subpopulation, which can be detected in patient samples and is associated with a mesenchymal-like phenotype. PDGFR marks immune-active yet ICB non-responding melanomas Lastly, we investigated whether this PDGFR + mesenchymal-like and cytokine-resistant melanoma population correlates with therapeutic response to ICB. The presence of T cells in the tumor microenvironment is required for ICB response and is associated with improved patient outcome 29 , 46 – 48 . An immune activity score (IAS low, intermediate, high) was generated based on the sum of averaged Z-score expression of each of the following gene signatures, namely: IFNg 35 , antigen presentation machinery (APM) 36 , T cells 37 , T cell reactivity 37 , TLS-chemokines 38 , TNF 39 and BATF3 40 , as well as antigenicity (i.e., total neoantigen counts) ( Fig. 6 A and S6 A ) . Two public RNA-seq melanoma datasets (i.e., Hugo 33 and Riaz 34 ) were combined for this analysis. We observed that both IAS-high and IAS-low patients were similarly distributed among responders ( HI_IAS/R and LO_IAS/R) and non-responders ( HI_IAS/NR and LO_IAS/NR, Fig. 6 A and 6 B ) . To determine whether the cytokine-resistant melanoma population marked by NGFR and PDGFR was associated with poor ICB response, expression of both PDGFR and NGFR was analyzed in the combined RNA-seq dataset ( Fig. 6 A and S6 B ) . We noticed a significantly higher expression of PDGFR in immune-active ICB non-responders (HI_IAS/NR) compared to both non-responders with low immune activation (LO_IAS/NR) and responders with a high immune activation score (HI_IAS/R) ( Fig. 6 C and S6 B ) . Supporting this observation, there was a similar significant enrichment of the MES signature in HI_IAS/NR tumors compared to both LO_IAS/NR and HI_IAS/R groups ( Fig. 6 D ) . This suggests that PDGFR + mesenchymal-like melanoma cells are significantly enriched in immune-active, yet ICB non-responding, melanoma patients. Lastly, to complement our RNA-seq findings, we investigated the presence of tumor- infiltrating T cells in PDGFRb + and PDGFRb − tumors by staining 33 melanoma samples for both PDGFRb/PAS and CD3 (with or without PAS). We identified 10 tumors harboring PDGFRb + melanoma cells, thereby confirming that PDGFRb is expressed in a subfraction of patient melanomas (Figure S6C) . Of these, nine tumors (90%) were infiltrated with T cells. Although PDGFRb − tumors showed a slightly higher average CD3 score (fold change 1.33) than PDGFRb + tumors, this difference was not statistically significant ( Fig. 6 E ) . Moreover, IHC analysis revealed that in most PDGFRb + melanomas, CD3 + T cells showed no clear association with either PDGFRb High or PDGFRb Low tumor regions (i.e., 7/9 infiltrated tumors (78%), Fig. 6 F ) . Together, these results show that PDGFR + mesenchymal-like tumor cells are enriched in T cell-infiltrated, immune-active, but ICB non-responding melanomas. DISCUSSION Melanoma is positioned at the extreme end of the mutational spectrum in cancer 2 , 50 . Combined with its high degree of phenotypic plasticity, this results in substantial heterogeneity, as we and others have shown previously 1 , 4 – 11 . For example, a small yet common population of melanoma cells marked by NGFR + expression is associated with reduced susceptibility to T cell killing 19 , 15 , NK-cell killing 24 , immune exclusion 19 , immunotherapy resistance 20 and phenotypic plasticity 18 , both in patients and mouse models. However, it is incompletely understood whether this NGFR + fraction represents a homogenous cell group or instead comprises yet additional heterogeneous cell populations, and whether this has consequences for therapy response. By integrating scRNA-seq, Visium spatial analysis, functional studies and analyses of clinical ICB cohorts, we find an unexpectedly high degree of heterogeneity even within these NGFR + cell fractions. This is manifested at the level of both expression and regulation of transcriptional networks, phenotypically, functionally and clinically. To begin dissecting NGFR heterogeneity in clinical samples, we employed IHC, confirming that NGFR + cells can be identified in most human melanomas. Previous studies have shown NGFR expression to be associated with an AXL high program, specifically marking the NCSC state in melanoma 7 , 8 , 11 , 19 , 51 . However, while we were able to confirm the positive association between NGFR and dedifferentiation marker AXL in clinical samples, we also observed expression of differentiation markers MITF and/or Melan-A in NGFR + tumor regions. This prompted us to investigate whether this heterogeneity is also seen at the single cell level. In vitro studies revealed that NGFR, after initial upregulation, is strongly downregulated when reaching a full undifferentiated/mesenchymal-like cell state upon long-term MAPK inhibition 17 . However, whether NGFR marks multiple melanoma transcriptional programs along the (de)differentiation scale in patient melanomas has not been addressed. By employing scRNA-seq, we observed that indeed, NGFR marks a spectrum of (de)differentiation in clinical melanoma samples. Specifically, NGFR + cell populations display differentiated, undifferentiated, transitory/neural crest-like, stress-like and mitotic signature categories, as well as combinations of these, within each patient, to different extents. This indicates that NFGR marks a transcriptomically diverse pool of melanoma cells in clinical melanoma samples. Considering the observed transcriptomic heterogeneity within the NGFR + cell fraction, we hypothesized that the different cell states marked by NGFR are governed by multiple gene regulatory networks (GRNs). To test this, SCENIC 30 , 45 was used to analyze the scRNA-seq patient data. At least nine different NGFR + cell states per patient were identified, each characterized by a unique GRN. It was recently shown that melanoma GRNs show diversity among individual patients 31 . Observing no common GRN across the patients studied here, our results demonstrate that this heterogeneity among patients even extends to the NGFR + subpopulations. Thus, NGFR marks patient-specific cell states characterized by unique GRNs. To capture the spatial distribution of NGFR + melanoma subpopulations, Visium 10X spatial transcriptomics was employed. Previous studies have shown the spatial distribution of NGFR + melanoma cells in skin- 20 , metastatic- 20 , 52 and patient-derived xenografts (PDX) lesions 53 . However, an in-depth spatial transcriptomic characterization of NGFR + tumor regions per patient was lacking. We show here that consistent with our scRNA-seq findings, there is a spatial spectrum of (de)differentiation within NGFR + patient tumor regions. Moreover, by combining scRNA-seq and Visium, we show that also the spatial distribution of NGFR + cell states is governed by unique GRNs. We observed that cell states are spatially organized in melanoma, often forming satellites in tumor areas dominated by other cell states. A potential force driving satellite formation may be stochastic variability in gene expression 54 with or without reinforcement of Lamarckian induction by microenvironmental stimuli and physical cues 3 , 55 . In other words, it is plausible that local microenvironmental changes are intertwined with the outgrowth of a particular cell state or affects cell state identity 54 , 55 . For instance, the stem-like murine melanoma phenotype was found to be associated with, and fueled by, endothelial cells 6 , 55 . Whether stochastic gene expression variation and Lamarckian induction are independently or cooperatively causing satellite formation of different NGFR + cell states in melanoma needs to be further investigated. From this study, we can conclude that while NGFR + cell states are mostly spatially organized, they also form satellites in tumor regions dominated by other cell states across different patients. Considering the transcriptomic and phenotypic diversity within the NGFR + melanoma cell fractions, we studied any functional consequence of this heterogeneity for immune cytokine sensitivity. We made use of a representative NGFR + cell line panel co-expressing additional phenotype switch markers. We observed that PDGFR expression was consistently associated with T cell cytokine resistance. No association between AXL expression and cytokine resistance was observed, which we had also noted in the context of T cell resistance 19 . By employing IHC, we identified NGFR + /PDGFR + cells often near the stroma/vessels. By leveraging an external scRNA-seq dataset, we validated the presence of NGFR + /PDGFR + cells in additional patient melanomas. These double-positive cells were enriched in the mesenchymal-like (MES) cell state tumor fraction 31 . Thus, while NGFR marks a relatively insensitive T cell cytokine-phenotype 19 , here we show that particularly PDGFR marks the most cytokine-resistant NGFR + mesenchymal-like subpopulations. Immunotherapy has improved overall survival in advanced melanoma, reaching a clinical benefit over 50% 56–58 . However, a significant proportion of patients either intrinsically harbors or develops acquired resistance to immunotherapy, due to diverse cell intrinsic or extrinsic mechanisms 32 , 59 , 60 . It is therefore imperative to explore novel treatment options 61 while simultaneously investing in biomarker identification to improve prediction of therapy resistance. This may also help circumventing unnecessary immune-related adverse events (irAEs) 62 . Thus, identifying upfront therapy resistance is a relevant clinical need 63 . Since PDGFR marks a mesenchymal-like cytokine-resistant melanoma population, we investigated whether this population was association with a therapeutic response to ICB. By analyzing two melanoma clinical datasets, we identified a subset of non-responders characterized by high immune activity, as reflected by the presence of active T cells, enriched with PDGFR and MES signature expression. Recently, it has been reported that an increase in the MES cell fraction, governed by TCF4, in early ICB on-treatment patient tumor biopsies is associated with ICB resistance 31 . In the current study, we show that the presence of PDGFR + MES cells at baseline is predictive of ICB resistance, particularly in the context of immune-active and T cell-infiltrated melanomas. Of note, TCF4 was not identified as a significant differentially active regulon in any of the studied scRNA-seq regulon clusters. This implies that TCF4 is not a key transcription factor driving any of the characterized NGFR + melanoma cell states, including the PDGFR + cell states. Given that NGFR + /PDGFR + melanoma cells pre-exist in untreated melanomas, a challenge for the future is to investigate which cell-autonomous or non-autonomous factors 53 – 55 contribute to their development and affect therapy. In conclusion, we show here a remarkable degree of heterogeneity in a heavily studied melanoma subpopulation, namely NGFR + cells. Several NGFR + melanoma subpopulations position along the (de)differentiation gray scale, driven by patient-specific GRNs and spatially organized as satellites in tumor regions dominated by other cell states. The co-expression of PDGFR, marking T cell cytokine resistance in immune-active ICB non-responding patients may serve as a biomarker for immunotherapy resistance, which may have diagnostic potential. Declarations Disclosure D.S.P. is co-founder, shareholder and advisor of Flindr Tx, which is unrelated to this study. ACKNOWLEDGEMENTS We thank all the members of the Peeper lab for their valuable input. Moreover, we thank our in-house flowcytometry, animal pathology, sequencing, and bioimaging facilities for their help and support. We would like to acknowledge Alexander van Akkooi, Winan van Houdt, John B. Haanen, Lisanne Zijlker and Max F. Madu for collecting patient samples. Moreover, we thank the NKI-AVL Core Facility Molecular Pathology & Biobanking (CFMPB) for supplying NKI-AVL Biobank material and performing IHC stainings. Lastly, we thank Joyce Sanders for helping with FFPE tumor block selection for Visium 10X. D.S.Peeper is funded by the Oncode Institute and by the Dutch Cancer Society KWF and this work was supported by the Koningin Wilhelmina Fonds (KWF; project no. 10425) and by MRA 681127 to D.S. Peeper. AUTHOR CONTRIBUTIONS: S.M.S and D.S.P: conceived study and designed experiments. S.M.S: performed experiments, investigation, visualization, (bioinformatic) data interpretation, writing–original draft and project administration. J.J.H.T: performed bioinformatic analyses and provided critical input. A.V , I.d.R , A.v.V , A.G and J.S.N : performed bioinformatic analyses. M.N: processed Visium 10X samples. S.B : Helped collecting clinical samples. S.M.S , J-Y.S, J.B, S.B, M.D.S.G and H.M.H : analyzed human pathology samples. I.H and L.d.V : performed IHC-, PAS- and HE stainings and prepared/sliced FFPE blocks for Visium 10X. A.E.D : imaging of Visium 10X (HE stained) patient samples. M.v.B : performed cell sorting. S.M.S and D.S.P: wrote the manuscript. The project was supervised by D.S.P . All authors reviewed and approved the manuscript. References Wouters, J., Kalender-Atak, Z., Minnoye, L., Spanier, K. I., De Waegeneer, M., Bravo González-Blas, C., Mauduit, D., Davie, K., Hulselmans, G., Najem, A., Dewaele, M., Pedri, D., Rambow, F., Makhzami, S., Christiaens, V., Ceyssens, F., Ghanem, G., Marine, J. C., Poovathingal, S., & Aerts, S. (2020). Robust gene expression programs underlie recurrent cell states and phenotype switching in melanoma. Nature cell biology , 22 (8), 986–998. https://doi.org/10.1038/s41556-020-0547-3 Alexandrov, L. B., Nik-Zainal, S., Wedge, D. C., Aparicio, S. A., Behjati, S., Biankin, A. V., Bignell, G. R., Bolli, N., Borg, A., Børresen-Dale, A. L., Boyault, S., Burkhardt, B., Butler, A. P., Caldas, C., Davies, H. R., Desmedt, C., Eils, R., Eyfjörd, J. E., Foekens, J. A., Greaves, M., … Stratton, M. R. (2013). Signatures of mutational processes in human cancer. Nature , 500 (7463), 415–421. https://doi.org/10.1038/nature12477 Marine, J. C., Dawson, S. J., & Dawson, M. A. (2020). Non-genetic mechanisms of therapeutic resistance in cancer. Nature reviews. Cancer , 20 (12), 743–756. https://doi.org/10.1038/s41568-020-00302-4 Hoek, K. S., Eichhoff, O. M., Schlegel, N. C., Döbbeling, U., Kobert, N., Schaerer, L., Hemmi, S., & Dummer, R. (2008). In vivo switching of human melanoma cells between proliferative and invasive states. Cancer research , 68 (3), 650–656. https://doi.org/10.1158/0008-5472.CAN-07-2491 Verfaillie, A., Imrichova, H., Atak, Z. K., Dewaele, M., Rambow, F., Hulselmans, G., Christiaens, V., Svetlichnyy, D., Luciani, F., Van den Mooter, L., Claerhout, S., Fiers, M., Journe, F., Ghanem, G. E., Herrmann, C., Halder, G., Marine, J. C., & Aerts, S. (2015). Decoding the regulatory landscape of melanoma reveals TEADS as regulators of the invasive cell state. Nature communications , 6 , 6683. https://doi.org/10.1038/ncomms7683 Karras, P., Bordeu, I., Pozniak, J., Nowosad, A., Pazzi, C., Van Raemdonck, N., Landeloos, E., Van Herck, Y., Pedri, D., Bervoets, G., Makhzami, S., Khoo, J. H., Pavie, B., Lamote, J., Marin-Bejar, O., Dewaele, M., Liang, H., Zhang, X., Hua, Y., Wouters, J., … Marine, J. C. (2022). A cellular hierarchy in melanoma uncouples growth and metastasis. Nature , 610 (7930), 190–198. https://doi.org/10.1038/s41586-022-05242-7 Rambow, F., Rogiers, A., Marin-Bejar, O., Aibar, S., Femel, J., Dewaele, M., Karras, P., Brown, D., Chang, Y. H., Debiec-Rychter, M., Adriaens, C., Radaelli, E., Wolter, P., Bechter, O., Dummer, R., Levesque, M., Piris, A., Frederick, D. T., Boland, G., Flaherty, K. T., … Marine, J. C. (2018). Toward Minimal Residual Disease-Directed Therapy in Melanoma. Cell , 174 (4), 843–855.e19. https://doi.org/10.1016/j.cell.2018.06.025 Tsoi, J., Robert, L., Paraiso, K., Galvan, C., Sheu, K. M., Lay, J., Wong, D. J. L., Atefi, M., Shirazi, R., Wang, X., Braas, D., Grasso, C. S., Palaskas, N., Ribas, A., & Graeber, T. G. (2018). Multi-stage Differentiation Defines Melanoma Subtypes with Differential Vulnerability to Drug-Induced Iron-Dependent Oxidative Stress. Cancer cell , 33 (5), 890–904.e5. https://doi.org/10.1016/j.ccell.2018.03.017 Müller, J., Krijgsman, O., Tsoi, J., Robert, L., Hugo, W., Song, C., Kong, X., Possik, P. A., Cornelissen-Steijger, P. D., Geukes Foppen, M. H., Kemper, K., Goding, C. R., McDermott, U., Blank, C., Haanen, J., Graeber, T. G., Ribas, A., Lo, R. S., & Peeper, D. S. (2014). Low MITF/AXL ratio predicts early resistance to multiple targeted drugs in melanoma. Nature communications , 5 , 5712. https://doi.org/10.1038/ncomms6712 Konieczkowski, D. J., Johannessen, C. M., Abudayyeh, O., Kim, J. W., Cooper, Z. A., Piris, A., Frederick, D. T., Barzily-Rokni, M., Straussman, R., Haq, R., Fisher, D. E., Mesirov, J. P., Hahn, W. C., Flaherty, K. T., Wargo, J. A., Tamayo, P., & Garraway, L. A. (2014). A melanoma cell state distinction influences sensitivity to MAPK pathway inhibitors. Cancer discovery , 4 (7), 816–827. https://doi.org/10.1158/2159-8290.CD-13-0424 Tirosh, I., Izar, B., Prakadan, S. M., Wadsworth, M. H., 2nd, Treacy, D., Trombetta, J. J., Rotem, A., Rodman, C., Lian, C., Murphy, G., Fallahi-Sichani, M., Dutton-Regester, K., Lin, J. R., Cohen, O., Shah, P., Lu, D., Genshaft, A. S., Hughes, T. K., Ziegler, C. G., Kazer, S. W., … Garraway, L. A. (2016). Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science (New York, N.Y.) , 352 (6282), 189–196. https://doi.org/10.1126/science.aad0501 Boiko, A. D., Razorenova, O. V., van de Rijn, M., Swetter, S. M., Johnson, D. L., Ly, D. P., Butler, P. D., Yang, G. P., Joshua, B., Kaplan, M. J., Longaker, M. T., & Weissman, I. L. (2010). Human melanoma-initiating cells express neural crest nerve growth factor receptor CD271. Nature , 466 (7302), 133–137. https://doi.org/10.1038/nature09161 Quintana, E., Shackleton, M., Foster, H. R., Fullen, D. R., Sabel, M. S., Johnson, T. M., & Morrison, S. J. (2010). Phenotypic heterogeneity among tumorigenic melanoma cells from patients that is reversible and not hierarchically organized. Cancer cell , 18 (5), 510–523. https://doi.org/10.1016/j.ccr.2010.10.012 Mehta, A., Kim, Y. J., Robert, L., Tsoi, J., Comin-Anduix, B., Berent-Maoz, B., Cochran, A. J., Economou, J. S., Tumeh, P. C., Puig-Saus, C., & Ribas, A. (2018). Immunotherapy Resistance by Inflammation-Induced Dedifferentiation. Cancer discovery , 8 (8), 935–943. https://doi.org/10.1158/2159-8290.CD-17-1178 Landsberg, J., Kohlmeyer, J., Renn, M., Bald, T., Rogava, M., Cron, M., Fatho, M., Lennerz, V., Wölfel, T., Hölzel, M., & Tüting, T. (2012). Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation. Nature , 490 (7420), 412–416. https://doi.org/10.1038/nature11538 Fallahi-Sichani, M., Becker, V., Izar, B., Baker, G. J., Lin, J. R., Boswell, S. A., Shah, P., Rotem, A., Garraway, L. A., & Sorger, P. K. (2017). Adaptive resistance of melanoma cells to RAF inhibition via reversible induction of a slowly dividing de-differentiated state. Molecular systems biology , 13 (1), 905. https://doi.org/10.15252/msb.20166796 Su, Y., Wei, W., Robert, L., Xue, M., Tsoi, J., Garcia-Diaz, A., Homet Moreno, B., Kim, J., Ng, R. H., Lee, J. W., Koya, R. C., Comin-Anduix, B., Graeber, T. G., Ribas, A., & Heath, J. R. (2017). Single-cell analysis resolves the cell state transition and signaling dynamics associated with melanoma drug-induced resistance. Proceedings of the National Academy of Sciences of the United States of America , 114 (52), 13679–13684. https://doi.org/10.1073/pnas.1712064115 Restivo, G., Diener, J., Cheng, P. F., Kiowski, G., Bonalli, M., Biedermann, T., Reichmann, E., Levesque, M. P., Dummer, R., & Sommer, L (2017) . The low affinity neurotrophin receptor CD271 regulates phenotype switching in melanoma. Nat Commun 8 , 1988. https://doi.org/10.1038/s41467-017-01573-6 Boshuizen, J., Vredevoogd, D. W., Krijgsman, O., Ligtenberg, M. A., Blankenstein, S., de Bruijn, B., Frederick, D. T., Kenski, J. C. N., Parren, M., Brüggemann, M., Madu, M. F., Rozeman, E. A., Song, J. Y., Horlings, H. M., Blank, C. U., van Akkooi, A. C. J., Flaherty, K. T., Boland, G. M., & Peeper, D. S. (2020) . Reversal of pre-existing NGFR-driven tumor and immune therapy resistance. Nat Commun 11 , 3946. https://doi.org/10.1038/s41467-020-17739-8 Liu, D., Lin, J. R., Robitschek, E. J., Kasumova, G. G., Heyde, A., Shi, A., Kraya, A., Zhang, G., Moll, T., Frederick, D. T., Chen, Y. A., Wang, S., Schapiro, D., Ho, L. L., Bi, K., Sahu, A., Mei, S., Miao, B., Sharova, T., Alvarez-Breckenridge, C., … Boland, G. M. (2021). Evolution of delayed resistance to immunotherapy in a melanoma responder. Nature medicine , 27 (6), 985–992. https://doi.org/10.1038/s41591-021-01331-8 Riesenberg, S., Groetchen, A., Siddaway, R., Bald, T., Reinhardt, J., Smorra, D., Kohlmeyer, J., Renn, M., Phung, B., Aymans, P., Schmidt, T., Hornung, V., Davidson, I., Goding, C. R., Jönsson, G., Landsberg, J., Tüting, T., & Hölzel, M. (2015). MITF and c-Jun antagonism interconnects melanoma dedifferentiation with pro-inflammatory cytokine responsiveness and myeloid cell recruitment. Nature communications , 6 , 8755. https://doi.org/10.1038/ncomms9755 Kim, Y. J., Sheu, K. M., Tsoi, J., Abril-Rodriguez, G., Medina, E., Grasso, C. S., Torrejon, D. Y., Champhekar, A. S., Litchfield, K., Swanton, C., Speiser, D. E., Scumpia, P. O., Hoffmann, A., Graeber, T. G., Puig-Saus, C., & Ribas, A. (2021). Melanoma dedifferentiation induced by IFN-γ epigenetic remodeling in response to anti-PD-1 therapy. The Journal of clinical investigation , 131 (12), e145859. https://doi.org/10.1172/JCI145859 Furuta, J., Inozume, T., Harada, K., & Shimada, S. (2014). CD271 on melanoma cell is an IFN-γ-inducible immunosuppressive factor that mediates downregulation of melanoma antigens. The Journal of investigative dermatology , 134 (5), 1369–1377. https://doi.org/10.1038/jid.2013.490 Lehmann, J., Caduff, N., Krzywińska, E., Stierli, S., Salas-Bastos, A., Loos, B., Levesque, M. P., Dummer, R., Stockmann, C., Münz, C., Diener, J., & Sommer, L. (2023). Escape from NK cell tumor surveillance by NGFR-induced lipid remodeling in melanoma. Science advances , 9 (2), eadc8825. https://doi.org/10.1126/sciadv.adc8825 Schieven, S. M., Traets, J. J. H., Vliet, A. V., Baalen, M. V., Song, J. Y., Guimaraes, M. D. S., Kuilman, T., & Peeper, D. S. (2023). The Elongin BC Complex Negatively Regulates AXL and Marks a Differentiated Phenotype in Melanoma. Molecular Cancer Research : MCR , 21 (5), 428–443. https://doi.org/10.1158/1541-7786.MCR-22-0648 Hao, Y., Hao, S., Andersen-Nissen, E., Mauck, W. M., 3rd, Zheng, S., Butler, A., Lee, M. J., Wilk, A. J., Darby, C., Zager, M., Hoffman, P., Stoeckius, M., Papalexi, E., Mimitou, E. P., Jain, J., Srivastava, A., Stuart, T., Fleming, L. M., Yeung, B., Rogers, A. J., … Satija, R. (2021). Integrated analysis of multimodal single-cell data. Cell , 184 (13), 3573–3587.e29. https://doi.org/10.1016/j.cell.2021.04.048 McGinnis, C. S., Murrow, L. M., & Gartner, Z. J. (2019). DoubletFinder: Doublet Detection in Single-Cell RNA Sequencing Data Using Artificial Nearest Neighbors. Cell systems , 8 (4), 329–337.e4. https://doi.org/10.1016/j.cels.2019.03.003 Choudhary, S., & Satija, R. (2022). Comparison and evaluation of statistical error models for scRNA-seq. Genome biology , 23 (1), 27. https://doi.org/10.1186/s13059-021-02584-9 Jerby-Arnon, L., Shah, P., Cuoco, M. S., Rodman, C., Su, M. J., Melms, J. C., Leeson, R., Kanodia, A., Mei, S., Lin, J. R., Wang, S., Rabasha, B., Liu, D., Zhang, G., Margolais, C., Ashenberg, O., Ott, P. A., Buchbinder, E. I., Haq, R., Hodi, F. S., … Regev, A. (2018). A Cancer Cell Program Promotes T Cell Exclusion and Resistance to Checkpoint Blockade. Cell , 175 (4), 984–997.e24. https://doi.org/10.1016/j.cell.2018.09.006 Aibar, S., González-Blas, C. B., Moerman, T., Huynh-Thu, V. A., Imrichova, H., Hulselmans, G., Rambow, F., Marine, J. C., Geurts, P., Aerts, J., van den Oord, J., Atak, Z. K., Wouters, J., & Aerts, S. (2017). SCENIC: single-cell regulatory network inference and clustering. Nature methods , 14 (11), 1083–1086. https://doi.org/10.1038/nmeth.4463 Pozniak, J., Pedri, D., Landeloos, E., Van Herck, Y., Antoranz, A., Vanwynsberghe, L., Nowosad, A., Roda, N., Makhzami, S., Bervoets, G., Maciel, L. F., Pulido-Vicuña, C. A., Pollaris, L., Seurinck, R., Zhao, F., Flem-Karlsen, K., Damsky, W., Chen, L., Karagianni, D., Cinque, S., … Marine, J. C. (2024). A TCF4-dependent gene regulatory network confers resistance to immunotherapy in melanoma. Cell , 187 (1), 166–183.e25. https://doi.org/10.1016/j.cell.2023.11.037 Widmer, D. S., Cheng, P. F., Eichhoff, O. M., Belloni, B. C., Zipser, M. C., Schlegel, N. C., Javelaud, D., Mauviel, A., Dummer, R., & Hoek, K. S. (2012). Systematic classification of melanoma cells by phenotype-specific gene expression mapping. Pigment cell & melanoma research , 25 (3), 343–353. https://doi.org/10.1111/j.1755-148X.2012.00986.x Hugo, W., Zaretsky, J. M., Sun, L., Song, C., Moreno, B. H., Hu-Lieskovan, S., Berent-Maoz, B., Pang, J., Chmielowski, B., Cherry, G., Seja, E., Lomeli, S., Kong, X., Kelley, M. C., Sosman, J. A., Johnson, D. B., Ribas, A., & Lo, R. S. (2016). Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. Cell , 165 (1), 35–44. https://doi.org/10.1016/j.cell.2016.02.065 Riaz, N., Havel, J. J., Makarov, V., Desrichard, A., Urba, W. J., Sims, J. S., Hodi, F. S., Martín-Algarra, S., Mandal, R., Sharfman, W. H., Bhatia, S., Hwu, W. J., Gajewski, T. F., Slingluff, C. L., Jr, Chowell, D., Kendall, S. M., Chang, H., Shah, R., Kuo, F., Morris, L. G. T., … Chan, T. A. (2017). Tumor and Microenvironment Evolution during Immunotherapy with Nivolumab. Cell , 171 (4), 934–949.e16. https://doi.org/10.1016/j.cell.2017.09.028 Ayers, M., Lunceford, J., Nebozhyn, M., Murphy, E., Loboda, A., Kaufman, D. R., Albright, A., Cheng, J. D., Kang, S. P., Shankaran, V., Piha-Paul, S. A., Yearley, J., Seiwert, T. Y., Ribas, A., & McClanahan, T. K. (2017). IFN-γ-related mRNA profile predicts clinical response to PD-1 blockade. The Journal of clinical investigation , 127 (8), 2930–2940. https://doi.org/10.1172/JCI91190 Thompson, J. C., Davis, C., Deshpande, C., Hwang, W. T., Jeffries, S., Huang, A., Mitchell, T. C., Langer, C. J., & Albelda, S. M. (2020). Gene signature of antigen processing and presentation machinery predicts response to checkpoint blockade in non-small cell lung cancer (NSCLC) and melanoma. Journal for immunotherapy of cancer , 8 (2), e000974. https://doi.org/10.1136/jitc-2020-000974 Chow, A., Uddin, F. Z., Liu, M., Dobrin, A., Nabet, B. Y., Mangarin, L., Lavin, Y., Rizvi, H., Tischfield, S. E., Quintanal-Villalonga, A., Chan, J. M., Shah, N., Allaj, V., Manoj, P., Mattar, M., Meneses, M., Landau, R., Ward, M., Kulick, A., Kwong, C., … Rudin, C. M. (2023). The ectonucleotidase CD39 identifies tumor-reactive CD8 + T cells predictive of immune checkpoint blockade efficacy in human lung cancer. Immunity , 56 (1), 93–106.e6. https://doi.org/10.1016/j.immuni.2022.12.001 Li, X., Wan, Z., Liu, X., Ou, K., & Yang, L. (2022). A 12-chemokine gene signature is associated with the enhanced immunogram scores and is relevant for precision immunotherapy. Medical oncology (Northwood, London, England) , 39 (4), 43. https://doi.org/10.1007/s12032-021-01635-2 Vredevoogd, D. W., Kuilman, T., Ligtenberg, M. A., Boshuizen, J., Stecker, K. E., de Bruijn, B., Krijgsman, O., Huang, X., Kenski, J. C. N., Lacroix, R., Mezzadra, R., Gomez-Eerland, R., Yildiz, M., Dagidir, I., Apriamashvili, G., Zandhuis, N., van der Noort, V., Visser, N. L., Blank, C. U., Altelaar, M., … Peeper, D. S. (2019). Augmenting Immunotherapy Impact by Lowering Tumor TNF Cytotoxicity Threshold. Cell , 178 (3), 585–599.e15. https://doi.org/10.1016/j.cell.2019.06.014 Hoefsmit, E. P., van Royen, P. T., Rao, D., Stunnenberg, J. A., Dimitriadis, P., Lieftink, C., Morris, B., Rozeman, E. A., Reijers, I. L. M., Lacroix, R., Shehwana, H., Ligtenberg, M. A., Beijersbergen, R. L., Peeper, D. S., & Blank, C. U. (2023). Inhibitor of Apoptosis Proteins Antagonist Induces T-cell Proliferation after Cross-Presentation by Dendritic Cells. Cancer immunology research , 11 (4), 450–465. https://doi.org/10.1158/2326-6066.CIR-22-0494 Najem, A., Wouters, J., Krayem, M., Rambow, F., Sabbah, M., Sales, F., Awada, A., Aerts, S., Journe, F., Marine, J. C., & Ghanem, G. E. (2021). Tyrosine-Dependent Phenotype Switching Occurs Early in Many Primary Melanoma Cultures Limiting Their Translational Value. Frontiers in oncology , 11 , 780654. https://doi.org/10.3389/fonc.2021.780654 Shaffer, S. M., Dunagin, M. C., Torborg, S. R., Torre, E. A., Emert, B., Krepler, C., Beqiri, M., Sproesser, K., Brafford, P. A., Xiao, M., Eggan, E., Anastopoulos, I. N., Vargas-Garcia, C. A., Singh, A., Nathanson, K. L., Herlyn, M., & Raj, A. (2017). Rare cell variability and drug-induced reprogramming as a mode of cancer drug resistance. Nature , 546 (7658), 431–435. https://doi.org/10.1038/nature22794 Kemper, K., Krijgsman, O., Kong, X., Cornelissen-Steijger, P., Shahrabi, A., Weeber, F., van der Velden, D. L., Bleijerveld, O. B., Kuilman, T., Kluin, R. J. C., Sun, C., Voest, E. E., Ju, Y. S., Schumacher, T. N. M., Altelaar, A. F. M., McDermott, U., Adams, D. J., Blank, C. U., Haanen, J. B., & Peeper, D. S. (2016). BRAF(V600E) Kinase Domain Duplication Identified in Therapy-Refractory Melanoma Patient-Derived Xenografts. Cell reports , 16 (1), 263–277. https://doi.org/10.1016/j.celrep.2016.05.064 Titz, B., Lomova, A., Le, A., Hugo, W., Kong, X., Ten Hoeve, J., Friedman, M., Shi, H., Moriceau, G., Song, C., Hong, A., Atefi, M., Li, R., Komisopoulou, E., Ribas, A., Lo, R. S., & Graeber, T. G. (2016). JUN dependency in distinct early and late BRAF inhibition adaptation states of melanoma. Cell discovery , 2 , 16028. https://doi.org/10.1038/celldisc.2016.28 Van de Sande, B., Flerin, C., Davie, K., De Waegeneer, M., Hulselmans, G., Aibar, S., Seurinck, R., Saelens, W., Cannoodt, R., Rouchon, Q., Verbeiren, T., De Maeyer, D., Reumers, J., Saeys, Y., & Aerts, S. (2020). A scalable SCENIC workflow for single-cell gene regulatory network analysis. Nature protocols , 15 (7), 2247–2276. https://doi.org/10.1038/s41596-020-0336-2 Tumeh, P. C., Harview, C. L., Yearley, J. H., Shintaku, I. P., Taylor, E. J., Robert, L., Chmielowski, B., Spasic, M., Henry, G., Ciobanu, V., West, A. N., Carmona, M., Kivork, C., Seja, E., Cherry, G., Gutierrez, A. J., Grogan, T. R., Mateus, C., Tomasic, G., Glaspy, J. A., … Ribas, A. (2014). PD-1 blockade induces responses by inhibiting adaptive immune resistance. Nature , 515 (7528), 568–571. https://doi.org/10.1038/nature13954 Litchfield, K., Reading, J. L., Puttick, C., Thakkar, K., Abbosh, C., Bentham, R., Watkins, T. B. K., Rosenthal, R., Biswas, D., Rowan, A., Lim, E., Al Bakir, M., Turati, V., Guerra-Assunção, J. A., Conde, L., Furness, A. J. S., Saini, S. K., Hadrup, S. R., Herrero, J., Lee, S. H., … Swanton, C. (2021). Meta-analysis of tumor- and T cell-intrinsic mechanisms of sensitization to checkpoint inhibition. Cell , 184 (3), 596–614.e14. https://doi.org/10.1016/j.cell.2021.01.002 Spranger, S., Bao, R., & Gajewski, T. F. (2015). Melanoma-intrinsic β-catenin signalling prevents anti-tumour immunity. Nature , 523 (7559), 231–235. https://doi.org/10.1038/nature14404 Ji, R. R., Chasalow, S. D., Wang, L., Hamid, O., Schmidt, H., Cogswell, J., Alaparthy, S., Berman, D., Jure-Kunkel, M., Siemers, N. O., Jackson, J. R., & Shahabi, V. (2012). An immune-active tumor microenvironment favors clinical response to ipilimumab. Cancer immunology, immunotherapy: CII , 61 (7), 1019–1031. https://doi.org/10.1007/s00262-011-1172-6 Sha, D., Jin, Z., Budczies, J., Kluck, K., Stenzinger, A., & Sinicrope, F. A. (2020). Tumor Mutational Burden as a Predictive Biomarker in Solid Tumors. Cancer discovery , 10 (12), 1808–1825. https://doi.org/10.1158/2159-8290.CD-20-0522 Marin-Bejar, O., Rogiers, A., Dewaele, M., Femel, J., Karras, P., Pozniak, J., Bervoets, G., Van Raemdonck, N., Pedri, D., Swings, T., Demeulemeester, J., Borght, S. V., Lehnert, S., Bosisio, F., van den Oord, J. J., Bempt, I. V., Lambrechts, D., Voet, T., Bechter, O., Rizos, H., … Marine, J. C. (2021). Evolutionary predictability of genetic versus nongenetic resistance to anticancer drugs in melanoma. Cancer cell , 39 (8), 1135–1149.e8. https://doi.org/10.1016/j.ccell.2021.05.015 Lauss, M., Phung, B., Borch, T. H., Harbst, K., Kaminska, K., Ebbesson, A., Hedenfalk, I., Yuan, J., Nielsen, K., Ingvar, C., Carneiro, A., Isaksson, K., Pietras, K., Svane, I. M., Donia, M., & Jönsson, G. (2024). Molecular patterns of resistance to immune checkpoint blockade in melanoma. Nature communications , 15 (1), 3075. https://doi.org/10.1038/s41467-024-47425-y Boe, R. H., Triandafillou, C. G., Lazcano, R., Wargo, J. A., & Raj, A. (2024). Spatial transcriptomics reveals influence of microenvironment on intrinsic fates in melanoma therapy resistance. bioRxiv: the preprint server for biology , 2024.06.30.601416. https://doi.org/10.1101/2024.06.30.601416 Bai, X., Fisher, D. E., & Flaherty, K. T. (2019). Cell-state dynamics and therapeutic resistance in melanoma from the perspective of MITF and IFNγ pathways. Nature reviews. Clinical oncology , 16 (9), 549–562. https://doi.org/10.1038/s41571-019-0204-6 Karras, P., Black, J. R. M., McGranahan, N., & Marine, J. C. (2024). Decoding the interplay between genetic and non-genetic drivers of metastasis. Nature , 629 (8012), 543–554. https://doi.org/10.1038/s41586-024-07302-6 Wolchok, J. D., Chiarion-Sileni, V., Gonzalez, R., Rutkowski, P., Grob, J. J., Cowey, C. L., Lao, C. D., Wagstaff, J., Schadendorf, D., Ferrucci, P. F., Smylie, M., Dummer, R., Hill, A., Hogg, D., Haanen, J., Carlino, M. S., Bechter, O., Maio, M., Marquez-Rodas, I., Guidoboni, M., … Larkin, J. (2017). Overall Survival with Combined Nivolumab and Ipilimumab in Advanced Melanoma. The New England journal of medicine , 377 (14), 1345–1356. https://doi.org/10.1056/NEJMoa1709684 57. Larkin, J., Chiarion-Sileni, V., Gonzalez, R., Grob, J. J., Cowey, C. L., Lao, C. D., Schadendorf, D., Dummer, R., Smylie, M., Rutkowski, P., Ferrucci, P. F., Hill, A., Wagstaff, J., Carlino, M. S., Haanen, J. B., Maio, M., Marquez-Rodas, I., McArthur, G. A., Ascierto, P. A., Long, G. V., … Wolchok, J. D. (2015). Combined Nivolumab and Ipilimumab or Monotherapy in Untreated Melanoma. The New England journal of medicine , 373 (1), 23–34. https://doi.org/10.1056/NEJMoa1504030 Hodi, F. S., O'Day, S. J., McDermott, D. F., Weber, R. W., Sosman, J. A., Haanen, J. B., Gonzalez, R., Robert, C., Schadendorf, D., Hassel, J. C., Akerley, W., van den Eertwegh, A. J., Lutzky, J., Lorigan, P., Vaubel, J. M., Linette, G. P., Hogg, D., Ottensmeier, C. H., Lebbé, C., Peschel, C., … Urba, W. J. (2010). Improved survival with ipilimumab in patients with metastatic melanoma. The New England journal of medicine , 363 (8), 711–723. https://doi.org/10.1056/NEJMoa1003466 Sharma, P., Hu-Lieskovan, S., Wargo, J. A., & Ribas, A. (2017). Primary, Adaptive, and Acquired Resistance to Cancer Immunotherapy. Cell , 168 (4), 707–723. https://doi.org/10.1016/j.cell.2017.01.017 Kalbasi, A., & Ribas, A. (2020). Tumour-intrinsic resistance to immune checkpoint blockade. Nature reviews. Immunology , 20 (1), 25–39. https://doi.org/10.1038/s41577-019-0218-4 Versluis, J. M., Thommen, D. S., & Blank, C. U. (2020). Rationalizing the pathway to personalized neoadjuvant immunotherapy: the Lombard Street Approach. Journal for immunotherapy of cancer , 8 (2), e001352. https://doi.org/10.1136/jitc-2020-001352 Rozeman, E. A., Hoefsmit, E. P., Reijers, I. L. M., Saw, R. P. M., Versluis, J. M., Krijgsman, O., Dimitriadis, P., Sikorska, K., van de Wiel, B. A., Eriksson, H., Gonzalez, M., Torres Acosta, A., Grijpink-Ongering, L. G., Shannon, K., Haanen, J. B. A. G., Stretch, J., Ch'ng, S., Nieweg, O. E., Mallo, H. A., Adriaansz, S., … Blank, C. U. (2021). Survival and biomarker analyses from the OpACIN-neo and OpACIN neoadjuvant immunotherapy trials in stage III melanoma. Nature medicine , 27 (2), 256–263. https://doi.org/10.1038/s41591-020-01211-7 Zila, N., Eichhoff, O. M., Steiner, I., Mohr, T., Bileck, A., Cheng, P. F., Leitner, A., Gillet, L., Sajic, T., Goetze, S., Friedrich, B., Bortel, P., Strobl, J., Reitermaier, R., Hogan, S. A., Martínez Gómez, J. M., Staeger, R., Tuchmann, F., Peters, S., Stary, G., … Paulitschke, V. (2024). Proteomic Profiling of Advanced Melanoma Patients to Predict Therapeutic Response to Anti-PD-1 Therapy. Clinical cancer research : an official journal of the American Association for Cancer Research , 30 (1), 159–175. https://doi.org/10.1158/1078-0432.CCR-23-0562 Additional Declarations There is NO Competing Interest. Supplementary Files SUPPLEMENTARYFIGURES.docx Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-6506453","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Article","associatedPublications":[],"authors":[{"id":451837754,"identity":"96fbdfdb-c33f-49c7-903b-56017e82b68f","order_by":0,"name":"Daniel Peeper","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAuklEQVRIiWNgGAWjYDCCA1Can4ENxLYgQYtkA1iLBAlaDA6wgSgitPAd4D344OMOm3zj222JBxh3EKFF8gBfsuHMM2mW2+4cO3CA8QwRWgwO8JhJ87YdNjC7kd5wgLGNOC3mv/+2/TcwnkGCFjNmxrYDBgYSaQeI0yJ5mC9ZsvdMsoHEjbSEA4nEaOE73nvww88ddgb8M9KMP3xssyGshYGZh4GBsQHKSSBCAxAgaxkFo2AUjIJRgA0AANNyO0/CBgfTAAAAAElFTkSuQmCC","orcid":"https://orcid.org/0000-0003-1293-3177","institution":"Netherlands Cancer Institute","correspondingAuthor":true,"prefix":"","firstName":"Daniel","middleName":"","lastName":"Peeper","suffix":""},{"id":451837755,"identity":"7bfb21f9-7650-4f38-98d8-35cb137d26ca","order_by":1,"name":"Sebastiaan Schieven","email":"","orcid":"","institution":"Antoni van Leeuwenhoek Hospital","correspondingAuthor":false,"prefix":"","firstName":"Sebastiaan","middleName":"","lastName":"Schieven","suffix":""},{"id":451837756,"identity":"2ab72a41-c9d3-44ea-ab00-e3bbbf0b52eb","order_by":2,"name":"Joleen Traets","email":"","orcid":"","institution":"Netherlands Cancer Institute","correspondingAuthor":false,"prefix":"","firstName":"Joleen","middleName":"","lastName":"Traets","suffix":""},{"id":451837757,"identity":"b2b5f0e2-6d8e-4a0e-acfa-2fd67dd745c1","order_by":3,"name":"Arno Velds","email":"","orcid":"","institution":"NKI","correspondingAuthor":false,"prefix":"","firstName":"Arno","middleName":"","lastName":"Velds","suffix":""},{"id":451837758,"identity":"53a84a99-2136-4407-bd19-a600292e792e","order_by":4,"name":"Iris de Rink","email":"","orcid":"","institution":"The Netherlands Cancer Institute","correspondingAuthor":false,"prefix":"","firstName":"Iris","middleName":"","lastName":"de Rink","suffix":""},{"id":451837759,"identity":"71ce7de0-bf22-4ff3-9598-16761ad6ef33","order_by":5,"name":"Juan Simon Nieto","email":"","orcid":"","institution":"NKI","correspondingAuthor":false,"prefix":"","firstName":"Juan","middleName":"Simon","lastName":"Nieto","suffix":""},{"id":451837760,"identity":"3b9d2bdb-e7f1-4a16-8806-df1a2b4d1c60","order_by":6,"name":"Ji-Ying Song","email":"","orcid":"","institution":"Netherlands Cancer Institute","correspondingAuthor":false,"prefix":"","firstName":"Ji-Ying","middleName":"","lastName":"Song","suffix":""},{"id":451837761,"identity":"a0f4e5b4-34a6-4262-9b0b-88ea7a8909fd","order_by":7,"name":"Alex Vliet","email":"","orcid":"","institution":"NKI","correspondingAuthor":false,"prefix":"","firstName":"Alex","middleName":"","lastName":"Vliet","suffix":""},{"id":451837762,"identity":"39eed3f8-a008-445e-895e-68de97dba1d7","order_by":8,"name":"Austin George","email":"","orcid":"","institution":"NKI","correspondingAuthor":false,"prefix":"","firstName":"Austin","middleName":"","lastName":"George","suffix":""},{"id":451837763,"identity":"9396abb4-34f3-4a39-9515-0656ca2289d3","order_by":9,"name":"Marja Nieuwland","email":"","orcid":"","institution":"NKI","correspondingAuthor":false,"prefix":"","firstName":"Marja","middleName":"","lastName":"Nieuwland","suffix":""},{"id":451837764,"identity":"e70df368-6208-4edd-8ddb-e99d2a5974f7","order_by":10,"name":"Ingrid Hofland","email":"","orcid":"","institution":"NKI","correspondingAuthor":false,"prefix":"","firstName":"Ingrid","middleName":"","lastName":"Hofland","suffix":""},{"id":451837765,"identity":"94692a35-29f7-4ec8-b763-e4e562b66255","order_by":11,"name":"Lex Vrije","email":"","orcid":"","institution":"NKI","correspondingAuthor":false,"prefix":"","firstName":"Lex","middleName":"","lastName":"Vrije","suffix":""},{"id":451837766,"identity":"3338ddcf-44fa-4a89-ab88-88e2d7698e4b","order_by":12,"name":"Stephanie Blankenstein","email":"","orcid":"","institution":"NKI","correspondingAuthor":false,"prefix":"","firstName":"Stephanie","middleName":"","lastName":"Blankenstein","suffix":""},{"id":451837767,"identity":"08d9e129-4ddc-4fbd-b219-3c888d9cb7b3","order_by":13,"name":"Julia Boshuizen","email":"","orcid":"","institution":"NKI","correspondingAuthor":false,"prefix":"","firstName":"Julia","middleName":"","lastName":"Boshuizen","suffix":""},{"id":451837768,"identity":"5dcd2aaf-1a59-46d0-bf09-24e21251dd3b","order_by":14,"name":"Marcos Da Silva Guimaraes","email":"","orcid":"","institution":"NKI","correspondingAuthor":false,"prefix":"","firstName":"Marcos","middleName":"Da Silva","lastName":"Guimaraes","suffix":""},{"id":451837769,"identity":"a9814a51-3de6-4140-a35a-d7b9c9417806","order_by":15,"name":"Hugo Horlings","email":"","orcid":"https://orcid.org/0000-0003-4782-8828","institution":"The Netherlands Cancer Institute","correspondingAuthor":false,"prefix":"","firstName":"Hugo","middleName":"","lastName":"Horlings","suffix":""},{"id":451837770,"identity":"a188f63c-37fc-404d-b05b-690ccf1df749","order_by":16,"name":"Amalie Dick","email":"","orcid":"","institution":"NKI","correspondingAuthor":false,"prefix":"","firstName":"Amalie","middleName":"","lastName":"Dick","suffix":""},{"id":451837771,"identity":"d3edd7da-3999-41f2-ab7f-74397609848c","order_by":17,"name":"Martijn van Baalen","email":"","orcid":"","institution":"NKI","correspondingAuthor":false,"prefix":"","firstName":"Martijn","middleName":"van","lastName":"Baalen","suffix":""}],"badges":[],"createdAt":"2025-04-22 17:26:43","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-6506453/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-6506453/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":82352055,"identity":"1125e4fe-052f-47f7-a3db-f01c950b991a","added_by":"auto","created_at":"2025-05-09 11:00:35","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":1314288,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eNGFR-expressing melanomas are highly heterogeneous for differentiation and dedifferentiation markers.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e(A)\u0026nbsp; MITF, Melan-A, AXL and NGFR IHC scores (percentage positive of total viable tumor cells) per patient. The maximum score per staining is a 100%. #: denotes patients 35, 37, 44 and 45.\u003c/p\u003e\n\u003cp\u003e(B)\u0026nbsp; Correlation analysis between percentage MITF and NGFR positive tumor cells across patient melanomas. Patient tumor numbers are depicted. Patients 35, 37, 44 and 45 are shown in grey, magenta, purple and cyan, respectively. Coefficient of determination (R\u003csup\u003e2\u003c/sup\u003e), \u003cem\u003eP\u003c/em\u003e value and Pearson correlation are depicted.\u003c/p\u003e\n\u003cp\u003e(C)\u0026nbsp; Correlation analysis between percentage Melan-A and NGFR positive tumor cells across patient melanomas. Patient tumor numbers are depicted. Patients 35, 37, 44 and 45 are shown in grey, magenta, purple and cyan, respectively.\u003cem\u003e \u003c/em\u003eCoefficient of determination (R\u003csup\u003e2\u003c/sup\u003e), \u003cem\u003eP\u003c/em\u003e value and Pearson correlation are depicted.\u003c/p\u003e\n\u003cp\u003e(D)\u0026nbsp; Correlation analysis between percentage AXL and NGFR positive tumor cells across patient melanomas. Patient tumor numbers are depicted. Patients 35, 37, 44 and 45 are shown in grey, magenta, purple and cyan, respectively. Coefficient of determination (R\u003csup\u003e2\u003c/sup\u003e), \u003cem\u003eP\u003c/em\u003e value and Pearson correlation are depicted.\u003c/p\u003e\n\u003cp\u003e(E)\u0026nbsp; Dedifferentiation trajectory is illustrated, beginning from a differentiated cell state characterized by MITF and Melan-A, transitioning to a NCSC state marked by NGFR, and ultimately reaching an undifferentiated state marked by AXL. Possible NGFR subphenotypes are also depicted.\u003c/p\u003e\n\u003cp\u003e(F)\u0026nbsp; NGFR, AXL, MITF and Melan-A IHC staining’s on patient tumor 35. Bar is 200 mm.\u003c/p\u003e\n\u003cp\u003e(G) NGFR, AXL, MITF and Melan-A IHC staining’s on patient tumor 37. Bar is 200 mm.\u003c/p\u003e\n\u003cp\u003e(H)\u0026nbsp; NGFR, AXL, MITF and Melan-A IHC staining’s on patient tumor 44. Bar is 500 mm.\u003c/p\u003e\n\u003cp\u003e(I)\u0026nbsp;\u0026nbsp;\u0026nbsp; NGFR, AXL, MITF and Melan-A IHC staining’s on patient tumor 45. Bar is 200 mm.\u003c/p\u003e\n\u003cp\u003e(J)\u0026nbsp;\u0026nbsp; Schematic overview of scRNA-seq and Visium 10X spatial transcriptomics pipelines employed on melanomas from patients 35, 37, 44 and 45.\u0026nbsp;\u003c/p\u003e","description":"","filename":"floatimage1.png","url":"https://assets-eu.researchsquare.com/files/rs-6506453/v1/4f28bcda4053a5758e0b73eb.png"},{"id":82353824,"identity":"60961099-d81a-4504-8e28-e13f65fba825","added_by":"auto","created_at":"2025-05-09 11:08:35","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":1080036,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eScRNA-seq identifies several transcriptomically distinct NGFR\u003c/strong\u003e\u003csup\u003e\u003cstrong\u003e+\u003c/strong\u003e\u003c/sup\u003e\u003cstrong\u003e melanoma phenotypes.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e(A)\u0026nbsp; CD45\u003csup\u003e-\u003c/sup\u003e/NGFR\u003csup\u003e+\u003c/sup\u003e sorting strategy of patients 35, 37, 44 and 45. Percentage CD45\u003csup\u003e-\u003c/sup\u003e/NGFR\u003csup\u003e+\u003c/sup\u003e cells from all CD45\u003csup\u003e-\u003c/sup\u003e/DAPI\u003csup\u003e- \u003c/sup\u003ecells is depicted per patient and FACS sorted for subsequent scRNA-seq.\u003c/p\u003e\n\u003cp\u003e(B)\u0026nbsp; Corresponding NGFR MFI (CD45\u003csup\u003e-\u003c/sup\u003e/DAPI\u003csup\u003e- \u003c/sup\u003efraction) of the four melanoma patients depicted in panel A.\u003c/p\u003e\n\u003cp\u003e(C)\u0026nbsp; UMAP of patients 35, 37, 44 and 45 NGFR\u003csup\u003e+\u003c/sup\u003e/CD45\u003csup\u003e-\u003c/sup\u003e melanoma cells. Only malignant cells, based on the Jerby-Arnon malignant signature\u003csup\u003e29\u003c/sup\u003e, are plotted.\u003c/p\u003e\n\u003cp\u003e(D)\u0026nbsp; Heatmap showing Z-scores of different reported melanoma signatures and phenotype switching marker genes across all identified scRNA-seq clusters of patients 35, 37, 44 and 45. Row scaling was performed per patient. For the bottom gene panel, row scaled expression values of depicted genes are shown per patient. Cluster numbers per patient are depicted at the bottom. Grey colored boxes: gene not detected.\u003c/p\u003e","description":"","filename":"floatimage2.png","url":"https://assets-eu.researchsquare.com/files/rs-6506453/v1/7150cef25466600ce853973f.png"},{"id":82356005,"identity":"bd484f40-ff84-4860-a561-669aab66789c","added_by":"auto","created_at":"2025-05-09 11:16:35","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":944917,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eNGFR\u003c/strong\u003e\u003csup\u003e\u003cstrong\u003e+\u003c/strong\u003e\u003c/sup\u003e\u003cstrong\u003e melanoma cells are characterized by patient-specific gene-regulatory networks.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e(A)\u0026nbsp; Regulon UMAP clusters are depicted for patients 35, 37, 44 and 45. Clustering was done based on differential regulon activity between the different identified clusters for each patient.\u003c/p\u003e\n\u003cp\u003e(B)\u0026nbsp; Heatmaps of patients 35, 37, 44 and 45 showing top three significant regulons identified per regulon cluster. Z-score of the differential regulon activity across clusters per patient is plotted.\u003c/p\u003e\n\u003cp\u003e(C)\u0026nbsp; Venn diagram showing the overlap of significantly differentially active regulons identified in the four patients.\u003c/p\u003e\n\u003cp\u003e(D)\u0026nbsp; Heatmap showing average reciprocal overlap of differentially active regulons for all four patients. Numbers of clusters corresponding to panel (A) are depicted.\u003c/p\u003e","description":"","filename":"floatimage3.png","url":"https://assets-eu.researchsquare.com/files/rs-6506453/v1/7264109bfff93a8f21e3bcca.png"},{"id":82353825,"identity":"0558ff6f-71cf-4575-82f8-4b1eb60ab235","added_by":"auto","created_at":"2025-05-09 11:08:35","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":1331426,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eSpatial organization of NGFR\u003c/strong\u003e\u003csup\u003e\u003cstrong\u003e+\u003c/strong\u003e\u003c/sup\u003e\u003cstrong\u003e subpopulations.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e(A)\u0026nbsp; NGFR IHC stainings of patients 35, 37, 44 and 45 (same as in S1B-E). Visium 10X quadrant area selection per patient is depicted. Bar is 5 mm for patients 35 and 37, 10 mm for patient 44 and 2 mm for patient 45.\u003c/p\u003e\n\u003cp\u003e(B)\u0026nbsp; NGFR staining on the same tumor region as the selected quadrant area in panel A for each Visium 10X experiment. Bar is 2 mm.\u003c/p\u003e\n\u003cp\u003e(C)\u0026nbsp; Normalized NGFR expression on the Visium 10X quadrant slices.\u003c/p\u003e\n\u003cp\u003e(D)\u0026nbsp; Seurat Louvain clustering using Visium 10X transcriptomic data as input. Clusters are depicted in the Visium images.\u003c/p\u003e\n\u003cp\u003e(E)\u0026nbsp; Heatmap showing Z-scores of different reported melanoma signatures and phenotype switching marker genes across all identified Visium gene expression clusters of patients 35, 37, 44 and 45. Row scaling was performed per patient. For the bottom gene panel, row scaled expression values of depicted genes are shown per patient. Cluster numbers per tumor are depicted at the bottom of the heatmaps.\u003c/p\u003e\n\u003cp\u003e(F)\u0026nbsp; Spatial location of scRNA-seq regulon clusters within the tumor tissue. ScRNA-seq regulon clusters that showed the most significant overlap with Visium regulon clusters (shown in Figure S4C) are depicted per patient.\u0026nbsp;\u003c/p\u003e","description":"","filename":"floatimage4.png","url":"https://assets-eu.researchsquare.com/files/rs-6506453/v1/a6920cbb276715089b832784.png"},{"id":82352058,"identity":"211d3bc4-0e0e-4d26-949c-7050a7b2b8e5","added_by":"auto","created_at":"2025-05-09 11:00:35","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":1071266,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003ePDGFR marks cytokine-resistant and mesenchymal-like NGFR\u003c/strong\u003e\u003csup\u003e\u003cstrong\u003e+\u003c/strong\u003e\u003c/sup\u003e\u003cstrong\u003e melanoma cells\u003c/strong\u003e.\u003c/p\u003e\n\u003cp\u003e(A)\u0026nbsp; Western blot analysis of different melanoma NGFR\u003csup\u003e+ \u003c/sup\u003ecell lines. Protein expression of PDGFRa, PDGFRb, AXL, NGFR and Melan-A was determined. Tubulin was used as loading control. The markings indicate the closest molecular weight marker. Two to three biological replicates were performed.\u003c/p\u003e\n\u003cp\u003e(B)\u0026nbsp; Quantification of cytotoxicity assay of IFNg and TNF on a panel of NGFR\u003csup\u003e+ \u003c/sup\u003emelanoma lines. Cells were treated for five days with 25ng/ml IFNg and 25ng/ml TNF before crystal violet staining. The relative viability to untreated control and to D10 cells is shown. Mean with SD of pooled technical replicates is shown. Two to five biological replicates were performed each with two or three technical replicates. Statistical test is One-way ANOVA (\u003cem\u003e***, P\u003c/em\u003e\u0026lt;0.001\u003cem\u003e (only the comparison of A875 with M026.X1.CL); ****, P\u003c/em\u003e\u0026lt;0.0001)\u003cem\u003e.\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003e(C)\u0026nbsp; Colony formation assay of two NGFR\u003csup\u003e+\u003c/sup\u003e/PDGFR\u003csup\u003e-\u003c/sup\u003e and two NGFR\u003csup\u003e+\u003c/sup\u003e/PDGFR\u003csup\u003e+\u003c/sup\u003e cell lines exposed to 0 or 25ng/ml IFNg and TNF. Data matched with figure 5B.\u003c/p\u003e\n\u003cp\u003e(D)\u0026nbsp; NGFR/PAS and PDGFRb/PAS stainings on patient tumors 7b, 36, 43 and 45. Pictures were made from the same tumor region. PAS staining (in pink) marks the vessel/stromal component. Zoom-ins on PDGFRb/PAS stainings are also depicted. Bar is 50 mm.\u003c/p\u003e\n\u003cp\u003e(E)\u0026nbsp; Z-scores for \u003cem\u003eNGFR\u003c/em\u003e, \u003cem\u003ePDGFRA\u003c/em\u003e, \u003cem\u003ePDGFRB\u003c/em\u003e genes across the different Pozniak melanoma phenotypes are depicted\u003csup\u003e31\u003c/sup\u003e. Only the malignant BT (i.e., before treatment) compartment was analyzed (n=19).\u0026nbsp;\u003c/p\u003e","description":"","filename":"floatimage6.png","url":"https://assets-eu.researchsquare.com/files/rs-6506453/v1/a25363329426368b94f5a0f6.png"},{"id":82352061,"identity":"43e2ab37-a35e-44db-b434-e5bad2629d58","added_by":"auto","created_at":"2025-05-09 11:00:35","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":873500,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cem\u003e\u003cstrong\u003ePDGFR\u003c/strong\u003e\u003c/em\u003e\u003cstrong\u003e marks immune-active yet ICB non-responding melanomas.\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e(A)\u0026nbsp; Heatmap of combined Hugo and Riaz datasets showing row Z-scores for genes and signatures per sample. Patient melanoma samples are annotated for dataset, response, calculated immune activity score (IAS) and \u003cem\u003eNGFR\u003c/em\u003e and \u003cem\u003ePDGFR \u003c/em\u003eexpression.\u003c/p\u003e\n\u003cp\u003e(B)\u0026nbsp; IAS score across the different IAS patient groups is shown. Data of the combined RNA-seq dataset was used. The upper and lower quantile are shown in the boxplots with the median indicated by the center line. The whiskers show the 1.5 interquartile range. All datapoints are shown. Statistical test is Wilcoxon rank test. \u003cem\u003eP\u003c/em\u003e values are depicted.\u003c/p\u003e\n\u003cp\u003e(C)\u0026nbsp; The average Z-score of \u003cem\u003ePDGFR \u003c/em\u003eis depicted across the different IAS patient groups. Data of the three combined RNA-seq dataset was used. The upper and lower quantile are shown in the boxplots with the median indicated by the center line. The whiskers show the 1.5 interquartile range. All datapoints are indicated. Statistical test is Wilcoxon rank test. \u003cem\u003eP\u003c/em\u003e values are depicted.\u003c/p\u003e\n\u003cp\u003e(D)\u0026nbsp; The average Z-score of the mesenchymal-like signature is depicted across the different IAS patient groups. Data of the combined RNA-seq dataset was used. The upper and lower quantile are shown in the boxplots with the median indicated by the center line. The whiskers show the 1.5 interquartile range. All datapoints are shown. Statistical test is Wilcoxon rank test. \u003cem\u003eP\u003c/em\u003e values are depicted.\u003c/p\u003e\n\u003cp\u003e(E)\u0026nbsp; CD3 score for PDGFRb\u003csup\u003e+\u003c/sup\u003e and PDGFRb\u003csup\u003e-\u003c/sup\u003e human melanoma samples. Mean with SD is shown (n=10 and n=23). Statistical test is unpaired \u003cem\u003eT\u003c/em\u003e test. \u003cem\u003eP\u003c/em\u003e value is depicted.\u003c/p\u003e\n\u003cp\u003e(F)\u0026nbsp; PDGFRb/PAS and CD3/PAS IHC staining’s of a PDGFRb\u003csup\u003e+\u003c/sup\u003e human melanoma. Zoom-ins are shown for PDGFRb\u003csup\u003eHigh\u003c/sup\u003e and PDGFRb\u003csup\u003eLow\u003c/sup\u003e regions. Bar indicates 100 mm.\u003c/p\u003e","description":"","filename":"floatimage8.png","url":"https://assets-eu.researchsquare.com/files/rs-6506453/v1/6179dfa1936319acd48feb81.png"},{"id":84575895,"identity":"2bd34283-263d-4b98-ac6a-32d02d3d75fa","added_by":"auto","created_at":"2025-06-13 16:49:23","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":7808288,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-6506453/v1/67fe7ce0-4602-4f8c-80fa-46223b74fe45.pdf"},{"id":82356008,"identity":"3d2283b5-cbd5-49c0-97b6-bea1742ffa6b","added_by":"auto","created_at":"2025-05-09 11:16:35","extension":"docx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":7112630,"visible":true,"origin":"","legend":"","description":"","filename":"SUPPLEMENTARYFIGURES.docx","url":"https://assets-eu.researchsquare.com/files/rs-6506453/v1/950f51a2ce8299734dc57ab6.docx"}],"financialInterests":"There is \u003cb\u003eNO\u003c/b\u003e Competing Interest.","formattedTitle":"\u003cp\u003eIntratumoral Heterogeneity in Ngfr\u003csup\u003e+\u003c/sup\u003e Melanoma Subpopulations Shapes Immune Evasion and Immunotherapy Resistance\u003c/p\u003e","fulltext":[{"header":"INTRODUCTION","content":"\u003cp\u003eMelanoma is a heterogeneous cancer type, in part owing to both its high mutational load and the presence of multiple transcriptional phenotypes\u003csup\u003e\u003cspan additionalcitationids=\"CR2\" citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e\u003c/sup\u003e. Originally, melanomas were classified as either proliferative or invasive\u003csup\u003e\u003cspan citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e,\u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u003c/sup\u003e, mostly as a function of the expression levels of the key melanocyte transcription factor MITF and its target genes. Later, this model was refined to include a higher phenotypic complexity at the transcriptomic level\u003csup\u003e\u003cspan additionalcitationids=\"CR7\" citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e\u003c/sup\u003e. The different melanoma phenotypes can transition from one to the other, for example upon inhibition of the MAPK pathway\u003csup\u003e\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e\u003c/sup\u003e. We and others have shown that such transitions are associated with cell state-specific gene signatures and marker genes, for example AXL and NGFR\u003csup\u003e\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e,\u003cspan additionalcitationids=\"CR10\" citationid=\"CR9\" class=\"CitationRef\"\u003e9\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e\u003c/sup\u003e. NGFR, a cell surface neurotrophic factor receptor marking neural crest stem cells (NCSCs)\u003csup\u003e\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e,\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e,\u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e13\u003c/span\u003e\u003c/sup\u003e stands out, as it can be co-expressed with either differentiation or dedifferentiation markers in melanoma cells\u003csup\u003e\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e,\u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e12\u003c/span\u003e,\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e,\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e\u003c/sup\u003e. Its expression can be induced through various mechanisms, including MAPK pathway inhibition\u003csup\u003e\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e,\u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e16\u003c/span\u003e,\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e\u003c/sup\u003e, TGF-β\u003csup\u003e\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e\u003c/sup\u003e, \u003cem\u003ein vitro\u003c/em\u003e T cell challenge\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e, adoptive cell transfer (ACT)\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e,\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e\u003c/sup\u003e, immunotherapy\u003csup\u003e\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e\u003c/sup\u003e and immune cytokines such as TNF\u003csup\u003e\u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e14\u003c/span\u003e,\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e,\u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e21\u003c/span\u003e\u003c/sup\u003e and IFNg\u003csup\u003e\u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e22\u003c/span\u003e,\u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e23\u003c/span\u003e\u003c/sup\u003e. This complex regulation of NGFR conceivably contributes to its heterogeneous expression pattern.\u003c/p\u003e \u003cp\u003eWe previously showed that NGFR\u003csup\u003e+\u003c/sup\u003e subpopulations commonly pre-exist in patient melanomas already before any treatment, and that they typically manifest in heterogeneous patterns\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e. Furthermore, we demonstrated that NGFR marks melanoma cells that are relatively insensitive to cytotoxic T cells, as well as their cytokines IFNg and TNF, and can be associated with T cell exclusion\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e. In line with these observations, others have revealed NGFR's involvement in rendering melanoma cells resistant to natural killer (NK) cells\u003csup\u003e\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e\u003c/sup\u003e. Additionally, in a patient case study, the proportion of NGFR-expressing cells increased following immune checkpoint blockade (ICB)\u003csup\u003e\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e\u003c/sup\u003e. Given that a substantial fraction of melanoma patients exhibits either intrinsic or acquired resistance to immunotherapy, potentially attributed to NGFR expression, we suggested previously that it may be beneficial to target NGFR to overcome resistance\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eThese observations together suggest that the NGFR\u003csup\u003e+\u003c/sup\u003e melanoma subpopulation can manifest in heterogeneous patterns, express different marker proteins, and is induced by distinct stimuli. However, it is unclear whether it reflects an intrinsically homogeneous cell group or that even this subpopulation is heterogeneous in nature. Clearly, this would have immediate consequences for responses to immune cells and therapy. Therefore, we dissected NGFR heterogeneity from clinical melanoma specimens using different technologies including single-cell RNA-sequencing (scRNA-seq) and spatial transcriptomics, in combination with functional analyses. Specifically, we investigated whether (i) this heterogeneity manifests only between cell groups or that even single cells can co-express different markers; (ii) this pattern of heterogeneity is consistent among patients; (iii) heterogeneous NGFR\u003csup\u003e+\u003c/sup\u003e cells are spatially organized; (iv) heterogeneous NGFR\u003csup\u003e+\u003c/sup\u003e subpopulations show differential immune responses; (v) these responses are reflected by differential clinical responses to immunotherapy.\u003c/p\u003e"},{"header":"MATERIAL AND METHODS","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003eCell lines and cell culture conditions\u003c/h2\u003e \u003cp\u003eAll melanoma cell lines are from the Peeper laboratory cell line stock. Melanoma cell lines were cultured in DMEM (Gibco) with fetal bovine serum (FBS, Sigma), 100 U/ml penicillin and 0.1 mg/ml streptomycin (both Gibco) under standard conditions. All cell lines were authenticated by short tandem repeat (STR) profiling (Promega) and regularly confirmed to be mycoplasma-free by PCR.\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eWestern blotting and antibodies\u003c/h3\u003e\n\u003cdiv class=\"Heading\"\u003eWestern blotting and antibodies\u003c/div\u003e \u003cp\u003eAs described previously\u003csup\u003e\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e\u003c/sup\u003e. Cell pellets were lysed in RIPA buffer (50 mM TRIS pH 8, 150 mM NaCl, 1% Nonidet P40, 0.5% sodium deoxycholate, 0.1% SDS) supplemented with HALT\u0026trade; protease and phosphatase inhibitor cocktail (100x) (Fisher Scientific, cat # 78444). Lysis was performed for 30 minutes and vortexed every 10 minutes. Samples were centrifuged at maximum speed for 10 minutes. Protein concentration was determined using a Bradford assay (Bio-Rad). Protein concentrations were normalized to each other and mixed with 4x LDS sample buffer (Fisher Scientific, 15484379) containing 10% β-Mercaptoethanol (Merck) (final concentration 2.5%) and incubated for 5 minutes at 95\u0026deg;C. Western blotting was performed with standard techniques using 4\u0026ndash;12% Bis-Tris polyacrylamide-SDS gels (NuPAGE, Life Technologies) and nitrocellulose membranes (Whatman, GE Healthcare). Blotting was performed using the iBlot dry blotting system from Invitrogen. Blots were blocked in 4% milk in PBS plus 0.2% Tween 100 and incubated with primary antibody: AXL (1:1,000, C89E7, Cell Signaling), Melan-A (1:10,000, M7196, Dako); PDGFRb (1:1,000, 3169, CST), PDGFRa (1:1,000, D1E1E, Cell Signaling), Tubulin (1:10,000, T9026, Sigma), and NGFR (1:1,000, #8238, Cell Signaling). The following secondary antibodies were used: goat anti-rabbit peroxidase conjugate (1:7,500, G21234) and goat anti-mouse (1:7,500, G21040), both purchased from Invitrogen. Immunoblots were incubated with Clarity\u0026trade; Western ECL Substrate (cat# 170\u0026ndash;5061, Biorad). Luminescence was captured by the Bio-Rad ChemiDoc imaging system. Both 8-bit tiff and 16-bit raw tiff images were used to make the figures.\u003c/p\u003e\n\u003ch3\u003eIHC of human melanomas\u003c/h3\u003e\n\u003cp\u003e \u003cb\u003e\u003c/b\u003e As described previously\u003csup\u003e\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e\u003c/sup\u003e. The collection and use of human tissue was approved by the Medical Ethical Review Board of the Antoni van Leeuwenhoek. Patients gave informed consent for secondary use of tumor tissue. The study was approved by the Institutional Review Board (IRB).\u003c/p\u003e \u003cp\u003eImmunohistochemistry of the FFPE tumor samples was performed on a BenchMark Ultra autostainer (Ventana Medical Systems). Briefly, paraffin sections were cut at 3 \u0026micro;m, heated at 75\u0026deg;C for 28 minutes and deparaffinized in the instrument with EZ prep solution (Ventana Medical Systems). Heat-induced antigen retrieval was carried out using Cell Conditioning 1 (CC1, Ventana Medical Systems) for 32 minutes at 95\u0026deg;C (MITF) and 40 minutes at 95\u0026deg;C (AXL). MITF was detected using clone C5/D5 (1:800 dilution, 32 minutes at 37\u0026deg;C., LSBio), AXL using clone C89E7 (1:100 dilution, 32 minutes at room temperature., Cell Signaling). For AXL, signal amplification was applied using the Optiview Amplification Kit (8 minutes, Ventana Medical Systems). Melan-A was detected using clone A103 (1:40 dilution, 32 minutes at 37\u0026deg;C, Agilent/Dako). NGFR was detected using clone D4B3 (1:400 dilution, 1 hour at room temperature, Cell signaling). For Melan-A, signal amplification was applied using the Optiview Amplification Kit (4 minutes, Ventana Medical Systems). \u003cspan fontcategory=\"NonProportional\" class=\"\" name=\"Emphasis\"\u003e\u003c/span\u003e CD3 was stained as described previously (RM-9107-S, Thermo Scientific)\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e. PDGFRb was detected using clone 28E1 (1:50 dilution, 60 minutes at room temperature, Cell Signaling). Stainings were developed using brown (OptiView DAB Detection Kit (Ventana Medical Systems) or red (UltraView Universal Alkaline Phosphatase Red Detection (Roche Diagnostics, Ventana)) visualization. Slides were counterstained with Hematoxylin and Bluing Reagent (Ventana Medical Systems).\u003c/p\u003e \u003cp\u003eFor a sequential Periodic Acid\u0026ndash;Schiff (PAS) staining, slides were removed from the BenchMark Ultra autostainer, rinsed in distilled water, incubated in 0.5% Periodic Acid (7 minutes, VWR) followed by Schiff\u0026rsquo;s Reagent (30 minutes, CellaVision / RAL Diagnostics) after washing steps in between. Slides were counterstained with Hematoxylin (KliniPath). A PANNORAMIC\u0026reg; 1000 scanner from 3DHISTECH was used to scan the slides at a 40x magnification.\u003c/p\u003e \u003cp\u003ePercentage positive tumor cells was quantified. Cytoplasmic and membrane staining were scored for AXL, NGFR, PDGFRβ and Melan-A. Only nuclear MITF staining was scored. Scoring was performed by a certified pathologist. A semi-quantitative scoring was applied to quantify CD3 infiltration. Tumor regions were scored as: not infiltrated (score of 0), lowly infiltrated (a score of 2), moderately infiltrated (a score of 4) and highly infiltrated (a score of 6). The average of all scores was taken if tumors received more than one score. Only intra-tumoral CD3\u003csup\u003e+\u003c/sup\u003e T cells were scored. The majority of the tumors studied were untreated/baseline and obtained from surgical resections or \u003cspan fontcategory=\"NonProportional\" class=\"\" name=\"Emphasis\"\u003e\u003c/span\u003e using a 14-gauge biopsy needle. One tumor was an anti-PD1 treatment relapse sample (patient tumor 44 in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA and \u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003eS1\u003c/span\u003eA).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e\n\u003ch3\u003eFlow cytometry\u003c/h3\u003e\n\u003cp\u003eMethod for cell surface staining as described previously\u003csup\u003e\u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e25\u003c/span\u003e\u003c/sup\u003e. Cells were stained with antibodies targeting surface molecules of interest according to manufacturer\u0026rsquo;s instructions and analyzed on a Fortessa flow cytometer LSR ( \u003cspan fontcategory=\"NonProportional\" class=\"\" name=\"Emphasis\"\u003e\u003c/span\u003e BD Bioscience). The following antibodies were used: AXL-PE conjugated antibody (1:200, FAB154P, R\u0026amp;D), PDGFRb-PE conjugated antibody (1:50, FAB1263P, R\u0026amp;D) and \u003cspan fontcategory=\"NonProportional\" class=\"\" name=\"Emphasis\"\u003e\u003c/span\u003e NGFR-APC (1:200, 345107, Biolegend) for 20 minutes at 4\u0026deg;C. Data depicted is mean fluorescence intensity (MFI) of sample - MFI of unstained sample, unless otherwise stated.\u003c/p\u003e\n\u003ch3\u003eColony formation assay\u003c/h3\u003e\n\u003cp\u003e100,000 cells per well were seeded in a 12 wells plate. The day after, cells were challenged with cytokines IFNg (300-02, Prepotech) and TNF (11343017, Immunotools). After 5 days, medium was removed and cells were \u003cspan fontcategory=\"NonProportional\" class=\"\" name=\"Emphasis\"\u003e\u003c/span\u003e stained with a crystal violet solution containing 0.1% crystal violet (Sigma) and 50% methanol (Honeywell) for 1 hour. \u003cspan fontcategory=\"NonProportional\" class=\"\" name=\"Emphasis\"\u003e\u003c/span\u003e For quantification, the crystal violet stain was dissolved in 10% acetic acid (Sigma). Absorbance of this solution was measured on an Infinite 200 Pro spectrophotometer (Tecan) at 595 nm.\u003c/p\u003e \u003cdiv id=\"Sec8\" class=\"Section2\"\u003e \u003ch2\u003eSingle-cell RNA-sequencing sample preparation\u003c/h2\u003e \u003cp\u003eCollected tumors from patients 35 (lymph node lesion), 37 (lymph node lesion), 44 (thigh lesion within fat tissue) and 45 (skin lesion, Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA) were digested in FBS free RPMI medium (Gibco) with 100 U/ml penicillin/0.1 mg/ml streptomycin, collagenase IV (1:50, 17104-019, ThermoFisher Scientific) and pulmozyme (1 mg/ml, Roche) for 30 minutes at 37\u0026deg;C on a rotator. Digest was filtered over a 100 \u0026micro;m filter and frozen down in 90% FBS and 10% DMSO (45-34943, Sigma-Aldrich). Next, digests were thawed in 3 ml cold RPMI medium with 10% FBS and 4 \u0026micro;l benzonase (1:1,000, 70746-3, VWR). Cells were spun for 8 minutes 14,000 rpm at 4\u0026deg;C. Supernatant was removed, and cells were resuspended in 4 ml cold RPMI with 4 \u0026micro;l benzonase (1:1,000). Samples were centrifuged for 8 minutes 14,000 rpm at 4\u0026deg;C followed by resuspension in cold PBS with 1% BSA. The following antibodies were used for staining: NGFR-APC (1:200, \u003cspan fontcategory=\"NonProportional\" class=\"\" name=\"Emphasis\"\u003e\u003c/span\u003e 345107, Biolegend) and CD45-Alexa Fluor 488 (1:400, 304019, Biolegend). Cells were incubated in staining solution (i.e., antibodies diluted in PBS\u0026thinsp;+\u0026thinsp;1% BSA) for 20\u0026ndash;30 minutes. Then, cells were washed twice with cold PBS with 1% BSA and centrifuged (8 minutes 14,000 rpm at 4\u0026deg;C). Cells were kept in PBS with 2% of FBS before sorting for NGFR\u003csup\u003e+\u003c/sup\u003e and CD45\u003csup\u003e\u0026minus;\u003c/sup\u003e cells. Cells were sorted by BD FACSAria Fusion. DAPI was used as a life dead marker. Cells were considered NGFR positive if they had a NGFR signal higher than unstained and/or fluorescence minus one (FMO) control. CD45\u003csup\u003e\u0026minus;\u003c/sup\u003e/NGFR\u003csup\u003e+\u003c/sup\u003e cells were used as input for scRNA-seq.\u003c/p\u003e \u003c/div\u003e\n\u003ch3\u003eSingle-cell RNA-sequencing\u003c/h3\u003e\n\u003cp\u003eFor each sample, the Chromium Controller platform of 10X Genomics was used for single-cell partitioning and barcoding. Each cell\u0026rsquo;s transcriptome was barcoded during reverse transcription, pooled cDNA was amplified and Single-Cell 3\u0026rsquo; Gene Expression were prepared according to the manufacturer\u0026rsquo;s protocol \u0026ldquo;Chromium NextGEM Single-Cell 3\u0026rsquo; Reagent Kits v3.1\u0026rdquo; (CG000315, 10X Genomics). Both Single-Cell 3\u0026rsquo; Gene Expression libraries were quantified on a 2100 Bioanalyzer Instrument following the manufacturer\u0026rsquo;s protocol \u0026ldquo;Agilent DNA 7500 kit\u0026rdquo; (G2938-90024, Agilent Technologies). These Single-Cell 3\u0026rsquo; Gene Expression libraries were combined to create one sequence library pool which was quantified by qPCR, according to manufacturer\u0026rsquo;s protocol \u0026ldquo;KAPA Library Quantification Kit Illumina\u0026reg; Platforms\u0026rdquo; (KR0405, KAPA Biosystems). A NovaSeq 6000 Illumina sequencing system was used for paired end sequencing of the Single-Cell 3\u0026rsquo; Gene Expression libraries at a sequencing depth of approximately 40,000-110,000 mean reads per cell. An estimated number of 5,000\u0026ndash;10,000 cells per sample were targeted. NovaSeq 6000 paired end sequencing was performed using 28 cycles for Read 1, 10 cycles for Read i7, 10 cycles for Read i5 and 90 cycles for Read 2, using NovaSeq SP Reagent Kit v1.5 (cat# 20028401, Illumina) and NovaSeq S2 Reagent Kit v1.5 (cat# 20028316, Illumina).\u003c/p\u003e\n\u003ch3\u003eSingle-cell RNA-sequencing data analysis\u003c/h3\u003e\n\u003cp\u003eScRNA-seq libraries were sequenced on an Illumina Novaseq 6000 using paired-end dual index reads. Demultiplexing and FastQ generation was performed using either bcl2fastq version 2.20.0 (Illumina) or BCLconvert version 3.9.3 (Illumina). FastQ data was aligned and quantified using the Cell Ranger package version 6.1.2 (10X Genomics). All analyses on the single-cell data were performed in R version 4.2.3 using Seurat version 4.3.0\u003csup\u003e26\u003c/sup\u003e. After loading the expression matrix DoubletFinder\u003csup\u003e\u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e27\u003c/span\u003e\u003c/sup\u003e was used to identify the heterotypic doublets. Cells labeled as doublets, having more than 25% mitochondrial reads or less than 1000 detected features were removed from the expression matrix which was subsequently normalized using SCTransform v2, regressing out cell cycle and mitochondrial fraction\u003csup\u003e\u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e28\u003c/span\u003e\u003c/sup\u003e. After dimensionality reduction and UMAP projection, the cells were clustered using the Louvain algorithm with a resolution of 1.0. ScopeLoomR was used to export and import the data for the SCENIC analysis. To identify malignant cells, each single-cell experiment was subset based on the AUCell score of Jerby-Arnon malignant gene set (\u0026gt;\u0026thinsp;0.12)\u003csup\u003e29\u003c/sup\u003e.\u003c/p\u003e \u003cdiv id=\"Sec11\" class=\"Section2\"\u003e \u003ch2\u003eSingle-cell regulatory network inference and clustering (SCENIC)\u003c/h2\u003e \u003cp\u003eSCENIC was run on the raw count matrix after mitochondrial read and doublet removal in addition to malignant cell identification as mentioned above. SCENIC analysis was performed using pySCENIC version 0.12.1. The co-expression modules procedure was run using GRNboost2, regulon prediction was run using the motif ranking databases hg38 500bp up/100bp down and 10kbp up/down (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://resources.aertslab.org\u003c/span\u003e\u003cspan address=\"https://resources.aertslab.org\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e). To limit potential undesirable effects of the stochastic nature of the gradient-boosting step within SCENIC, the complete pipeline was run ten times on the four patient samples individually (i.e., patients 35, 37, 44 and 45 in Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA). The set of predicted regulons were imported into R and filtered based on recurrence of both the detected regulons (6 out of the 10 runs) as well as the predicted target genes (6 out of 10 runs). The filtered regulons were then used as input for AUCell\u003csup\u003e\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section2\"\u003e \u003ch2\u003eRegulon activity and differentially activated regulons\u003c/h2\u003e \u003cp\u003eThe activity matrix generated by AUCell was added as an assay to the Seurat object. After scaling, the AUCell matrix was used for dimensionality reduction, UMAP projection and regulon clustering (Louvain). Differentially activated regulons were identified using Seurat\u0026rsquo;s FindAllMarkers using the Wilcoxon rank sum test requiring an adjusted P value\u0026thinsp;\u0026lt;\u0026thinsp;0.05, a log fold-change threshold of 0.01 and a minimum positive fraction of 0.95.\u003c/p\u003e \u003cp\u003eTo calculate the reciprocal overlap of regulons between two clusters across samples, the ratio of the overlapping regulons within each cluster was averaged.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec13\" class=\"Section2\"\u003e \u003ch2\u003eReactome analysis\u003c/h2\u003e \u003cp\u003eDifferentially active regulon intersect of patients 35, 37, 44 and 45 clusters was introduced to Reactome to perform pathway overrepresentation analysis. Negative log10 FDR values were used to plot the results, as they include a corrected over-representation probability.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec14\" class=\"Section2\"\u003e \u003ch2\u003eSpatial transcriptomics\u003c/h2\u003e \u003cp\u003eTo perform spatial transcriptomics, we made use of the Visium Spatial Gene Expression for FFPE platform of 10X Genomics. HE staining, tissue adhesion test quality check and all other procedures were performed as explained by Visium Spatial Gene Expression for FFPE from 10X Genomics ( \u003cspan fontcategory=\"NonProportional\" class=\"\" name=\"Emphasis\"\u003e\u003c/span\u003e CG000407, \u003cspan fontcategory=\"NonProportional\" class=\"\" name=\"Emphasis\"\u003e\u003c/span\u003e CG000408 and CG000409). Patient tumor material in FFPE was used as input material for the Visium. Region of interest was marked and slices were taken for DV200 estimation (Agilent TapeStation). FFPE blocks were chosen for Visium if the DV200 value was equal to or higher than 50% and had NGFR protein expression. Selected tumor quadrants were sliced and placed on the Visium slide. An extra slice of the same quadrant region was taken along for an additional NGFR staining.\u003c/p\u003e \u003cp\u003eVisium Spatial Gene Expression libraries were prepared according to the Visium Spatial Gene Expression User Guide for FFPE (CG000407). For each experiment, the Visium Spatial Gene Expression libraries were pooled to create one sequencing pool. The sequencing pool was quantified by qPCR, according to the KAPA Library Quantification Kit Illumina\u0026reg; Platforms protocol (KR0405, KAPA Biosystems). Paired end sequencing was performed on NextSeq 550 and NovaSeq 6000 Systems (Illumina) using Illumina sequencing reagent kits (cat no. 20024906, cat no. 20028401, Illumina), at a sequencing depth of approximately 30,000\u0026ndash;85,000 read pairs per tissue covered spot.\u003c/p\u003e \u003cp\u003eVisium libraries were sequenced on an Illumina Novaseq 6000 using paired-end dual index reads. Demultiplexing and FastQ generation was performed using Illumina BCLconvert version 3.9.3 (Illumina). FastQ data was quantified and spatially resolved using the Space Ranger package version 1.3.0 (10X Genomics). All analyses on the spatial data were performed in R version 4.2.3 using Seurat version 4.3.0\u003csup\u003e26\u003c/sup\u003e. The expression matrix was normalized using SCTransform v2, clustered using the Louvain algorithm with a resolution of 1.0.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec15\" class=\"Section2\"\u003e \u003ch2\u003eImaging Visium slide\u003c/h2\u003e \u003cp\u003eThe images have been acquired on an ZeissAxiover 200M microscope equipped with a 20x /0.75 NA Plan-Achromat objective, a 0.55 NA condenser and an Axiocam 512 color camera (pixel size 0.246 x 0.246 \u0026micro;m) in the ZEN 2.3 software. Images were stitched and exported as merged RGB color tiff.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec16\" class=\"Section2\"\u003e \u003ch2\u003eGene signatures scoring\u003c/h2\u003e \u003cp\u003eAUCell was used to score the enrichment of the different melanoma gene signatures\u003csup\u003e\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e,\u003cspan additionalcitationids=\"CR6 CR7\" citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e,\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e,\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e\u003c/sup\u003e for both the single-cell and the Visium data. For the Pozniak signatures, the top 100 most significant genes were selected for each signature. Of note, we excluded the Pozniak melanocytic signature for the Visium analyses, since it largely consisted of ribosomal genes which are not measured with the Visium FFPE probes set. The AUC value matrix was transformed to a Z-score for each gene-set. Then, the values were averaged for every Louvain cluster and the resulting matrix was shown as a heatmap using an additional scaling step for each row (gene signature) to emphasize the differences between the clusters of each patient tumor.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec17\" class=\"Section2\"\u003e \u003ch2\u003eSpatial mapping of single-cell regulon clusters\u003c/h2\u003e \u003cp\u003eSCENIC derived regulons from patient tumors 35, 37, 44 and 45 were used for scoring regulon activity, Louvain clustering and determining the differentially active regulons (identical to the single-cell method) on the matched Visiums. Next, the differentially active regulons of the clusters from both the scRNA-seq and the Visium 10X platforms were compared to determine the degree of similarity. When a scRNA-seq regulon scored the highest significant overlap with a Visium regulon cluster, the scRNA-seq regulon cluster number was transferred to the Visium regulon cluster.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec18\" class=\"Section2\"\u003e \u003ch2\u003eAnalysis of patient RNA-sequencing datasets\u003c/h2\u003e \u003cp\u003eRaw transcript counts from all datasets (Hugo\u003csup\u003e\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e\u003c/sup\u003e and Riaz\u003csup\u003e34)\u003c/sup\u003e were compiled using Kallisto 0.48.0. Samples with abnormal counts, as detected in principal component analysis (PCA) or \u0026le;5% of the median total transcript counts, were excluded. Similarly, transcripts with 0 counts or total counts \u0026le;0.1% of the median were also left out of the analysis. CombatSeq function from the SVA package (3.42.0) was used for dataset-related batch effect removal. Transcript counts were then transformed into TPMs considering the Kallisto produced length matrix. TPMs were then aggregated by gene and scaled for downstream analysis.\u003c/p\u003e \u003cp\u003eTo limit non-melanoma cell noise, samples were filtered based on purity and tumor markers. The ESTIMATE deconvolution algorithm in the immunedeconv package (2.1.0) was used to predict tumor purity. Samples that had lower than 0.4 estimated tumor purity score and were below 50% median TPM counts for the sum of the expression of \u003cem\u003eS100B\u003c/em\u003e, \u003cem\u003eSOX10\u003c/em\u003e, \u003cem\u003eMLANA\u003c/em\u003e, \u003cem\u003eMITF\u003c/em\u003e genes were filtered out.\u003c/p\u003e \u003cp\u003eNext, tumor samples were scored for the following signatures: the IFNg axis was based on the 28 gene signature by Ayers and colleagues\u003csup\u003e\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e\u003c/sup\u003e. The APM score was calculated using the eight selected genes by Thompson and colleagues \u003csup\u003e\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e\u003c/sup\u003e. T cell reactivity and T cell signatures were taken from Chow and colleagues\u003csup\u003e\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e\u003c/sup\u003e. The tertiary lymphoid structure (TLS) chemokines signature was derived from Li and colleagues\u003csup\u003e\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e\u003c/sup\u003e. For the TNF signature, we used the \u003cspan fontcategory=\"NonProportional\" class=\"\" name=\"Emphasis\"\u003e\u003c/span\u003e PID_TNF_PATHWAY\u003csup\u003e\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e\u003c/sup\u003e. The BATF3 signature was used as described by Hoefsmit and colleagues\u003csup\u003e\u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e40\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eSince antigenicity was calculated using different approaches in the considered datasets, a different calculation method was applied. The antigenicity predictor was scaled for each dataset independently, to avoid biases when employing different methods. For the Riaz and Hugo datasets we took their neopeptides that had an affinity \u0026le;500 nM and/or a ranking percentage \u0026le;2%. If expression data was available, only expressed neopeptides were considered. Final Z-scores where then scaled across all patients.\u003c/p\u003e \u003cp\u003eThe scaled TPM expression of each signature gene is averaged, then the resulting average is scaled across all patients. The Immune Activity Score (IAS) is calculated as the average of Z-scores for all the different signatures, including antigenicity, and aims to give an indication of the overall T cell presence, activity and priming status in the biopsied sample. Based on the median IAS score, we separated patient tumors in high and low immune active respectively (HI_IAS and LO_IAS). Each cohort contained responders and non-responders (R and NR) creating a final classification of four patient groups. All analyses on the combined RNA-seq data were performed in R version 4.1.3\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec19\" class=\"Section2\"\u003e \u003ch2\u003eSingle-cell analysis of publicly available datasets\u003c/h2\u003e \u003cp\u003eSingle-cell RNA-sequencing data was used from a study by Pozniak and colleagues\u003csup\u003e\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e\u003c/sup\u003e. Key marker genes for the mesenchymal-like cell state were calculated using Seurat\u0026rsquo;s (4.4.0) FindMarkers function across all available before treatment (BT) samples and on the SCT data. Those markers with an adjusted \u003cem\u003eP\u003c/em\u003e value\u0026thinsp;\u0026lt;\u0026thinsp;0.01, an average log2FoldChange\u0026thinsp;\u0026gt;\u0026thinsp;1.5 and pct.1\u0026thinsp;\u0026gt;\u0026thinsp;0.2 of the cluster\u0026rsquo;s cells were considered for further enrichment analysis of the combined RNA-seq dataset. Analyses were performed in R version 4.1.3.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec20\" class=\"Section2\"\u003e \u003ch2\u003eStatistics\u003c/h2\u003e \u003cp\u003eTo compare two means, a two-tailed Student \u003cem\u003eT\u003c/em\u003e test was used. One-way ANOVA test for more than 2 comparisons. For bulk RNA-seq, all signatures, scores and gene expression comparisons were statistically assessed using Wilcoxon rank tests. Statistics were performed by Prism (Graphpad Software Inc., version 9.0) or in R (4.1.3.). A \u003cem\u003eP\u003c/em\u003e value of lower than 0.05 was regarded as being statistically significant.\u003c/p\u003e \u003c/div\u003e"},{"header":"Results","content":"\u003cdiv id=\"Sec22\" class=\"Section2\"\u003e \u003ch2\u003eNGFR-expressing melanomas are highly heterogeneous for differentiation and dedifferentiation markers\u003c/h2\u003e \u003cp\u003eTo determine the degree of heterogeneity within NGFR\u003csup\u003e+\u003c/sup\u003e human melanomas, we first determined the expression of differentiation and dedifferentiation markers. We stained a panel of 65 clinical human melanoma samples derived from 46 patients for the differentiation markers MITF and Melan-A and the dedifferentiation markers AXL and NGFR\u003csup\u003e\u003cspan additionalcitationids=\"CR8 CR9 CR10\" citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e\u003c/sup\u003e \u003cb\u003e(Figure S1A)\u003c/b\u003e. Because we wished to study NGFR\u003csup\u003e+\u003c/sup\u003e melanoma subpopulations, we focused on 41 (63% of all stained patient tumors) of those, displaying a range of 1-100% NGFR positivity on viable tumor cells \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA\u003cb\u003e)\u003c/b\u003e. We observed that across these 41 tumors, NGFR correlated inversely with both Melan-A and MITF expression, while positively correlating with AXL, in agreement with previous observations \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eB-D\u003cb\u003e)\u003c/b\u003e\u003csup\u003e\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e,\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e,\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e\u003c/sup\u003e.\u003c/p\u003e \u003cp\u003eMITF and Melan-A mark a differentiated melanocytic phenotype, whereas NGFR marks a neural crest stem cell (NCSC) phenotype and AXL an undifferentiated mesenchymal phenotype \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eE\u003cb\u003e)\u003c/b\u003e\u003csup\u003e\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e,\u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e\u003c/sup\u003e. However, we noted that \u0026gt;\u0026thinsp;50% of NGFR\u003csup\u003e+\u003c/sup\u003e melanomas (22 out of 41) also comprised regions positive for differentiation (Melan-A and/or MITF) and dedifferentiation (AXL) markers \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA-I and \u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003eS1\u003c/span\u003eB-E). This included patient melanomas #35, 37, 44 and 45, which showed different degrees and patterns of marker heterogeneity and for which sufficient high-quality material was available for several downstream analyses \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eA-D, \u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eF-I and \u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003eS1\u003c/span\u003eB-E\u003cb\u003e)\u003c/b\u003e.\u003c/p\u003e \u003cp\u003eThus, the NGFR\u003csup\u003e+\u003c/sup\u003e fraction does not represent a single population but instead is heterogeneous, comprising a phenotypically diverse pool of differentiated and undifferentiated melanoma cells, and possibly intermediate cell states. This observation led us to investigate NGFR\u003csup\u003e+\u003c/sup\u003e heterogeneity in more detail, focusing on single cells, inter-patient variation, spatial organization, functional consequences and clinical responses. Therefore, we set up a pipeline to determine both transcriptional and spatial heterogeneity in NGFR\u003csup\u003e+\u003c/sup\u003e patient melanomas, combining scRNA-seq on FACS-sorted NGFR\u003csup\u003e+\u003c/sup\u003e viable cells with spatial transcriptomics on matched formalin-fixed paraffin-embedded (FFPE) samples \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003eJ\u003cb\u003e)\u003c/b\u003e, which was complemented with functional analyses and clinical responses to immunotherapy.\u003c/p\u003e \u003cdiv id=\"Sec23\" class=\"Section3\"\u003e \u003ch2\u003e​​​ ​\u003c/h2\u003e \u003cdiv id=\"Sec24\" class=\"Section4\"\u003e \u003ch2\u003eScRNA-seq identifies several transcriptomically distinct NGFR\u003csup\u003e+\u003c/sup\u003e melanoma phenotypes\u003c/h2\u003e \u003cp\u003eFirst, we set out to investigate whether the NGFR\u003csup\u003e+\u003c/sup\u003e fraction shows heterogeneity in gene expression at the single-cell level. ScRNA-seq was performed on FACS-sorted NGFR\u003csup\u003e+\u003c/sup\u003e melanoma cells from the above-mentioned four patient melanomas, from whom a total of 24,953 NGFR\u003csup\u003e+\u003c/sup\u003e/CD45\u003csup\u003e\u0026minus;\u003c/sup\u003e cells were isolated, expressing a wide range of cell surface NGFR expression levels \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;2A-B and \u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003eS2\u003c/span\u003eA\u003cb\u003e)\u003c/b\u003e. Using Uniform Manifold Approximation and Projection (UMAP) on the gene expression as a dimension reduction approach and Louvain clustering, we observed that NGFR\u003csup\u003e+\u003c/sup\u003e cells were represented by several different gene expression clusters. Melanoma cells from patients 35, 37 and 44 were grouped into 13 clusters and patient 45 in 17 clusters \u003cb\u003e(Figure S2B)\u003c/b\u003e. Expression of several markers, namely \u003cem\u003eSOX10\u003c/em\u003e, \u003cem\u003eS100B\u003c/em\u003e, \u003cem\u003eMITF\u003c/em\u003e and \u003cem\u003eMLANA\u003c/em\u003e, demonstrated that most sorted NGFR\u003csup\u003e+\u003c/sup\u003e cells were positive for markers from the melanocytic lineage \u003cb\u003e(Figure S2C).\u003c/b\u003e Furthermore, almost all sorted cells scored positively for a malignant melanoma signature\u003csup\u003e\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e\u003c/sup\u003e, confirming that we had selected cancer cells \u003cb\u003e(Figure S2D).\u003c/b\u003e Those melanoma cells were then subjected to a second round of UMAP analysis and Louvain clustering, which revealed 12, 12, 10 and 15 clusters for patients 35, 37, 44 and 45, respectively \u003cb\u003e(Fig.\u0026nbsp;2C)\u003c/b\u003e.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eTo determine the phenotypic identities of the NGFR\u003csup\u003e+\u003c/sup\u003e single cells within each patient tumor, the expression of key melanoma markers\u003csup\u003e\u003cspan additionalcitationids=\"CR8 CR9 CR10\" citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e,\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e,\u003cspan additionalcitationids=\"CR42\" citationid=\"CR41\" class=\"CitationRef\"\u003e41\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e43\u003c/span\u003e\u003c/sup\u003e was assessed and a series of melanoma signatures that were previously reported and commonly used were scored\u003csup\u003e\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e,\u003cspan additionalcitationids=\"CR6 CR7\" citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e,\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e,\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e\u003c/sup\u003e. This was done using AUCell\u003csup\u003e\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e\u003c/sup\u003e, an algorithm that quantifies the enrichment of input signature genes within the expressed genes of each single cell \u003cb\u003e(Fig.\u0026nbsp;2D)\u003c/b\u003e. All clusters showed absence of \u003cem\u003ePTPRC\u003c/em\u003e expression (encoding CD45), confirming the exclusion of immune cells. The signatures were grouped into five main categories: differentiated, undifferentiated, transitory/neural crest-like, stress-like and mitotic. The gene expression clusters from all four patients showed great diversity for these signatures, implying that cells that are positive for the neural crest/stem cell marker NGFR can exist in several different states along the (de)differentiation trajectory. Gene expression analysis of melanoma phenotype switch markers (i.e., \u003cem\u003eMITF\u003c/em\u003e, \u003cem\u003eMelan-A\u003c/em\u003e, \u003cem\u003eNGFR\u003c/em\u003e, \u003cem\u003eAXL\u003c/em\u003e, \u003cem\u003ePDGFRA\u003c/em\u003e, \u003cem\u003ePDGFRB\u003c/em\u003e and \u003cem\u003eEGFR\u003c/em\u003e) confirmed the presence of both differentiated and undifferentiated NGFR\u003csup\u003e+\u003c/sup\u003e subpopulations \u003cb\u003e(Fig.\u0026nbsp;2D)\u003c/b\u003e. Strikingly, several clusters scored positively for multiple signature categories. For example, patient 35 cluster 9 scored highly positive for undifferentiated, transitory/neural crest-like and stress-like categories, while patient 45 cluster 7, scored positively for differentiated, transitory/neural crest-like and stress-like groups. Together, these results indicate that individual patients harbor an NGFR\u003csup\u003e+\u003c/sup\u003e fraction comprising subpopulations that exist at different positions along the (de)differentiation trajectory.\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv id=\"Sec25\" class=\"Section3\"\u003e \u003ch2\u003eNGFR\u003csup\u003e+\u003c/sup\u003e melanoma cells are characterized by patient-specific gene-regulatory networks\u003c/h2\u003e \u003cp\u003eWe next sought to find a mechanistic explanation for the observed great inter-patient diversity in gene expression within the NGFR\u003csup\u003e+\u003c/sup\u003e cell fractions. We considered the possibility that NGFR marks different melanoma cell states governed by unique master regulators. To investigate this, we employed motif based \u003cspan fontcategory=\"NonProportional\" class=\"\" name=\"Emphasis\"\u003e\u003c/span\u003e single-cell regulatory network inference and clustering (SCENIC) on the scRNA-seq data followed by UMAP analysis and Louvain clustering. SCENIC is a computational tool for gene-regulatory network (GRN) inference and cell state discovery\u003csup\u003e\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e\u003c/sup\u003e, taking a three-step approach: inference of candidate target genes through co-expression inference, followed by motif-based filtering and AUC transcription factor (TF) activity scoring\u003csup\u003e\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e\u003c/sup\u003e. This allowed us to first identify key underlying GRNs of relative cell states within patients and, second, to interrogate whether these GRNs are shared between patients. SCENIC revealed 9\u0026ndash;13 independent cell states or regulon clusters within the different patients \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e3\u003c/span\u003eA\u003cb\u003e).\u003c/b\u003e Most regulon clusters, including cluster 4 in patient 37 and cluster 3 in patient 44, were characterized by dozens of differentially active TFs \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e3\u003c/span\u003eB and \u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003eS3\u003c/span\u003eA). Of note, an overlay of gene expression and regulon clusters revealed that some gene expression clusters were associated with the same regulons. This was evident for regulon cluster 2 in patient 35, which is composed of gene expression clusters 2 and 6; for regulon cluster 4 in patient 37 (gene expression clusters 8, and 9 and part of 1) and cluster 7 in patient 45 (gene expression clusters 12 and part of 0) \u003cb\u003e(Figure S3B)\u003c/b\u003e. These results suggest that NFGR marks several unique cell states per patient that are characterized by highly diverse GRNs.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eNext, the overlap between regulons across the four patients was determined. We found 103\u0026ndash;145 regulons per patient, 50 of which were shared among the four patients. Using Reactome enrichment analysis we observed a common enrichment across the four samples for Nerve Growth Factor (NGF)-stimulated transcription, consistent with our NGFR-focused scRNA-seq approach \u003cb\u003e(Figure S3C\u003c/b\u003e). Many other regulons were unique to individual patients, raising the possibility that cell states are established in a patient-dependent fashion. To investigate this further, we determined the average reciprocal overlap of differentially active regulons between clusters across patients. Only significant regulons per cluster were considered as cluster-specific cell state signatures. Little overlap was observed between clusters of the patients investigated \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e3\u003c/span\u003eD\u003cb\u003e)\u003c/b\u003e. Together, these results indicate that NGFR\u003csup\u003e+\u003c/sup\u003e melanoma cells are characterized by patient-specific GRNs.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec26\" class=\"Section3\"\u003e \u003ch2\u003eSpatial organization of NGFR\u003csup\u003e+\u003c/sup\u003e subpopulations\u003c/h2\u003e \u003cp\u003eWhile highlighting the diversity of melanoma GRNs, the analyses above do not consider how the corresponding cells are spatially organized. We considered this a key parameter to include in our study, given the great intrinsic diversity of NGFR\u003csup\u003e+\u003c/sup\u003e melanoma fractions as revealed by our immunohistochemical \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e\u003cb\u003e)\u003c/b\u003e and scRNA-seq analyses \u003cb\u003e(\u003c/b\u003eFigs.\u0026nbsp;2 and \u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e3\u003c/span\u003e\u003cb\u003e)\u003c/b\u003e. Therefore, we used Visium 10X spatial transcriptomics to dissect the spatial and transcriptomic NGFR heterogeneity in the same patient melanomas. For this analysis, selection of tumor quadrant areas was guided by the presence of NGFR protein expression \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e4\u003c/span\u003eA\u003cb\u003e).\u003c/b\u003e Two slices were taken of the quadrant area of each single FFPE block. One was used for H\u0026amp;E staining on the Visium 10X slide, and the other one was stained for NGFR \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e4\u003c/span\u003eB and \u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003eS4\u003c/span\u003eA\u003cb\u003e)\u003c/b\u003e. In agreement with the NGFR IHC staining results \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e4\u003c/span\u003eB\u003cb\u003e)\u003c/b\u003e, \u003cem\u003eNGFR\u003c/em\u003e expression was detected in the Visium 10X quadrants, again showing highly heterogeneous patterns \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e4\u003c/span\u003eC\u003cb\u003e)\u003c/b\u003e. Furthermore, all samples showed large areas scoring positively for the malignant melanoma signature\u003csup\u003e\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e\u003c/sup\u003e, indicative of the presence of melanoma cells \u003cb\u003e(Figure S4B)\u003c/b\u003e.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eTo investigate the spatial gene expression patterns of the NGFR\u003csup\u003e+\u003c/sup\u003e tumor regions, gene expression Louvain clustering was performed, followed by melanoma signature scoring and phenotype switch marker gene expression analysis. Among the patients, we identified 10\u0026ndash;13 unique gene expression clusters that were spatially organized \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e4\u003c/span\u003eD\u003cb\u003e)\u003c/b\u003e. Next, we scored these clusters for a series of common signatures representative of different melanoma cell states\u003csup\u003e\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e,\u003cspan additionalcitationids=\"CR6 CR7\" citationid=\"CR5\" class=\"CitationRef\"\u003e5\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e,\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e,\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e\u003c/sup\u003e. We observed again a high degree of variation in melanoma cell state signatures and phenotype switch marker genes \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e4\u003c/span\u003eE\u003cb\u003e)\u003c/b\u003e. This result supports our scRNA-seq analyses showing that NGFR\u003csup\u003e+\u003c/sup\u003e melanoma cells can exist in several different states along the (de)differentiation trajectory \u003cb\u003e(Fig.\u0026nbsp;2D)\u003c/b\u003e.\u003c/p\u003e \u003cp\u003eTo complement this analysis with spatial annotation, we determined the location of the scRNA-seq regulon clusters. We first generated Visium regulon clusters for each patient. The single cell regulon list was used as input to score the Visium spots, employing AUCell followed by Louvain clustering \u003cb\u003e(Figure S4C)\u003c/b\u003e. Next, we highlighted scRNA-seq regulon clusters that most significantly overlapped with Visium regulon clusters in the Visium spots. Whereas most scRNA-seq regulon clusters were spatially organized, there were clear patterns of mutual cell spreading to other compartments, which was apparent for all patients. For example, in patient 44, single cell regulon clusters 2 and 5 formed spatially different compartments, but each formed several satellites in the other cluster \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e4\u003c/span\u003eF\u003cb\u003e)\u003c/b\u003e. Similarly, satellite formation in other regulon clusters was apparent for clusters 0, 5 and 10 in patient 45, and for clusters 2, 6 and 8 in patient 35. Thus, we observe a pattern, common among melanoma patients, of NGFR\u003csup\u003e+\u003c/sup\u003e cell states that are mostly spatially organized yet also form satellites in spatial regions dominated by other cell states.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec27\" class=\"Section3\"\u003e \u003ch2\u003ePDGFR marks cytokine-resistant and mesenchymal-like NGFR\u003csup\u003e+\u003c/sup\u003e melanoma cells\u003c/h2\u003e \u003cp\u003eWe previously reported that NGFR\u003csup\u003e+\u003c/sup\u003e melanoma populations are less sensitive to T cell cytokines\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e. The findings above indicate that even within the melanoma NGFR\u003csup\u003e+\u003c/sup\u003e cell fraction there is great transcriptional and phenotypic heterogeneity. This raises the possibility that this diversity is associated with functional differences. We thus investigated any functional consequences of NGFR heterogeneity for cytokine sensitivity. A panel of representative NGFR\u003csup\u003e+\u003c/sup\u003e melanoma cell lines was selected using the established phenotype switch markers Melan-A, AXL, PDGFRa and PDGFRb \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e5\u003c/span\u003eA and \u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003eS5\u003c/span\u003eA-C\u003cb\u003e)\u003c/b\u003e. Next, we examined the cytotoxic potential of T cell cytokines IFNg and TNF across this cell line panel. We observed that specifically the cell lines marked by PDGFRa and PDGFRb expression (\u0026ldquo;PDGFR\u0026rdquo;, i.e., A875, M080.X1.CL and M019R.X1.CL), showed the highest cytokine resistance \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e5\u003c/span\u003eB, C and \u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003eS5\u003c/span\u003eD\u003cb\u003e)\u003c/b\u003e. These results show that while NGFR\u003csup\u003e+\u003c/sup\u003e melanoma fractions already show reduced sensitivity to cytokines\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e, NGFR\u003csup\u003e+\u003c/sup\u003e/PDGFR\u003csup\u003e+\u003c/sup\u003e double-positive tumor cells show an even higher degree of cytokine resistance.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eTo complement these functional and our Visium spatial analyses, we next examined the location of NGFR\u003csup\u003e+\u003c/sup\u003e/PDGFR\u003csup\u003e+\u003c/sup\u003e double-positive cells in clinical melanoma samples. We performed NGFR and PDGFRb IHC staining on consecutive FFPE slices combined with Periodic Acid\u0026ndash;Schiff (PAS) to distinguish the tumor from the stromal component (i.e., vessels and collagen). PDGFRb expression was detected both on melanoma cells in NGFR\u003csup\u003e+\u003c/sup\u003e tumor regions and in the vascular area, the latter being indicative of pericytes \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e5\u003c/span\u003eD\u003cb\u003e)\u003c/b\u003e. To validate these observations, we analyzed a scRNA-seq melanoma cohort\u003csup\u003e\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e\u003c/sup\u003e. We found that the \u003cem\u003eNGFR\u003c/em\u003e\u003csup\u003e\u003cem\u003e+\u003c/em\u003e\u003c/sup\u003e\u003cem\u003e/PDGFR\u003c/em\u003e\u003csup\u003e\u003cem\u003e+\u003c/em\u003e\u003c/sup\u003e melanoma cell fraction was highly enriched in the mesenchymal-like cell fraction (\u0026ldquo;MES\u0026rdquo;) in baseline melanomas \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e5\u003c/span\u003eE\u003cb\u003e)\u003c/b\u003e. In the other three patients studied we were unable to convincingly identify NGFR\u003csup\u003e+\u003c/sup\u003e/PDGFRb\u003csup\u003e+\u003c/sup\u003e melanoma cells \u003cb\u003e(Figure S5E)\u003c/b\u003e, confirming the low frequency of this double-positive population seen in the scRNA-seq analyses \u003cb\u003e(Fig.\u0026nbsp;2D)\u003c/b\u003e. Therefore, we extended the panel of melanoma samples with three additional patients (7b, 36, 43). This analysis validated our finding of NGFR\u003csup\u003e+\u003c/sup\u003e/PDGFRb\u003csup\u003e+\u003c/sup\u003e melanoma cells, which were again often located near the stromal component including vessels \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e5\u003c/span\u003eD\u003cb\u003e)\u003c/b\u003e. Together, PDGFR marks a cytokine-resistant melanoma subpopulation, which can be detected in patient samples and is associated with a mesenchymal-like phenotype.\u003c/p\u003e \u003cp\u003e \u003cb\u003ePDGFR\u003c/b\u003e \u003cb\u003emarks immune-active yet ICB non-responding melanomas\u003c/b\u003e\u003c/p\u003e \u003cp\u003eLastly, we investigated whether this PDGFR\u003csup\u003e+\u003c/sup\u003e mesenchymal-like and cytokine-resistant melanoma population correlates with therapeutic response to ICB. The presence of T cells in the tumor microenvironment is required for ICB response and is associated with improved patient outcome\u003csup\u003e\u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e29\u003c/span\u003e, \u003cspan additionalcitationids=\"CR47\" citationid=\"CR46\" class=\"CitationRef\"\u003e46\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e48\u003c/span\u003e\u003c/sup\u003e. An immune activity score (IAS low, intermediate, high) was generated based on the sum of averaged Z-score expression of each of the following gene signatures, namely: IFNg\u003csup\u003e\u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e35\u003c/span\u003e\u003c/sup\u003e, antigen presentation machinery (APM)\u003csup\u003e\u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e36\u003c/span\u003e\u003c/sup\u003e, T cells\u003csup\u003e\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e\u003c/sup\u003e, T cell reactivity\u003csup\u003e\u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e37\u003c/span\u003e\u003c/sup\u003e, TLS-chemokines\u003csup\u003e\u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e38\u003c/span\u003e\u003c/sup\u003e, TNF\u003csup\u003e\u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e39\u003c/span\u003e\u003c/sup\u003e and BATF3\u003csup\u003e40\u003c/sup\u003e, as well as antigenicity (i.e., total neoantigen counts) \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e6\u003c/span\u003eA and \u003cspan refid=\"Fig11\" class=\"InternalRef\"\u003eS6\u003c/span\u003eA\u003cb\u003e)\u003c/b\u003e. Two public RNA-seq melanoma datasets (i.e., Hugo\u003csup\u003e\u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e33\u003c/span\u003e\u003c/sup\u003e and Riaz\u003csup\u003e\u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e34\u003c/span\u003e\u003c/sup\u003e) were combined for this analysis.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eWe observed that both IAS-high and IAS-low patients were similarly distributed among responders \u003cb\u003e(\u003c/b\u003eHI_IAS/R and LO_IAS/R) and non-responders \u003cb\u003e(\u003c/b\u003eHI_IAS/NR and LO_IAS/NR, Fig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e6\u003c/span\u003eA and \u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e6\u003c/span\u003eB\u003cb\u003e)\u003c/b\u003e. To determine whether the cytokine-resistant melanoma population marked by NGFR and PDGFR was associated with poor ICB response, expression of both \u003cem\u003ePDGFR and NGFR\u003c/em\u003e was analyzed in the combined RNA-seq dataset \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e6\u003c/span\u003eA and \u003cspan refid=\"Fig11\" class=\"InternalRef\"\u003eS6\u003c/span\u003eB\u003cb\u003e)\u003c/b\u003e. We noticed a significantly higher expression of \u003cem\u003ePDGFR\u003c/em\u003e in immune-active ICB non-responders (HI_IAS/NR) compared to both non-responders with low immune activation (LO_IAS/NR) and responders with a high immune activation score (HI_IAS/R) \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e6\u003c/span\u003eC and \u003cspan refid=\"Fig11\" class=\"InternalRef\"\u003eS6\u003c/span\u003eB\u003cb\u003e)\u003c/b\u003e. Supporting this observation, there was a similar significant enrichment of the MES signature in HI_IAS/NR tumors compared to both LO_IAS/NR and HI_IAS/R groups \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e6\u003c/span\u003eD\u003cb\u003e)\u003c/b\u003e. This suggests that \u003cem\u003ePDGFR\u003c/em\u003e\u003csup\u003e+\u003c/sup\u003e mesenchymal-like melanoma cells are significantly enriched in immune-active, yet ICB non-responding, melanoma patients.\u003c/p\u003e \u003cp\u003eLastly, to complement our RNA-seq findings, we investigated the presence of tumor- infiltrating T cells in PDGFRb\u003csup\u003e+\u003c/sup\u003e and PDGFRb\u003csup\u003e\u0026minus;\u003c/sup\u003e tumors by staining 33 melanoma samples for both PDGFRb/PAS and CD3 (with or without PAS). We identified 10 tumors harboring PDGFRb\u003csup\u003e+\u003c/sup\u003e melanoma cells, thereby confirming that PDGFRb is expressed in a subfraction of patient melanomas \u003cb\u003e(Figure S6C)\u003c/b\u003e. Of these, nine tumors (90%) were infiltrated with T cells. Although PDGFRb\u003csup\u003e\u0026minus;\u003c/sup\u003e tumors showed a slightly higher average CD3 score (fold change 1.33) than PDGFRb\u003csup\u003e+\u003c/sup\u003e tumors, this difference was not statistically significant \u003cb\u003e(\u003c/b\u003eFig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e6\u003c/span\u003eE\u003cb\u003e)\u003c/b\u003e. Moreover, IHC analysis revealed that in most PDGFRb\u003csup\u003e+\u003c/sup\u003e melanomas, CD3\u003csup\u003e+\u003c/sup\u003e T cells showed no clear association with either PDGFRb\u003csup\u003eHigh\u003c/sup\u003e or PDGFRb\u003csup\u003eLow\u003c/sup\u003e tumor regions (i.e., 7/9 infiltrated tumors (78%), Fig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e6\u003c/span\u003eF\u003cb\u003e)\u003c/b\u003e. Together, these results show that \u003cem\u003ePDGFR\u003c/em\u003e\u003csup\u003e+\u003c/sup\u003e mesenchymal-like tumor cells are enriched in T cell-infiltrated, immune-active, but ICB non-responding melanomas.\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e"},{"header":"DISCUSSION","content":"\u003cp\u003eMelanoma is positioned at the extreme end of the mutational spectrum in cancer\u003csup\u003e\u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2\u003c/span\u003e,\u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e50\u003c/span\u003e\u003c/sup\u003e. Combined with its high degree of phenotypic plasticity, this results in substantial heterogeneity, as we and others have shown previously\u003csup\u003e\u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e1\u003c/span\u003e,\u003cspan additionalcitationids=\"CR5 CR6 CR7 CR8 CR9 CR10\" citationid=\"CR4\" class=\"CitationRef\"\u003e4\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e\u003c/sup\u003e. For example, a small yet common population of melanoma cells marked by NGFR\u003csup\u003e+\u003c/sup\u003e expression is associated with reduced susceptibility to T cell killing\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e,\u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e15\u003c/span\u003e\u003c/sup\u003e, NK-cell killing\u003csup\u003e\u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e24\u003c/span\u003e\u003c/sup\u003e, immune exclusion\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e, immunotherapy resistance\u003csup\u003e\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e\u003c/sup\u003e and phenotypic plasticity\u003csup\u003e\u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e18\u003c/span\u003e\u003c/sup\u003e, both in patients and mouse models. However, it is incompletely understood whether this NGFR\u003csup\u003e+\u003c/sup\u003e fraction represents a homogenous cell group or instead comprises yet additional heterogeneous cell populations, and whether this has consequences for therapy response. By integrating scRNA-seq, Visium spatial analysis, functional studies and analyses of clinical ICB cohorts, we find an unexpectedly high degree of heterogeneity even within these NGFR\u003csup\u003e+\u003c/sup\u003e cell fractions. This is manifested at the level of both expression and regulation of transcriptional networks, phenotypically, functionally and clinically.\u003c/p\u003e \u003cp\u003eTo begin dissecting NGFR heterogeneity in clinical samples, we employed IHC, confirming that NGFR\u003csup\u003e+\u003c/sup\u003e cells can be identified in most human melanomas. Previous studies have shown NGFR expression to be associated with an AXL\u003csup\u003ehigh\u003c/sup\u003e program, specifically marking the NCSC state in melanoma\u003csup\u003e\u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e7\u003c/span\u003e,\u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e8\u003c/span\u003e,\u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e11\u003c/span\u003e,\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e,\u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e51\u003c/span\u003e\u003c/sup\u003e. However, while we were able to confirm the positive association between NGFR and dedifferentiation marker AXL in clinical samples, we also observed expression of differentiation markers MITF and/or Melan-A in NGFR\u003csup\u003e+\u003c/sup\u003e tumor regions. This prompted us to investigate whether this heterogeneity is also seen at the single cell level.\u003c/p\u003e \u003cp\u003e \u003cem\u003eIn vitro\u003c/em\u003e studies revealed that NGFR, after initial upregulation, is strongly downregulated when reaching a full undifferentiated/mesenchymal-like cell state upon long-term MAPK inhibition\u003csup\u003e\u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e17\u003c/span\u003e\u003c/sup\u003e. However, whether NGFR marks multiple melanoma transcriptional programs along the (de)differentiation scale in patient melanomas has not been addressed. By employing scRNA-seq, we observed that indeed, NGFR marks a spectrum of (de)differentiation in clinical melanoma samples. Specifically, NGFR\u003csup\u003e+\u003c/sup\u003e cell populations display differentiated, undifferentiated, transitory/neural crest-like, stress-like and mitotic signature categories, as well as combinations of these, within each patient, to different extents. This indicates that NFGR marks a transcriptomically diverse pool of melanoma cells in clinical melanoma samples.\u003c/p\u003e \u003cp\u003eConsidering the observed transcriptomic heterogeneity within the NGFR\u003csup\u003e+\u003c/sup\u003e cell fraction, we hypothesized that the different cell states marked by NGFR are governed by multiple gene regulatory networks (GRNs). To test this, SCENIC\u003csup\u003e\u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e30\u003c/span\u003e,\u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e45\u003c/span\u003e\u003c/sup\u003e was used to analyze the scRNA-seq patient data. At least nine different NGFR\u003csup\u003e+\u003c/sup\u003e cell states per patient were identified, each characterized by a unique GRN. It was recently shown that melanoma GRNs show diversity among individual patients\u003csup\u003e\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e\u003c/sup\u003e. Observing no common GRN across the patients studied here, our results demonstrate that this heterogeneity among patients even extends to the NGFR\u003csup\u003e+\u003c/sup\u003e subpopulations. Thus, NGFR marks patient-specific cell states characterized by unique GRNs.\u003c/p\u003e \u003cp\u003eTo capture the spatial distribution of NGFR\u003csup\u003e+\u003c/sup\u003e melanoma subpopulations, Visium 10X spatial transcriptomics was employed. Previous studies have shown the spatial distribution of NGFR\u003csup\u003e+\u003c/sup\u003e melanoma cells in skin-\u003csup\u003e\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e\u003c/sup\u003e, metastatic-\u003csup\u003e\u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e20\u003c/span\u003e,\u003cspan citationid=\"CR52\" class=\"CitationRef\"\u003e52\u003c/span\u003e\u003c/sup\u003e and patient-derived xenografts (PDX) lesions\u003csup\u003e\u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e53\u003c/span\u003e\u003c/sup\u003e. However, an in-depth spatial transcriptomic characterization of NGFR\u003csup\u003e+\u003c/sup\u003e tumor regions per patient was lacking. We show here that consistent with our scRNA-seq findings, there is a spatial spectrum of (de)differentiation within NGFR\u003csup\u003e+\u003c/sup\u003e patient tumor regions. Moreover, by combining scRNA-seq and Visium, we show that also the spatial distribution of NGFR\u003csup\u003e+\u003c/sup\u003e cell states is governed by unique GRNs. We observed that cell states are spatially organized in melanoma, often forming satellites in tumor areas dominated by other cell states. A potential force driving satellite formation may be stochastic variability in gene expression\u003csup\u003e\u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e54\u003c/span\u003e\u003c/sup\u003e with or without reinforcement of Lamarckian induction by microenvironmental stimuli and physical cues\u003csup\u003e\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e3\u003c/span\u003e,\u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e55\u003c/span\u003e\u003c/sup\u003e. In other words, it is plausible that local microenvironmental changes are intertwined with the outgrowth of a particular cell state or affects cell state identity\u003csup\u003e\u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e54\u003c/span\u003e,\u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e55\u003c/span\u003e\u003c/sup\u003e. For instance, the stem-like murine melanoma phenotype was found to be associated with, and fueled by, endothelial cells\u003csup\u003e\u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e6\u003c/span\u003e,\u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e55\u003c/span\u003e\u003c/sup\u003e. Whether stochastic gene expression variation and Lamarckian induction are independently or cooperatively causing satellite formation of different NGFR\u003csup\u003e+\u003c/sup\u003e cell states in melanoma needs to be further investigated. From this study, we can conclude that while NGFR\u003csup\u003e+\u003c/sup\u003e cell states are mostly spatially organized, they also form satellites in tumor regions dominated by other cell states across different patients.\u003c/p\u003e \u003cp\u003eConsidering the transcriptomic and phenotypic diversity within the NGFR\u003csup\u003e+\u003c/sup\u003e melanoma cell fractions, we studied any functional consequence of this heterogeneity for immune cytokine sensitivity. We made use of a representative NGFR\u003csup\u003e+\u003c/sup\u003e cell line panel co-expressing additional phenotype switch markers. We observed that PDGFR expression was consistently associated with T cell cytokine resistance. No association between AXL expression and cytokine resistance was observed, which we had also noted in the context of T cell resistance\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e. By employing IHC, we identified NGFR\u003csup\u003e+\u003c/sup\u003e/PDGFR\u003csup\u003e+\u003c/sup\u003e cells often near the stroma/vessels. By leveraging an external scRNA-seq dataset, we validated the presence of NGFR\u003csup\u003e+\u003c/sup\u003e/PDGFR\u003csup\u003e+\u003c/sup\u003e cells in additional patient melanomas. These double-positive cells were enriched in the mesenchymal-like (MES) cell state tumor fraction\u003csup\u003e\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e\u003c/sup\u003e. Thus, while NGFR marks a relatively insensitive T cell cytokine-phenotype\u003csup\u003e\u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e19\u003c/span\u003e\u003c/sup\u003e, here we show that particularly PDGFR marks the most cytokine-resistant NGFR\u003csup\u003e+\u003c/sup\u003e mesenchymal-like subpopulations.\u003c/p\u003e \u003cp\u003eImmunotherapy has improved overall survival in advanced melanoma, reaching a clinical benefit over 50%\u003csup\u003e56\u0026ndash;58\u003c/sup\u003e. However, a significant proportion of patients either intrinsically harbors or develops acquired resistance to immunotherapy, due to diverse cell intrinsic or extrinsic mechanisms\u003csup\u003e\u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e32\u003c/span\u003e,\u003cspan citationid=\"CR59\" class=\"CitationRef\"\u003e59\u003c/span\u003e,\u003cspan citationid=\"CR60\" class=\"CitationRef\"\u003e60\u003c/span\u003e\u003c/sup\u003e. It is therefore imperative to explore novel treatment options\u003csup\u003e61\u003c/sup\u003e while simultaneously investing in biomarker identification to improve prediction of therapy resistance. This may also help circumventing unnecessary immune-related adverse events (irAEs)\u003csup\u003e\u003cspan citationid=\"CR61\" class=\"CitationRef\"\u003e62\u003c/span\u003e\u003c/sup\u003e. Thus, identifying upfront therapy resistance is a relevant clinical need\u003csup\u003e\u003cspan citationid=\"CR62\" class=\"CitationRef\"\u003e63\u003c/span\u003e\u003c/sup\u003e. Since PDGFR marks a mesenchymal-like cytokine-resistant melanoma population, we investigated whether this population was association with a therapeutic response to ICB. By analyzing two melanoma clinical datasets, we identified a subset of non-responders characterized by high immune activity, as reflected by the presence of active T cells, enriched with \u003cem\u003ePDGFR\u003c/em\u003e and MES signature expression.\u003c/p\u003e \u003cp\u003eRecently, it has been reported that an increase in the MES cell fraction, governed by TCF4, in early ICB on-treatment patient tumor biopsies is associated with ICB resistance\u003csup\u003e\u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e31\u003c/span\u003e\u003c/sup\u003e. In the current study, we show that the presence of \u003cem\u003ePDGFR\u003c/em\u003e\u003csup\u003e+\u003c/sup\u003e MES cells at baseline is predictive of ICB resistance, particularly in the context of immune-active and T cell-infiltrated melanomas. Of note, TCF4 was not identified as a significant differentially active regulon in any of the studied scRNA-seq regulon clusters. This implies that TCF4 is not a key transcription factor driving any of the characterized NGFR\u003csup\u003e+\u003c/sup\u003e melanoma cell states, including the \u003cem\u003ePDGFR\u003c/em\u003e\u003csup\u003e\u003cem\u003e+\u003c/em\u003e\u003c/sup\u003e cell states. Given that NGFR\u003csup\u003e+\u003c/sup\u003e/PDGFR\u003csup\u003e+\u003c/sup\u003e melanoma cells pre-exist in untreated melanomas, a challenge for the future is to investigate which cell-autonomous or non-autonomous factors\u003csup\u003e\u003cspan additionalcitationids=\"CR54\" citationid=\"CR53\" class=\"CitationRef\"\u003e53\u003c/span\u003e\u0026ndash;\u003cspan citationid=\"CR55\" class=\"CitationRef\"\u003e55\u003c/span\u003e\u003c/sup\u003e contribute to their development and affect therapy.\u003c/p\u003e \u003cp\u003eIn conclusion, we show here a remarkable degree of heterogeneity in a heavily studied melanoma subpopulation, namely NGFR\u003csup\u003e+\u003c/sup\u003e cells. Several NGFR\u003csup\u003e+\u003c/sup\u003e melanoma subpopulations position along the (de)differentiation gray scale, driven by patient-specific GRNs and spatially organized as satellites in tumor regions dominated by other cell states. The co-expression of PDGFR, marking T cell cytokine resistance in immune-active ICB non-responding patients may serve as a biomarker for immunotherapy resistance, which may have diagnostic potential.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eDisclosure\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eD.S.P. is co-founder, shareholder and advisor of Flindr Tx, which is unrelated to this study.\u0026nbsp;\u003c/p\u003e\u003cp\u003e\u003cstrong\u003eACKNOWLEDGEMENTS\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eWe thank all the members of the Peeper lab for their valuable input. Moreover, we thank our in-house flowcytometry, animal pathology, sequencing, and bioimaging facilities for their help and support. We would like to acknowledge Alexander van Akkooi, Winan van Houdt, John B. Haanen, Lisanne Zijlker and Max F. Madu for collecting patient samples. Moreover, we thank the NKI-AVL Core Facility Molecular Pathology \u0026amp; Biobanking (CFMPB) for supplying NKI-AVL Biobank material and performing IHC stainings. \u0026nbsp;Lastly, we thank Joyce Sanders for helping with FFPE tumor block selection for Visium 10X. D.S.Peeper is funded by the Oncode Institute and by the Dutch Cancer Society KWF and this work was supported by the Koningin Wilhelmina Fonds (KWF; project no. 10425) and by MRA 681127 to D.S. Peeper.\u003c/p\u003e\n\u003cp\u003e\u0026nbsp;\u003cstrong\u003eAUTHOR CONTRIBUTIONS:\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eS.M.S and D.S.P:\u0026nbsp;\u003c/strong\u003econceived study and\u003cstrong\u003e\u0026nbsp;\u003c/strong\u003edesigned experiments.\u003cstrong\u003e\u0026nbsp;S.M.S:\u003c/strong\u003e performed experiments, investigation, visualization, (bioinformatic) data interpretation, writing\u0026ndash;original draft and project administration. \u003cstrong\u003eJ.J.H.T:\u0026nbsp;\u003c/strong\u003eperformed bioinformatic analyses and provided critical input. \u003cstrong\u003eA.V\u003c/strong\u003e, \u003cstrong\u003eI.d.R\u003c/strong\u003e, \u003cstrong\u003eA.v.V\u003c/strong\u003e, \u003cstrong\u003eA.G\u0026nbsp;\u003c/strong\u003eand\u003cstrong\u003e\u0026nbsp;J.S.N\u003c/strong\u003e: performed bioinformatic analyses. \u003cstrong\u003eM.N:\u0026nbsp;\u003c/strong\u003eprocessed Visium 10X samples. \u003cstrong\u003eS.B\u003c/strong\u003e:\u003cstrong\u003e\u0026nbsp;\u003c/strong\u003eHelped collecting clinical samples. \u003cstrong\u003eS.M.S\u003c/strong\u003e, \u003cstrong\u003eJ-Y.S, J.B, S.B, M.D.S.G and H.M.H\u003c/strong\u003e: analyzed human pathology samples.\u003cstrong\u003e\u0026nbsp;I.H and L.d.V\u003c/strong\u003e: performed IHC-, PAS- and HE stainings and prepared/sliced FFPE blocks for Visium 10X. \u003cstrong\u003eA.E.D\u003c/strong\u003e: imaging of Visium 10X (HE stained) patient samples. \u003cstrong\u003eM.v.B\u003c/strong\u003e: performed cell sorting. \u003cstrong\u003eS.M.S and D.S.P:\u0026nbsp;\u003c/strong\u003ewrote the manuscript. The project was supervised by \u003cstrong\u003eD.S.P\u003c/strong\u003e. All authors reviewed and approved the manuscript.\u0026nbsp;\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eWouters, J., Kalender-Atak, Z., Minnoye, L., Spanier, K. I., De Waegeneer, M., Bravo Gonz\u0026aacute;lez-Blas, C., Mauduit, D., Davie, K., Hulselmans, G., Najem, A., Dewaele, M., Pedri, D., Rambow, F., Makhzami, S., Christiaens, V., Ceyssens, F., Ghanem, G., Marine, J. C., Poovathingal, S., \u0026amp; Aerts, S. (2020). Robust gene expression programs underlie recurrent cell states and phenotype switching in melanoma. \u003cem\u003eNature cell biology\u003c/em\u003e, \u003cem\u003e22\u003c/em\u003e(8), 986\u0026ndash;998. https://doi.org/10.1038/s41556-020-0547-3\u003c/li\u003e\n\u003cli\u003eAlexandrov, L. B., Nik-Zainal, S., Wedge, D. C., Aparicio, S. A., Behjati, S., Biankin, A. V., Bignell, G. R., Bolli, N., Borg, A., B\u0026oslash;rresen-Dale, A. L., Boyault, S., Burkhardt, B., Butler, A. P., Caldas, C., Davies, H. R., Desmedt, C., Eils, R., Eyfj\u0026ouml;rd, J. E., Foekens, J. A., Greaves, M., \u0026hellip; Stratton, M. R. (2013). Signatures of mutational processes in human cancer. \u003cem\u003eNature\u003c/em\u003e, \u003cem\u003e500\u003c/em\u003e(7463), 415\u0026ndash;421. https://doi.org/10.1038/nature12477\u003c/li\u003e\n\u003cli\u003eMarine, J. C., Dawson, S. J., \u0026amp; Dawson, M. A. (2020). Non-genetic mechanisms of therapeutic resistance in cancer. \u003cem\u003eNature reviews. Cancer\u003c/em\u003e, \u003cem\u003e20\u003c/em\u003e(12), 743\u0026ndash;756. https://doi.org/10.1038/s41568-020-00302-4\u003c/li\u003e\n\u003cli\u003eHoek, K. S., Eichhoff, O. M., Schlegel, N. C., D\u0026ouml;bbeling, U., Kobert, N., Schaerer, L., Hemmi, S., \u0026amp; Dummer, R. (2008). In vivo switching of human melanoma cells between proliferative and invasive states. \u003cem\u003eCancer research\u003c/em\u003e, \u003cem\u003e68\u003c/em\u003e(3), 650\u0026ndash;656. https://doi.org/10.1158/0008-5472.CAN-07-2491 \u003c/li\u003e\n\u003cli\u003eVerfaillie, A., Imrichova, H., Atak, Z. K., Dewaele, M., Rambow, F., Hulselmans, G., Christiaens, V., Svetlichnyy, D., Luciani, F., Van den Mooter, L., Claerhout, S., Fiers, M., Journe, F., Ghanem, G. E., Herrmann, C., Halder, G., Marine, J. C., \u0026amp; Aerts, S. (2015). Decoding the regulatory landscape of melanoma reveals TEADS as regulators of the invasive cell state. \u003cem\u003eNature communications\u003c/em\u003e, \u003cem\u003e6\u003c/em\u003e, 6683. https://doi.org/10.1038/ncomms7683 \u003c/li\u003e\n\u003cli\u003eKarras, P., Bordeu, I., Pozniak, J., Nowosad, A., Pazzi, C., Van Raemdonck, N., Landeloos, E., Van Herck, Y., Pedri, D., Bervoets, G., Makhzami, S., Khoo, J. H., Pavie, B., Lamote, J., Marin-Bejar, O., Dewaele, M., Liang, H., Zhang, X., Hua, Y., Wouters, J., \u0026hellip; Marine, J. C. (2022). A cellular hierarchy in melanoma uncouples growth and metastasis. \u003cem\u003eNature\u003c/em\u003e, \u003cem\u003e610\u003c/em\u003e(7930), 190\u0026ndash;198. https://doi.org/10.1038/s41586-022-05242-7 \u003c/li\u003e\n\u003cli\u003eRambow, F., Rogiers, A., Marin-Bejar, O., Aibar, S., Femel, J., Dewaele, M., Karras, P., Brown, D., Chang, Y. H., Debiec-Rychter, M., Adriaens, C., Radaelli, E., Wolter, P., Bechter, O., Dummer, R., Levesque, M., Piris, A., Frederick, D. T., Boland, G., Flaherty, K. T., \u0026hellip; Marine, J. C. (2018). Toward Minimal Residual Disease-Directed Therapy in Melanoma. \u003cem\u003eCell\u003c/em\u003e, \u003cem\u003e174\u003c/em\u003e(4), 843\u0026ndash;855.e19. https://doi.org/10.1016/j.cell.2018.06.025 \u003c/li\u003e\n\u003cli\u003eTsoi, J., Robert, L., Paraiso, K., Galvan, C., Sheu, K. M., Lay, J., Wong, D. J. L., Atefi, M., Shirazi, R., Wang, X., Braas, D., Grasso, C. S., Palaskas, N., Ribas, A., \u0026amp; Graeber, T. G. (2018). Multi-stage Differentiation Defines Melanoma Subtypes with Differential Vulnerability to Drug-Induced Iron-Dependent Oxidative Stress. \u003cem\u003eCancer cell\u003c/em\u003e, \u003cem\u003e33\u003c/em\u003e(5), 890\u0026ndash;904.e5. https://doi.org/10.1016/j.ccell.2018.03.017\u003c/li\u003e\n\u003cli\u003eM\u0026uuml;ller, J., Krijgsman, O., Tsoi, J., Robert, L., Hugo, W., Song, C., Kong, X., Possik, P. A., Cornelissen-Steijger, P. D., Geukes Foppen, M. H., Kemper, K., Goding, C. R., McDermott, U., Blank, C., Haanen, J., Graeber, T. G., Ribas, A., Lo, R. S., \u0026amp; Peeper, D. S. (2014). Low MITF/AXL ratio predicts early resistance to multiple targeted drugs in melanoma. \u003cem\u003eNature communications\u003c/em\u003e, \u003cem\u003e5\u003c/em\u003e, 5712. https://doi.org/10.1038/ncomms6712\u003c/li\u003e\n\u003cli\u003eKonieczkowski, D. J., Johannessen, C. M., Abudayyeh, O., Kim, J. W., Cooper, Z. A., Piris, A., Frederick, D. T., Barzily-Rokni, M., Straussman, R., Haq, R., Fisher, D. E., Mesirov, J. P., Hahn, W. C., Flaherty, K. T., Wargo, J. A., Tamayo, P., \u0026amp; Garraway, L. A. (2014). A melanoma cell state distinction influences sensitivity to MAPK pathway inhibitors. \u003cem\u003eCancer discovery\u003c/em\u003e, \u003cem\u003e4\u003c/em\u003e(7), 816\u0026ndash;827. https://doi.org/10.1158/2159-8290.CD-13-0424\u003c/li\u003e\n\u003cli\u003eTirosh, I., Izar, B., Prakadan, S. M., Wadsworth, M. H., 2nd, Treacy, D., Trombetta, J. J., Rotem, A., Rodman, C., Lian, C., Murphy, G., Fallahi-Sichani, M., Dutton-Regester, K., Lin, J. R., Cohen, O., Shah, P., Lu, D., Genshaft, A. S., Hughes, T. K., Ziegler, C. G., Kazer, S. W., \u0026hellip; Garraway, L. A. (2016). Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. \u003cem\u003eScience (New York, N.Y.)\u003c/em\u003e, \u003cem\u003e352\u003c/em\u003e(6282), 189\u0026ndash;196. https://doi.org/10.1126/science.aad0501\u003c/li\u003e\n\u003cli\u003eBoiko, A. D., Razorenova, O. V., van de Rijn, M., Swetter, S. M., Johnson, D. L., Ly, D. P., Butler, P. D., Yang, G. P., Joshua, B., Kaplan, M. J., Longaker, M. T., \u0026amp; Weissman, I. L. (2010). Human melanoma-initiating cells express neural crest nerve growth factor receptor CD271. \u003cem\u003eNature\u003c/em\u003e, \u003cem\u003e466\u003c/em\u003e(7302), 133\u0026ndash;137. https://doi.org/10.1038/nature09161\u003c/li\u003e\n\u003cli\u003eQuintana, E., Shackleton, M., Foster, H. R., Fullen, D. R., Sabel, M. S., Johnson, T. M., \u0026amp; Morrison, S. J. (2010). Phenotypic heterogeneity among tumorigenic melanoma cells from patients that is reversible and not hierarchically organized. \u003cem\u003eCancer cell\u003c/em\u003e, \u003cem\u003e18\u003c/em\u003e(5), 510\u0026ndash;523. https://doi.org/10.1016/j.ccr.2010.10.012\u003c/li\u003e\n\u003cli\u003eMehta, A., Kim, Y. J., Robert, L., Tsoi, J., Comin-Anduix, B., Berent-Maoz, B., Cochran, A. J., Economou, J. S., Tumeh, P. C., Puig-Saus, C., \u0026amp; Ribas, A. (2018). Immunotherapy Resistance by Inflammation-Induced Dedifferentiation. \u003cem\u003eCancer discovery\u003c/em\u003e, \u003cem\u003e8\u003c/em\u003e(8), 935\u0026ndash;943. https://doi.org/10.1158/2159-8290.CD-17-1178\u003c/li\u003e\n\u003cli\u003eLandsberg, J., Kohlmeyer, J., Renn, M., Bald, T., Rogava, M., Cron, M., Fatho, M., Lennerz, V., W\u0026ouml;lfel, T., H\u0026ouml;lzel, M., \u0026amp; T\u0026uuml;ting, T. (2012). Melanomas resist T-cell therapy through inflammation-induced reversible dedifferentiation. \u003cem\u003eNature\u003c/em\u003e, \u003cem\u003e490\u003c/em\u003e(7420), 412\u0026ndash;416. https://doi.org/10.1038/nature11538\u003c/li\u003e\n\u003cli\u003eFallahi-Sichani, M., Becker, V., Izar, B., Baker, G. J., Lin, J. R., Boswell, S. A., Shah, P., Rotem, A., Garraway, L. A., \u0026amp; Sorger, P. K. (2017). Adaptive resistance of melanoma cells to RAF inhibition via reversible induction of a slowly dividing de-differentiated state. \u003cem\u003eMolecular systems biology\u003c/em\u003e, \u003cem\u003e13\u003c/em\u003e(1), 905. https://doi.org/10.15252/msb.20166796\u003c/li\u003e\n\u003cli\u003eSu, Y., Wei, W., Robert, L., Xue, M., Tsoi, J., Garcia-Diaz, A., Homet Moreno, B., Kim, J., Ng, R. H., Lee, J. W., Koya, R. C., Comin-Anduix, B., Graeber, T. G., Ribas, A., \u0026amp; Heath, J. R. (2017). Single-cell analysis resolves the cell state transition and signaling dynamics associated with melanoma drug-induced resistance. \u003cem\u003eProceedings of the National Academy of Sciences of the United States of America\u003c/em\u003e, \u003cem\u003e114\u003c/em\u003e(52), 13679\u0026ndash;13684. https://doi.org/10.1073/pnas.1712064115\u003c/li\u003e\n\u003cli\u003eRestivo, G., Diener, J., Cheng, P. F., Kiowski, G., Bonalli, M., Biedermann, T., Reichmann, E., Levesque, M. P., Dummer, R., \u0026amp; Sommer, L (2017)\u003cem\u003e.\u003c/em\u003e The low affinity neurotrophin receptor CD271 regulates phenotype switching in melanoma. \u003cem\u003eNat Commun\u003c/em\u003e \u003cem\u003e8\u003c/em\u003e, 1988. https://doi.org/10.1038/s41467-017-01573-6 \u003c/li\u003e\n\u003cli\u003eBoshuizen, J., Vredevoogd, D. W., Krijgsman, O., Ligtenberg, M. A., Blankenstein, S., de Bruijn, B., Frederick, D. T., Kenski, J. C. N., Parren, M., Br\u0026uuml;ggemann, M., Madu, M. F., Rozeman, E. A., Song, J. Y., Horlings, H. M., Blank, C. U., van Akkooi, A. C. J., Flaherty, K. T., Boland, G. M., \u0026amp; Peeper, D. S. (2020)\u003cem\u003e.\u003c/em\u003e Reversal of pre-existing NGFR-driven tumor and immune therapy resistance. \u003cem\u003eNat Commun\u003c/em\u003e \u003cem\u003e11\u003c/em\u003e, 3946. https://doi.org/10.1038/s41467-020-17739-8 \u003c/li\u003e\n\u003cli\u003eLiu, D., Lin, J. R., Robitschek, E. J., Kasumova, G. G., Heyde, A., Shi, A., Kraya, A., Zhang, G., Moll, T., Frederick, D. T., Chen, Y. A., Wang, S., Schapiro, D., Ho, L. L., Bi, K., Sahu, A., Mei, S., Miao, B., Sharova, T., Alvarez-Breckenridge, C., \u0026hellip; Boland, G. M. (2021). Evolution of delayed resistance to immunotherapy in a melanoma responder. \u003cem\u003eNature medicine\u003c/em\u003e, \u003cem\u003e27\u003c/em\u003e(6), 985\u0026ndash;992. https://doi.org/10.1038/s41591-021-01331-8\u003c/li\u003e\n\u003cli\u003eRiesenberg, S., Groetchen, A., Siddaway, R., Bald, T., Reinhardt, J., Smorra, D., Kohlmeyer, J., Renn, M., Phung, B., Aymans, P., Schmidt, T., Hornung, V., Davidson, I., Goding, C. R., J\u0026ouml;nsson, G., Landsberg, J., T\u0026uuml;ting, T., \u0026amp; H\u0026ouml;lzel, M. (2015). MITF and c-Jun antagonism interconnects melanoma dedifferentiation with pro-inflammatory cytokine responsiveness and myeloid cell recruitment. \u003cem\u003eNature communications\u003c/em\u003e, \u003cem\u003e6\u003c/em\u003e, 8755. https://doi.org/10.1038/ncomms9755\u003c/li\u003e\n\u003cli\u003eKim, Y. J., Sheu, K. M., Tsoi, J., Abril-Rodriguez, G., Medina, E., Grasso, C. S., Torrejon, D. Y., Champhekar, A. S., Litchfield, K., Swanton, C., Speiser, D. E., Scumpia, P. O., Hoffmann, A., Graeber, T. G., Puig-Saus, C., \u0026amp; Ribas, A. (2021). Melanoma dedifferentiation induced by IFN-\u0026gamma; epigenetic remodeling in response to anti-PD-1 therapy. \u003cem\u003eThe Journal of clinical investigation\u003c/em\u003e, \u003cem\u003e131\u003c/em\u003e(12), e145859. https://doi.org/10.1172/JCI145859\u003c/li\u003e\n\u003cli\u003eFuruta, J., Inozume, T., Harada, K., \u0026amp; Shimada, S. (2014). CD271 on melanoma cell is an IFN-\u0026gamma;-inducible immunosuppressive factor that mediates downregulation of melanoma antigens. \u003cem\u003eThe Journal of investigative dermatology\u003c/em\u003e, \u003cem\u003e134\u003c/em\u003e(5), 1369\u0026ndash;1377. https://doi.org/10.1038/jid.2013.490\u003c/li\u003e\n\u003cli\u003eLehmann, J., Caduff, N., Krzywińska, E., Stierli, S., Salas-Bastos, A., Loos, B., Levesque, M. P., Dummer, R., Stockmann, C., M\u0026uuml;nz, C., Diener, J., \u0026amp; Sommer, L. (2023). Escape from NK cell tumor surveillance by NGFR-induced lipid remodeling in melanoma. \u003cem\u003eScience advances\u003c/em\u003e, \u003cem\u003e9\u003c/em\u003e(2), eadc8825. https://doi.org/10.1126/sciadv.adc8825\u003c/li\u003e\n\u003cli\u003eSchieven, S. M., Traets, J. J. H., Vliet, A. V., Baalen, M. V., Song, J. Y., Guimaraes, M. D. S., Kuilman, T., \u0026amp; Peeper, D. S. (2023). The Elongin BC Complex Negatively Regulates AXL and Marks a Differentiated Phenotype in Melanoma. \u003cem\u003eMolecular Cancer Research : MCR\u003c/em\u003e, \u003cem\u003e21\u003c/em\u003e(5), 428\u0026ndash;443. https://doi.org/10.1158/1541-7786.MCR-22-0648\u003c/li\u003e\n\u003cli\u003eHao, Y., Hao, S., Andersen-Nissen, E., Mauck, W. M., 3rd, Zheng, S., Butler, A., Lee, M. J., Wilk, A. J., Darby, C., Zager, M., Hoffman, P., Stoeckius, M., Papalexi, E., Mimitou, E. P., Jain, J., Srivastava, A., Stuart, T., Fleming, L. M., Yeung, B., Rogers, A. J., \u0026hellip; Satija, R. (2021). Integrated analysis of multimodal single-cell data. \u003cem\u003eCell\u003c/em\u003e, \u003cem\u003e184\u003c/em\u003e(13), 3573\u0026ndash;3587.e29. https://doi.org/10.1016/j.cell.2021.04.048\u003c/li\u003e\n\u003cli\u003eMcGinnis, C. S., Murrow, L. M., \u0026amp; Gartner, Z. J. (2019). DoubletFinder: Doublet Detection in Single-Cell RNA Sequencing Data Using Artificial Nearest Neighbors. \u003cem\u003eCell systems\u003c/em\u003e, \u003cem\u003e8\u003c/em\u003e(4), 329\u0026ndash;337.e4. https://doi.org/10.1016/j.cels.2019.03.003\u003c/li\u003e\n\u003cli\u003eChoudhary, S., \u0026amp; Satija, R. (2022). Comparison and evaluation of statistical error models for scRNA-seq. \u003cem\u003eGenome biology\u003c/em\u003e, \u003cem\u003e23\u003c/em\u003e(1), 27. https://doi.org/10.1186/s13059-021-02584-9\u003c/li\u003e\n\u003cli\u003eJerby-Arnon, L., Shah, P., Cuoco, M. S., Rodman, C., Su, M. J., Melms, J. C., Leeson, R., Kanodia, A., Mei, S., Lin, J. R., Wang, S., Rabasha, B., Liu, D., Zhang, G., Margolais, C., Ashenberg, O., Ott, P. A., Buchbinder, E. I., Haq, R., Hodi, F. S., \u0026hellip; Regev, A. (2018). A Cancer Cell Program Promotes T Cell Exclusion and Resistance to Checkpoint Blockade. \u003cem\u003eCell\u003c/em\u003e, \u003cem\u003e175\u003c/em\u003e(4), 984\u0026ndash;997.e24. https://doi.org/10.1016/j.cell.2018.09.006 \u003c/li\u003e\n\u003cli\u003eAibar, S., Gonz\u0026aacute;lez-Blas, C. B., Moerman, T., Huynh-Thu, V. A., Imrichova, H., Hulselmans, G., Rambow, F., Marine, J. C., Geurts, P., Aerts, J., van den Oord, J., Atak, Z. K., Wouters, J., \u0026amp; Aerts, S. (2017). SCENIC: single-cell regulatory network inference and clustering. \u003cem\u003eNature methods\u003c/em\u003e, \u003cem\u003e14\u003c/em\u003e(11), 1083\u0026ndash;1086. https://doi.org/10.1038/nmeth.4463\u003c/li\u003e\n\u003cli\u003ePozniak, J., Pedri, D., Landeloos, E., Van Herck, Y., Antoranz, A., Vanwynsberghe, L., Nowosad, A., Roda, N., Makhzami, S., Bervoets, G., Maciel, L. F., Pulido-Vicu\u0026ntilde;a, C. A., Pollaris, L., Seurinck, R., Zhao, F., Flem-Karlsen, K., Damsky, W., Chen, L., Karagianni, D., Cinque, S., \u0026hellip; Marine, J. C. (2024). A TCF4-dependent gene regulatory network confers resistance to immunotherapy in melanoma. \u003cem\u003eCell\u003c/em\u003e, \u003cem\u003e187\u003c/em\u003e(1), 166\u0026ndash;183.e25. https://doi.org/10.1016/j.cell.2023.11.037\u003c/li\u003e\n\u003cli\u003eWidmer, D. S., Cheng, P. F., Eichhoff, O. M., Belloni, B. C., Zipser, M. C., Schlegel, N. C., Javelaud, D., Mauviel, A., Dummer, R., \u0026amp; Hoek, K. S. (2012). Systematic classification of melanoma cells by phenotype-specific gene expression mapping. \u003cem\u003ePigment cell \u0026amp; melanoma research\u003c/em\u003e, \u003cem\u003e25\u003c/em\u003e(3), 343\u0026ndash;353. https://doi.org/10.1111/j.1755-148X.2012.00986.x\u003c/li\u003e\n\u003cli\u003eHugo, W., Zaretsky, J. M., Sun, L., Song, C., Moreno, B. H., Hu-Lieskovan, S., Berent-Maoz, B., Pang, J., Chmielowski, B., Cherry, G., Seja, E., Lomeli, S., Kong, X., Kelley, M. C., Sosman, J. A., Johnson, D. B., Ribas, A., \u0026amp; Lo, R. S. (2016). Genomic and Transcriptomic Features of Response to Anti-PD-1 Therapy in Metastatic Melanoma. \u003cem\u003eCell\u003c/em\u003e, \u003cem\u003e165\u003c/em\u003e(1), 35\u0026ndash;44. https://doi.org/10.1016/j.cell.2016.02.065\u003c/li\u003e\n\u003cli\u003eRiaz, N., Havel, J. J., Makarov, V., Desrichard, A., Urba, W. J., Sims, J. S., Hodi, F. S., Mart\u0026iacute;n-Algarra, S., Mandal, R., Sharfman, W. H., Bhatia, S., Hwu, W. J., Gajewski, T. F., Slingluff, C. L., Jr, Chowell, D., Kendall, S. M., Chang, H., Shah, R., Kuo, F., Morris, L. G. T., \u0026hellip; Chan, T. A. (2017). Tumor and Microenvironment Evolution during Immunotherapy with Nivolumab. \u003cem\u003eCell\u003c/em\u003e, \u003cem\u003e171\u003c/em\u003e(4), 934\u0026ndash;949.e16. https://doi.org/10.1016/j.cell.2017.09.028\u003c/li\u003e\n\u003cli\u003eAyers, M., Lunceford, J., Nebozhyn, M., Murphy, E., Loboda, A., Kaufman, D. R., Albright, A., Cheng, J. D., Kang, S. P., Shankaran, V., Piha-Paul, S. A., Yearley, J., Seiwert, T. Y., Ribas, A., \u0026amp; McClanahan, T. K. (2017). IFN-\u0026gamma;-related mRNA profile predicts clinical response to PD-1 blockade. \u003cem\u003eThe Journal of clinical investigation\u003c/em\u003e, \u003cem\u003e127\u003c/em\u003e(8), 2930\u0026ndash;2940. https://doi.org/10.1172/JCI91190\u003c/li\u003e\n\u003cli\u003eThompson, J. C., Davis, C., Deshpande, C., Hwang, W. T., Jeffries, S., Huang, A., Mitchell, T. C., Langer, C. J., \u0026amp; Albelda, S. M. (2020). Gene signature of antigen processing and presentation machinery predicts response to checkpoint blockade in non-small cell lung cancer (NSCLC) and melanoma. \u003cem\u003eJournal for immunotherapy of cancer\u003c/em\u003e, \u003cem\u003e8\u003c/em\u003e(2), e000974. https://doi.org/10.1136/jitc-2020-000974\u003c/li\u003e\n\u003cli\u003eChow, A., Uddin, F. Z., Liu, M., Dobrin, A., Nabet, B. Y., Mangarin, L., Lavin, Y., Rizvi, H., Tischfield, S. E., Quintanal-Villalonga, A., Chan, J. M., Shah, N., Allaj, V., Manoj, P., Mattar, M., Meneses, M., Landau, R., Ward, M., Kulick, A., Kwong, C., \u0026hellip; Rudin, C. M. (2023). The ectonucleotidase CD39 identifies tumor-reactive CD8\u003csup\u003e+\u003c/sup\u003e T cells predictive of immune checkpoint blockade efficacy in human lung cancer. \u003cem\u003eImmunity\u003c/em\u003e, \u003cem\u003e56\u003c/em\u003e(1), 93\u0026ndash;106.e6. https://doi.org/10.1016/j.immuni.2022.12.001 \u003c/li\u003e\n\u003cli\u003eLi, X., Wan, Z., Liu, X., Ou, K., \u0026amp; Yang, L. (2022). A 12-chemokine gene signature is associated with the enhanced immunogram scores and is relevant for precision immunotherapy. \u003cem\u003eMedical oncology (Northwood, London, England)\u003c/em\u003e, \u003cem\u003e39\u003c/em\u003e(4), 43. https://doi.org/10.1007/s12032-021-01635-2\u003c/li\u003e\n\u003cli\u003eVredevoogd, D. W., Kuilman, T., Ligtenberg, M. A., Boshuizen, J., Stecker, K. E., de Bruijn, B., Krijgsman, O., Huang, X., Kenski, J. C. N., Lacroix, R., Mezzadra, R., Gomez-Eerland, R., Yildiz, M., Dagidir, I., Apriamashvili, G., Zandhuis, N., van der Noort, V., Visser, N. L., Blank, C. U., Altelaar, M., \u0026hellip; Peeper, D. S. (2019). Augmenting Immunotherapy Impact by Lowering Tumor TNF Cytotoxicity Threshold. \u003cem\u003eCell\u003c/em\u003e, \u003cem\u003e178\u003c/em\u003e(3), 585\u0026ndash;599.e15. https://doi.org/10.1016/j.cell.2019.06.014\u003c/li\u003e\n\u003cli\u003eHoefsmit, E. P., van Royen, P. T., Rao, D., Stunnenberg, J. A., Dimitriadis, P., Lieftink, C., Morris, B., Rozeman, E. A., Reijers, I. L. M., Lacroix, R., Shehwana, H., Ligtenberg, M. A., Beijersbergen, R. L., Peeper, D. S., \u0026amp; Blank, C. U. (2023). Inhibitor of Apoptosis Proteins Antagonist Induces T-cell Proliferation after Cross-Presentation by Dendritic Cells. \u003cem\u003eCancer immunology research\u003c/em\u003e, \u003cem\u003e11\u003c/em\u003e(4), 450\u0026ndash;465. https://doi.org/10.1158/2326-6066.CIR-22-0494\u003c/li\u003e\n\u003cli\u003eNajem, A., Wouters, J., Krayem, M., Rambow, F., Sabbah, M., Sales, F., Awada, A., Aerts, S., Journe, F., Marine, J. C., \u0026amp; Ghanem, G. E. (2021). Tyrosine-Dependent Phenotype Switching Occurs Early in Many Primary Melanoma Cultures Limiting Their Translational Value. \u003cem\u003eFrontiers in oncology\u003c/em\u003e, \u003cem\u003e11\u003c/em\u003e, 780654. https://doi.org/10.3389/fonc.2021.780654\u003c/li\u003e\n\u003cli\u003eShaffer, S. M., Dunagin, M. C., Torborg, S. R., Torre, E. A., Emert, B., Krepler, C., Beqiri, M., Sproesser, K., Brafford, P. A., Xiao, M., Eggan, E., Anastopoulos, I. N., Vargas-Garcia, C. A., Singh, A., Nathanson, K. L., Herlyn, M., \u0026amp; Raj, A. (2017). Rare cell variability and drug-induced reprogramming as a mode of cancer drug resistance. \u003cem\u003eNature\u003c/em\u003e, \u003cem\u003e546\u003c/em\u003e(7658), 431\u0026ndash;435. https://doi.org/10.1038/nature22794\u003c/li\u003e\n\u003cli\u003eKemper, K., Krijgsman, O., Kong, X., Cornelissen-Steijger, P., Shahrabi, A., Weeber, F., van der Velden, D. L., Bleijerveld, O. B., Kuilman, T., Kluin, R. J. C., Sun, C., Voest, E. E., Ju, Y. S., Schumacher, T. N. M., Altelaar, A. F. M., McDermott, U., Adams, D. J., Blank, C. U., Haanen, J. B., \u0026amp; Peeper, D. S. (2016). BRAF(V600E) Kinase Domain Duplication Identified in Therapy-Refractory Melanoma Patient-Derived Xenografts. \u003cem\u003eCell reports\u003c/em\u003e, \u003cem\u003e16\u003c/em\u003e(1), 263\u0026ndash;277. https://doi.org/10.1016/j.celrep.2016.05.064\u003c/li\u003e\n\u003cli\u003eTitz, B., Lomova, A., Le, A., Hugo, W., Kong, X., Ten Hoeve, J., Friedman, M., Shi, H., Moriceau, G., Song, C., Hong, A., Atefi, M., Li, R., Komisopoulou, E., Ribas, A., Lo, R. S., \u0026amp; Graeber, T. G. (2016). JUN dependency in distinct early and late BRAF inhibition adaptation states of melanoma. \u003cem\u003eCell discovery\u003c/em\u003e, \u003cem\u003e2\u003c/em\u003e, 16028. https://doi.org/10.1038/celldisc.2016.28\u003c/li\u003e\n\u003cli\u003eVan de Sande, B., Flerin, C., Davie, K., De Waegeneer, M., Hulselmans, G., Aibar, S., Seurinck, R., Saelens, W., Cannoodt, R., Rouchon, Q., Verbeiren, T., De Maeyer, D., Reumers, J., Saeys, Y., \u0026amp; Aerts, S. (2020). A scalable SCENIC workflow for single-cell gene regulatory network analysis. \u003cem\u003eNature protocols\u003c/em\u003e, \u003cem\u003e15\u003c/em\u003e(7), 2247\u0026ndash;2276. https://doi.org/10.1038/s41596-020-0336-2 \u003c/li\u003e\n\u003cli\u003eTumeh, P. C., Harview, C. L., Yearley, J. H., Shintaku, I. P., Taylor, E. J., Robert, L., Chmielowski, B., Spasic, M., Henry, G., Ciobanu, V., West, A. N., Carmona, M., Kivork, C., Seja, E., Cherry, G., Gutierrez, A. J., Grogan, T. R., Mateus, C., Tomasic, G., Glaspy, J. A., \u0026hellip; Ribas, A. (2014). PD-1 blockade induces responses by inhibiting adaptive immune resistance. \u003cem\u003eNature\u003c/em\u003e, \u003cem\u003e515\u003c/em\u003e(7528), 568\u0026ndash;571. https://doi.org/10.1038/nature13954 \u003c/li\u003e\n\u003cli\u003eLitchfield, K., Reading, J. L., Puttick, C., Thakkar, K., Abbosh, C., Bentham, R., Watkins, T. B. K., Rosenthal, R., Biswas, D., Rowan, A., Lim, E., Al Bakir, M., Turati, V., Guerra-Assun\u0026ccedil;\u0026atilde;o, J. A., Conde, L., Furness, A. J. S., Saini, S. K., Hadrup, S. R., Herrero, J., Lee, S. H., \u0026hellip; Swanton, C. (2021). Meta-analysis of tumor- and T cell-intrinsic mechanisms of sensitization to checkpoint inhibition. \u003cem\u003eCell\u003c/em\u003e, \u003cem\u003e184\u003c/em\u003e(3), 596\u0026ndash;614.e14. https://doi.org/10.1016/j.cell.2021.01.002\u003c/li\u003e\n\u003cli\u003eSpranger, S., Bao, R., \u0026amp; Gajewski, T. F. (2015). Melanoma-intrinsic \u0026beta;-catenin signalling prevents anti-tumour immunity. \u003cem\u003eNature\u003c/em\u003e, \u003cem\u003e523\u003c/em\u003e(7559), 231\u0026ndash;235. https://doi.org/10.1038/nature14404\u003c/li\u003e\n\u003cli\u003eJi, R. R., Chasalow, S. D., Wang, L., Hamid, O., Schmidt, H., Cogswell, J., Alaparthy, S., Berman, D., Jure-Kunkel, M., Siemers, N. O., Jackson, J. R., \u0026amp; Shahabi, V. (2012). An immune-active tumor microenvironment favors clinical response to ipilimumab. \u003cem\u003eCancer immunology, immunotherapy: CII\u003c/em\u003e, \u003cem\u003e61\u003c/em\u003e(7), 1019\u0026ndash;1031. https://doi.org/10.1007/s00262-011-1172-6\u003c/li\u003e\n\u003cli\u003eSha, D., Jin, Z., Budczies, J., Kluck, K., Stenzinger, A., \u0026amp; Sinicrope, F. A. (2020). Tumor Mutational Burden as a Predictive Biomarker in Solid Tumors. \u003cem\u003eCancer discovery\u003c/em\u003e, \u003cem\u003e10\u003c/em\u003e(12), 1808\u0026ndash;1825. https://doi.org/10.1158/2159-8290.CD-20-0522\u003c/li\u003e\n\u003cli\u003eMarin-Bejar, O., Rogiers, A., Dewaele, M., Femel, J., Karras, P., Pozniak, J., Bervoets, G., Van Raemdonck, N., Pedri, D., Swings, T., Demeulemeester, J., Borght, S. V., Lehnert, S., Bosisio, F., van den Oord, J. J., Bempt, I. V., Lambrechts, D., Voet, T., Bechter, O., Rizos, H., \u0026hellip; Marine, J. C. (2021). Evolutionary predictability of genetic versus nongenetic resistance to anticancer drugs in melanoma. \u003cem\u003eCancer cell\u003c/em\u003e, \u003cem\u003e39\u003c/em\u003e(8), 1135\u0026ndash;1149.e8. https://doi.org/10.1016/j.ccell.2021.05.015\u003c/li\u003e\n\u003cli\u003eLauss, M., Phung, B., Borch, T. H., Harbst, K., Kaminska, K., Ebbesson, A., Hedenfalk, I., Yuan, J., Nielsen, K., Ingvar, C., Carneiro, A., Isaksson, K., Pietras, K., Svane, I. M., Donia, M., \u0026amp; J\u0026ouml;nsson, G. (2024). Molecular patterns of resistance to immune checkpoint blockade in melanoma. \u003cem\u003eNature communications\u003c/em\u003e, \u003cem\u003e15\u003c/em\u003e(1), 3075. https://doi.org/10.1038/s41467-024-47425-y\u003c/li\u003e\n\u003cli\u003eBoe, R. H., Triandafillou, C. G., Lazcano, R., Wargo, J. A., \u0026amp; Raj, A. (2024). Spatial transcriptomics reveals influence of microenvironment on intrinsic fates in melanoma therapy resistance. \u003cem\u003ebioRxiv: the preprint server for biology\u003c/em\u003e, 2024.06.30.601416. https://doi.org/10.1101/2024.06.30.601416\u003c/li\u003e\n\u003cli\u003eBai, X., Fisher, D. E., \u0026amp; Flaherty, K. T. (2019). Cell-state dynamics and therapeutic resistance in melanoma from the perspective of MITF and IFN\u0026gamma; pathways. \u003cem\u003eNature reviews. Clinical oncology\u003c/em\u003e, \u003cem\u003e16\u003c/em\u003e(9), 549\u0026ndash;562. https://doi.org/10.1038/s41571-019-0204-6\u003c/li\u003e\n\u003cli\u003eKarras, P., Black, J. R. M., McGranahan, N., \u0026amp; Marine, J. C. (2024). Decoding the interplay between genetic and non-genetic drivers of metastasis. \u003cem\u003eNature\u003c/em\u003e, \u003cem\u003e629\u003c/em\u003e(8012), 543\u0026ndash;554. https://doi.org/10.1038/s41586-024-07302-6\u003c/li\u003e\n\u003cli\u003eWolchok, J. D., Chiarion-Sileni, V., Gonzalez, R., Rutkowski, P., Grob, J. J., Cowey, C. L., Lao, C. D., Wagstaff, J., Schadendorf, D., Ferrucci, P. F., Smylie, M., Dummer, R., Hill, A., Hogg, D., Haanen, J., Carlino, M. S., Bechter, O., Maio, M., Marquez-Rodas, I., Guidoboni, M., \u0026hellip; Larkin, J. (2017). Overall Survival with Combined Nivolumab and Ipilimumab in Advanced Melanoma. \u003cem\u003eThe New England journal of medicine\u003c/em\u003e, \u003cem\u003e377\u003c/em\u003e(14), 1345\u0026ndash;1356. https://doi.org/10.1056/NEJMoa1709684 \u003c/li\u003e\n\u003cli\u003e57. Larkin, J., Chiarion-Sileni, V., Gonzalez, R., Grob, J. J., Cowey, C. L., Lao, C. D., Schadendorf, D., Dummer, R., Smylie, M., Rutkowski, P., Ferrucci, P. F., Hill, A., Wagstaff, J., Carlino, M. S., Haanen, J. B., Maio, M., Marquez-Rodas, I., McArthur, G. A., Ascierto, P. A., Long, G. V., \u0026hellip; Wolchok, J. D. (2015). Combined Nivolumab and Ipilimumab or Monotherapy in Untreated Melanoma. \u003cem\u003eThe New England journal of medicine\u003c/em\u003e, \u003cem\u003e373\u003c/em\u003e(1), 23\u0026ndash;34. https://doi.org/10.1056/NEJMoa1504030 \u003c/li\u003e\n\u003cli\u003eHodi, F. S., O\u0026apos;Day, S. J., McDermott, D. F., Weber, R. W., Sosman, J. A., Haanen, J. B., Gonzalez, R., Robert, C., Schadendorf, D., Hassel, J. C., Akerley, W., van den Eertwegh, A. J., Lutzky, J., Lorigan, P., Vaubel, J. M., Linette, G. P., Hogg, D., Ottensmeier, C. H., Lebb\u0026eacute;, C., Peschel, C., \u0026hellip; Urba, W. J. (2010). Improved survival with ipilimumab in patients with metastatic melanoma. \u003cem\u003eThe New England journal of medicine\u003c/em\u003e, \u003cem\u003e363\u003c/em\u003e(8), 711\u0026ndash;723. https://doi.org/10.1056/NEJMoa1003466\u003c/li\u003e\n\u003cli\u003eSharma, P., Hu-Lieskovan, S., Wargo, J. A., \u0026amp; Ribas, A. (2017). Primary, Adaptive, and Acquired Resistance to Cancer Immunotherapy. \u003cem\u003eCell\u003c/em\u003e, \u003cem\u003e168\u003c/em\u003e(4), 707\u0026ndash;723. https://doi.org/10.1016/j.cell.2017.01.017\u003c/li\u003e\n\u003cli\u003eKalbasi, A., \u0026amp; Ribas, A. (2020). Tumour-intrinsic resistance to immune checkpoint blockade. \u003cem\u003eNature reviews. Immunology\u003c/em\u003e, \u003cem\u003e20\u003c/em\u003e(1), 25\u0026ndash;39. https://doi.org/10.1038/s41577-019-0218-4\u003c/li\u003e\n\u003cli\u003eVersluis, J. M., Thommen, D. S., \u0026amp; Blank, C. U. (2020). Rationalizing the pathway to personalized neoadjuvant immunotherapy: the Lombard Street Approach. \u003cem\u003eJournal for immunotherapy of cancer\u003c/em\u003e, \u003cem\u003e8\u003c/em\u003e(2), e001352. https://doi.org/10.1136/jitc-2020-001352\u003c/li\u003e\n\u003cli\u003eRozeman, E. A., Hoefsmit, E. P., Reijers, I. L. M., Saw, R. P. M., Versluis, J. M., Krijgsman, O., Dimitriadis, P., Sikorska, K., van de Wiel, B. A., Eriksson, H., Gonzalez, M., Torres Acosta, A., Grijpink-Ongering, L. G., Shannon, K., Haanen, J. B. A. G., Stretch, J., Ch\u0026apos;ng, S., Nieweg, O. E., Mallo, H. A., Adriaansz, S., \u0026hellip; Blank, C. U. (2021). Survival and biomarker analyses from the OpACIN-neo and OpACIN neoadjuvant immunotherapy trials in stage III melanoma. \u003cem\u003eNature medicine\u003c/em\u003e, \u003cem\u003e27\u003c/em\u003e(2), 256\u0026ndash;263. https://doi.org/10.1038/s41591-020-01211-7\u003c/li\u003e\n\u003cli\u003eZila, N., Eichhoff, O. M., Steiner, I., Mohr, T., Bileck, A., Cheng, P. F., Leitner, A., Gillet, L., Sajic, T., Goetze, S., Friedrich, B., Bortel, P., Strobl, J., Reitermaier, R., Hogan, S. A., Mart\u0026iacute;nez G\u0026oacute;mez, J. M., Staeger, R., Tuchmann, F., Peters, S., Stary, G., \u0026hellip; Paulitschke, V. (2024). Proteomic Profiling of Advanced Melanoma Patients to Predict Therapeutic Response to Anti-PD-1 Therapy. \u003cem\u003eClinical cancer research : an official journal of the American Association for Cancer Research\u003c/em\u003e, \u003cem\u003e30\u003c/em\u003e(1), 159\u0026ndash;175. https://doi.org/10.1158/1078-0432.CCR-23-0562\u003cem\u003e\u003cbr\u003e\u003c/em\u003e\u003c/li\u003e\n\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":true,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"NGFR heterogeneity, spatial transcriptomics, PDGFR, cytokine resistance and immune checkpoint blockade resistance","lastPublishedDoi":"10.21203/rs.3.rs-6506453/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-6506453/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eHuman melanomas dedifferentiate into a neural crest-like cell state when exposed to T cell cytokines or MAPK pathway inhibitors. This transformation is associated with cellular heterogeneity and the emergence of small, therapy-resistant melanoma populations characterized by elevated nerve growth factor receptor (NGFR) expression. However, the extent of this heterogeneity and its impact on immunotherapy response remain unclear. By dissecting intratumor heterogeneity in patient melanomas, we show here that even within NGFR\u003csup\u003e+\u003c/sup\u003e tumor subpopulations, remarkable phenotypic and functional diversity exists. Combined single-cell RNA sequencing (scRNA-seq) of NGFR\u003csup\u003e+\u003c/sup\u003e fractions and spatial transcriptomics uncovered pronounced diversity among single-cell clusters, characterized by patient-specific gene regulatory networks (GRNs) and distinct spatial organization. Furthermore, we identify an NGFR\u003csup\u003e+\u003c/sup\u003e subpopulation marked by co-expression of Platelet-Derived Growth Factor Receptor (PDGFR), which is associated with increased resistance to the T cell cytokines IFNg and TNF. Clinically corroborating these findings, we observed that \u003cem\u003eNGFR\u003c/em\u003e\u003csup\u003e\u003cem\u003e+\u003c/em\u003e\u003c/sup\u003e\u003cem\u003e/PDGFR\u003c/em\u003e\u003csup\u003e\u003cem\u003e+\u003c/em\u003e\u003c/sup\u003e mesenchymal-like cells are enriched in melanomas infiltrated with active T cells yet failing to respond to immune checkpoint blockade treatment. Our results highlight extreme heterogeneity within human melanoma, which is spatially organized and regulated by patient-specific GRNs, and harboring a distinct subfraction linked to immunotherapy resistance.\u003c/p\u003e","manuscriptTitle":"Intratumoral Heterogeneity in Ngfr+ Melanoma Subpopulations Shapes Immune Evasion and Immunotherapy Resistance","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2025-05-09 11:00:30","doi":"10.21203/rs.3.rs-6506453/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"ab4c11fb-28cd-4a9c-8a35-cfe1bebf6b94","owner":[],"postedDate":"May 9th, 2025","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[{"id":48132806,"name":"Biological sciences/Cancer/Skin cancer/Melanoma"},{"id":48132807,"name":"Biological sciences/Computational biology and bioinformatics/Gene regulatory networks"},{"id":48132808,"name":"Biological sciences/Immunology/Immune evasion"}],"tags":[],"updatedAt":"2025-06-13T16:41:13+00:00","versionOfRecord":[],"versionCreatedAt":"2025-05-09 11:00:30","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-6506453","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-6506453","identity":"rs-6506453","version":["v1"]},"buildId":"8U1c8b4HqxoKbykW_rLl7","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}