MicroRNA Isoforms Contribution to Melanoma Pathogenesis

Cutaneous melanoma (CM) is the most lethal tumor among skin cancers, and its incidence is constantly increasing. A deeper understanding of the molecular processes guiding melanoma pathogenesis could improve diagnosis, treatment and prognosis. MicroRNAs play a key role in melanoma biology. Recently, next generation sequencing (NGS) experiments, designed to assess small-RNA expression, revealed the existence of microRNA variants with different length and sequence. These microRNA isoforms are known as isomiRs and provide an additional layer to the complex non-coding RNA world. Here, we collected data from NGS experiments to provide a comprehensive characterization of miRNA and isomiR dysregulation in benign nevi (BN) and early-stage melanomas. We observed that melanoma and BN express different and specific isomiRs and have a different isomiR abundance distribution. Moreover, isomiRs from the same microRNA can have opposite expression trends between groups. Using The Cancer Genome Atlas (TCGA) dataset of skin cancers, we analyzed isomiR expression in primary melanoma and melanoma metastasis and tested their association with NF1, BRAF and NRAS mutations. IsomiRs differentially expressed were identified and catalogued with reference to the canonical form. The reported non-random dysregulation of specific isomiRs contributes to the understanding of the complex melanoma pathogenesis and serves as the basis for further functional studies.


Introduction
Malignant cutaneous melanoma (CM) results from transformation of pigments-producing cells known as melanocytes. Worldwide incidence of melanoma has progressively increased over the last several decades, reaching an age-standardized rate of 3.4 per 100,000. However, in fair-skinned population regions, including Australia-New Zealand, Northern Europe and Northern America, this rate exceeds 15 per 100,000 (35.8, 17.8 and 16.1, respectively) (source IARC 2020).
Although malignant melanoma accounts only 3-5% of all skin cancers, it is the most lethal form, and determines up to 65% of the deaths [1]. The pathogenesis of melanoma is complex, and the onset of melanoma has been associated with several risk factors that include genetic, environmental, and phenotypic factors such as ultraviolet (UV) exposure, fair phototypes, multiple dysplastic nevi and a positive family/personal history of cutaneous melanoma [2].
The gain-of function mutations of BRAF and NRAS lead to constitutive activation of the signaling of the mitogen-activated protein kinase (MAPK) pathway and thus to tumor proliferation and progression [10]. On the other hand, the loss-of function of NF1 leads to hyperactivation of RAS proteins and to promotion of MAPK pathway signaling [11]. The discovery that BRAF and NRAS mutations are two of the main oncogenic drivers of melanoma proliferation and survival has resulted in important therapeutic implications and changed the management of cutaneous melanoma patients with the development of specific targeted therapies based on MAPK pathway inhibitors [12].
In the past decade, the characterization of the tissue miRNA profile has taken advantage of the spreading of deep sequencing techniques to investigate the sequence and abundance of small noncoding RNAs. Most small RNA sequencing studies focused on the identification of canonical miRNAs, whose sequence is reported in the miRBase database [18]. Recently, many deep sequencing experiments in human tissues evidenced the existence of alterations in miRNA length and sequence [19], suggesting the importance of including the analysis of miRNA variants or isoforms (isomiRs) in miRNA studies. The isomiR biogenesis has not been completely clarified yet, but it seems that isomiRs could derive from alternative Drosha and or Dicer cleavage of the precursor miRNAs [20], or from post-transcriptional modifications made by nucleotidyl transferase, which adds nucleotides to the pre-miRNA or mature miRNA ends [21][22][23]. Stability and abundance can vary between the canonical miRNA and its isoforms. The proposed role of isomiRs is to reinforce the regulatory activity of canonical miRNAs [24], though there is still much to discover. Furthermore, isomiR expression can discriminate human cancers [25].
Recently, we investigated the global canonical miRNA profile of formalin-fixed paraffinembedded (FFPE) samples of benign nevi, single and multiple primary melanomas [26]. In addition, we analyzed the smallRNA sequencing data at the isomiR level. Combining all the samples together, we identified a panel of highly expressed isomiRs, with an expression up to 10 times higher than the canonical miRNA.
This prompted us to further investigate the contribution of isomiRs in melanoma progression. Specifically, in this study we identify and classify isomiRs detectable in FFPE benign nevi and melanoma samples using proprietary small RNA sequencing data and assess isomiR expression in relation to the corresponding canonical miRNA. In addition, we compare the isomiR expression profile of benign nevus and melanoma. Finally, using publicly available smallRNA sequencing data from fresh-frozen skin cutaneous melanoma (SKCM) of The Cancer Genome Atlas (TCGA) database, we identify and compare the isomiR expression profile of primary melanoma and metastasis, and identify isomiRs whose expression is associated with NF1, BRAF and NRAS mutations.
In what follows, the term "canonical" has been used to indicate the sequence of a microRNA as reported in miRBase. On the other hand, "isomiRs", "miRNA variants" and "miRNA isoforms" indicates all the variant forms that differentiate from the canonical miRNA. Finally, the set of canonical miRNAs and their isomiRs has been grouped and called "mature microRNAs".

Mature microRNA Profile and miRNA Variants Characterization in Early-Stage Melanoma Samples
We analyzed the global mature microRNA (canonical and isomiR) profile of 23 FFPE samples: three benign nevi (BN) and 20 early-stage primary cutaneous melanoma (CM) samples [26]. Using a dedicated bioinformatics pipeline called isoMiRmap [27], we identified 494 mature microRNAs, including 130 canonical miRNAs and 364 microRNA isoforms that passed quality trimming and filtering. Specifically, we identified 90 canonical miRNAs with at least one isomiR, for a total of 324 isomiRs with a detectable canonical form. About 28% of these 90 canonical miRNAs had only one isomiR, 30% had two or three isomiRs and the remained 42% consisted in miRNAs that had a number of isomiRs that ranged from four to thirteen (Table 1, Figure 1A). The characterization of all mature miRNAs evidenced that 40 expressed canonical miRNAs did not have any isomiR in our data. As a counterpart, canonical miRNA sequences were not detected in our dataset for about 11% of the identified isomiRs. In this case, isomiRs are called "orphan" isomiRs (Supplementary Table S1). Since isomiR sequences differentiate from canonical miRNAs in different ways, we classified them into five categories: start-site isomiR, end-site isomiR, 3 non-templated addition isomiR, shifted isomiR and the remaining "mixed" isomiR. The first and second groups present nucleotides addition or deletion respectively at the 5 and at the 3 end of the canonical sequence. In these two groups, added nucleotides "match" with the sequence of the stem loop. On the contrary, when the added nucleotides in 5 or 3 end are different compared to those in the stem loop sequence, we defined 5 non-templated addition and 3 non-templated addition isomiRs, respectively. In our dataset, we did not identify any 5 non-templated addition isomiR. Conversely, several isomiRs present 3 non-templated addition and most of them consist in uridylation, namely addition of one or more uracils. Shifted isomiRs are long as the canonical counterpart, however their start-sites and, consequently, their end-sites are shifted to the left (in 5 direction) or to the right (in 3 direction) to one or more positions. Finally, the mixed isomiR group includes isomiRs that manifest at least two of these differences compared to the canonical form. An example of isomiR classification is reported in Figure 1B.

End-Site isomiRs Are the Most Abundant and Expressed isomiRs in FFPE Samples
Most of the detected isomiRs belonged to the end-site class of isomiRs (61.4%), while the other four classes were less represented. In addition, looking at the mixed group, we observed a prevalence of isomiRs that present the combination of different end-site modifications (deletions or additions of "matched" nucleotides) and 3 non-templated additions when compared to canonical miRNA (Table 1).
Looking at the normalized expression of the 364 isomiRs in the 20 FFPE early-stage melanoma samples, we observed that they can vary between 5-6 to over 2300. We selected the most expressed isomiRs in these 20 CM samples (normalized expression >20) and obtained a list of 126 isomiRs that are reported in Supplementary Table S2. We found that the majority of expressed isomiRs (about 75%) belong to the end-site group. Less frequent are the isomiRs belonging to the other groups: 9% mixed group, 9% 3 non-template addition group, 4% start-site group and 3% shifted group.

Identification of isomiRs More Expressed Than the Canonical Form in Early-Stage Melanoma
To identify the most relevant isomiRs in early-stage melanoma, we focused on isomiRs that showed a higher or an equal expression level than the corresponding canonical mRNA form. For this purpose, we calculated the expression ratio between isomiRs and canonical miRNA. For this analysis, we excluded miRNAs without isomiRs and orphan isomiRs. We considered the average expression of the selected 412 mature miRNAs in 20 early-stage melanomas and calculated the expression ratio for each isomiR. We obtained a wide set of ratios: some isomiRs are expressed up to 3000 times less than the canonical miRNAs, as expected, but some isomiRs are upregulated up to 50 times more than their canonical counterparts (Supplementary Table S3).

Classification of isomiRs with a Potentially Relevant Role in Early-Stage Melanoma
We decided to further investigate the expression of mature miRNAs in FFPE samples showing two characteristics: (i) at least one isoform is differentially expressed between early-stage CM and BN and (ii) at least one isoform is more expressed than the canonical form. We obtained a list of 10 miRNAs, which were further divided into four classes: (a) isomiRs with a similar trend in early-stage CM vs. BN and similar relative abundance; (b) isomiRs with a similar trend in early-stage CM vs. BN and different relative abundance; (c) isomiRs with opposite trend in early-stage CM vs. BN and similar relative abundance; and (d) isomiRs with opposite trend in early-stage CM vs. BN and different relative abundance.

IsomiRs with a Similar Trend in CM vs. BN and Different Relative Abundance
The second class is formed by miRNAs with the same trend of variation for all their isomiRs in CM vs. BN, but with different abundance distribution between CM and BN. miR-101-3p and miR-27b-3p belong to this class. Almost all variants and canonical miR-101-3p have a significantly lower expression in CM ( Figure 4A,B). Mature microRNAs of miR-101-3p are not distributed in the same way in CM and BN, though. As represented in Figure 4C, we can observe that the major expressed isomiR is miR-101-3p|0|0(+1U) in CM and miR-101-3p|−1|−1 in BN. All 9 isomiRs and canonical form of miR-27b-3p are downregulated in CM compared to BN with a significant difference in most of the cases ( Figure 4A,B). Focusing on abundance distribution of the different variants and canonical miRNA, we observed that miR-27b-3p|0|−1 is the most expressed form in CM, while the canonical miRNA is the most expressed mature miRNA in BN ( Figure 4B,C).

IsomiRs with Opposite Trend in CM vs. BN and Similar Relative Abundance
The only miRNA of the third class is miR-141-3p. This isomiR is downregulated in CM compared to BN (miR-141-3p|0|−1, p value = 0.0023), but the canonical miRNA and miR-141-3p|0|+1 have a comparable expression in CM and BN. miR-141-3p|0|+2 is only detectable in CM samples ( Figure 5A,B). The abundance distribution is similar between CM and BN, in fact miR-141-3p|0|−1 isomiR is the most expressed mature microRNA in both groups ( Figure 5C).

IsomiRs with Opposite Trend in CM vs. BN and Different Relative Abundance
The last class comprises three miRNAs, namely miR-30a-5p, miR-30d-5p and miR-203a-3p, which show different isomiRs abundance distribution in CM and BN and a different expression trend between isomiRs and the canonical form ( Figure 6A-C).

IsomiR Classification in Fresh-Frozen Primary Melanoma and Metastasis from TCGA Database
In order to identify isomiRs potentially associated with tumor progression, we interrogated the SKCM cohort of The Cancer Genome Atlas (TCGA) database. The TCGA database collects clinical and molecular data from fresh-frozen tissues, including smallRNA seq data. We explored the contribution of isomiRs in 323 melanoma patients, including 63 primary cutaneous melanomas (CM), and 260 melanoma metastases (MM), and assessed the association of isomiR expression with melanoma mutated genes. We detected a total of 3690 mature microRNAs, including 474 canonical microRNAs and 3216 isomiRs (Supplementary Table S4). Among the canonical miRNAs, only 85 miRNAs did not have any detected isomiR, while 389 canonical microRNAs had at least one co-expressed isomiR. Again, we observed the presence of orphan isomiRs (n = 258). However, for the majority of isomiRs, the canonical forms are detectable (n = 2958). The classification in isomiR types is reported in Table 4. The most frequent isomiR group is the end-site type. Interestingly, in this isomiRs list we showed the presence of a 5 non-templated addition group, namely isomiRs with the addition of one or more nucleotides in the 5 end, which did not match with the stem-loop sequence. In addition, we found four novel miRNAs and three isomiRs of novel miRNAs. These novel miRNAs have been previously detected by examining small RNA-seq samples [28] (Supplementary Table S4).

Identification of the isomiRs More Expressed Than the Canonical Form in Fresh Primary Melanoma Samples from TCGA
To identify the most abundant isomiRs in fresh-frozen melanomas we investigated isomiRs that showed a higher or an equal expression level than the corresponding canonical mRNA form in primary melanomas from the TCGA SKCM cohort, which includes mostly stage II and stage III fresh primary melanoma.
We considered the average expression of mature microRNAs in 63 CM and we calculated the expression ratio between 2958 isomiRs and the corresponding canonical miRNA. We obtained a wide set of ratios that are reported in Supplementary Table S5, and 1373 of them showed a ratio higher than 1. Among these 1373 isomiRs, there were also 21 of the isomiRs identified in early-stage FFPE melanomas (mostly stage I melanomas) reported in Table 2, confirming the identification of these isomiRs also in fresh melanoma tissue and late-stage melanomas (Table 5).

Identification of isomiRs Associated with Driver GeneMutations in TCGA Samples
We obtained the mutation status of BRAF, NRAS and NF1 for the melanoma samples from the TCGA database. We reported the mutation status of the three genes (when this information was available), in Supplementary Table S7. The most frequent mutation in BRAF is the missense mutation V600E, which is present in 145 samples. Among these 145 samples, 21 samples have also the mutation V600M. NRAS mutations are mostly missense mutation in Q61-position. Specifically, we found 33 samples with Q61R, 28 samples with Q61K, 11 samples with Q61L and five samples with Q61H.
NF1 presents a wide variety of different mutations and different types of mutations, including missense, stop codon, frameshift, and splice mutations. We analyzed the isomiR differential expression in mutated NF1, BRAF and NRAS compared to wild type samples.

Name
Sequence Type

Discussion
The miRBase database provides a single mature sequence for each miRNA, which is usually the sequence with the highest coverage reported in small RNA sequencing experiments. However, we discovered that in both FFPE and fresh-frozen melanoma samples several isomiRs are more abundant than the canonical forms, being up to 50 times more expressed in FFPE early-stage melanoma samples. We observed that end-site isomiRs are the most represented group in FFPE and fresh melanomas. Similarly, in colon cancer it was observed that 3 isomiRs are the most frequent microRNA variants [29].
In the literature, there is a concern that sequencing and/or mapping artefacts may result in an overrepresentation of isomiRs [30]. However, previous studies have shown that isomiRs represent actual molecules rather than sequencing artifacts. First, by using aggressive filtering, most categories of isomiRs are detected at levels above thresholds comparable to their canonical miRNA sequence. In many cases, the isoform is present at increased levels as compared to the canonical sequences [19]. Furthermore, various studies have shown that isomiR profiles across tissues differ. In other words, the most abundant isoform does differ across tissues and disease states. This suggests that isomiR biogenesis may be a regulated process [31]. Furthermore, sequence library preparation could cause miRNA degradation at the miRNA ends, resulting in the identification of isomiRs. If true, these effects cannot explain the bias of the heterogeneity between the two miRNA ends in which the 3 -ends are usually more diverse than the 5 -ends [32]. Additionally, a significant number of isomiRs are also longer and more abundant than the canonical miRNA. These observations suggest that isomiRs are not experimental artifacts but rather true biological events. Our sequence analysis pipeline uses isoMiRmap [27]. In a deterministic and exhaustive manner, isoMiRmap comprehensively reports all isomiRs whose sequences exist within known miRNA sequences. Additionally, this approach identifies 3 -non-templated isomiRs.
Similarly to canonical miRNAs, isomiR expression can discriminate human cancers, [24]. In addition, recent studies have demonstrated that isomiRs could be used as prognostic and diagnostic biomarkers in various cancers and related subtypes; in fact, they can differentiate between healthy and non-healthy individuals [33]. To analyze the contribution of isomiRs to melanoma pathogenesis, we compared the isomiR profile of benign nevi and early-stage primary cutaneous melanoma FFPE samples, and of primary cutaneous melanoma and melanoma metastasis of fresh-frozen samples.
A differential expression of isomiRs was observed between benign nevi and earlystage melanoma FFPE samples; specifically, 55 mature microRNAs (including the canonical form and all isoforms) were differently expressed. Interestingly, for some miRNAs only the isomiR(s) was differently expressed and not the canonical miRNA. Among these 55 mature microRNAs, we identified 16 isomiRs that were more abundant than their canonical counterpart. We focused on these 16 isomiRs, which were classified in four groups according to the expression trend in tumors and miRNA/isomiR relative abundance distribution.
Among the mature microRNAs with similar expression trends and similar abundance distribution in nevi and early-stage melanomas, we found miR-19b-3p, miR-27a-3p, miR-29a-3p and miR-222-3p. Despite our data that showed that miR-19b-3p and its isomiR are downregulated in early-stage FFPE melanoma, an in vitro study showed that miR-19b-3p expression is higher in most melanoma cell lines compared with normal melanocytes, where it promotes cell proliferation [34]. The literature suggests that melanoma tissues upregulate miR-27a-3p and its expression is positively correlated with tumor stage and lymph node metastasis of melanoma tissues. The downregulation of miR-27a-3p induces autophagy and apoptosis of melanoma cells [35]. No data exist about the differential expression of this miRNA between nevi and melanoma. We found that miR-29a-3p expression is lower in early-stage melanoma. This data is in agreement with the literature; in fact it was observed that the overexpression of miR-29a-3p inhibits proliferation, migration, and invasion, and promotes apoptosis by blocking Wnt/β-catenin and NF-κB pathways [36]. The higher expression of miR-222-3p and its isomiRs in early-stage melanoma suggests an oncogenic role, in fact miR-222-3p belongs to the miR-221 family, which is an important oncogenic microRNA family involved in the progression of melanoma and have a key role in EMT [37]. For this class of isomiRs, we hypothesized a cooperative effect of the isomiRs on the target genes.
In the second class we identified miR-101-3p and miR-27b-3p mature microRNAs, which have the same expression trend but different abundance distribution in benign nevi and early-stage melanomas. All their mature microRNAs are downregulated in melanoma.
These results agree with previous studies describing these two canonical miRNAs as tumor suppressor miRNAs. miR-101-3p directly targets MITF and EZH2, resulting in inhibition of invasion and proliferation. In addition, low expression of miR-101-3p correlated with a poor survival in stage IV melanoma patients [38]. Canonical miR-27b-3p inhibits the development of melanoma by targeting MYC [39].
We identified a specific isomiR of miR-141-3p, miR-141-3p|0|−1, which is downregulated in early-stage melanoma, while the canonical miRNA is not differentially expressed. In both groups, namely benign nevi and early-stage melanoma, this specific isomiR is the most abundant variant. This seems to suggest that only this isomiR could have a protective role, and not the canonical form or the other variants. The canonical miRNA was found to inhibit the proliferation of melanoma cells [40].
Our analysis revealed that miR-30a-5p, miR-203a-3p and miR-30d-5p showed a different expression trend for at least two isoforms, and a different mature microRNA prevalence in benign nevi and early-stage melanomas. The canonical miRNAs of these isomiRs have been extensively studied in melanoma. MiR-30a-5p is considered a tumor suppressor mi-croRNA in melanoma. In fact, miR-30a-5p is normally downregulated in melanoma tissues and cell lines. The overexpression of miR-30a-5p leads to a significant inhibition of proliferation, migration, and invasion of melanoma cells. Moreover, miR-30a-5p in vivo delays tumor growth [41] and inhibits metastasis [42]. In our analysis, the canonical miR-30a-5p is more expressed in benign nevi than early-stage melanomas, although the difference is not statistically significant. We can speculate that the tumor suppressor effects derive from the cumulative activity of all the isoforms.
MiR-203a-5p is an important tumor suppressor miRNA in melanoma, involved in the regulation of proliferation, cell cycle, migration, and invasion [43,44]. The expression of canonical miR-203a-3p is downregulated in melanoma, especially in metastatic melanoma. Low miR-203a-3p expression is associated with poor overall survival [45]. MiR-30d-5p acts as an oncomiR in melanoma. High expression of miR-30d-5p in melanoma is associated with stage, metastatic potential, shorter time to recurrence, and reduced overall survival. miR-30d-5p induces metastatic behavior of melanoma cells [46]. In this analysis, BN and melanomas show similar expression levels of miR-203a-3p and miR-30d-5p canonical miRNAs, which are probably linked with the early stages of the melanomas, but several isoforms of miR-30 are significantly downregulated in tumors.
Mutations in key genes of the MAPK/ERK pathway, namely NF1, NRAS and BRAF, identify three main molecular subtypes of melanoma. Several microRNAs are known to affect or to be affected by the altered activation of this pathway in melanoma, and some of them are responsible for drug resistance [47]. Since the study of the expression of isomiRs in relation to the mutation status of these proteins is missing in the literature, we focused on the different isomiR profiles associated with the different mutation status of NF1, BRAF and NRAS.
We collected information on the mutation status of BRAF, NRAS and NF1 genes for SKCM TCGA samples and investigated the association between isomiR expression and the most clinically relevant mutations in melanoma.
Several canonical miRNAs of isomiRs that are differently expressed in mutated BRAF samples were already studied in relation to BRAF. It was observed that Vemurafenib, a BRAF inhibitor, alters the miRNA expression in melanoma cells, leading to an upregulation of miR-181a-2-3p [49]. MiR-181c-5p was found to be associated with the BRAF mutation in thyroid follicular adenomas. Indeed, its expression is upregulated in RAS-or BRAFmutated follicular adenomas compared to wild-type tumors [50]. However, its role in melanoma has yet to be established.
It has been described that B-Raf/MKK/ERK can upregulate miRNAs in melanoma cells, specifically the miR-17-92 microRNA family [57]. We observed a remarkable upregulation of several isomiRs from the miR-17-92 cluster in NRAS mutant melanomas.
IsomiR expression validation with methods other than NGS is still an open issue, given the lack of commercial kits or largely recognized approaches. The most promising technique to quantify miRNA variants is Dumbbell PCR. Honda et al. developed this method to quantify specific small RNA variants with a single nucleotide resolution at terminal sequences (5 and 3 ) [58]. In our previous work, we combined two different commercial kits to quantify the sum of 5 and 3 isoforms and the canonical miRNA [26]. The development of reliable assays will be of the utmost importance for further functional studies involving isomiRs.

FFPE Samples and Small-RNA Sequencing Data
Patients characteristics and experimental procedure are described in our previous work [26]. In this paper, we used small-RNA sequencing data from three benign nevi (BN) and 20 melanomas (CM).
Samples were obtained from 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). Histopathologic specimens were evaluated by two dermato-pathologists and classified into benign nevi and melanomas. BN includes patients with no prior diagnosis of CM or non-melanoma skin cancer and a follow-up of at least 10 years. CM includes 19 patients with superficial cutaneous melanoma (stage I, Breslow range = 0.3-.), and one patient with nodular cutaneous melanoma (stage II, Breslow = 2.5). Patient specimens were formalin-fixed and paraffin-embedded (FFPE), which were used to RNA extraction using the miRNeasy FFPE kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. The 23 small-RNA libraries were generated using TruSeq Small RNA Library PrepKit v2 (Illumina, San Diego, CA, USA, RS-200-0012/24/36/48), according to the manufacturer's directions. Raw base-call data were demultiplexed using an Illumina BaseSpace Sequence Hub and converted to FASTQ format [26]. The small RNA sequencing raw data (FASTQ format) are available at the European Nucleotide Archive (ENA) with the following accession number: PRJEB35819.

IsomiR Analysis of TCGA SKCM Cohort
We collected isomiR data and mutation status of BRAF, NRAS and NF1 from 323 melanoma fresh-frozen samples from The Cancer Genome Atlas (TCGA), including 63 primary cutaneous melanomas (CM) and 260 melanoma metastases. In the TCGA database, "metastatic melanoma" refers to the metastatic tissue and not to the primary melanoma that have metastasized. To avoid misunderstanding, we preferred to call this group melanoma metastasis (MM), instead of metastatic melanoma. CMs are mostly II or III stage (

IsomiR Identification
The mapping of sequence reads and isomiR quantification was performed using the isoMiRmap tools [27]. isoMiRmap identifies and quantifies all isomiRs through the direct processing of short RNA-sequencing datasets. Briefly, for each short RNA-sequence dataset, the sequence reads were quality trimmed using the Cutadapt tool [59], keeping only those reads between 18 and 26 nts in length. Each resulting read is subsequently compared to a "lookup" table containing all possible isomiRs from miRbase release 22 hairpin sequences. All sequences are required to have exact matches to the isomiR lookup table, after which they are then counted and quantified in read per million (RPM). This approach identifies 5 -isomiRs, 3 -isomiRs, 5 -and-3 mixed isomiRs, as well as 3 nontemplated post-transcriptional additions. This approach was used for both those datasets generated here and those samples represented in TCGA.

IsomiR Labeling and Classification
The annotation system developed by Loher et al. was used to label isomiRs [32]. Briefly, nomenclature indicates the name of the canonical miRNA, 5 end (start site), 3 end (end-site) and the eventual insertion of uracil of the isomiR compared to the canonical miRNA sequence in miRBase database. Specifically, to define the start site and end site, zero indicates the same terminus of the canonical miRNA, while positive sign (+) or negative sign (-) followed by the number indicates the isomiR nucleotide shift in the 3 or 5 direction, respectively, when compared with miRBase microRNA sequence. IsomiR were classified in six different groups depending on the change from the canonical miRNA: start site isomiRs (those isomiRs that varies in the start site with deletion or addition of nucleotides that match with stem loop), end site isomiRs (those isomiRs that varies in the end site, with deletion or addition of nucleotides that match with stem loop), 5 non-templated addition isomiRs (isomiRs with addition of nucleotides at the 5 end, which do not match with stem loop), 3 non-templated addition isomiRs (isomiRs with addition of nucleotides at the 3 that do not match with stem loop sequence), shifted isomiRs (start site and end site are both shifted to 5 or to 3 direction, but the length of the isomiR is identical to the canonical miRNA), and mixed isomiRs (isomiRs with more than one variation).

Statistical Analysis
Normalized data were analyzed with GeneSpring GX software (Agilent Technologies). Differentially expressed mature miRNAs (canonical miRNAs + isomiRs) were identified using a fold change > 2 filter and moderated t-test (FDR 5% with Benjamini-Hochberg correction) in benign nevi vs. melanoma comparison (small-RNA seq data). For TCGA SKCM data, we applied a fold-change > 1.5 and t-test unpaired (FDR 5% with Benjamini-Hochberg correction) in primary melanoma vs. melanoma metastasis comparison. To analyze the association between isomiR and the melanoma-linked mutation, we performed an unpaired t-test with unequal variance (Welch correction) (FDR 5% with Benjamini-Hochberg correction). Graphpad Prism 6 (GraphPad Software, San Diego, CA, USA) was used to perform specific isomiR comparisons using the Mann-Whitney non-parametric test.

Conclusions
Tumor characterization at the small RNA level, when focusing only on canonical miRNAs, has the risk to overlook and lose important information that could improve our understanding of cancer pathogenesis. Here, we provided a comprehensive characterization of isomiR dysregulation in benign nevi, melanoma at different stages and sources, and melanoma metastasis. We showed that isomiRs can be more expressed than their canonical miRNA and also that different isoforms are modulated in different and sometimes opposite ways in normal and tumor samples, a fact that was not expected. We reported that specific isomiRs are associated with tumor mutation status, sometimes diverging from the trend of the canonical miRNA. Our results and observations about a non-random dysregulation of specific isomiRs in melanoma serve as the basis for further functional studies.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/ncrna7040063/s1. Table S1: List of mature microRNAs detected in the small-RNA sequencing of 23 samples. Table S2: List of isomiRs that are most expressed in 20 early-stage melanoma samples. Table S3: List of isomiR/canonical miRNA ratios in 20 early-stage melanoma samples. Table S4: List of mature microRNAs identified in TCGA SKCM data. Table S5: List of isomiR/canonical miRNA ratios in 63 fresh-frozen melanoma samples from TCGA. Table S6: List of mature microRNAs differentially expressed in melanoma metastasis compared to primary melanomas from TCGA SKCM dataset. Table S7: List of TCGA melanoma samples with the mutation status of BRAF, NRAS and NF1. Figure S1: List of 332 mature microRNAs differently expressed in melanoma metastasis compared to primary melanoma. Figure S2: Opposite expression trend of mature microRNAs in melanoma metastasis and primary melanoma.  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The small RNA sequencing raw data (FASTQ format) are available at the European Nucleotide Archive (ENA) with the following accession number: PRJEB35819.

Conflicts of Interest:
The authors have no conflict of interest to declare. All co-authors agree with the contents of the manuscript and have no financial interests to report.