The Transcriptional Landscape of BRAF Wild Type Metastatic Melanoma: A Pilot Study

Melanoma is a relatively rare disease worldwide; nevertheless, it has a great relevance in some countries, such as in Europe. In order to shed some light upon the transcriptional profile of skin melanoma, we compared the gene expression of six independent tumours (all progressed towards metastatic disease and with wild type BRAF) to the expression profile of non-dysplastic melanocytes (considered as a healthy control) in a pilot study. Paraffin-embedded samples were manually micro-dissected to obtain enriched samples, and then, RNA was extracted and analysed through a microarray-based approach. An exhaustive bioinformatics analysis was performed to identify differentially expressed transcripts between the two groups, as well as enriched functional terms. Overall, 50 up- and 19 downregulated transcripts were found to be significantly changed in the tumour compared to the control tissue. Among the upregulated transcripts, the majority belonged to the immune response group and to the proteasome, while most of the downregulated genes were related to cytosolic ribosomes. A Gene Set Enrichment Analysis (GSEA), along with the RNA-Seq data retrieved from the TCGA/GTEx databases, confirmed the general trend of downregulation affecting cytoribosome proteins. In contrast, transcripts coding for mitoribosome proteins showed the opposite trend.


Introduction
According to the definition reported by the Dictionary of Cancer Terms (https:// www.cancer.gov/publications/dictionaries/cancer-terms, accessed on 22 February 2022), melanoma is "A form of cancer that begins in melanocytes (cells that make the pigment melanin). It may begin in a mole (skin melanoma), but can also begin in other pigmented tissues, such as in the eye or in the intestines".
Overall, analysing the global incidence and mortality worldwide, melanoma is a relatively rare disease that affected 325,000 people in 2020, and 57,000 died because of the disease (source: Globocan 2020, https://gco.iarc.fr/today/online-analysis-table, accessed on 18 December 2021). Nevertheless, skin melanoma has a great relevance in certain countries, such as Australia and Europe, where the incidence rate is higher. In Europe, skin melanoma represents the seventh-most frequent malignancy and accounts for roughly 46% of the incident cases in the whole world (https://gco.iarc.fr/today/data/factsheets/ cancers/16-Melanoma-of-skin-fact-sheet, accessed on 18 December 2021). Moreover, it has been shown that some complex diseases, such as diabetes mellitus, might represent a risk factor for melanoma occurrence [1], and Endothelial Progenitor Cells (EPCs) are known to be selectively recruited within the tumour [2] and might represent a potential tool for therapy in this group of patients, especially the younger ones in which the EPC levels are higher [3].
Melanocytes can give rise to benign lesions called melanocytic naevi that can progress towards malignant lesions termed melanomas. Melanomas are classified according to the TNM staging system (AJCC staging manual 8th edition, issued in 2016 and updated in 2018) [4], and a global melanoma database has been released [5]. The data obtained from the clinical and pathological evaluations of melanoma are combined to divide patients into staging groups with different outcomes [6]. In addition, specific classification systems for melanoma were defined by Clark [7] and Breslow [8] long ago. The expression profiles and somatic mutations of advanced lesions and metastases have been defined and are reported in the Cancer Genome Atlas Network [9], while less is known about the initial phases of melanoma progression [10]. The most frequent genetic alteration described in melanoma is BRAF mutation, present in roughly 50% of the patients, according to the COSMIC database (Catalogue Of Somatic Mutations In Cancer) [11]. The vast majority of BRAF mutations are represented by a missense mutation named V600E, leading to the substitution of glutamic acid with a valine in codon 600 [9,12]. The final phenotype is characterised by the constitutive activation of the mitogen-activated protein kinase (MAPK) pathway, sustaining cell proliferation and preventing apoptosis. BRAF inhibitors (such as Vemurafenib and Dabrafenib) have been approved for the treatment of metastatic melanoma [13], since they significantly improve progression-free and overall survival, although resistance is rapidly acquired [14]. Currently, for patients not carrying BRAF mutations, no target therapy is available; therefore, a significant effort is needed to better define their molecular profile in order to identify potential molecular markers and targets for therapy. Therapy for this group of patients mainly relies on immunotherapy and checkpoint inhibitors [15][16][17], and searching for potential predictive biomarkers in response to biological agents is warranted in melanoma, as well as in other metastatic tumours [18].
The aim of this pilot study was to analyse the transcriptomic profile of patients suffering from metastatic melanoma without BRAF mutations in order to evaluate the possible differences in the expression profiles between advanced melanoma cells and non-dysplastic melanocytes serving as the healthy control.

RNA Extraction and Array Hybridisation
In order to obtain a comparative transcriptomic profile of melanoma cells relative to the healthy tissue (i.e., non-dysplastic naevi composed of healthy melanocytes representing the normal counterpart of melanoma cells), a pilot study was designed, and we performed a manual microdissection of paraffin-embedded surgical samples [19] of both normal tissue and BRAF wild type melanoma from patients whose clinicopathological characteristics are in Table 1. The samples then enriched in the melanocytic population ( Figure 1) were processed for RNA extraction and then hybridised on Agilent arrays after RNA validation (see the details in Section 4, Supplementary Figure S1).

Differential Expression Analysis
The gene expression data from the microarray experiments were pre-processed according to a standard pipeline, as described in Section 4. The so-obtained normalised log 2 expression data were then filtered and subjected to a differential expression analysis (DEA) using the rank product statistics (see Section 4 for details). In particular, n = 6 independent samples of BRAF wild type metastatic melanoma were compared to a reference array representing a pool of healthy tissues enriched with non-dysplastic melanocytes (biological average of n = 4 independent samples). Transcripts featuring a BH-FDR q-value < 0.05 and a |log 2 FC| > 0.5 were deemed as differentially expressed genes (DEGs). Overall, the DEA returned a list of 84 statistically significant probes, but only 69 of them could be annotated. Specifically, 50 up-and 19 downregulated transcripts were found to be significantly changed in the tumour compared to the control tissue, as reported in Table 2 and Table 3, respectively.     In order to further strengthen the transcriptomic analysis results, the differential expression of six DEGs (APOE, GAPDH, ACTB, RNA28S1N5, RPL31, and RPS17), selected on the basis of the most relevant fold changes and p-values emerging from the microarray experiments, was further assessed by RT-qPCR in a small subset of melanoma samples from the microarray analysis cohort (as described in the relative Section 4). Considering a p-value ≤ 0.05 to assess the significance, four of the six selected DEGs (ACTB, RPS17, RPL31, and RN28S1N5) showed statistically significant differences between the melanoma samples and healthy melanocytes, with fold change directions (expressed as relative expression fold changes performed by the ∆∆CT method) consistent with the microarray screening results (Supplementary Figure S2). GAPDH and APOE differential expression failed to reach statistical significance when analysed by RT-qPCR in the same sample subset but still maintained the fold change consistency with the microarray results (Supplementary Figure S2).

Genes Related to Antigen Processing and Presentation Are Upregulated in Tumour vs. Control
DEGs that were found to be upregulated in metastatic melanoma samples compared to non-dysplastic controls were tested for functional enrichment using the ToppFun web tool (https://toppgene.cchmc.org/, accessed on 22 February 2022). The full table of the statistically significant terms retrieved from such a query can be found as Table S1 in the Supplementary Materials section. Inspecting the results, it is noticeable how the top-most ranked functional terms and pathways were almost all related to the positive regulation of some features of the immune response, involving 32 out of the 50 upregulated DEGs resulting from the DEA. For example, the most relevant GO terms referring to biological processes (BPs) were innate immune response, defense response to other organism, regulation of immune system process, cell activation, response to external biotic stimulus, leukocyte mediated immunity, and antigen processing and presentation of exogenous peptide antigen via MHC class I (BH-FDR < 1.7 × 10 −4 ). Accordingly, the most significant KEGG pathway [20,21] was antigen processing and presentation (BH-FDR = 6.4 × 10 −4 ), accounting for five DEGs having a central role in the MHC class I pathway: HLA-A, HLA-C, CALR, PSME2, and HSP90AA1 ( Figure 2, genes in magenta). In cyan are the elements of the pathway additionally detected by the GSEA leading edge analysis. The GSEA also revealed a significant involvement of the proteasome complex (in green). Using the same colour code, the Venn diagram in the upper inset shows the complete lists of the official gene symbols found to be upregulated within the two KEGG pathways. The KEGG pathway map is hsa04612-antigen processing and presentation-Homo sapiens (human), modified and published with permission from Kanehisa Laboratories as the copyright holder.
Such a finding was confirmed by the GSEA (see Section 4), according to which, the gene set corresponding to this pathway was positively enriched (NES = 2.03, FDR q-value = 0.064, Figure 3A). More in detail, the leading edge analysis identified 13 main genes involved in the upregulation of both the MHC class I and class II pathways, thus extending the previous set of five DEGs detected on the basis of gene-wise hypothesis testing (Figure 2, genes in magenta plus genes in cyan). In addition, the GSEA pointed at a significant positive regulation of the proteasome complex (NES = 2.07, FDR q-value = 0.069, Figure 3B), another KEGG pathway term deeply connected to the previous one, with a leading edge featuring nine genes coding for different proteasome subunits, including pivotal proteasome activator subunit 2 (PSME2) as a linker between the two gene sets (Figure 2, green ellipse and Venn diagram).

Cytosolic Ribosome Proteins Are Downregulated in Tumour vs. Control
Strikingly, the vast majority of the transcripts found to be downregulated in metastatic melanoma compared to the control reference (i.e., non-dysplastic melanocytes) were related to cytosolic ribosomes (see Table 3). More in detail, 13 out of the 19 downregulated DEGs corresponded to ribosomal proteins (rProteins) of either the large (60S) or the small (40S) cytosolic ribosome subunit. In addition, two different probes targeting the product of the RNA28SN5 gene-the ribosomal RNA, giving rise to the 28S subunit-were among the 19 DEGs featured by the list of downregulated genes.
In order to confirm and extend these results, we ran a GSEA, testing the whole spectrum of rProteins of both cytosolic and mitochondrial origin. To do this, we took advantage of the already available Ribosomal Protein Gene Set (RPGS), which is the complete list of all human gene symbols related to ribosomes we assembled for a recent work in order to answer a similar scientific question [19]. Such an analysis confirmed a significant downregulation of the structural constituents of both the large (60S) and the small (40S) cytosolic ribosome subunits ( Figure 4A-C). On the contrary, and most interestingly, mitochondrial rProteins did not show any significant downregulation but, rather, an opposite trend ( Figure 4D-F).

External Validation through TCGA vs. GTEx Cohorts
In order to rule out any technical artefact related to microarray hybridisation or the sample origin, we decided to externally validate these findings using the UCSC Xena Browser (University of California, Santa Cruz, CA, USA, http://xena.ucsc.edu/, accessed on 24 February 2022) that provides a convenient way to access gene expression data stored in TCGA database for the comparative analysis of tumour samples with the normal analogies available from GTEx database (https://gtexportal.org/home/, accessed on 24 February 2022) [22,23]. TCGA samples were thus filtered based on cancer type (Skin Cutaneous Melanoma, SKCM), stage (metastatic), and genomic subtype (BRAF wild type). The so-obtained cohort featured 179 SKCM samples that were compared with the corresponding healthy GTEx cohort of normal skin tissue made out of 557 samples for a total sample size of n = 736. A thorough validation was carried out for the following sets of genes resulting from the corresponding GSEA leading edge analysis shown in Figures 3 and 4: MHC pathway (14 genes), proteasome (9 genes), cytosolic rProteins (14 genes), and mitochondrial rProteins (14 genes). Notably, almost all differential expressions we tested could be confirmed by TCGA/GTEx RNA-Seq data in terms of both the change direction (log 2 FC sign) and statistical significance ( Figure 5). The detailed validation scores are given in Table 4. normalised counts, gene-wise mean-centred, and shown as boxplots for the 14 cytosolic (magenta and teal) and the 14 mitochondrial (red and green) rProteins that emerged from the leading edge analysis of microarray data. Overall, 24 out of 28 differential expressions were confirmed by TCGA data, supporting the evidence of a generalised downregulation of the cytosolic rProteins and an overexpression of the mitochondrial ones. The four unconfirmed comparisons are marked with the symbol Ø (see Table 2 for more details). To control for the possible effects from the age and sex of the patients, these two covariates were also considered after downloading the specific metadata from the consortium portals (see Section 4). The results of such an analysis are presented as dot plots in Figure 6 and numerically in Supplementary Table S2 for the four gene sets of interest separately. Even if the overall dysregulation patterns could be substantially confirmed in all sub-cohorts, age appeared to be an important exacerbating factor (compare the Old Patients with Young Patients rows in the four panels of Figure 6). On the contrary, sex did not seem to be a discriminating factor, except, perhaps, in the downregulation of cytosolic rProteins, which were almost unaltered in the sub-cohort of Young Females (upper-left panel in Figure 6).

Discussion
In the pilot study reported in the present paper, the RNA extracted from paraffinembedded samples was hybridised on Agilent microarrays to assess the transcriptomic profile of BRAF wild type metastatic melanoma compared to the transcriptional reference of non-dysplastic melanocytes serving as the healthy control. Beyond the canonical DEA, followed by the functional enrichment analysis of the resulting DEGs, the GSEA was extensively used both to deepen the involvement of some relevant pathways of interest (i.e., MHC and proteasome) and to quantify the overall dysregulation of the whole rProtein spectrum. In addition, because of the limited sample size of this pilot study, and due to the particular nature of the RNA starting material, all our findings were validated querying the TCGA/GTEx gene expression databases to consider larger cohorts of patients, and at the same time, the hallmark gene sets from MSigDB were tested through a dedicated GSEA to check the consistency and the reliability of the expression levels, as measured by our microarray experiment (see below and Supplementary Table S3).
The GSEA and ToppFun functional enrichment analysis of the upregulated DEGs showed a consistent involvement of the immune system, with 70% of the overexpressed genes coherently annotated to some immune response-related process. This is in line with the well-known immunogenic nature of melanoma and the recent literature pointing at MHC-I/II protein expression as a powerful prognostic marker to predict the effectiveness of anti-CTLA-4 and anti-PD-1 immunotherapy in metastatic melanoma and other cancer types [24][25][26][27][28]. Even though most of these papers agree on the fact that a transcriptional downregulation of the MHC-I and MHC-II genes is a common feature of advanced untreated melanomas, in our study, the opposite seems to be true. Notably, this cannot be ascribed to some spurious effects induced by drugs-such as antibodies targeting the immune checkpoints-since all the samples we used for RNA extraction were excised from the patients before any therapeutic schedule. Moreover, microarray technical validation performed via RT-qPCR in a sample subset (Supplementary Figure S2 Figure S3). Rather, fold-change signs could be dependent on the particular stage at which melanoma samples were collected. Indeed, the MHC gene expression profile has already been reported to be heavily dependent on tumour progression, and its gradual loss is likely to facilitate the evasion of cancer cells from immune surveillance [29].
The other gene set we found to be upregulated in our cohort of metastatic melanoma patients compared to healthy controls was related to proteasomal function. Beyond its increased expression, the proteasome complex in melanoma cells may also be overactive because of the overexpressed PSME2 gene and the proteasome activator complex subunit 2 (aka PA28B), thus contributing, in turn, to the increased antigen presentation by the MHC class I pathway discussed above (see Figure 2). These data agree with the notion that melanoma cells heavily rely on proteasomal function to survive, so that selective proteasome inhibitors have already been used as new attractive therapeutics for this type of cancer [30][31][32].
On the other hand, the DEGs downregulated in the melanoma samples compared to the healthy controls were mostly related to cytosolic ribosomal proteins (rProteins). This was not completely unexpected given the accumulating evidence that relates cancer onset and progression with alterations of cell translational machinery [33]. Specifically, both enhanced and reduced ribosome biogenesis and protein synthesis have been reported to be associated with cancer in mammals, depending on the particular type of tissue and stage taken into account [34][35][36][37][38]. For this reason, rProteins configuration in metastatic melanoma was evaluated more in depth by running a GSEA of all the structural constituents of both cytosolic and mitochondrial ribosomes. Interestingly, the two ribosome types showed opposite patterns of deregulation: while cytosolic rProteins tended to be underexpressed in metastatic melanoma, the mitochondrial ones were sharply upregulated.
Such a finding is of great interest considering, in particular, the data we recently published in another paper addressing transcriptional alterations in colorectal carcinoma [19]. As in the present case, even in that study, we were able to find a consistent change in the rProtein expression but with a different FC sign, in that the upregulation of rProteins concerned cytosolic ribosomes and not the mitochondrial ones. Importantly, in both studies, all the ribosomal transcriptional alterations we reported found confirmation in the TCGA/GTEx large cohorts of patients.
As a final note on the DEG lists emerged from our analysis, it is worth noting that the set of significantly dysregulated genes was obtained through the hard thresholding (cut-off 0.05) of the whole transcriptome genes sorted by increasing q-values, with an additional cut-off on the fold change (|log 2 FC| > 0.5), as reported in Section 4. While this conservative approach is effective in reducing the number of false-positive hits (i.e., controlling for type I error), the statistical power may be affected, resulting in an increased number of false negatives, especially in the case of small sample sizes and RNA partial degradation (as in the case of paraffin-embedded samples). For this reason, it is not surprising that the DEG lists lack some genes whose dysregulation is expected in metastatic melanoma. Nevertheless, when the GSEA was performed testing the "Hallmark gene sets" (the H-collection of the MSigDB) and the weighted contribution from all the log 2 FCranked genes of the array was taken into account, many significant hallmark gene sets emerged related to: proliferation processes (MYC_TARGETS_V1 and P53_PATHWAY); cancer (mTORC1_SIGNALING, EPITHELIAL_MESENCHIMAL_TRANSITION, and UN-FOLDED_PROTEIN_RESPONSE); the immune system response (COMPLEMENT, ALLO-GRAFT_REJECTION, and INTERFERON_GAMMA_RESPONSE); and alteration of the metabolism (GLYCOLYSIS and OXIDATIVE_PHOSPHORYLATION). This is in excellent agreement with the results of other GSEAs performed on several different melanoma datasets (see, e.g., [39]), ultimately confirming the reliability of the gene expression profile returned by our microarray experiments (see Supplementary Table S3 for a complete list of the significant Hallmark Gene Sets, as returned by the GSEA).
As for the energy metabolic pathways involved in the early phases of melanoma pathogenesis, the key process is represented by glycolysis, and after the occurrence of BRAF mutations, the stimulation of transcription factors acting as key regulators of such process makes it even more effective [40,41]. Moreover, in BRAF-mutated cells, the Oxidative Phosphorylation (OXPHOS) is inhibited [42]. It is well-known that, between these two metabolic phenotypes, a dynamic switch occurs, and plasticity plays a key role in melanoma [43][44][45], leading to metabolic reprogramming of the cells. To make the picture more complex, it has been shown that some melanomas are able to exploit diverse nutrients and energy sources to adapt to different extracellular conditions, thus showing a "hybrid" glycolysis/OXPHOS metabolic phenotype [46][47][48]. Finally, the so-called "Reverse Warburg" effect has been described in melanoma cells [45,47,49]. This effect relies on the stimulation of cancer-associated fibroblasts (CAFs) that increase their glucose upload and lactate secretion through Monocarboxylate Transporter (MCT) family proteins [50]. Moreover, lactate can be internalised by cancer cells via MCT and conveyed into the Krebs cycle, thus fuelling OXPHOS. In this view, immune cells can deregulate metabolic pathways, representing a link between the deregulated pathways that emerged in this paper. As a further confirmation, the transcript of SLC66A1, encoding an MCT, was upregulated in our cohort (see Table 2).
All the proteins encoded by mitochondrial DNA are involved in the assembly and functioning of the respiratory complexes, along with the proteins encoded by nuclear DNA. For this reason, the OXPHOS biogenesis is subjected to a synchronised regulation of the mitochondrial and cytoplasmic ribosomes. Considering that-in contrast to the mitochondrial rRNA-mitochondrial rProteins are synthesised in the cytosol after the translation of mRNA of nuclear origin, the interplay between the nuclear and mitochon-drial components for ribosome production and the consequent synthesis of the various proteins involved in glycolysis, and OXPHOS is extremely complex. Since our data derive from a transcriptomic approach, they cannot give insights on the expression and function of the glycolytic and OXPHOS enzymes; therefore, no robust hypothesis on functional significance can be proposed. Nevertheless, these data could pave the road for further evaluations in which biochemical and physiological assays, together with proteomic and metabolomic approaches, can be used to define the activity and expression levels of the key glycolytic/OXPHOS enzymes.
The results reported in this paper might be relevant for two main reasons: (i) the central role of protein synthesis and energy metabolism in cancer and (ii) the fact that, despite the many recent reports about cytosolic ribosome aberrant function in cancer, there are still few data about the 55S mitochondrial counterparts and their functional interplay with 80S ribosomes. For example, the evidence of possible mitoribosome onco-patterns could provide a new rationale for the design (or repurposing) of novel antibiotics specific for cancer treatment, a still-debated clinical practice [51].
Taken together, our data point at a complete and deep remodelling of protein synthesis and degradation in metastatic melanoma that suggests, alongside biopsy genotyping, a more integrated evaluation of specific gene expression patterns-in particular, those related to MHC, proteasomes, and rProteins-as a practice that could help in choosing the most effective treatment in a context of personalised medicine.

Patients
Six patients (2 females, 4 males with mean age at diagnosis of 60.3 years, range 46-70) suffering from metastatic melanoma not harbouring BRAF mutations were enrolled for the study between April 2016 and October 2018 within the OMITERC study coordinated by Medical Oncology Unit, Azienda Ospedaliero-Universitaria Careggi (Florence). All the patients provided informed written consent, and the study was approved by the local Ethical Committee of Azienda Ospedaliero-Universitaria Careggi (BIO. 16 (Table 1). Moreover, 4 non-dysplastic naevi were also collected from the same institution as above.

Sample Preparation
In order to obtain tumour-enriched samples, paraffin-embedded specimens were manually micro-dissected, applying the same protocol as in [15]. Briefly, 20-µm-thick sections were cut, put on no positively charged slides, and counterstained with Meyer's Haematoxylin following the standard protocol. In order to achieve the enrichment of the tumour component of the metastatic melanoma samples, tumour areas were identified by an experienced operator (EL), collected through a sterile needle, and transferred to a fresh Eppendorf tube for further processing.

RNA Extraction and Quality Control
The total RNA was extracted from the enriched samples with the AllPrepDNA/RNA FFPE kit (Qiagen, Hilden, Germany), according to the manufacturer's protocols. The extracted RNA was then checked for its quality and integrity by the Agilent 2100 Bioanalyzer with the RNA 6000 Nano kit (Agilent Technologies, Santa Clara, CA, USA). The RNA concentration was also measured by a Nanodrop ND-1000 (Thermo Scientific, Waltham, MA, USA).

Microarray Hybridisation
A one-color microarray-based gene expression analysis was applied to analyse the RNA samples on the Agilent-026652 Whole Human Genome Microarray 4 × 44 K v2 platform (Agilent Technologies, Santa Clara, CA, USA), according to the manufacturer's protocols. To scan the microarrays, an Agilent G49000 DA SureScan Microarray scanner (Agilent Technologies, Santa Clara, CA, USA) was used, and subsequently, the data were extracted by Agilent Feature Extraction (Agilent Technologies, Santa Clara, CA, USA).

Differential Expression Analysis
Raw data obtained from microarray scanning were processed using Bioconductor software packages in the R environment. Briefly, fluorescence intensities were backgroundsubtracted, log 2 -transformed, and quantile-quantile normalised to get the gene expression. Based on the results of hierarchical clustering and a PCA on the samples, one array (ID Melanoma_5) was excluded from the subsequent steps of the analysis. Low-intensity probes (featuring a log 2 expression below 6.3 in more than one melanoma sample) were filtered out of the expression matrix as probes targeting unexpressed genes. Overall, 13,455 probes out of 34,127 (~40%) were retained at the end of the filtering procedure, and their log 2 expression values were tested for differential expression using rank product statistics. In particular, n = 5 melanoma biological replicates were compared against the single reference represented by the healthy biological mRNA pool of n = 4 non-dysplastic naevi (RankProd v3.18.0 Bioconductor package, one-sample design) [52][53][54][55][56]. p-values were adjusted for multiple comparisons, and all genes with a q-value (Benjamini-Hochberg False Discovery Rate, BH-FDR) < 0.05 were deemed as differentially expressed genes (DEGs) [57]. Finally, an additional cut-off on the fold changes (FCs) was applied to expunge from the DEG lists genes with a |log 2 FC| < 0.5.

Enrichment Analysis
The ToppFun web tool (by ToppGene Suite, https://toppgene.cchmc.org/, accessed on 22 February 2022) was used to analyse the DEG lists for functional enrichment through a hypergeometric hypothesis test [58]. All terms with a BH-FDR q-value < 0.05 were considered statistically significant. A Gene Set Enrichment Analysis was performed using GSEA v4.2.2 with the MSigDB database v7.5.1 (updated January 2022) [58,59]. Expression data from the microarray experiments were provided in the form of a pre-ranked list of genes (log 2 FC metric). Probes were collapsed into unique gene symbols before the analysis, and a standard (weighted) enrichment statistic was chosen. The Normalized Enriched Score (NES) and BH-FDR q-values are reported in the main text or Supplementary Materials for each gene set of interest. Within the context of the GSEA, the threshold of the q-value for a gene set to be considered statistically significant was set at 0.25. To evaluate the global transcriptional alterations affecting ribosomal proteins (rProteins), a custom gene set including all rProtein and rRNA genes was used. Details about such a custom Ribosomal Protein Gene Set (RPGS) have already been provided elsewhere [19].

TCGA/GTEx Validation
To provide an external validation of our main findings regarding the MHC pathway, proteasome, and rProteins, we used the UCSC Xena Browser (University of California, Santa Cruz, http://xena.ucsc.edu/, accessed on 24 February 2022) [22], which allows the direct comparison of tumour expression data stored in the database of The Cancer Genome Atlas (TCGA) consortium with healthy samples from the Genotype-Tissue Expression (GTEx) project database (https://gtexportal.org/home/, accessed on 24 February 2022) [23]. Specifically, we filtered TCGA data in order to keep samples only from the Skin Cutaneous Melanoma (SKCM) study of the metastatic type (TM, excised from patients not harbouring any BRAF mutation. As for the control group, all the normal skin samples retrieved from GTEx could be used. This led to a final comparison between n = 179 tumour samples and n = 557 normal tissues. This final cohort of n = 736 patients was then further characterised to take into account the contribution of age and sex as covariates. Namely, metadata for the TCGA consortium was downloaded with an ad hoc pipeline (https://github.com/MrHed mad/Edmund, accessed on 21 May 2022), while GTEx data was retrieved directly from the project's data portal. Using these metadata, the samples were divided into nine categories: all samples (regardless of metadata), biological male patients, biological female patients, old patients (defined as greater than 50 years of age at the cancer diagnosis or sample acquisition for the TCGA and GTEx samples, respectively), young patients (defined as not old), young male patients, young female patients, old male patients, and old female patients. For each sub-cohort, the expression of the genes of interest (see Table 4 and Figure 6) was compared between the cancer (TCGA) and healthy (GTEx) samples. A comparison was performed with the Student's t-test, double-sided, and the groupwise family error rate was corrected with the Benjamini-Hochberg procedure.

Real-Time PCR Validation
In order to further validate the transcriptomic data, 6 DEGs were selected according to the log 2 FC and p-values (see Tables 2 and 3) for further validation by RT-qPCR. The selected gene expressions were assessed by using commercially available KiCqStart SYBR Green Primer pairs (Merck Millipore, Burlington, MA, USA), following the manufacturer's protocols. Relative expression quantification was performed by the ∆∆CT method [60] using the gene expression normalisation approach identified in [61] for the melanoma samples. Gene expression was assessed in a small melanoma sample subset (n = 3) from the microarray cohort compared with the healthy control primary epidermal melanocytes.

Acknowledgments:
The authors thank tLuisa Di Cerbo, Lorenzo Antonuzzo, Lorenzo Leoncini, Michele Maio, and Felice Arcuri for coordinating the sample collection. The authors also thank Barbara Stecca for the kind gift of RNA from the primary melanocytes.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of the data; in the writing of the manuscript; or in the decision to publish the results.