High-Plex Spatial RNA Proﬁling Reveals Cell Type ‒ Speciﬁc Biomarker Expression during Melanoma Development

6 Early diagnosis of melanoma is critical for improved survival. However, the biomarkers of early melanoma evolution and their origin within the tumor and its microenvironment, including the keratinocytes, are poorly deﬁned. To address this, we used spatial transcript proﬁling that maintains the morphological tumor context to measure the expression of > 1,000 RNAs in situ in patient-derived formalin-ﬁxed, parafﬁn-embedded tissue sections in primary melanoma and melanocytic nevi. We proﬁled 134 regions of interest (each 200 m m in diameter) enriched in melanocytes, neighboring keratinocytes, or immune cells. This approach captured distinct expression patterns across cell types and tumor types during melanoma development. Unexpectedly, we discovered that S100A8 is expressed by keratinocytes within the tumor microenvironment during melanoma growth. Immunohistochemistry of 252 tumors showed prominent keratinocyte-derived S100A8 expression in melanoma but not in benign tumors and conﬁrmed the same pattern for S100A8’s binding partner S100A9, suggesting that injury to the epidermis may be an early and readily detectable indicator of melanoma development. Together, our results establish a framework for high-plex, spatial, and cell type ‒ speciﬁc resolution of gene expression in archival tissue applicable to the development of biomarkers and characterization of tumor microenvironment interactions in tumor evolution.


INTRODUCTION
Melanoma, the deadliest of the common skin cancers, is curable with early diagnosis and treatment (Gershenwald et al., 2017). However, histopathologic diagnosis of melanoma can be complicated by morphological mimicry, especially in its early forms, by a subset of melanocytic nevi. Because the development of melanoma is a stepwise process in which melanocytes accrue mutations and escape environmental controls on proliferation (Villanueva and Herlyn, 2008), understanding the interaction of melanocytes with neighboring cell types is crucial to the development of diagnostic tools and effective treatments.
Many melanoma-associated genes have been identified (Bastian, 2014;Charbel et al., 2014;Roh et al., 2015;Shain and Bastian, 2016;Shain et al., 2015), and molecular tests for diagnosis and prognosis of melanoma are gradually being introduced (Clarke et al., 2015;Gerami et al., 2015aGerami et al., , 2015b, but markers of early melanoma development, including within the tumor microenvironment, remain lacking. In addition, although the treatment of metastatic melanoma has changed drastically since the development of immune checkpoint inhibitor therapies (Khair et al., 2019), biomarkers predicting durable treatment response are largely unknown. Given the heterogeneity, low cellularity, and spatial context of immune and microenvironment responses (Finotello and Eduati, 2018), spatially resolved techniques are likely to outperform bulk molecular profiling for the discovery of early-stage and predictive biomarkers.
Previous studies have revealed the importance of keratinocyte (KC)-derived growth factors and cell adhesion molecules in limiting melanocyte proliferation in normal skin and elucidated the mechanisms by which malignant melanocytes escape this regulation (Haass et al., 2005;Villanueva and Herlyn, 2008). However, these experiments relied on the use of co-culture systems or heterologous expression of KC-derived genes in melanocytes, neither of which capture the spatial element of melanocyte-KC interactions in situ. Furthermore, single-cell RNA sequencing (RNA-seq) studies on melanoma have largely focused on melanoma metastases, overlooking the KC microenvironment of primary melanomas (Jerby-Arnon et al., 2018;Tirosh et al., 2016). Because single-cell RNA-seq relies on fresh tissue, studies on benign melanocytic tumors in humans are also lacking.
To address the challenges outlined above, we took a discovery-based approach to studying tumor microenvironment interactions during melanoma evolution within the native morphological context of the tumor, including within the KC microenvironment. Specifically, we examined the expression of over 1,000 genes (Supplementary Table S1) in 134 regions of interest (ROIs) enriched for melanocytes, neighboring KCs, or immune cells in patient-derived formalin-fixed, paraffin-embedded (FFPE) tissue sections from 12 melanocytic tumors, ranging from benign to malignant, using the NanoString GeoMx Digital Spatial Profiler (DSP) (Beechem, 2020) ( Figure 1). We discovered that the damage-associated molecular pattern (DAMP) S100A8, which is a known melanoma marker (Clarke et al., 2015) thought to be expressed by immune cells (Wagner et al., 2019), is KC derived in melanoma and confirmed this finding by immunohistochemistry (IHC) in a cohort of 252 melanomas.

Cell type and case type both influence the expression profile of each ROI
We validated the technical performance of DSP on one FFPE tissue from each of four case types (common nevi [CNs], dysplastic nevi [DNs], melanoma in situ [MIS], and invasive melanoma [MM]), observing high reproducibility across replicate experiments and expected gene expression patterns across ROI types using a 108 gene test panel (1,048 probes, 10 probes per gene for most genes) (Supplementary Figure S1). Next, we used a 1,412 gene panel (4,998 probes, 3 probes per gene for most genes) to profile three tumors per case type. Two of the three melanoma cases were stage pT1a, and the third was stage pT2b. We selected 6-16 circular ROIs (200 mm diameter) in each tissue, with 134 ROIs in total (Supplementary Figure S2). One KC-rich ROI on the periphery of each tissue was chosen to serve as normal control ( Figure 2a). To facilitate comparison across tissues, we assigned all other ROIs to one of five categories (immune rich, melanocyte rich, immune/melanocyte rich, KC/melanocyte rich, or mixed) ( Figure 2a) and evaluated the cell type composition of each ROI (Figure 2b). Raw counts were normalized to the upper quartile (Bullard et al., 2010) for each ROI to enable comparison across ROIs.
To further validate the performance of DSP, we examined the expression of known melanomagenesis-associated genes across the cohort. We observed that the melanoma biomarker PRAME (preferentially expressed antigen in melanoma) (Clarke et al., 2015;Gerami et al., 2017) and the melanomaassociated developmental marker PMEL (Sarantou et al., 1997) were significantly and specifically elevated in melanocyte-containing ROIs in the three MM cases (Figure 2c). The keratin gene KRT14 and the chemokine gene CXCL14 (Riker et al., 2008) were elevated in KC-containing ROIs relative to melanocyte-rich or immune-rich ROIs in all the four tumor types and were further upregulated in the melanomas (Figure 2c). Finally, the pan-leukocyte marker PTPRC (CD45) and the T-cell chemoattractant gene CXCL9, which are included in a gene expression-based diagnostic test for melanoma (Clarke et al., 2015), were detected specifically in immune-rich ROIs, particularly in melanomas ( Figure 2c). Together, these results validated previous data and also revealed cell type-specific expression of known melanoma biomarkers.
Unbiased clustering of ROIs based on pairwise correlation coefficients revealed that cell type and tumor type both affect the similarity between ROIs (Figure 2d). The ROIs clustered into five groups, two of which consisted of melanocyte-rich ROIs (M1 and M2; Figure 2d), two consisted of immunerich ROIs (I1 and I2; Figure 2d), and one consisted of KC-rich ROIs (K; Figure 2d). For the M and I clusters, the differences between groups 1 and 2 were based on case type: M1 and I1 contained only melanoma or MIS ROIs, whereas M2 and I2 had no melanoma ROIs (Figure 2d). The K cluster contained all the 12 control KC ROIs and most of the mixed KC/melanocyte ROIs, suggesting that KC-specific genes drive the clustering of these ROIs.

A framework for identification of cell type-specific gene expression
Because the I1 and M1 clusters were composed entirely of ROIs from the melanoma or MIS cases, the genes expressed only in these clusters have biomarker potential. Importantly, because ROIs were selected on the basis of enrichment of certain cell populations, we were able to perform a cell typespecific analysis. We performed linear regression to identify the genes that were (i) significantly enriched in M1 or I1 ROIs compared with those in other ROIs containing the same cell type(s) and (ii) not detected in any ROIs from CN. For M1, the known biomarker PRAME was by far the most significantly enriched gene (Figure 3a). The other top hits LEF1, CD276, BCL2A1, and SLC7A5 have been shown to be upregulated or amplified in melanoma but are also detectable in benign nevi or normal melanocytes (Haq et al., 2013;Tekle et al., 2012;Wang et al., 2014;Xu et al., 2015). For the immune ROIs in I1, the two most highly enriched genes in melanoma that were also not detected in nevi were PTPRC and CXCL9 (Figure 3b), followed by genes broadly related to leukocyte biology (Figure 3b), including major histocompatibility complex class II antigen presentation (HLA-D, CTSS).
Enrichment analysis does not provide information on the relationship of genes with each other, so to further leverage the spatial resolution of DSP data, we used the dimensionality reduction technique Uniform Manifold Approximation and Projection (UMAP) to visualize the relative spatial expression profiles of all genes (Figure 3c and d). This approach can be used to determine which genes tend to be expressed in the same location (or not), a property that is lost in bulk measurements such as RNA-seq. Importantly, marker genes with different spatial expression profiles can provide orthogonal information, whereas genes with the same expression profile are more likely to be reporting on similar biology. We chose UMAP over other methods such as tdistributed stochastic neighbor embedding (t-SNE) because it preserves both local and global data structure and captures meaningful biological relationships (Becht et al., 2019). At the global level, UMAP analysis produced three gene clusters on the left side of the plot space, and we validated the UMAP output by confirming that gene clustering correlated with cell Panel content is approximately 35% immune related, 40% tumor-related, and 20% microenvironment related, with 1% housekeeping genes and 3% negative probes. (c) Schematic of probe design for DSP with NGS readout. Each probe contains an antisense sequence that hybridizes to target mRNA (green), a photocleavable linker (purple), an RNA ID that identifies the mRNA target (red), and a unique molecular identifier (blue) to allow removal of PCR duplicates when converting reads to digital counts. DSP probe pools target each gene with 1-10 probes that hybridize to different sequences along the mRNA transcript and contain >80 negative probes that target scrambled or nonhuman sequences. (d) Workflow for DSP with NGS readout. Collected oligos are PCR amplified using indexing primers to preserve ROI identity, pooled, purified, and sequenced. (e) Illustration of ROI selection process. Top images: ROIs selected by a pathologist on the basis of enrichment for melanocytes, keratinocytes, or immune cells in H&E section. Bottom images: ROIs collected from a serial section during DSP. Fluorescent antibodies to melanocyte markers S100B and PMEL, T-cell marker CD3, lymphocyte marker CD45, and DNA stain SYTO 13 were used as visualization markers during DSP to guide the matching of ROIs to the H&E sections. Bar ¼ 0.2 mm. CN, common nevus; CT, cancer/testis; DN, dysplastic nevus; DSP, Digital Spatial Profiler; ECM, extracellular matrix; EMT, epithelial to mesenchymal transition; FFPE, formalin-fixed, paraffin-embedded; ID, identifier; MIS, melanoma in situ; MM, invasive melanoma; NGS, next-generation sequencing; PI3K, phosphatidylinositol 3 kinase; ROI, region of interest; TLR, toll-like receptor; UMI, unique molecular identifier. Notably, the expression profiles of PMEL, CTNNB1, LDHB, and CDK2 were more similar to that of PRAME compared with the expression profile of the genes enriched in melanoma melanocytes (Figure 3c and d), but of these genes, only CDK2 was not detected in nevi (Figure 3e). Genes enriched in melanoma immune ROIs were not tightly UMAP adjacent ( Figure 3c); instead, they clustered next to genes expressed by the same cell type or within the same pathway ( Figure 3d): for example, PTPRC and LCP1 are expressed by lymphocytes, and HLA-DMA and HLA-DQA1/2 are expressed by antigenpresenting cells, while CXCL9 clustered with other IFN-gstimulated genes such as GBP1. All of these genes were primarily detected in MIS or MM (Figure 3f). The gene signatures shown in Figure 3e and f were both predictive of cell type and case type (Supplementary Figure S3). Together, these results established a framework for the identification of cell type-specific gene expression during melanoma evolution.

Analysis of intermediate ROI clusters reveals the components of the KC microenvironment in melanoma
Because the K1/M2 clusters encompass most epidermal KCcontaining ROIs (Figure 2d), we reasoned that K1/M2 might reveal the markers of MIS, which grows within the epidermis. We performed linear regression comparing the seven malignant ROIs in K1 and M2 with all the other ROIs in K1 and M2. The most highly enriched gene in this analysis that was also not detected in CN was the S100 calcium-binding protein family member and DAMP S100A8 (Figure 4a).
We determined by IHC that S100A8 is expressed by KCs rather than by melanocytes ( Figure 4b); yet, S100A8 was most strongly expressed in ROIs containing >50% melanocytes rather than in those containing >50% KCs (Figure 4c), suggesting that KCs may express S100A8 in response to consumption of the epidermis by malignant melanocytes. Indeed, we observed that the ROIs with the highest S100A8 expression had melanocytes scattered throughout the epidermis, a histopathologic feature of MIS ( Figure 4d). Expression of S100A8 was not strongly correlated with KC content (Figure 4d), further indicating that high S100A8 expression is not merely a function of KC abundance. KRT17 and KRT6 were the closest points to S100A8 in UMAP space, and indeed, these genes were enriched in the same ROIs ( Figure 4f). This profile was not seen for other keratins genes such as KRT14 (Figure 4f). KRT17 and KRT6 are known to be upregulated in wounded skin (Zhang et al., 2019); thus, enrichment of these genes in the same ROIs as S100A8 supports the notion that the growth of MIS within the epidermis elicits an injury response within the KCs. S100A8/A9 are detected in the KC microenvironment of melanoma To further validate the data from spatial transcript profiling, we performed S100A8 IHC on an independent cohort of 252 melanocytic tumors (68 CN, 66 dysplastic nevi, 69 MIS, 49 MM) (Table 1). S100A8 expression was scored on the basis of the percentage of epidermis expressing S100A8 that was directly associated with the tumor (epidermis containing tumor and/or epidermis overlying intradermal tumor; a schematic of the scoring is shown in Supplementary Figure S4). Notably, the KC microenvironment of melanoma (Figure 5a and b and Supplementary Figure S5) and many cases of MIS (Figure 5c and d) showed prominent expression of S100A8, whereas the KC microenvironment of most dysplastic nevi (Figure 5e and f) and CN (Figure 5g and h and Supplementary Figure S5) lacked or had only limited staining (Table 2) (P < 0.001; Figure 5i; area under the curve ¼ 0.83) (the representative images of scores 1-5 are shown in Supplementary Figure S6). A binary logistic regression model showed that increased S100A8 IHC score was significantly associated with MM tumor type (OR ¼ 2.49, 95% confidence interval ¼ 1.93-3.21) and remained significant after adjusting for sex, anatomic site, and age (OR ¼ 2.54, 95% confidence interval ¼ 1.92-3.36). S100A8 expression was also apparent in scattered immune cells within the dermis of both nevi and melanoma as well as within the normal hair follicle epithelium ( Figure 5 and Supplementary Figure S5), but this expression was not taken into account when scoring S100A8 staining. Remarkably, if the tumor showed skip areas, that is, areas of uninvolved skin, the epidermis lacked S100A8 expression in these foci ( Figure 5j). S100A8 is generally coexpressed with S100A9, forming a complex known as calprotectin (Gebhardt et al., 2006). Because the DSP probe pool used in our experiments did not target S100A9 and our S100A8 antibody (CF-145) is not cross-reactive, we analyzed the expression of S100A9 in a subset of tumors by IHC. Similar to S100A8, S100A9 was expressed by the KCs associated with melanoma but not with nevi (Supplementary Figure S7 and Supplementary Table S2). Together, these data show that epidermal KCs express S100A8/A9 in response to melanoma growth, revealing the cellular origin of a melanoma biomarker within the tumor microenvironment and emphasizing the interaction between melanoma and the KC microenvironment.

DISCUSSION
Melanoma is the fifth most common cancer type in the United States and causes the vast majority of skin cancer deaths (Surveillance, Epidemiology, and End Results Program, 2019). Despite advances in immunotherapy, markers predicting sustained treatment response are inadequate. Understanding the interplay between melanocytes and neighboring KCs and immune cells will be crucial to the development of improved diagnostic and prognostic tools and therapeutic targets. We aimed to study tumor microenvironment interactions in melanoma within their native morphological context using GeoMx DSP.
After validating the performance and reproducibility of the DSP assay on our FFPE tissue sections, we profiled 134 ROIs in a cohort of 12 common benign and malignant melanocytic tumors with a panel that targets 1,400þ immune-oncologyrelated genes. We detected over 900 targets with high confidence in a single experiment, which would traditionally = normalized counts for the 923 genes detected in the experiment (489 genes below the detection threshold in all the 134 ROIs were removed before clustering). Tumor type and ROI type are indicated by colored bars. The five largest clusters (hclust method) are boxed and named according to their predominant cell type (M1-2 for melanocyte rich, I1-2 for immune rich, K1 for keratinocyte rich). DSP, Digital Spatial Profiler; ROI, region of interest. ROIs classified as keratinocyte rich, melanocyte rich, keratinocyte/melanocyte, or mixed. Significance (-log 10 of P-value) was determined by linear regression with a term for random effects from intertissue variation. (b) Volcano plot comparing gene expression in I1 ROIs with that in all the other ROIs classified as immune or immune/melanocyte. Significance (-log 10 of P-value) was determined by linear regression with a term for random effects from intertissue variation. In a and b, genes were only considered if they were above the detection threshold in at least three ROIs, and gene names are only shown if the gene was below the detection threshold in all common nevus ROIs. (c, d) UMAP analysis comparing the spatial expression profiles of all 923 genes detected in at least one ROI. require tedious microdissection or be limited to a few genes (Wang et al., 2012). Other technologies for transcriptomescale spatial profiling have recently been applied to melanoma, but our experiment detected approximately 10 times more transcripts per square micron than a recent study (Thrane et al., 2018). Importantly, the compatibility of GeoMx DSP with FFPE enables the profiling of archival tissues (all specimens in our study were at least aged 2 years), which has been especially difficult for skin samples (Kwong et al., 2018). This is particularly important for the study of melanoma evolution in patient-derived benign and malignant primary tumors and the KC microenvironment, which have been mainly overlooked in single-cell RNA-seq studies of melanoma. (Jerby-Arnon et al., 2018;Tirosh et al., 2016) A cell type-enriched ROI selection strategy enabled us to directly compare similar cell populations across tumors when looking for gene enrichment in melanomas (Figures 2 and 3). This approach, coupled with a dimensionality reductionbased approach to identify coexpressed genes, successfully identified known melanoma-associated markers and showed their specificity to melanocytes, immune infiltrates, or the epidermal (KC) microenvironment (Figures 3 and 4). Interestingly, genes enriched in melanocytes (PRAME), KCs (S100A8), or immune cells (PTPRC, CXCL9) in our study correspond to the three components of a commercially available diagnostic test; our results may explain why the three components do not correlate well with each other and are all required to generate the most predictive score (Clarke et al., 2015).
Importantly, a published bulk RNA-seq dataset comparing 57 melanomas and 23 nevi (Kunz et al., 2018) validates our data and highlights the advantages of a spatially resolved, cell type-specific analysis, especially regarding the tumor microenvironment. Although the melanoma markers specific to melanocyte ROIs were also melanoma enriched in the bulk RNA-seq data, markers of the immune and KC microenvironment were not as strongly melanoma associated (Supplementary Figure S8), likely reflecting the masking of less abundant cell types in bulk measurements. By contrast, enrichment of these gene sets in melanoma was of similar magnitude for melanocyte-associated versus for immuneassociated genes in the DSP data ( Figure 3a and b), suggesting that this technology may be more sensitive for identifying gene products that originate from nontumor cells in the tumor microenvironment. In fact, a recent study using GeoMx DSP to study protein expression in melanoma found that expression of PD-L1 in macrophages rather than in tumor cells was predictive of response to immunotherapy (Toki et al., 2019). Because the understanding of tumor microenvironment interactions in melanoma is much needed in the era of melanoma immunotherapy, additional studies are warranted to validate the genes enriched in melanoma immune infiltrates in our DSP cohort. S100A8 is a DAMP that has multiple roles in promoting immune responses and inflammation (Nukui et al., 2008;Srikrishna, 2012). It is most well-known as part of the S100A8/A9 complex (calprotectin), which is canonically expressed and secreted by neutrophils, monocytes, and macrophages (Gebhardt et al., 2006). S100A8/A9 is also upregulated in a number of inflammatory disorders such as psoriasis and cystic fibrosis (Gebhardt et al., 2006;Nukui et al., 2008) and is expressed in a variety of epithelial tumor types (Gebhardt et al., 2006), where it is associated with invasiveness and poor prognosis (Arai et al., 2008). Although the upregulation of S100A8/A9 in melanoma has been established using bulk methods (Clarke et al., 2015;Kunz et al., 2018;Wagner et al., 2019), previous studies have assumed that it is expressed by immune cells (Clarke et al., 2015;Wagner et al., 2019). Our analysis reveals the KCderived origin of S100A8/A9 during melanoma development (Figures 4a-c and 5 and Table 2 and Supplementary  Figures S5-7). This observation may explain why S100A8/ A9 is strongly detected in primary melanomas but not in metastases (Xiong et al., 2019). Furthermore, it emphasizes the importance of understanding a biomarker's cellular origin because noncutaneous metastasis and microdissected tumors may lack KCs and thus produce falsely low expression levels of S100A8/A9 in the current commercially available tests. In addition, expression of S100A8/A9 specifically at the skin surface in early melanoma could potentially be exploited to increase the sensitivity of an adhesive patch biopsy assay, similar to that already available for PRAME (Gerami et al., 2017). Finally, the results suggest a potential role for S100A8 IHC as an ancillary test for the diagnosis of melanoma, especially because IHC-based testing can be readily adopted for use in most pathology laboratories.
Owing to the correlation of S100A8 expression with melanocyte growth within the epidermis (Figure 4b-d) and expression of the wound-associated keratins 6 and 17 (Figure 4e), S100A8 expression in this context may be a response to inflammation or destruction of the epidermis by the melanocytes. This notion is supported by the previous literature on S100A8 being expressed in epithelial cells in response to stress and inflammation (Kerkhoff et al., 2012). Cytokines secreted by nearby tumor cells likely play a role as well because KCs overlying melanoma purely within the dermis also strongly expressed S100A8 ( Figure 5), and multiple cytokines are known to induce S100A8 in normal KCs (Nukui et al., 2008). Because S100A8/A9 is a chemoattractant for melanocytes that express certain cell adhesion molecules such as MCAM, ALCAM, and RAGE (Ruma et al., 2016), induction of S100A8 in the epidermis may stimulate melanocyte migration; indeed, S100A8/A9 has been implicated in metastasis of melanoma and other tumor types (Arai et al., 2008;Ruma et al., 2016;Saha et al., 2010). Additional  . Analysis of K1/M2 clusters reveals enrichment of S100A8 in keratinocytes during melanomagenesis. (a) Volcano plot comparing gene expression in the subset of ROIs classified as melanoma in situ by a pathologist with that in all other ROIs in K1/M2. Significance (-log 10 of P-value) was determined by linear regression with a term for random effects from intertissue variation. Genes were only considered if they were above the detection threshold in at least three ROIs, and gene names are only shown if the gene was not detected in any CN ROIs. (b) Representative IHC image (left panel) and corresponding H&E image (right panel) showing that S100A8 is expressed by keratinocytes (arrowhead) rather than by melanocytes (arrow). (c) Ternary plot of S100A8 expression in all the ROIs, with a zoomed-in view of keratinocyte/melanocyte ROIs at left. (d) H&E images of selected ROIs containing a mixture of keratinocytes and melanocytes. Norm S100A8 counts are plotted against keratinocyte score and colored by tumor type. (e) Ridgeline plot of Norm counts for selected genes in all the ROIs, organized by tumor type. The blue line indicates LOQ in each ROI. Bar ¼ 0.05 mm. CN, common nevus; DN, dysplastic nevus; FC, fold change; IHC, immunohistochemistry; LOQ, limit of quantitation; MIS, melanoma in situ; MM, invasive melanoma; Norm, normalized; Q, quartile; ROI, region of interest.
DSP studies profiling a larger number of patients and ROIs are warranted to further resolve the interplay between KCs and melanocytes during melanomagenesis.

Cases
This study was approved, and a waiver of the informed consent requirement was granted by the institutional review board of the University of California Davis (Sacramento, CA) because specimens were deidentified. Pathology archives were searched for CN, dysplastic nevi, MIS, and MM, and archived H&E-stained sections were reviewed by a dermatopathologist (MK) to confirm diagnoses.
Serial sections of 5-mm thickness derived from FFPE tissue blocks were cut. One section was stained with H&E, and two unstained sections were mounted on positively charged histology slides for in situ hybridization and DSP. Whole-slide imaging of the H&E sections was performed on a Nikon TE2000-E microscope (Nikon, Tokyo, Japan), and ROIs representative of the tumor and its microenvironment were selected by the pathologist.

In situ hybridization
See Supplementary Materials and Methods for details.

Digital spatial profiling
Slides were blocked for 30 minutes as described (Merritt et al., 2020) and then incubated with 300 nM SYTO 13 and fluorescently conjugated antibodies to CD3, S100B, PMEL, and CD45 for 1 hour and washed in Â2 saline-sodium citrate. Slides were loaded onto a DSP instrument and submerged/washed in PBS with 0.1% Tween 20 as described (Merritt et al., 2020). After Â20 fluorescence scanning to obtain a high-resolution image of the tissue, ROIs were selected by matching to the pathologist-selected regions on the H&E-stained serial section. Between 6 and 16 ROIs were chosen per tissue (all were 200 mm diameter circles). Indexing oligos were released from each ROI by exposure to UV light as described (Merritt et al., 2020), and 10 ml of liquid from above the ROI was collected by a microcapillary tip and deposited in a 96-well plate.

Library preparation and sequencing
Indexing oligos from each ROI were PCR amplified using primers that (i) hybridize to constant regions and (ii) contain unique dualindexed barcoding sequences to preserve ROI identity. PCR products were pooled and purified twice with Ampure XP beads (Beckman Coulter, Brea, CA). Library concentration and purity were measured using a high-sensitivity DNA Bioanalyzer chip (Agilent Technologies, Santa Clara, CA). Paired-end (2 Â 75 bp reads) sequencing was performed on an Illumina MiSeq instrument (pilot

Data analysis
A high-confidence detection threshold was set at the geometric mean plus 2.5 SDs of the negative probes. A total of 923 of 1,412 genes (65%) in the full panel were above the detection threshold in at least one ROI. The 489 genes below the Figure 5. S100A8 is detected in the keratinocyte microenvironment of melanoma. (a-h) Representative images of S100A8 IHC (upper panel) and corresponding H&E staining (lower panel) of (a, b) invasive melanoma, (c, d) melanoma in situ, (e, f) dysplastic nevus, and (g, h) common nevus. S100A8 (brown) is expressed by keratinocytes in (a) invasive melanoma and (c) melanoma in situ but not in (e) dysplastic nevus or (g) common nevus. (i) S100A8 ROC curve for in situ or invasive melanoma. (j) Low power image of S100A8 IHC of melanoma with in situ and invasive components and an area of uninvolved skin (upper panel). S100A8 (brown) is expressed by keratinocytes in the epidermis directly associated with melanoma (insets 1 and 3) but not in the uninvolved epidermis (inset 2). Letter e indicates epidermal layer containing keratinocytes, and letter t indicates tumor. Melanin pigment (brown) is present in (e) dysplastic nevus and in association with invasive melanoma (*) in j. Bar ¼ 0.05 mm for a--h. Bar ¼ 0.5 mm for j (upper panel). Bar ¼ 0.1 mm for j (insets). AUC, area under the curve; IHC, immunohistochemistry; ROC, receiver operating characteristic.
detection threshold in all ROIs were excluded from further analysis.
IHC S100A8 mouse mAb (CF-145) (catalog number 14-9745-82, eBioscience, San Diego, CA) and S100A9 mouse mAb (Clone #747315) (catalog number MAB5578-SP, R&D Systems, Minneapolis, MN) were used. Staining was scored with 100% consensus agreement by two observers (MK and SW) on the basis of the percentage of epidermis stained that is directly associated with the tumor, using a previously described scoring system with appropriate modifications (Lezcano et al. Figure S3).

Data availability statement
Datasets related to this article can be found at https://www.ncbi.nlm. nih.gov/geo/ using the accession number GSE168013.  Percentage of epidermis associated with tumor stained: M Kiuru et al. Pathology archives were searched for common nevi, dysplastic nevi, melanomas in situ, and invasive melanomas, and archived H&E-stained sections were reviewed by a dermatopathologist (MK) to confirm diagnoses. One common nevus, one dysplastic nevus, one melanoma in situ, and one invasive melanoma were included in the first phase of digital spatial profiling. Three common nevi, three dysplastic nevi, three melanomas in situ, and three invasive melanomas from trunk or extremity of male and female patients aged 49-69 years were included in the second phase of digital spatial profiling. A total of 12 common nevi, 12 dysplastic nevi, 11 melanomas in situ, and 12 invasive melanomas were included in the immunohistochemistry (IHC) study (as shown later).
Serial sections of 5-mm thickness derived from formalinfixed, paraffin-embedded tissue blocks were cut. One section was stained with H&E, and two unstained sections were mounted on positively charged histology slides for in situ hybridization and digital spatial profiling. Whole-slide imaging of the H&E sections was performed on a Nikon TE2000-E microscope (Nikon, Tokyo, Japan), and regions of interest (ROIs) representative of the tumor and its microenvironment were selected by the pathologist.

In situ hybridization
The pilot experiment probe pool contained probes that target 108 mRNA transcripts and 80 negative probes (1,048 probes total, a median of 10 probes per target). The full panel included probes that target 1,412 mRNA transcripts and 130 negative probes (4,998 probes total, median of 3 probes per target).
To expose epitopes, unstained formalin-fixed, paraffinembedded sections were deparaffinized, heated in ER2 solution (Leica, Wetzlar, Germany) at 100 C for 20 minutes, and treated with 1 mg/ml proteinase K (Ambion, Austin, TX) at 37 C for 15 minutes on a BOND RXm autostainer (Leica). An overnight in situ hybridization was performed as described (Merritt et al., 2020) with a final probe concentration of 4 nM per probe. The pilot experiment probe pool contained probes that target 108 mRNA transcripts as well as 80 negative probes (1,048 probes total, a median of 10 probes per target). The full panel included probes that target 1,412 mRNA transcripts as well as 130 negative probes (4,998 probes total, a median of 3 probes per target). Slides were washed twice at 37 C for 25 minutes with 50% formamide/Â2 saline-sodium citrate buffer to remove unbound probes.

Digital spatial profiling
Slides were blocked for 30 minutes as described (Merritt et al., 2020) and then were incubated with 300 nM SYTO 13, and fluorescently conjugated antibodies to CD3, S100B, PMEL, and CD45 were added for 1 hour and washed in Â2 saline-sodium citrate. Slides were loaded onto a digital spatial profiler instrument and submerged/washed in PBS with 0.1% Tween 20 as described (Merritt et al., 2020). After Â20 fluorescence scanning to obtain a high-resolution image of the tissue, ROIs were selected by matching to the pathologist-selected regions on the H&E-stained serial section (as discussed earlier). Between 6 and 16 ROIs were chosen per tissue (all were 200 mm diameter circles). Indexing oligos were released from each ROI by exposure to UV light as described (Merritt et al., 2020), and 10 ml of liquid from above the ROI was collected by a microcapillary tip and deposited in a 96-well plate.

Library preparation and sequencing
Indexing oligos from each ROI were PCR amplified using primers that (i) hybridize to constant regions and (ii) contain unique dual-indexed barcoding sequences to preserve ROI identity. PCR products were pooled and purified twice with AMPure XP beads (Beckman Coulter, Brea, CA). Library concentration and purity were measured using a highsensitivity DNA Bioanalyzer chip (Agilent Technologies, Santa Clara, CA). Paired-end (2 Â 75 base pair reads) sequencing was performed on an Illumina MiSeq instrument (pilot experiment) or Illumina HiSeq 2500 instrument (full panel experiment) (Illumina, San Diego, CA).

Data processing
After sequencing, reads were trimmed, merged, and aligned to a list of indexing oligos to identify the source probe. The unique molecular identifier region of each read was used to remove PCR duplicates and duplicate reads, thus converting reads into digital counts. A total of 40% of reads were unique in the pilot experiment, and 2% of reads were unique in the full panel experiment, indicating that in both cases the sequencing depth was sufficient to sample the population. For each gene in each sample, the reported count value is the mean of the individual probe counts after removal of outlier probes.

Data analysis
A high-confidence detection threshold was set at the geometric mean plus 2.5 SDs of the negative probes. A total of 923 of 1,412 genes (65%) in the full panel were above the detection threshold in at least one ROI. The 489 genes below the detection threshold in all ROIs were excluded from further analysis.

Statistics
IHC data were summarized with descriptive statistics. A chisquare test was implemented to examine whether the frequency of S100A8 scores or S100A9 scores were significantly different between the tumor types. Chi-square tests and ANOVA tests were utilized to test for associations between S100A8 score and sex, age, location of tumor, nevus growth pattern, invasive melanoma tumor thickness, invasive melanoma ulceration, and primary tumor stage. To evaluate the association of S100A8 IHC score, sex, location, and age with the occurrence of invasive melanoma tumor, we used multivariable binary logistic regression to calculate adjusted ORs and 95% confidence interval. All analyses were conducted using IBM SPSS Statistics for Windows, version 26 software (IBM, Armonk, NY).

Study approval
This study was approved, and a waiver of the informed consent requirement was granted by the institutional review board of the University of California Davis (Sacramento).

Spatial RNA Profiling in Melanoma
Supplementary Figure S1. DSP yields highly reproducible gene counts that are consistent with sample biology. (a) Scatterplot comparing raw counts from replicate library prep PCRs performed from the same DSP aspirate. Each dot represents a gene count from two measurements of the same ROI. (b) Scatterplot comparing raw counts from replicate DSP experiments performed on serial sections from the same tissue. Each dot represents a gene count from replicate ROIs on serial sections. (c) Raw probe counts (open circles) and gene counts (filled circles) for representative melanocyte-rich and immune-rich ROIs (ROIs 13 and 14 from pilot melanoma case). Counts for probes targeting the same gene tracked together. The geometric mean of the negative probes for each ROI is indicated by a red dot, and the detection threshold (geometric mean plus 2.5 SDs of the negative probes) is indicated by a red line. Examples of genes expected to be enriched in melanocytes or immune cells are indicated by green or pink vertical bars, respectively. (d) Counts for the genes used as visualization markers correlate qualitatively with the corresponding fluorescent signals during DSP. Left: H&E image of ROIs 12-17 from the pilot melanoma case. Center: DSP IHC image of ROIs 12-17 from the pilot melanoma case. Right: Bar plot of counts from ROIs shown for genes corresponding to the DSP visualization markers. The filled color represents the IHC channel; S100B and PMEL antibodies both fluoresced in the green channel. Counts were scaled by gene to allow for the plotting of different genes on the same axis. Black lines represent two (*) detection thresholds because two genes are plotted per bar. DSP, digital spatial profiling; IHC, immunohistochemistry; ROI, region of interest.

Spatial RNA Profiling in Melanoma
Supplementary Figure S6. Representative images of S100A8 staining scores 1--5. (a-d) Score 1: Representative images of (c, d) S100A8 IHC staining score 1 and (a, b) corresponding H&E staining in a common nevus. (e--h) Score 2: Representative images of (g, h) S100A8 IHC staining score 2 and (e, f) corresponding H&E staining in a dysplastic nevus. (i-k) Score 3: Representative images of (k, l) S100A8 IHC staining score 3 and (i, j) corresponding H&E staining in melanoma in situ; an associated melanocytic nevus is present. (m--p) Score 4: Representative images of (o, p) S100A8 IHC staining score 4 and (m, n) corresponding H&E Supplementary Figure S7. S100A9 is expressed in the epidermis associated with melanoma. (a-d) Representative images of S100A9 IHC (right panel) and corresponding H&E staining (left panel) of (a, b) invasive melanoma and (c, d) common nevus. S100A9 (brown) is expressed by keratinocytes in (b) invasive melanoma but not in (d) common nevus. Letter e indicates epidermis and letter t indicates tumor. IHC, immunohistochemistry.

Supplementary Table 1. Continued
Gene Description