A digital mRNA expression signature to classify challenging Spitzoid melanocytic neoplasms

Spitzoid neoplasms are a challenging group of cutaneous melanocytic proliferations. They are characterized by epithelioid and/or spindle‐shaped melanocytes and classified as benign Spitz nevi (SN), atypical Spitz tumors (AST), or malignant Spitz tumors (MST). The intermediate AST category represents a diagnostically challenging group since on purely histopathological grounds, their benign or malignant character remains unpredictable. This results in uncertainties in patient treatment and prognosis. The molecular properties of Spitzoid lesions, especially their transcriptomic landscape, remain poorly understood, and genomic alterations in melanoma‐associated oncogenes are typically absent. The aim of this study was to characterize their transcriptome with digital mRNA expression profiling. Formalin‐fixed paraffin‐embedded samples (including 27 SN, 10 AST, and 14 MST) were analyzed using the NanoString nCounter PanCancer Pathways Panel. The number of significantly differentially expressed genes in SN vs. MST, SN vs. AST, and AST vs. MST was 68, 167, and 18, respectively. Gene set enrichment analysis revealed upregulation of pathways related to epithelial–mesenchymal transition and immunomodulatory‐, angiogenesis‐, hormonal‐, and myogenesis‐associated processes in AST and MST. A molecular signature of SN vs. MST was discovered based on the top‐ranked most informative genes: NRAS, NF1, BMP2, EIF2B4, IFNA17, and FZD9. The AST samples showed intermediate levels of the identified signature. This implies that the gene signature can potentially be used to distinguish high‐grade from low‐grade AST with a larger study cohort in the future. This combined histopathological and transcriptomic methodology is promising for prospective diagnostics of Spitzoid neoplasms and patient management in dermatological oncology.

Spitzoid neoplasms are a challenging group of cutaneous melanocytic proliferations. They are characterized by epithelioid and/or spindle-shaped melanocytes and classified as benign Spitz nevi (SN), atypical Spitz tumors (AST), or malignant Spitz tumors (MST). The intermediate AST category represents a diagnostically challenging group since on purely histopathological grounds, their benign or malignant character remains unpredictable. This results in uncertainties in patient treatment and prognosis. The molecular properties of Spitzoid lesions, especially their transcriptomic landscape, remain poorly understood, and genomic alterations in melanoma-associated oncogenes are typically absent. The aim of this study was to characterize their transcriptome with digital mRNA expression profiling. Formalin-fixed paraffin-embedded samples (including 27 SN, 10 AST, and 14 MST) were analyzed using the NanoString nCounter PanCancer Pathways Panel. The number of significantly differentially expressed genes in SN vs. MST, SN vs. AST, and AST vs. MST was 68, 167, and 18, respectively. Gene set enrichment analysis revealed upregulation of pathways related to epithelial-mesenchymal transition and immunomodulatory-, angiogenesis-, hormonal-, and myogenesis-associated processes in AST and MST. A molecular signature of SN vs. MST was discovered based on the top-ranked most informative genes: NRAS, NF1, BMP2, EIF2B4, IFNA17, and FZD9. The AST samples showed intermediate levels of the identified signature. This implies that the gene signature can potentially be used to distinguish high-grade from low-grade AST with a larger study cohort in the future. This combined histopathological and transcriptomic methodology is promising for prospective diagnostics of Spitzoid neoplasms and patient management in dermatological oncology. Spitzoid melanocytic neoplasms comprise uncommon melanocytic skin lesions usually arising in children and adolescents, but may also occur in older individuals [1]. Key challenges at the time of primary histopathological evaluation of Spitzoid melanocytic neoplasms are the distinction between benign or malignant processes and the choice of optimal therapy in case of melanocytic atypia or malignancy. Spitzoid neoplasms were first described by the American pathologist Sophie Spitz in 1948 as 'juvenile melanomas' or as 'melanomas of childhood' [2] and are composed of large epithelioid and/ or spindle-shaped melanocytes with large nuclei that contain vesicular chromatin and often prominent nucleoli [3]. They may be purely junctional, compound, or dermal, and have a background with variable amounts of lymphocytes, blood vessels, and sclerosis.
According to the recently launched fourth edition of the World Health Organization classification of skin tumors, these lesions range from benign Spitz nevi (SN) to highly proliferative and pleomorphic malignant Spitz tumor (MST) or Spitzoid melanoma [4]. Complicating the diagnostic issue, the major challenge lies in a subset of lesions that fall in between the two extreme types, that is, the atypical Spitz tumor (AST) or Spitzoid tumor of uncertain malignant potential (STUMP, [4,5]). The biological potential of AST/ STUMP is currently impossible to predict on histopathological grounds alone [6][7][8].
The to-date identified molecular alterations in these lesions have not yet resulted in one stand-alone or a clearly defined distinct set of ancillary biomarkers that can reliably predict whether an ambiguous Spitzoid neoplasm has benign or malignant potential (for a review, see Ref. [9][10][11]). Consequently, potentially benign AST may be misdiagnosed as MST, thereby exposing the patient to aggressive surgical and oncological treatment. Until now, only few transcriptomic studies on Spitzoid melanocytic lesions have been performed [12][13][14]. Advancing in that direction, we recently identified a distinct molecular expression profile differentiating SN from common nevocellular nevi (NCN), thereby suggesting a differential pathogenesis of the two groups of melanocytic neoplasms [15].
Based on these findings, we wanted to know whether such a differential gene expression profile can be found along the whole spectrum of Spitzoid melanocytic neoplasms (SN vs. AST vs. MST), and whether such a profile might be of help in the differential diagnosis in daily histopathology practice, clinical treatment decisions, and patient management. The aim of this study was thus to characterize the transcriptomic landscape of SN, AST, and MST with digital mRNA expression profiling. Our research question was threefold: (a) Are there differentially expressed genes in the three categories of Spitzoid melanocytic neoplasms (AST vs. MST, SN vs. AST, and SN vs. MST)? (b) Do possibly identified mRNA transcripts belong to distinct gene pathways?; And (c) can a distinct gene expression signature of top-ranked genes be identified to discriminate benign from malignant lesions?

Patients and tissues
Formalin-fixed paraffin-embedded (FFPE) excision specimens of patients harboring a SN, an AST, or a MST were retrieved from the archives of the Department of Pathology at the University of Leuven, KUL Belgium, the Maastricht Pathology Tissue Collection from the Department of Pathology, Maastricht University Medical Center, (MUMC+), the Netherlands, and the Department of Dermatology, Ruhr-University Bochum, Bochum, Germany. The patients kindly provided their written informed consent at hospital admission to the processing of tissue and personal data. All samples had been excised for diagnostic and therapeutic reasons. All use of tissue and patient data was in agreement with the Dutch Code of Conduct for Observational Research with Personal Data (2004) and Tissue (2001) and in accordance with Ethical Principles for Medical Research Involving Human Subjects (World Medical Association Declaration of Helsinki). Histopathological diagnoses were previously defined in routine diagnostics and were confirmed by three experienced dermatopathologists (VW, JVO, and LMH).
A total of 51 cases were included in the analysis. Of these, 27 were histologically diagnosed as banal SN, 10 as AST, and 14 as MST. SN were defined as benign melanocytic nevi composed of large epithelioid, oval, or spindled melanocytes arranged in nests and/ or fascicles without significant cytonuclear atypia. The diagnosis of an AST was made if there was an increase in at least one worrisome histological feature, that is, ulceration, size > 5 mm, increase in cell density with confluent growth, infiltrative growth into subcutaneous tissue with 'pushing' margins, increased cytonuclear atypia, > 2 dermal mitoses, absence of junctional clefts, few or no Kamino bodies, and more extensive pagetoid extension. The diagnosis of MST was based on the following histopathological criteria: ulceration, asymmetrical architecture, infiltrative growth, severe and/or confluent cytonuclear atypia, frequent dermal mitoses, pushing borders, epidermal effacement, and pagetoid extension [4]. Figure 1 illustrates the morphology of a prototypic SN (Fig. 1A,D), AST (Fig. 1B,E), and MST (Fig. 1C,F). The study was approved by the Maastricht Ethics Committee of the University of Maastricht, the Netherlands, and by the Institutional Review Board of the University Hospitals of Leuven, Belgium (project number S 59659).

Patient cohort
The SN group consisted of 27 patients with 15 females and 11 males ranging in age between 2 and 62 years (mean 31.3 years, median 31.5 years; Table S1); from one case (SN23), clinical data were not available. Histologically, nine SN were diagnosed as dermal, five SN as junctional, and 13 as compound lesions. SN4 showed dysplastic features with mild cytonuclear atypia and variation in cellular morphology but did not fulfill criteria for the diagnosis of an AST. The AST group consisted of 10 patients, with seven female and three male patients showing a similar age distribution ranging from 9 to 55 years (mean 31.3 years, median 29 years; Table S1). The AST group consisted of two dermal and eight compound lesions. Like the SN group, the AST group showed localization on the extremities, trunk, and head and neck area. A sentinel node procedure was not performed in the AST cohort in accordance with the wishes of the patients.
The MST group consisted of 14 cases with nine females and five males ranging in age from 20 to 66 years (mean 53.1 years, median 26.5 years; Table S1). The MST group showed a mean Breslow thickness of 1.39 mm (range from 0.3 to 2.98 mm, median 1.27 mm). Three patients showed regional lymph node metastasis (MST2, MST4, and MST9) with in two of them development of distant metastasis during follow-up (MST2 and MST4).
Follow-up was available for all patients except for cases SN23, AST1, AST2, and AST10 (Table S1).

Tissue preparation and mRNA expression analysis
Serial sections (5 µm) were cut from FFPE tissue blocks. One section per sample was mounted onto a glass slide and

Data normalization
Gene expression counts were normalized based on total library size (total number of counts per sample; Fig. S1 and Table S3). As such, both the housekeeping and endogenous genes were used for normalization. This method is commonly used when many genes are analyzed simultaneously (NanoString recommends > 300). The studied samples showed a mean library size of 307 914 counts [median 292 412; range 789 054 (sample AST4) to 19 042 (sample MST4)]. All analyses were performed using R and Bioconductor packages.

Identification of differentially expressed genes
Marker discovery and validation of differentially expressed genes from the nCounter data was done using the limma and voom method in R [16]. Three conditions were studied: (a) SN vs. MST, (bb) SN vs. AST, and (c) AST vs. MST. The limma procedure borrows information across all genes, which makes the analysis robust even when the sample size is small [17]. The voom procedure was applied to nonparametrically estimate the mean-variance relationship of the log-transformed gene counts, which is important when dealing with expression count data [18]. This procedure takes as input the raw counts, and total library size (total number of counts) per sample is used as a scaling factor to normalize the counts (Table S3). The voom procedure generates precision weights, which are incorporated into the limma linear modeling procedure, and this removes the mean-variance relationship in the log-counts. The normalized log-counts were used in subsequent analyses (Fig. S2).
To control for multiple testing, an adjusted P-value [i.e., false discovery rate (FDR) q-values] threshold of 25% was used for statistical significance. A heatmap of the expression data was generated using pheatmap in R. The data were scaled and centered before plotting using scale in R.

Identification of a gene signature
The Least Absolute Shrinkage and Selection Operator (LASSO) procedure was applied in R to create a gene expression signature for classifying SN in contrast to MST that was based on a small subset of the most informative gene expression markers. To prevent overfitting, all the genes were used as input for the analyses. A 10-fold cross-validation (CV) and the binomial deviance criterion were used to identify the optimal value for the tuning parameter that resulted in the smallest binomial deviance for classifying. In order to generate a robust gene expression signature, the LASSO procedure was run 500 times with each time using a randomly selected CV split of the data. Genes that were selected in at least 250 of the 500 LASSO individual models were included in the final signature. Each gene's average model coefficient (or relative weight in the signature) across all repetitions was calculated. Using this method, a model of six top-ranked gene expression markers was identified where every marker had a corresponding LASSO coefficient to distinguish SN from MST. A sum of the expression weights of the six top-ranked gene expression markers was calculated per sample and evaluated in relation to clinical and histological parameters. A genetic dataset for skin cutaneous melanoma (SKCM) provided by The Cancer Genome Atlas (TCGA, [19]) was downloaded from University of Santa Cruz (UCSC) Xena genome browser (https://xenabrowser.net/datapages/?co hort=TCGA%20Melanoma%20(SKCM)&removeHub= https%3A%2F%2Fxena.tre), providing mutational data for 473 melanoma samples. Mean expression values and Pvalues (> 0.05) were calculated using a t-test for each of the top-ranked signature genes in primary SKCM cases and metastatic samples.

Pathway analysis
A gene set enrichment analysis (GSEA) was done using the camera method (limma in R) to identify pathways that are significantly up-or downregulated in the three studied contrasts [20]. A matrix of gene expression levels (No. genes = 770) in the different samples was used as input for the analysis. An advantage of this method is that it accounts for intergenomic correlation. The gene expression data were compared to the HALLMARK Molecular Signatures Database (MSigDB, Broad Institute; No. = 50). The gene set was downloaded from: http://bioinf.wehi.edu.au/ software/MSigDB/.  Table S2 for the three studied conditions. For each identified gene, the log FC value, the average expression of the normalized log 2 intensity value from limma, the P-value, the t-value, and the adjusted P-value (or FDR q-value) used to correct for multiple testing are given. Volcano plots of the results from the three studied contrasts are shown in Fig. 2A-C. The blue dots indicate the significantly differentially expressed gene transcripts with the top five genes labeled with their gene names.

Differential gene expression in SN vs. MST
Sixty-eight out of the 770 genes showed significantly differential expression in the SN vs. the MST group (P-value < 0.05; Table S2.1). In the SN group, 27 transcripts had a higher mean expression level with a positive log FC and 41 transcripts showed a lower mean expression level with a negative log FC when compared to the MST group ( Fig. 2A). The gene transcript with the lowest log FC (log FC = À1.17) was WIF1. The transcript SPP1 showed the highest log FC (log FC = 1.59). An absolute log FC of at least 1 was identified in the four gene transcripts LAMB4, SPP1, PLA1A, and WIF1 (Table S2.1).

Differential gene expression in SN vs. AST
In the group of SN vs. AST, 167 genes out of the 770 genes showed differential expression. From these, 36 genes showed a positive log FC and 131 had a negative log FC. Figure 2B shows the volcano plots with the identified differentially expressed genes in SN vs. AST. The lowest log FC (log FC = À1.57) was seen for the transcript ZBTB16, whereas the highest log FC (log FC = 1.36) was identified for RXRG (Table S2.2). Out of the 167 genes, 15 gene transcripts showed an absolute log FC of 1. The gene transcripts ZBTB16, FGF14, SGK2, WT1, and CNTFR had a log FC < À1, and the gene transcripts ITGB3, PLA2G2A, SHC4, FOSL1, NGFR, SOCS3, IL13RA2, COMP, NTRK1, and RXRG showed a log FC > 1 (Table S2.2).

Differential gene expression in AST vs. MST
In the AST vs. MST cohort, there were 18 differentially expressed genes with 10 transcripts showing a positive log FC and 8 transcripts with a negative log FC (range À1.72 to 0.96; Fig. 2C, Table S2.3). NGFR, PLA2G2A, COL11A1, ITGA3, ITGB8, and COL6A6 were downregulated with a log FC < À1. The FGF19 transcript had the highest log FC (log FC = 0.96), and NGFR showed the lowest log FC (log FC = À1.72). None of the 18 identified transcripts showed a log FC> 1.

Pathway analysis
Gene set enrichment analysis was done to identify upand downregulated pathways in the three studied groups of Spitzoid melanocytic lesions using the camera method. Out of the 50 investigated MSigDB gene sets from Hallmark, a top of ten gene pathways was identified. The MSigDB gene sets are included in the following broad categories: (a) cell cycle/ proliferation; (b) cytokine/ immune/ inflammatory; (c) matrix/ adhesion; (d) hormone/ receptor/ signal transduction; and (e) transport and other. Detailed characteristics of the identified gene pathways including number of corresponding genes, direction of regulation, P-values, and FDR q-values are given in the supplemental data (Table S4) and Fig. 2D-F. A stringent FDR q-value threshold of 0.05 was used for statistical significance. As such, the number of significant pathways for SN vs. MST was two with a FDR q-value of 0.0071 for the Hallmark G2M checkpoint pathway and a FDR q-value of 0.0124 for the E2F targets pathway. There was downregulation of the Hallmark estrogen response early pathway and upregulation of the Hallmark spermatogenesis pathway in SN vs. MST (Fig. 2D). In the group of SN vs. AST, there was significant upregulation of the Hallmark epithelial-mesenchymal transition pathway (FDR q-value = 0.0001) and the Hallmark coagulation pathway (P = 0.0001; Fig. 2E). For the AST vs. MST condition, there were four significant pathways including the Hallmark epithelial-mesenchymal transition pathway (FDR q-value = 0.0001), the Hallmark coagulation pathway (FDR qvalue = 0.0209), the Hallmark G2M checkpoint pathway (FDR q-value = 0.0406), and the Hallmark E2F targets pathway (FDR q-value = 0.0468; Fig. 2F).
Interestingly, in SN vs. AST there was significant upregulation of the Hallmark epithelial-mesenchymal transition pathway and the Hallmark coagulation pathway, whereas in AST vs. MST, these pathways were significantly downregulated. The same trend was seen for the Hallmark myogenesis pathway, the Hallmark angiogenesis pathway, and the Hallmark complement pathway with upregulation in the group of SN  Table S4. vs. AST and downregulation in the group of AST vs. MST. However, the different regulation of these pathways failed to be statistically significant.

Identification of a molecular gene expression signature
The statistical learning method LASSO was used to create a molecular profile of the SN vs. MST group. From the analysis, a final gene expression signature was derived. A 10-fold CV identified the value for the tuning parameter that resulted in the lowest binomial deviance for classifying (Fig. 3A). The gene signature included six top-ranked gene transcripts: NRAS, EIF2B4, NF1, BMP2, FZD9, and IFNA17 (Fig. 3C). Detailed characteristics of the identified signature genes are summarized in Table 1 with chromosomal location, target sequence, subcellular localization, function of the target protein, and cellular pathway mapping. The signature levels of each sample are given as a multiplication from the gene expression level of the six genes with the associated LASSO coefficients ( Table 1, last column). The gene transcripts NRAS, IFNA17, and FZD9 were upregulated, and EIF2B4, BMP2, and NF1 were downregulated in the SN group compared to the MST group.
The identified gene transcripts showed genome-wide localization (including Chr. 1, 2, 7, 9, 17, and 20). The distribution of the differentially expressed gene transcripts across the genome is highlighted in a Manhattan plot (Fig. S3). Boxplot diagrams of the signature levels in the SN, AST, and MST are shown in Fig. 3B. Of note, the top-ranked signature gene IFNA17 identified with the LASSO procedure is not within the group of the 68 differentially expressed genes in SN vs. MST (Table S2.1).

Molecular signature in relation to histopathology and clinical features
In the following, the gene expression weight of the obtained molecular signature was calculated for all samples. The expression values of the defined signature genes NRAS, EIF2B, NF1, BMP2, FZD9, and IFNA17 were multiplied with the LASSO coefficients per sample. Transformed proportions of the signature expression weights are given in Table S1 (last column).
The SN samples had a 'low' expression value (< 0.4) with a mean value of 0.24 (range 0-0.54). Samples SN2, SN10, and SN17 formed the exception (> 0.4). These samples showed a compound histopathological architecture. SN2 had been excised from the arm of a 9-year-old boy, SN10 from the cervical area of a 39-year-old female patient, and SN17 from the shoulder area of a 32-year-old female patient.
Interestingly, the AST samples had intermediate levels of the signature in between the group of SN and MST (mean value 0.43, range 0.19-0.81). The samples AST1-3 and AST7 showed a 'high' gene expression weight (> 0.4) with a compound histopathological architecture and a trend toward a more central localization on the body (upper leg, gluteal area, back, and abdomen; Table S1). The samples AST4-6 and AST8-10 with a 'low' gene expression weight (< 0.4) showed localization on the extremities and central localization on the body. The AST samples with a 'low' gene expression weight had a dermal or compound histopathological architecture (Table S1). There was no significant difference in histopathological features such as mitotic index, or grade of nuclear atypia in between the 'low' and 'high' AST gene expression groups. All AST samples showed maturation within the dermal depth (data not shown).
The MST cases had a 'high' gene expression value (> 0.4) with a mean expression weight of 0.66 with exception of case MST4 (=0.38). There was a trend toward a higher Breslow thickness (mm) in association with the expression weight of the molecular signature; however, this failed to be statistically significant.

Molecular signature in relation to public datasets
In Table S5, mean expression values and P-values of the identified signature genes are shown for melanoma samples from the TCGA dataset [19]. The signature genes were significantly differentially expressed (P > 0.05) in primary vs. metastatic melanomas. The only exception was the gene transcript BMP2 with a P-value of 0.079. The IFNA17 transcript showed no expression in the primary melanoma group, and little expression in the metastatic samples when compared to absolute expression values from the other transcripts. The FZD9 transcript showed the highest range (1.10) in primary melanoma samples in relation to metastatic samples.

Discussion
In the present study, we identified genes that are differentially expressed in the three categories of Spitzoid melanocytic neoplasms (SN, AST, and MST). Although some molecular studies have been performed on Spitzoid tumors, the transcriptome of these distinct lesions remains largely unknown to date [12,13,15]. Therefore, the aim of this study was to identify a discriminating expression profile of SN, AST, and MST that would allow to discriminate these lesions from each other, and to define a list of genes that would add new insights into the underlying molecular basis of these lesions.
We performed the analysis using the PanCancer Pathways Panel from NanoString nCounter with a probe set of 770 genes related to cancer, growth, and invasion. NanoString nCounter gene expression profiling has been used to identify RNA expression patterns to accurately diagnose and classify different tumor types from various body sites [21][22][23]. This platform does not require reverse transcription or amplification, and the measurement of gene expression at the RNA level yields information on the actual functional state of a cell. Moreover, it allows robust gene expression profiling from FFPE tissue. However, despite the large number of more than 700 investigated genes the gene selection represents a systematic bias. The PanCancer Pathways Panel is focused on intrinsic tumoral characteristics. Thus, (micro)environmental tumoral features such as inflammation and immune response might be underexplored in our study. It might be worthwhile to  address Spitzoid melanocytic lesions with the PanCancer Immune Profiling Panel in future studies. The three studied groups showed differential gene expression with 68 significant genes in the group SN vs. MST, 167 significant genes when SN were compared with AST, and 18 significant genes in AST vs. MST. The reason why the number of differentially expressed genes differs so much in the three studied conditions might be due to the applied statistical methodology. The group AST vs. MST showed the lowest sample size (n = 10 for AST and n = 14 for MST), and thus, differentially expressed genes may have failed to overcome the applied level of statistical significance (FDR q-value with a threshold of 25%). However, the limma and voom methods were applied in this study to test for differential gene expression, which makes the analysis robust even when the sample size is small [17]. Another more likely explanation might be a higher genetic similarity (and thus a lower level of difference in the genetic make-up) between AST and MST than between AST and banal SN. This explanation also reflects the experience from histopathological observation where the morphological differentiation between high-grade AST and MST remains extremely challenging [5,7,8].
We identified a molecular gene expression signature including the six top-ranked gene transcripts NRAS, NF1, EIF2B4, IFNA17, BMP2, and FZD9 that could distinguish SN from MST. The gene transcripts NF1 and NRAS are associated with the Ras/ Raf/ MEK/ ERK (MAPK) pathway, which is well known to be deregulated in melanoma development and pathogenesis of Spitzoid neoplasms [19,[24][25][26]. Additionally, NRAS also plays an important role in the PI3/AKT (AKT) pathway [27,28]. We detected upregulation of NRAS transcripts, whereas NF1 transcripts showed to be downregulated in SN vs. MST. The observed reciprocal expression is reasonable from a functional point of view because NF1 encodes for neurofibromin, which stimulates the GTPase activity of NRAS, thereby inhibiting MAPK signaling [29][30][31]. The oncogenic activation of NRAS has been reported in 15-25% of classical cutaneous melanomas [19] and has also been found to be alternated in the subgroup of Spitzoid melanocytic neoplasms [10,32,33]. NF1 alterations are found in approximately 10 % of cutaneous melanomas [19] and have also been reported in subgroups of Spitzoid melanocytic neoplasms [25,34].
The identified signature gene EIF2B4 encodes for eukaryotic initiation factor 2 (EIF2), which is composed of five different subunits and catalyzes the exchange of EIF2-bound GDP for GTP [35]. The protein encoded by this gene is the fourth, or delta, subunit and has been associated with leukodystrophic neurodegenerative disease [36][37][38][39]. To our knowledge, there are no reports in the literature yet which have linked EIF2B4 to melanocyte pathology.
Although IFNA17 was not within the group of the 68 differentially expressed genes in SN vs. MST identified by the LASSO procedure, in combination with the other transcripts NRAS, NF1, BMP2, EIF2B4, and FZD9, it qualified as one of the top 6-ranked signature genes. IFNA17 encodes for the interferon (IFN)-alpha 17 molecule and plays a role in the JAK/ STAT/ phosphatidylinositol 3'-kinase (PI3K) pathway [40], which has also been reported to be deregulated in melanocytic neoplasms [41,42]. At the molecular level, it is well known that the IFNs have potent immune modulating and antitumor impacts on The BMP2 transcript encodes for bone morphogenetic proteins (BMPs), a group of signaling molecules that are part of TGF-b signaling [54], which has been reported to stimulate tyrosine gene expression and melanogenesis [55]. One may speculate that downregulation of BMP2 may relate to the low levels of pigmentation, often observed in SN. To our knowledge, there are no reports in the literature yet that have linked BMP2 to distinct categories of Spitzoid neoplasms yet.
Finally, the FZD9 gene also qualified as one of the top-ranked six transcripts for the gene signature and was upregulated in SN compared to MST. The transcript encodes for a frizzled class receptor (FZD) transmembrane receptor family member, which are part of the Wnt/ beta-catenin signaling pathway [56,57]. The Wnt/ beta-catenin pathway is well known to play a role in both development of melanocytes and their transformation to melanoma [58][59][60]. The major effect of Wnt ligand binding to its receptor is the stabilization of cytoplasmic beta-catenin through inhibition of the beta-catenin degradation complex. Betacatenin is then free to enter the nucleus and activate Wnt-regulated genes [61,62]. In line with FZD9 upregulation, the WIF1 transcript that encodes for Wnt inhibitory factor 1 (WIF1) was significantly downregulated with the lowest log FC of À1.17 from all 68 significantly differentially expressed genes in SN vs. MST. WIF1 interacts with the Wnt protein and therefore negatively regulates receptor binding to the FZD receptor family that initiate Wnt/ beta-catenin signaling [63]. When comparing our signature genes with the TCGA expression data, the FZD9 transcript showed a higher mean expression value in primary cutaneous melanoma samples in relation to metastatic samples. A possible explanation is that FZD9 may be a tumor suppressor in the presence of Wnt ligands in different kinds of cancer [64]. Our data on the FZD9 and WIF1 transcripts suggest an important role for the Wnt/ beta-catenin pathway in Spitzoid lesions and deserve further study in a larger cohort of lesions.
The AST samples showed levels of the identified signature intermediate between the group of SN and the group of MST (Fig. 3B). This implies that the six gene expression signature can potentially be used to segregate low-grade from high-grade AST (those that share features of SN vs. those that are more similar to MST). Of note, histopathological features that are commonly taken as indicators of melanocytic malignancy, such as nuclear atypia, absence of maturation in the dermis, deep extension, and deep dermal mitoses, most often fail to differentiate benign from malignant potential in Spitzoid lesions [14,65]. Our findings are in line with this inconsistency. Although the identified signature could segregate 'low' and 'high' gene expression groups in the gray-zone AST category, these two groups showed no significant differences in worrisome histopathological features when re-examining cases. Thus, patients diagnosed with a gray-zone AST category on histopathological grounds on forehand might receive more adopted treatment and ongoing follow-up in the future in line with the underlying molecular biology of the lesion.
Our results indicate that the identified signature genes from our cohort are of relevance in several altered signal pathways in cutaneous melanoma such as reported by the TCGA group [19]. However, one should be careful to draw too meaningful conclusions from our comparative analysis with the dataset from the TCGA. First, despite the large number of melanoma samples in the TCGA SKCM dataset the majority of tissue for molecular analysis originated from metastatic sites (~80%, [19]). In contrast, in our series only the tissue from case MST4 was derived from a lymph node metastasis. All other tissue for molecular analysis yielded from the primary MST samples. Second, in the TCGA SKCM dataset 89% of the primary melanoma cases were obtained from advanced-stage T4 melanomas with a Breslow thickness > 4 mm [19]. The authors reported that lower grade primary melanoma cases often had to be excluded from the investigation because there was not enough tissue available for further molecular analysis. In our series, the majority of MST represented lower stage melanomas (86%, n = 14, Breslow thickness between =/< 1.0 mm (T1) or 1.0-2.0 mm (T2); Table S1). Third, unfortunately the TCGA did not provide information about the histopathological subtypes of the investigated melanomas. Thus, a selection from the TCGA SKCM dataset of only Spitzoid melanoma subtypes for a comparative analysis with our cases was not possible.
In our study, the GSEA revealed downregulation of the Hallmark estrogen response early pathway and upregulation of the Hallmark spermatogenesis pathway in SN vs. MST (Fig. 2D). Indeed, eruptive SN have been documented during pregnancy and puberty, suggesting a potential role of hormonal deregulation in the pathogenesis of Spitz tumors [66]. Other case reports of disseminated SN also give hints to an endocrine pathogenesis with an affected patient suffering from Addison's disease [67]. Thus, our data provide further evidence that hormonal alterations are worthwhile to be considered in the pathogenesis of Spitzoid melanocytic neoplasms.
The performed GSEA also revealed significant alterations in the Hallmark epithelial-mesenchymal transition pathway and the Hallmark coagulation pathway in the SN vs. AST and AST vs. MST conditions. Interestingly, in SN vs. AST there was significant upregulation of both pathways whereas in the group of AST vs. MST there was significant downregulation of these pathways. The same trend was seen for the Hallmark myogenesis pathway, the Hallmark angiogenesis pathway, and the Hallmark complement pathway with upregulation in the group of SN vs. AST and downregulation in the group of AST vs. MST. The opposing regulations of these gene sets again point to the intermediate state of AST between SN and MST.
Although the performed GSEA revealed no differential expression of gene pathways related to immunomodulatory and inflammatory processes, there were several differentially expressed immune-related genes identifiable between the subgroups of Spitzoid melanocytic lesions (i.e., SN vs. AST vs. MST). These included IL11RA, IL10, IL12RB2, IL11RB, IL13RA2, CD40, IL20RRB, IL7, IL20RA, LEF1, and IL1R2 (Table S2). In an earlier study, we have compared the gene expression profile of SN in relation to conventional NCN and observed significant upregulation of gene pathways with increased expression of transcripts related to immunomodulatory and inflammatory processes in the SN group [15]. Whereas NCN are rarely inflamed or hosted by immune cells, a reasonable subset of Spitzoid melanocytic lesions, independent from the SN, AST, or MST group, is accompanied by a prominent inflammatory infiltrate of immune cells [1,68,69]. However, a more detailed characterization of this inflammatory microenvironment has not been performed yet. Interestingly, it is well known that similar to MST, AST frequently metastasize to locoregional lymph nodes, but in contrast to MST, rarely if ever spread to distant organs [6,70]. The mechanism behind this paradoxical behavior remains entirely unknown. As we have shown that there are only little molecular genetic differences in the tumor genome of the ambiguous AST category in comparison with the group of highly aggressive MST, one may speculate that the environmental immunologic neighborhood hosting the Spitzoid melanocytic neoplasms might define the lesion's biological behavior with either arrest of metastatic potential in the locoregional lymph node or development of distant metastasis.
We are currently aiming to validate the identified top-ranked gene transcripts of SN with immunohistochemistry to find out whether the signature is also discriminative on protein level. In the future, we also aim to validate this signature by Nano-String nCounter gene expression analysis in a larger group of SN, AST, and MST to find out whether it can contribute to diagnosis and ongoing treatment in daily surgical pathology practice and clinical patient management.

Conclusions
This is the first study employing transcriptomic Nano-String nCounter analyses as a novel tool in the evaluation of challenging Spitzoid melanocytic neoplasms. We identified differentially expressed genes in the group of SN compared to MST and AST. A molecular gene expression signature comprising the six topranked differentially expressed genes NRAS, EIF2B4, NF1, BMP2, FZD9, and IFNA17 could differentiate in between SN and MST. The identified gene transcripts belong to important pathways involved in the development of melanocytic tumors and melanoma: the MAPK pathway, the AKT pathway, and the Wnt/ beta-catenin pathway. Notably, AST samples showed intermediate levels of the signature, thereby implying that the six gene expression signature can potentially be used to distinguish more from less aggressive AST, thereby improving treatment decisions, and optimize clinical patient management in the future.
We also want to thank the reviewers and the editors for their expert input with very useful, constructive, and critical comments, which substantially helped us to revise the manuscript.
This work was partially supported by the academic incentive fund from the Maastricht University Hospital Center, MUMC+, the Netherlands.

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article.    Table S3. Library size of the studied cases (SN, AST and MST). Table S4. Characteristics of Identified Gene Pathways. Table S5. Gene expression values of the signature genes in the Skin Cutaneous Melanoma (SKCM) dataset provided by the Cancer Genome Atlas (TCGA, [19]).