Unraveling the Role of microRNA and isomiRNA Networks in Multiple Primary Melanoma Pathogenesis

Background Malignant cutaneous melanoma (CM) is a potentially lethal form of skin cancer whose worldwide incidence has been constantly increasing over the past decades. During their lifetime about 8% of patients with CM will develop multiple primary melanomas (MPM). Patients affected by MPM could have a genetically determined susceptibility, though germline mutations in hereditary melanoma genes are rarely detected. Methods To better characterize the biology of this subset of melanomas, we explored the miRNome of 24 single and multiple primary melanomas, including multiple tumors from the same patient, using a smallRNA sequencing approach and bioinformatic detection of miRNA isoforms. The differential expression of specic miRNAs/isomiRs was obtained using quantitative PCR. Results From a supervised analysis, 22 miRNAs were differentially expressed in MPM compared to single CM, including key miRNAs involved in epithelial-mesenchymal transition (EMT). Moreover, the rst and second melanoma from the same patient presented a different miRNA prole. Ten miRNAs, including miR-25-3p, 149-5p, 92b-3p, 211-5p, 125a-5p, 125b-5p, 205-5p, 200b-3p, 21-5p and 146a-5p, were further validated in a larger cohort of single and multiple melanoma samples (N=47). Overall, the Pathway Enrichment Analysis revealed a more differentiated and less invasive status of MPMs. Analyzing our smallRNA seq data, we detected a panel of melanoma-specic miRNA isoforms (isomiRs), which were validated in The Cancer Genome Atlas SKCM cohort. Specically, we identied hsa-miR-125a-5p|0|-2 isoform as 10-fold over-represented in melanoma and differentially expressed in MPMs. IsomiR-specic target analysis revealed that the miRNA shortening confers a novel pattern of target gene regulations, including genes implicated in melanocyte differentiation and cell adhesion. Conclusions Overall we provided a comprehensive characterization of the miRNA/isomiRNA regulatory network of multiple primary melanomas, highlighting mechanisms of tumor development and correlating in an average time of 33 months (range 3–98). MPM patients were tested for germline genetic alterations in CDKN2A, CDK4 and MITF gene 15 and only one patient was found to have a germline CDKN2A mutation (c.249C > A p.His83Gln exon 2) of unknown clinical signicance. Melanoma specimens were examined by two dermato-pathologists.

Regarding genetic factors, somatic mutations of BRAF gene have been found in almost 40-50% of sporadic CMs located in body sites with intermittent UV exposure; 15-20% of the other cases are associated to NRAS mutations and correlated with chronic UV exposure 24 . A small portion of melanomas occurs in acral or mucosal locations and a subset of them are related to KIT and GNAQ mutations 7 . These ndings have brought important therapeutic implications and changed the management of CM patients with the development of speci c target therapies. Germline mutations instead, can be found in multiple or familial cases of CM. The most frequently described germline mutation is in CDKN2A (cyclin-dependent kinase inhibitor 2A) gene occurring in 8-15% of subjects diagnosed with multiple primary melanomas (MPMs) without familial history and up to 40% of patients with hereditary CM 6,25,40,44 . Mutations in other susceptibility genes such as CDK4 (cyclin-dependent kinase 4), MITF (microphtalmia-associated transcription factor) and POT1 (protection of telomeres 1) are less frequently detected 4,14 . During their lifetime about 8% of patients with cutaneous melanoma will develop multiple primary melanomas, usually at a young age and within 3 years from the rst tumor/diagnosis 18 .
The occurrence of MPMs in the same patient is thought to be related to a personal genetic susceptibility in association with environmental factors. These patients may represent a model of high-risk CM occurrence. As a matter of fact, it is estimated that a personal history of CM is a strong risk factor for the development of a subsequent primary CM 18,43 . The excision of a prior CM determines a risk up to 8.5% to develop another CM and the frequency of MPMs is reported to be between 0.2 and 10% 13,25,29 . The above-reported rates may underestimate the lifetime rates due to limited series of patients and different follow-up periods. Variability may also arise due to differences in environmental factors such as ultraviolet radiation exposure across geographical regions. Among the cases of MPMs, 13-40% of patients are diagnosed with synchronous lesions (i.e. a subsequent primary CM diagnosed within 3 months from the prior diagnosis), while the remainder develop metachronous lesions 1,29,41 . The risk of a subsequent CM is highest in the rst year following the diagnosis of the primary CM; however, this risk remains increased for at least 20 years 1  Moreover, the frequency of germline mutations in melanoma susceptibility genes (CDKN2A, CDK4, MITF, POT1/ACD/TERF2IP, TERT, BAP1) is lower than expected in MPM patients 5,6,9,25 . Therefore, a better characterization of MPM pathogenesis and biological features is of the outmost importance.
The dysregulation of small noncoding RNAs, speci cally microRNAs (miRNAs, 18-22 nucleotides in length), plays a signi cant role in tumorigenesis, including melanoma onset and progression 46 . MiRNAs regulate multiple and speci c target genes, determining an oncogenic or tumor-suppressive function, being implicated in the proliferation, apoptosis and tumor progression. Moreover, miRNA global expression pro le faithfully re ects the overall expression pro le of normal and pathological cells and tissues, with the advantage to be feasible also from formalin-xed and para n-embedded (FFPE) tissues.
In this study, we investigated the global miRNA and isomiRNA expression pro le of multiple primary melanomas using an unbiased smallRNA sequencing approach. A comparison of familial/non familial MPM vs. single primary melanoma miRNome was established in order to investigate the possible similarities. Moreover, the evolution of MPM miRNA pro le was assessed matching multiple tumors from the same patient.

Clinical samples
A retrospective series of 47 samples from 29 patients was collected. Patients were selected among those referring to the melanoma center of the Dermatology Unit at Bologna University Hospital. The study was approved by Comitato Etico Indipendente di Area Vasta Emilia Centro -CE-AVEC, Emilia-Romagna Region (number 417/2018/Sper/AOUBo). Before study entry, all the patients provided written and voluntary informed consent for inclusion, collection and use of clinical-pathological data and samples and data privacy.
The specimens were classi ed into three groups: benign nevi, single primary cutaneous melanoma (CM) and multiple primary melanoma (MPM). Group 1 (n = 3), benign nevi of 3 patients with no prior diagnosis of CM or non-melanoma skin cancer and follow up of at least 10 years. Group 2 (n = 35), MPM samples from 17 patients with prior diagnosis of ≥ 2 CMs. 3 out of 17 patients had positive family history of CM (FAM). MPM patients were tested for CDKN2A, MITF and CDK4 genetic alterations and only 1 patient had a mutation in CDKN2A gene (c.249C > A p.His83Gln). Group 3 (n = 9), 9 samples from CM patients with no history of prior CMs and a follow-up of at least 10 years.
Tumor and nevi samples were formalin-xed and para n-embedded (FFPE). For each sample, 5/6 tissue sections on glass slides were obtained. One section was stained with hematoxylin-eosin (HE) and examined by an expert pathologist to select the tumor/nevus area, which was grossly dissected before RNA extraction.

RNA extraction
RNA was isolated from 10 µm-thick FFPE sections using miRNeasy FFPE kit (Qiagen) according to the manufacturer's instructions. Depara nization was performed with xylene followed by an ethanol wash. RNA was eluted in 30 µL of RNAse-free water and quanti ed by absorbance at 260 and 280 nm.

SmallRNA sequencing
We analyzed 3 benign nevi, 4 single CM, 17 multiple primary or familial melanomas from 8 different patients. The 24 smallRNA libraries were generated using TruSeq Small RNA Library PrepKit v2 (Illumina, RS-200-0012/24/36/48), according to manufacturer's indications. Brie y, 35 ng of puri ed RNA was linked to RNA 3' and 5' adapters, converted into cDNA, and ampli ed using Illumina primers containing unique indexes for each sample. High Sensitivity DNA kit was adopted for libraries quanti cation using Agilent Bioanalyzer (Agilent Technologies, California, USA5067-4626) and the 24 DNA libraries were combined in equal amount to generate a libraries pool.
Pooled libraries underwent to size selection employing magnetic beads (Agencourt) and amplicons with a length in the 130-160 bp range, were recovered.
Finally, 20pM of pooled libraries, quanti ed using the HS-DNA Kit (Agilent) were denatured, neutralized and combined with a Phix control library (standard library normalizator). A 1.8 pM nal concentration of pooled libraries (obtained by dilution with a dedicated buffer as described in Illumina protocol guidelines) was obtained and sequenced using NextSeq 500/550 High Output Kit v2 (75 cycles) (Illumina, FC-404-2005) on the Illumina NextSeq500 platform.
Raw base-call data were demultiplexed using Illumina BaseSpace Sequence Hub and converted to FASTQ format. After a quality check with FastQC tool, the adapter sequences were trimmed using Cutadapt, which was also used to remove sequences shorter than 16 nucleotides and longer than 30 nucleotides. Reads were mapped using the STAR algorithm. Only reads that mapped unambiguously to the genome (at least 16 nucleotides aligned, with a 10% mismatch allowed) were used for the analyses. The reference genome consisted in human miRNA sequences from the miRbase 21 database. Raw counts from mapped reads were obtained using the htseq-count script from the HTSeq tools 3 . Counts were normalized using DESeq2 bioconductor package 37

Quanti cation of isomiRs
IsomiRs were identi ed in our NGS dataset of 24 samples as described in Loher et al. 36 . Brie y, sequence reads were quality trimmed using the cutadapt tool, and mapped unambiguously using SHRIMP2 (PMID: 21278192) to the human genome assembly GRCh38. During the mapping, no insertions or deletions, and at most one mismatch was permitted. IsomiRs were identi ed as done previously 36,52 .
For TCGA isomiR analysis, short RNA-seq Aligned BAM les were downloaded from the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/) for all 32 cancer types. IsomiR pro les were generated using the same approach as described in Loher et al. 36 To simpli ed the labeling of the isomiRs, we used the annotation system developed by Loher et al. This nomenclature speci es the name of the canonical miRNA, the start site (5' end) of the isomiR compared to the canonical miRNA sequence in miRBase, the end-site (3' end) and the eventual insertion of uracil. In particular, to annotate the start and end site of an isomiR, a negative (-) or positive sign (+) followed by the number of nucleotides is used to indicate how many nucleotides the isomiRs terminus has, when compared to the canonical miRNA sequence. Zero indicate the same terminus of the canonical miRNA sequence.
We quanti ed isomiR abundances in reads per million (RPM). Only reads that passed quality trimming and ltering and could be aligned exactly to miRNA arms were used in the denominator of this calculation.
IsomiR targets were predicted using the RNA22 algorithm 35 and targets were allowed to be present in the 5´UTR, CDS, and 3´UTR of the candidate mRNA. We selected only those targets that had a p-value < 0.01 and a predicted binding energy <-16 while also allowing G:U wobbles and bulge's within the seed region.

Statistical analysis
Normalized sequencing data were imported and analyzed in Genespring GX software (Agilent Technologies). Differentially expressed miRNAs were identi ed using a fold change > 1.5 lter and moderated t-test (FDR 5% with Benjamini-Hochberg correction) in CM vs. MPM comparison and foldchange > 1.2 and paired t-test (p < 0.05) in 1st vs. 2nd MPM comparison. Cluster Analysis, was performed using Manhattan correlation as a similarity measure. Principal Component Analysis was performed on 24 samples using all human miRNAs detected by NGS analysis (n = 1629).
Graphpad Prism 6 (GraphPad Software) was used for statistical analyses. Group comparison was performed using unpaired t-test, when data had a normal distribution, with or without Welch's correction according to the signi cance of the variance test. Data that did not present a normal distribution were compared using Mann-Whitney non-parametric test.
Association of gene expression with overall survival in TCGA SKCM cohort, was obtained using Oncolnc website (http://www.oncolnc.org), logrank test was used to calculate the p-value.

Pathway Analysis
Pathway and network analysis of differentially expressed miRNAs, miR-125a-5p isoforms and their targets was investigated using the web-based software MetaCore (GeneGo, Thomson Reuters). A p value of 0.05 was used as a cut off to determine signi cant enrichment.

Patient Characteristics
Demographic, clinical and pathological features of 29 patients are summarized in Table 1. A total number of 16 males and 13 females were included, with a mean age at rst diagnosis of 59 years for single primary melanomas and 53 years for multiple primary melanomas. Nine patients had single cutaneous melanoma and 10 years of follow up; 17 developed more than one primary melanoma in an average time of 33 months (range 3-98). MPM patients were tested for germline genetic alterations in CDKN2A, CDK4 and MITF gene 15 and only one patient was found to have a germline CDKN2A mutation (c.249C > A p.His83Gln exon 2) of unknown clinical signi cance. Melanoma specimens were examined by two dermato-pathologists.
The microRNA pro le of multiple primary melanoma The global miRNA pro le of 17 multiple primary melanomas, obtained from 8 patients, was analyzed using a smallRNA sequencing approach. For each MPM patient we analyzed the rst and second primary tumor, and for one case also a third one. Three patients had a family history of melanoma. We compared the global miRNA pro le of MPM toward that of 4 single melanomas and 3 benign nevi. From the smallRNA sequencing data, we identi ed 1629 mature miRNAs expressed in melanoma and nevus cells.
The unsupervised Principal Component Analysis (PCA) of all miRNAs and all samples (n = 24) revealed that familial and non-familiar multiple primary melanomas have a greatly overlapping miRNA pro le (Fig. 1A), which is different from single cutaneous melanoma (CM) and benign nevi (BN). Indeed, a statistical comparison between familial and non-familial MPMs did not provide any signi cant result. Therefore, we considered familial and non-familial melanomas as a unique group in all subsequent analyses. From the PCA we can already observe that MPMs displayed a miRNA pro le more similar to benign nevi than CMs. When we compared multiple and single melanoma tumors, we obtained a markedly different miRNA expression pro le and a list of 22 miRNAs differentially expressed (adjusted p < 0.05, Table 2), which are represented with a Volcano plot in Fig. 1B. Cluster analysis of these samples based on the expression of the 22 differentially expressed miRNAs con rmed the separation between single and multiple tumors Fig. 1C. The MPM group was constituted by the paired rst and second tumors (and one additional tumor in one case) developed by the same patient over the years (Table 1).
Comparing the miRNA pro le of these two groups using a paired statistical analysis, miRNAs that characterize the second tumors, usually thinner and less aggressive than the rst melanoma given their early diagnosis, were identi ed. Despite the similarities between the two matching MPMs, a variation in miRNA expression was observed (Fig. 1B). Speci cally, thirty-seven miRNAs were differentially expressed between the rst and second MPM (paired t-test, p < 0.05, Table 3) and a signi cant separation was obtained applying the cluster analysis (Fig. 1D).  Validation of microRNA differential expression in single and multiple primary melanomas and paired primary tumors from the same patient Nine miRNAs were selected for an independent technical validation using quantitative RT-PCR in 47 novel samples including BN, CM and 1st and 2nd MPM. Speci cally, we included the miRNAs differentially expressed between CM and MPM (miR-21-5p, miR-25-3p, miR-125b-5p, miR-146a-5p, miR-205-5p, miR-149-5p) and others between the rst (MPM 1st ) and second (MPM 2nd ) melanoma within the same MPM patient (miR-149-5p, miR-92b-3p, miR-200b-3p, miR-125a-5p).
According to the smallRNA NGS results, an upregulated expression of miR-21-5p, miR-25-5p, miR-146a-5p, and a downregulated expression of miR-125b-5p, miR-149-5p and miR-205-5p in CM compared to MPM were expected. In MPM samples, all selected miRNAs are upregulated in the MPM 2nd compared to MPM 1st. In the validation experiment, we included also miR-211-5p, considering that it is a melanocyte speci c miRNA, whose genetic locus is located inside melastatin gene and whose expression is particularly high in nevi. The expression of this miRNA was higher in BN, with borderline statistical signi cance when compared to CM or MPM in our NGS data.
The validation was performed in a cohort of 29 patients, as described in Table 1, using miR-16-5p as a reference gene due to its invariant expression in NGS data. Expression distributions of selected miRNAs in benign nevi, cutaneous melanoma and multiple primary melanoma samples are represented in Fig. 2 upregulation of miR-200b-3p and miR-205-5p was observed in MPM. A similar trend can be observed for miR-149-5p. As expected, the melanocyte speci c, MITF-regulated miR-211-5p is progressively downregulated in multiple and single melanomas (Fig. 2).
Functional annotation of multiple primary melanoma miRNA signature The list of 22 miRNAs differentially expressed in multiple vs. single melanomas, was uploaded in Metacore software (Clarivate Analytics) to identify both the pathways that are signi cantly regulated by these miRNAs (Supplementary Table 1, Additional le 1) and the most signi cant miRNAs/targets networks ( Supplementary Fig. 1A, Additional le 2).
Multiple primary melanomas were found to have a higher expression or miR-200 family, miR-205-5p and miR-149-5p compared to single CM and even nevi ( Fig. 2 and Fig. 4). These microRNAs target ZEB1/TCF8 and ZEB2/SIP1 genes, and by doing so they inhibit the epithelial-mesenchymal transition (EMT) pathway. This pathway therefore appears to be speci cally activated in single melanomas (Fig. 4).
From MetaCore network analysis, three hub genes (TLR4, ITGA6 and BTG2) were identi ed as targeted by multiple miRNAs, either up-or down-regulated in multiple melanomas. When we assessed the association of TLR4, ITGA6 and BTG2 gene expression with melanoma prognosis, we observed that their higher expression (median cutoff) was signi cantly associated with a worse overall survival in TCGA SKCM cohort of 458 samples ( Supplementary Fig. 1B, Additional le 2).
IsomiRNA analysis revealed that miR-125a-5p isoforms are dysregulated in multiple primary melanoma Interestingly, miR-125a-5p differential expression in MPM was not con rmed by qPCR technology and we wondered about a possible explanation. We observed that the reads generated by the smallRNA sequencing experiment and attributed to mature miR-125a-5p following the standard matching pipeline were actually shorter by 1, or most frequently 2 nucleotides (lack of GA at the 3' end) in all samples (Supplementary Fig. 2A, Additional le 3). Although miRBase database reports a unique mature sequence for each miRNA, the so called canonical form, many evidences from deep sequencing experiments suggest that miRNAs have frequent modi cations in length and sequence in human tissues. These miRNA isoforms are called isomiRs 11 . We analyzed the isomiR expression level in all single and multiple primary CMs from our NGS experiments. We found 90 miRNAs with sequence and length heterogeneity, generating 324 different isomiRs, and 40 canonical microRNAs without any isomiR. In addition, we found 40 isomiRs named "orphan", because their canonical miRNA sequences could not be detected. For each isomiR, we calculated the average expression in melanoma samples and the ratio between each isomiR and its canonical miRNA. Finally, we obtained a panel of 17 miRNAs whose isoforms are 3-to 10-fold more abundant in melanoma than their canonical form (Table 4). Among them, hsa-miR-125a-5p|0|-2 isoform was differentially expressed in multiple vs. single primary melanomas and between the rst and second tumor of the same patient (paired t-test P = 0.0006). Unusually, miR-125a-5p canonical and 3' shorter isoforms show an opposite expression trend in nevi, single and multiple primary melanomas (Fig. 5A). We studied two different technical approaches for miRNA quanti cation based on qPCR (miRCURY LNA and miSCRIPT, both by Qiagen), to selectively quantify miR-125a-5p isoforms in all samples and validate the NGS data. Speci cally, we used miR-125a-5p miRCURY LNA assay (Exiqon/Qiagen) for the quanti cation of the canonical, 24nt-long isoform ( Supplementary Fig. 2B, Additional le 3). Results revealed a lack of variation of this mature isoform between single and multiple melanomas, and a higher expression in the rst vs. second melanoma (Fig. 5B,C). To quantify the miR-125a-5p 22nt-long isoform, we selected the miSCRIPT assay by Qiagen. The assay can quantify both the long and short isoforms of miR-125a-5p due to the use of a universal 3' primer for miRNA ampli cation. Given the high predominance of the short isoform in our NGS data, we assumed this assay could provide a bona de quanti cation of the short 22nt-long isoform ( Supplementary Fig. 2B, Additional le 3). As expected, an increase in miR-125a-5p levels in MPMs vs. CMs and in the second tumor from the same patient was observed (Fig. 5B,C). We examined the expression of hsa-miR-125a-5p|0|-2 and 0|0 (WT) isoforms across TCGA tumor types and discovered an overall higher expression of the shorter form in human cancers and a speci cally altered ratio of the two forms in SKCM (cutaneous melanoma cohort), which shows the largest variation (Fig. 6).

Discussion
The risk of melanoma development is in uenced by environmental and genetic factors. Families with history of melanoma could have a germline mutation that confers hereditary susceptibility, and this is particularly demonstrated in families where more members develop multiple primary melanomas. In 1968, Lynch and Krush described the familial atypical multiple mole-melanoma (FAMMM syndrome) which encompasses an association between pancreatic cancer, multiple nevi, and melanoma 38 . In the 70 s Clark described a similar phenotype, the B-K mole syndrome, consisting of familial melanoma in the setting of numerous atypical nevi. In the early 1990′s, germline mutations in the cell cycle gene, p16 (CDKN2A), were reported among a subset of FAMMM kindreds. Nowadays, most studies report a very low prevalence of CDKN2A /CDK4 in familial or multiple melanoma patients, especially in the Southern Europe countries 9 . Though MPM patients often report similar sun exposure experiences, the high percentage of atypical nevi in these patients and their family members, the frequent family history of melanoma, as well as the early onset of melanoma (young age) suggest that predisposing factors for the development of multiple melanomas are involved. Regardless of family history, they are reported also cases of multiple primary melanoma in individuals without familial history of melanoma. In these cases, germline mutations in melanoma predisposing genes are rarely detected. Therefore, it is evident that some other genetic or epigenetic factor is active in multiple primary melanoma to fuel multiple events of melanocytic transformation.
In this study, we provide the rst comprehensive molecular characterization of MPMs by assessing their miRNome with a smallRNA sequencing approach. The global microRNA expression re ects the mRNA expression of cells and tissues, with the advantage of being assessable in FFPE tissues. This analysis revealed a speci c expression pattern of multiple melanoma tumors when compared to single cutaneous melanoma. MPM miRNome is more similar to benign nevi, thus suggesting a less aggressive and more differentiated phenotype. We validated a panel of microRNAs in additional samples, including also multiple tumors from the same patient, obtaining a panel of microRNAs differentially expressed in tumors from the same patient. We provide here evidence that MPMs, from a biological point of view, have a less invasive phenotype as pointed out by the main regulatory pathways activated in these tumors, thus providing further elements of discussion to support MPM less aggressive evolution. It is worth mentioning that microRNAs known to inhibit epithelial-mesenchymal transition (e.g. miR-200 family, miR-205, miR-149) are more expressed in multiple primary melanoma compared to single melanoma. Tumor cells promote EMT to escape from the microenvironment and migrate to a new location to develop metastasis 22 . The acquisition of a mesenchymal phenotype promotes the production of extracellular matrix proteins, the resistance to apoptosis, the invasiveness and the migration 30 . EMT results from the loss of cell-to cell junctions, induced by the loss of E-cadherin; the process is mediated by transcription factors, including SNAIL, SLUG, SIP1, and E2A, and affected by regulatory proteins such as TGFβ, EGF, PDGF, ERK/MAPK, PI3K/AKT, SMADS, RHOB, β-catenin, LEF, RAS, C-FOS, integrins β4 and integrin α5 53 ., EMT has been reported in melanoma cells, despite their origin from neural crest-derived melanocytes. In fact, EMT promotes the metastatic phenotype of malignant melanocytes 2 8 . Moreover, melanocytes express Ecadherin, which mediates the adhesion between melanocytes and keratinocytes 51 . Many studies described the loss of E-cadherin in melanoma 32,49 , and CDH1 ectopic expression was associated with the downregulation of adhesion receptors, such as MCAM/MUC18 and β3 integrin subunit, resulting in suppression of melanoma cells invasion 26 . Hao et al. observed a switch from E-cadherin to N-cadherin expression in melanoma progression, a process regulated by PI3K and PTEN through TWIST and SNAIL 23 .
Consistently, we examined the main cellular hubs regulated by MPM speci c miRNAs and discovered that they are centered in TLR4, ITGA6 and BTG2 proteins. MicroRNAs regulating these hubs are mostly downregulated in MPMs and high expression of these three genes is associated with a more favorable prognosis in TCGA SKCM cohort.
Integrin α6 (ITGA6), also known as CD49f, is a transmembrane glycoprotein adhesion receptor that mediates cell-matrix and cell-cell interactions. ITGA6 was identi ed and described as an important stem cell biomarker. Indeed, it is the only common gene expressed in embryonic stem cells, neural stem cells and hematopoietic stem cells 28,45 . It is also expressed in more than 30 stem cell populations, including cancer stem cells 31 . ITGA6 can combine with other integrins such as integrin beta 1 and integrin beta 4 to form respectively integrin VLA-6 and TSP180. The role of ITGA6 in melanoma is not clear but our observation point toward its upregulation in MPMs upon miR-25 and 29 downregulation.
BTG2 is part of the anti-proliferative BTG/TOB family and its expression is p53-dependent 47 . This protein is involved in several cellular processes, including cell cycle regulation, DNA damage repair, cell differentiation, proliferation and apoptosis. However, its role is often cell-type dependent 39 . In fact, BTG2 inhibits proliferation and migration, acting as a tumor suppressor protein, in gastric cancer cells 59 and in lung cancer cells 57 , while in bladder cancer it promotes cancer cell migration 55 . In B16 melanoma cells it was shown that miR-21 promotes a metastatic behavior through the downregulation of many tumor suppressor proteins, including PTEN, PDCD4 and BTG2 58 . In MPMs, we observe the downregulation of several miRNAs targeting BTG2, including miR-21-5p, 146a-5p, 132-3p, 15a-5p. Therefore, an upregulation of BTG2 is to be expected.
Toll-like Receptor 4 (TLR4) belongs to TLR family and plays an important role in in ammation and cancer. TLR4 protein is expressed at very low levels in melanoma cells in vivo (Human protein atlas) but its activation has been reported to promote an in ammatory microenvironment and tumor progression in vitro 20 . In addition, TLR4 is associated with induction of proliferation and migration of melanoma cells 50 . TLR4 plays an important role in melanoma also because it interacts with TRIM44, a negative prognostic factor in melanoma. In particular, TRIM44 binds and stabilizes TLR4 leading to the activation of AKT/mTOR signaling, which results in EMT promotion 56 . This biological role for TLR4 in melanoma is partially in contrast with our observation of a better survival in melanoma patients with higher TLR4 levels.
Finally, we extended our molecular investigation to miRNA isoforms that were most abundant in our samples. According to the recent observation that miRNA isoforms can discriminate human cancers 52 , we detected a relevant number of miRNA variants in our dataset of single and multiple melanomas. A speci c isoform of miR-125a-5p, lacking 2 nucleotides at the 3'end, was detected as differentially expressed in MPMs. This isoform is highly abundant in melanoma, as we con rmed by analyzing its levels across 32 tumor types from TCGA database; and the ratio between miR-125a-5p isoform and canonical form is the broader in TCGA SKCM tumors (range 0.1-1100 times) and 2-6 logs more abundant in nevi and melanomas in our study. Moreover, we detected a speci c dysregulation of the isoform, but not the canonical form, in multiple melanomas. Bioinformatic analyses revealed that miR-125a-5p shorted isoform loses the ability to target and regulate a group of genes speci cally involved in cell adhesion and cell differentiation. Particularly relevant seems to be the lack of regulation of genes involved in neuronal differentiation. Indeed, miR-125a is the human ortholog of lin-4, the very rst miRNA identi ed in C. Elegans in 1993 34 . In mammals, miR-125 is expressed in embryonic stem cells and promote cell differentiation. Speci cally, miR-125 has a speci c role in adult nervous system development and neuronal differentiation 10,33 . The imbalance between major miR-125 isoforms in melanocytes could re ect a major role for miR-125 in melanocyte development and differentiation from the neural crest 42 , differentiating this lineage from other common ancestor cells. A role that is consequently re ected in melanoma development and progression.

Conclusions
Overall, we provide here a comprehensive characterization of microRNA/isomiRNA dysregulation and regulatory network in single and multiple primary melanomas. The pattern of miRNA alterations supports a less aggressive phenotype of multiple primary melanomas, whilst isomiR-125a-5p levels proved to be        Differential microRNA expression in 1st vs. 2nd melanoma from the same patient. Before-after plot showing the paired expression of 4 selected microRNAs in 17 multiple primary melanoma (MPM) patients. miR-92b-3p, miR-205-5p, miR-200b-5p and miR-149-5p are signi cantly downregulated in the 1st melanoma compared to the 2nd melanoma. Each miRNA was tested in triplicate by quantitative RT-PCR.
Relative miRNA expression was normalized on invariant miR-16-5p. Paired P-value is reported.

Figure 3
Differential microRNA expression in 1st vs. 2nd melanoma from the same patient. Before-after plot showing the paired expression of 4 selected microRNAs in 17 multiple primary melanoma (MPM) patients. miR-92b-3p, miR-205-5p, miR-200b-5p and miR-149-5p are signi cantly downregulated in the 1st melanoma compared to the 2nd melanoma. Each miRNA was tested in triplicate by quantitative RT-PCR.
Relative miRNA expression was normalized on invariant miR-16-5p. Paired P-value is reported.