The Circulating Transcriptome as a Source of Biomarkers for Melanoma

The circulating transcriptome is a valuable source of cancer biomarkers, which, with the exception of microRNAs (miRNAs), remains relatively unexplored. To elucidate which RNAs are present in plasma from melanoma patients and which could be used to distinguish cancer patients from healthy individuals, we used next generation sequencing (NGS), and validation was carried out by qPCR and/or ddPCR. We identified 442 different microRNAs in samples, eleven of which were differentially expressed (p < 0.05). Levels of miR-134-5p and miR-320a-3p were significantly down-regulated (p < 0.001) in melanoma samples (n = 96) compared to healthy controls (n = 28). Differentially expressed protein-encoding mRNA 5′-fragments were enriched for the angiopoietin, p21-activated kinase (PAK), and EIF2 pathways. Levels of ATM1, AMFR, SOS1, and CD109 gene fragments were up-regulated (p < 0.001) in melanoma samples (n = 144) compared to healthy controls (n = 41) (AUC = 0.825). Over 40% of mapped reads were YRNAs, a class of non-coding RNAs that to date has been little explored. Expression levels of RNY3P1, RNY4P1, and RNY4P25 were significantly higher in patients with stage 0 disease than either healthy controls or more advanced stage disease (p < 0.001). In conclusion, we have identified a number of novel RNA biomarkers, which, most importantly, we validated in multi-center retrospective and prospective cohorts, suggesting potential diagnostic use of these RNA species.


Introduction
Although malignant melanoma accounts for~4% of skin cancer cases, it accounts for~75% of all associated mortalities. In the US alone, it is estimated there were 87,000 new cases and nearly 10,000 deaths due to melanoma in 2017 [1]. Furthermore, the incidence of melanoma has been increasing at faster rate than any other cancer type, having doubled since 1973 [2]. The clinical case for early diagnosis of melanoma is compelling, as if detected early enough, (stage I) 5 year survival is >95%, whereas in advanced melanoma (stage IV) survival is just 10-20% [3].
Non-invasive diagnostics, or liquid biopsies, represent a major advance towards earlier diagnosis and disease monitoring of cancer patients, including those with melanoma. As a consequence, there has been a great deal of interest in recent years in the potential of circulating nucleic acids, and in particular circulating microRNAs (miRNAs) [4,5]. In melanoma, several studies have implemented circulating miRNA in metastasis and risk of recurrence [6,7]. However, outside of miRNAs there has been little research on other cell-free (cf) RNA species in the circulating transcriptome. Part of the reason for this paucity of knowledge is the presence of high levels of RNase activity in blood, which typically results in fragmentation of longer RNA species such as mRNA [8,9]. This makes detection of these molecules particularly challenging. The advent of next-generation sequencing (NGS) technology, however, presents new opportunities for the field, as sequences can be elucidated in a 'bottom-up' manner without the need of a priori probe sequence knowledge. We used next generation sequencing (NGS) to characterize and compare the circulating transcriptomes of plasma from melanoma patients with different stage disease, along with sex-/age-matched healthy individuals, in order to identify novel biomarkers for this cancer.
We identified not only circulating miRNAs with biomarker potential, but also mRNA gene fragments and non-coding YRNAs. YRNAs (Ro-associated Y) are poorly characterized, non-coding RNAs, which were initially identified in the blood of rheumatic autoimmune disease patients [10]. They are a family based around four highly conserved sequences (RNY1, RNY3, RNY3, RNY5) involved in Ro60 inhibition, DNA replication, and quality control of non-coding RNAs [11][12][13]. Our possible biomarkers were validated in an independent cohort of 327 plasma samples from melanoma patients, collected retrospectively and prospectively. This study provides evidence that the largely unexplored circulating transcriptome could provide a valuable source of liquid biopsy biomarkers for melanoma in particular, and cancer in general.

Sequencing the Circulating Transcriptome of Melanoma Patients
We optimized library construction by comparing several protocols, including different ribodepletion methodologies, using plasma from healthy controls in a pilot study, before settling on the protocol described in materials and methods [14]. It should be noted, however, that this protocol might result in some ligation-based bias [15]. Using this protocol, we sequenced cfRNA derived from plasma of melanoma patients, along with age-/sex-matched controls (Table 1).
Due to the low quantity of cfRNA in plasma, we decided to pool plasma samples according to clinical stage, as depicted in Table 1. For each pooled sample, 40-50 million reads were obtained, with an average Phred score of 37.9 (Supplementary Table S1). Between 50-55% of reads were mapped to the human genome (hg19), and approximately half of these sequences were considered small RNA (sRNA), representing sequences between 18-43 nt (Supplementary Table S1; Figure 1a). By far the largest category by frequency of reads was YRNAs, accounting for an average of 40.6% of reads (range 34-48%; Figure 1). Indeed, this category accounted for >95% of reads that composed the major peak at 32 nt seen in the size profiles of samples (typical example shown in Figure 1a).  Figure 1a). By far the largest category by frequency of reads was YRNAs, accounting for an average of 40.6% of reads (range 34-48%; Figure 1). Indeed, this category accounted for >95% of reads that composed the major peak at 32 nt seen in the size profiles of samples (typical example shown in Figure 1a).

mRNA Fragment Expression
Because mRNA is degraded by RNase activity in the blood into fragments with an average size of 200 bases [17], rather than mapping reads to whole genes or exons, we mapped them to short annotated probe sets (<500 bases). There were 3672 different probe sets (median length of 244 nt), with greater than 25 reads per probe set mapping to 13,641 different transcripts. For each probe set we analyzed their relative position within the respective transcripts as a percentage of the entire transcript length. By far the largest proportion (5350/13641 (39%)) of probe sets mapped to the first decile (10%) (i.e., 5′-end) of the transcript, compared to an average of 6% for the other nine deciles of the transcript length ( Figure 3a). (a,b) Levels of miR-320a and miR-134 were measured in a cohort of 124 plasma samples (28 controls and 96 melanoma patients). Expression levels were compared using the Kruskal-Wallis multiple comparison test (both miR-134 and miR-320a had p-values of <0.0001) and the Mann-Whitney independent t-test to carry out a pairwise comparison between groups (* p < 0.5; *** p < 0.001) (c,d) ROC analysis of miRNA probe expression levels as diagnostic biomarker (i.e., control vs. all melanoma patients, regardless of stage).

mRNA Fragment Expression
Because mRNA is degraded by RNase activity in the blood into fragments with an average size of 200 bases [8], rather than mapping reads to whole genes or exons, we mapped them to short annotated probe sets (<500 bases). There were 3672 different probe sets (median length of 244 nt), with greater than 25 reads per probe set mapping to 13,641 different transcripts. For each probe set we analyzed their relative position within the respective transcripts as a percentage of the entire transcript length. By far the largest proportion (5350/13641 (39%)) of probe sets mapped to the first decile (10%) (i.e., 5 -end) of the transcript, compared to an average of 6% for the other nine deciles of the transcript length ( Figure 3a).  Expression levels were compared using the Kruskal-Wallis multiple comparison test (all the mRNA fragments had p-values of <0.0001), and the Mann-Whitney independent t-test to carry out a pairwise comparison between groups (* p < 0.5; *** p < 0.001) (f-i) ROC analysis of mRNA probe expression levels as diagnostic biomarker (i.e., control vs. melanoma patient). (j) Panel (ATM, AMFR, and SOS1) performance shown with black line, and SOS1 is shown as a gray line for comparison.
Seventy of the probe sets were identified as being differentially expressed between samples (p < 0.05; Table A1). We designed custom Taqman probes to detect fragments from the five most differentially expressed probe sets (i.e., corresponding to ATM, ARHGAP, AMFR, CD109, and SOS1 genes) (Supplementary Table S2; Figure S1). We were unable to design a specific probe for the ARHGAP probe set due to a high level of repetitive sequences. None of the mRNAs validated are predicted targets of miRNAs studied (i.e., miR-134-5p and miR-320a) using the predictive algorithms TargetScan and miRDB [17]. The four probe sets were measured by qRT-PCR in 47 control and 173 melanoma patient samples ( Table 1).
Levels of ATM, AMFR, CD109, and SOS1 were all significantly higher (p < 0.001) in plasma samples from either stage 0, stage I/II, or stage III/IV melanoma patients than in healthy controls (Figure 3b-e, respectively), consistent with the NGS data (Table A1). Surprisingly, levels of AMFR and CD109 were higher in plasma from stage 0 patients than samples with more advanced disease. We carried out ROC analysis to determine the diagnostic ability of the mRNA fragments (Table 3; Figure 3f-i), and to discriminate between different disease stages ( Table 3; Supplementary Figure S2a-d). We looked at combinations of these biomarkers using the PanelomiX ROC comparison algorithm [18]. A combination of ATM, SOS1, and AMFR with cut-off values of 2.13, 2.96, and 2.26 respectively, gave the best diagnostic accuracy (AUC = 0.825) (Table 3, Figure 3j).

YRNA Expression
As nearly half of all mapped sRNA reads were identified as YRNA sequences (Figure 1), we examined this class of non-coding RNAs further. There were 322 different YRNA and YRNA-associated sequences identified in our samples, consisting of three YRNA sequences (RNY1, RNY3, and RNY4) representing an average of 26.1% of reads, 30 YRNA pseudogenes representing an average of 48.4% of reads, 69 7SK sequences (average 0.05% of reads), and 194 Rfam predicted YRNA sequences representing an average of 25.5% of reads (Supplementary Table S3). The vast majority of reads were represented by RNY4 and RNY4P sequences, accounting for >98% of their respective YRNA class (Figure 4a).
We compared the expression of YRNAs between low-stage disease (i.e., stage 0 and I/II) and high-stage disease (i.e., stage III and IV), and identified five differentially expressed YRNA fragments (p < 0.05; Table 4). We designed custom Taqman probes to measure three of these YRNAs (RNY3P1, RNY4P1, and RNY4P25), selected on the basis of fold-change and read count. We measured levels of the YRNAs in a validation cohort of 80 samples (22 controls and 58 melanoma patients, Table 1). Levels of all three YRNAs were significantly higher in stage 0 samples than control samples or stage I/II samples (Figure 4b-d).  Values are shown as absolute copies per µ l. Expression levels were compared using the Kruskal-Wallis multiple comparison test (all the YRNA fragments had p-values of < 0.001), and the Mann-Whitney independent t-test to carry out a pairwise comparison between groups, (* p < 0.5; **p < 0.01*** p < 0.001).

Discussion
The presence and relative stability of cfRNA in biological fluids has led to a great deal of interest in their use as 'liquid biopsies' for disease, in particular for cancer. However, with the exception of miRNAs, the circulating transcriptome remains largely unexplored. While NGS offers researchers the ability to elucidate the circulating transcriptome in its entirety, and therefore to identify novel biomarkers of disease, the application of RNAseq to biofluids such as plasma poses many challenges, not least the low quantity and quality of RNA present in these samples. As a consequence, studies to date have focused on the technical optimization of these techniques [20][21][22]. However, very few studies to date have sought to assess the potential usefulness of their findings through validation in independent cohorts.
In order to fully explore the complexity and biomarker potential of the melanoma circulating transcriptome, we pooled samples to maximize the starting quantity of cfRNA. As a result, we were able to obtain 40-50 million reads per pooled sample, an order of magnitude higher than comparable studies [21,23]. In contrast to exosomal cfRNA [23], we found that miRNAs only represented a minor component (<3%) of the whole plasma circulating transcriptome, levels similar to other plasma NGS Values are shown as absolute copies per µl. Expression levels were compared using the Kruskal-Wallis multiple comparison test (all the YRNA fragments had p-values of < 0.001), and the Mann-Whitney independent t-test to carry out a pairwise comparison between groups, (* p < 0.5; **p < 0.01*** p < 0.001).

Discussion
The presence and relative stability of cfRNA in biological fluids has led to a great deal of interest in their use as 'liquid biopsies' for disease, in particular for cancer. However, with the exception of miRNAs, the circulating transcriptome remains largely unexplored. While NGS offers researchers the ability to elucidate the circulating transcriptome in its entirety, and therefore to identify novel biomarkers of disease, the application of RNAseq to biofluids such as plasma poses many challenges, not least the low quantity and quality of RNA present in these samples. As a consequence, studies to date have focused on the technical optimization of these techniques [19][20][21]. However, very few studies to date have sought to assess the potential usefulness of their findings through validation in independent cohorts.
In order to fully explore the complexity and biomarker potential of the melanoma circulating transcriptome, we pooled samples to maximize the starting quantity of cfRNA. As a result, we were able to obtain 40-50 million reads per pooled sample, an order of magnitude higher than comparable studies [20,22]. In contrast to exosomal cfRNA [22], we found that miRNAs only represented a minor component (<3%) of the whole plasma circulating transcriptome, levels similar to other plasma NGS studies [19][20][21]. We identified 442 different miRNAs in our samples, somewhat higher than that reported in comparable studies [23,24], probably as a result of the pooled design and the higher quantity of RNA that we used. Consistent with other studies, we found that let-7b, miR-423, and miR-320a-3p were the most highly expressed miRNAs in our plasma samples [20]. We identified eleven miRNAs that were differentially expressed between healthy controls and the different clinical stages of melanoma (Table 2). This included miR-21, which has previously been shown to be upregulated in melanoma plasma samples [25], and miR-92b and miR-628, both of which are more highly expressed in plasma from monosomy 3 uveal melanoma patients [26].
Based on our sequencing results, we measured the expression of miR-320a-3p and miR-134-5p, the two miRNAs most differentially expressed between samples, in a validation cohort of 96 melanoma patients and 28 controls. Both miRNAs were significantly down-regulated (p < 0.0001) in plasma from all stages of melanoma patients when compared to samples from healthy controls. MiR-320a has also been found to be down-regulated in melanoma tumor cells when compared to heathy skin samples [27]. Furthermore, this miRNA was shown to function as an inhibitor of cell proliferation. The down-regulation of miR-320a has been observed in the blood of several cancers including colorectal cancer [28], gastric cancer [29], and retinoblastoma [30]. Moreover, miR-320a is up-regulated in melanoma cells after treatment with bevacizumab or rapamycin + bevacizumab [31]. Furthermore, miR-134 has been characterized as a tumor suppressor, acting to regulate proliferation, apoptosis, and invasion and migration in a wide range of cancer types, including melanoma [32,33]. ROC analysis of miRNA expression gave AUC values of 0.798 and 0.788 respectively. The miR-320a miRNA had a sensitivity of 90%, whereas miR-134 had a specificity of 96%, suggesting these two miRNAs in combination could be useful biomarkers for melanoma.
Even though circulating extracellular mRNA was first detected in 1999 (in melanoma) [34], as the vast majority of circulating mRNA is degraded by blood RNase activity [35], this potential source of biomarkers has largely been overlooked, even though mRNA fragments can represent up to 75% of total cfRNA [19]. In our study, just over 5% of mapped reads corresponded to protein-encoding mRNA fragments. We detected 3672 probes (<500 bases in length) that had at least 25 mapped reads in our samples. Nearly 40% of the probes mRNA fragments mapped to the 5´-terminus (i.e., first 10%) of their respective gene transcripts, probably reflecting the 3´to 5´cleavage activity of RNase A, the major RNase species in blood [26]. We did not notice a corresponding shift in the length profile between healthy and melanoma patient samples [36].
Pathway analysis of the genes corresponding to differentially-expressed mRNA fragments showed significant enrichment in the angiopoietin, p21-activated kinase (PAK) and Eukaryotic Initiation Factor 2 (EIF2) pathways. It has been previously reported that circulating level of Angiopoietin-2 (Ang-2) protein in melanoma patient sera closely correlates with disease progression [37]. Similarly, amplification of the PAK (p21-activated kinase) pathway is characteristic of BRAF-wild type melanoma [38], while in BRAF-mutant melanoma it is responsible for resistance to MAPK-inhibitor treatment [39]. Interestingly, both SOS1 and ATM1, which were the third and fourth most differentially expressed probe sets in our analysis, form part of the angiopoietin, PAK, and EIF2 pathways. We measured levels of SOS1, ATM1, CD109, and AMFR mRNA fragments in plasma samples from 173 melanoma patients and 47 healthy controls. With the exception of CD109, these mRNA fragments mapped to regions corresponding to the 5 -terminus of the reference gene transcript and included the initiation codon. Consistent with the NGS data (Table A1), levels of all these mRNA fragments were up-regulated in melanoma patient samples compared to samples from healthy controls. Particularly intriguing was the up-regulation of CD109 and AMFR in stage 0 samples compared to samples from more advanced stage melanoma, suggesting that these mRNA fragments could be used for early diagnosis of melanoma, although we do not have data on how many of these patients went on to develop advanced disease. (Figure 3b-e). CD109 has been identified as an important regulator of the Epithelial-mesenchymal transition (EMT pathway), and has also been found to be down-regulated in more advanced stage hepatocellular carcinoma [40]. The product of the AMFR gene, gp78, also regulates EMT, and increasing levels of AMFR are associated with metastatic melanoma [41,42]. Intriguingly, CD109 is a predicted target of miR-134 by the Targetscan algorithm; we are currently carrying out experiments to confirm this. ATM1 is a serine/threonine kinase induced by DNA damage and associated with risk in many cancer types [43]. SOS1 is a guanine nucleotide exchange factor for RAS proteins frequently mutated in melanoma [44]. Interestingly, all four of these gene fragments were more highly expressed not only in advanced stage disease, but also stage 0 disease; indeed, levels of CD109 and AMFR were higher in plasma from stage 0 disease than more advanced stage disease, suggesting that these biomarkers maybe non-tumoral in origin. Consistent with this hypothesis, the release of CD109 by bone marrow mesenchymal stem cells has recently been shown to attenuate EMT in skin squamous cell carcinoma [45], and AMFR plays an important role in regulation of the anti-cancer immune Stimulator of interferon genes (STING) pathway [46].
To test the potential diagnostic ability of these biomarkers, we carried out ROC analysis, however the results from individual mRNA fragments were poor (AUC range 0.722 (SOS1) to 0.767 (ATM1). In contrast, a combination of ATM1, SOS1, and AMFR resulted in an AUC value of 0.825 with a sensitivity of 75% and specificity of 92%. Although these findings need to be confirmed independently, this combination compares very favorably with existing sera markers such as LDH and S100B with reported sensitivities/specificities of 41.6/84.2% and 36.3/96.5%, respectively [47].
By far the largest class of circulating cfRNA that we identified in the samples corresponded to YRNA sequences, accounting for close to 50% of sRNA mapped reads. Remarkably, despite the prevalence of cfYRNAs in the blood, there is very little known about this class of ncRNA. YRNAs are short 80-110 nt ncRNAs, first identified in the early 1980s as an RNA component of the soluble Ro60 ribonucleoprotein particle found in the blood of patients with autoimmune diseases [48]. The function of YRNAs is still poorly understood; they appear to be essential for DNA replication [49] and are up-regulated in cancers [50], presumably as a result of their association with apoptosis [51]. The first description of circulating cfYRNAs came in 2013 from Dhahbi et al., who observed that 33% of mapped reads from sera of healthy individuals were YRNA sequences [52]. The same group later reported that YRNA accounts for 38% of cfRNA in sera from breast cancer patients [53], and subsequently, in the sera of head and neck squamous cell carcinoma patients [54]. More recently, a study of 183 plasma samples from healthy individuals found that YRNAs accounted for 63% of cfRNA [55]. As far as we are aware, apart from a recent study that measured YRNAs in the sera of 30 renal carcinoma patients [56], this the first study to look at the biomarker potential of YRNAs in cancer patients. Interestingly, we found that levels of YRNAs were significantly higher in samples from patients with stage 0 disease, maybe pointing to increased levels of tumor-associated apoptosis [51] even despite the small tumor sizes compared to more advanced disease stages.
In summary, we have elucidated the circulating transcriptome of plasma samples from melanoma patients and found a number of novel RNA biomarker species that we validated independently using qRT-PCR and ddRT-PCR in combined retrospective and prospective cohorts, detected from only 1mL of serum. These findings have potential clinical utility as new tools for early detection of melanoma, particularly as our results suggest that these biomarkers can detect disease much earlier than current diagnostic techniques. Furthermore, as blood-based biomarkers, there is potential for screening of non-symptomatic individuals. Whilst it is clear that much further validation is required, this study provides strong evidence that the circulating transcriptome holds much promise as a source of liquid biopsies for melanoma that surely merits further exploration.

Patient Cohorts
Patient plasma samples were collected both retrospectively (n = 119) and prospectively (n = 289). Retrospective samples were obtained from the John Radcliffe Hospital, Oxford (Oxford cohort; n = 30), and the AVAST-M multi-center phase 3 clinical trial (AVAST-M cohort; n = 89) [57]. Samples collected prospectively came either from the Hospital 12 de Octubre in Madrid (Madrid cohort; n = 67) or the Onkologikoa Cancer Hospital and Donostia University Hospital of San Sebastián (San Sebastián cohort; n = 102): a total of 327 melanoma patients ( Table 5). The clinical details of the individual patients included in these cohorts are given in Supplementary Tables S4-S7. Samples were collected at the time of diagnosis and prior to any treatment. Unfortunately, no information was available regarding the sequence of CDKN2A. Plasma samples from age/sex matched healthy controls (n = 99) were obtained from the Basque Biobank for Research O+EHUN. Plasma preparation was carried out within 1 h of phlebotomy, with blood collected in EDTA-coated tubes followed by centrifugation for 1000× g for 15 min at 4 • C. With the exception of the 37 prospectively collected samples used for NGS (Table 1), only a limited volume (1-2 mL) of plasma was available for validation studies. Therefore, we divided samples into three separate validation cohorts; a miRNA cohort of 96 melanoma patients and 28 controls; an mRNA cohort of 173 melanoma patients and 47 controls; and a YRNA cohort of 58 melanoma patients and 22 controls. The clinical details of the patients used are summarized in Table 1, and details of individual cohorts are provided in Table 4 and Supplementary Tables S4-S7. Written informed consent was obtained from patients for the inclusion of their samples in this study, and samples were collected in accordance with the Declaration of Helsinki and approved by local ethics committees (CEIC-Euskadi approval number PI2015024).
For the plasma samples used for NGS (5-8 mL), cfRNA was purified using the plasma RNA purification kit from Norgen Biotek (Ontario, Canada), and for validation studies (1 mL samples) the miRCURY™ RNA Isolation Kit Biofluids from Exiqon (Vedbaek, Denmark) was used.

Library Construction and Next-Generation Sequencing
The samples used for NGS were pooled according to disease stage, as shown in Table 1. Ribosomal RNA (rRNA) was removed from total cfRNA using the Ribozero Magnetic Human/Mouse/Rat kit (Epicentre (Maidon, WI, USA), #MRZH116), according to the low input protocol recommended by the manufacturer. Phosphatase and T4 polynucleotide kinase (PNK) treatments were carried out on the ribo-depleted RNA, and Illumina small RNA adapters ligated. Libraries were amplified using 15 cycles of PCR of barcoded primers [58]. Sequencing was performed on an Illumina HiSeq 2500 as 50 PE in rapid mode.

Bioinformatic Analysis
Sequencing reads were quality filtered using the fastx_artifacts_filter tool, and ligation adapters were removed using the AdRec.jar program (seqBuster suite of programs (omicX (Le-Petit-Quevilly, Le Petit-Quevilly, France))). Reads were mapped to the GRCh37 build of the human genome using the Bowtie 2.0 algorithm. A custom annotated probe set was built by combining probes from GENCODE version 8 [59], supplemented with rRNA and repeat annotations from RepeatMasker GRCh37, and snoRNA annotations from the UCSC table browser [60]. Expression of miRNA was calculated using the miraligner algorithm from the seqBuster suite, and YRNA expression was calculated using the HTseq-count algorithm. Differential expression analysis was carried out using the DESeq bioconductor package [61].
4.4. qRT-PCR (mRNA and miRNA) and ddPCR(yRNAs) mRNA was reverse transcribed (RT) using random primers with the High Capacity cDNA Reverse Transcription Kit from Applied Biosystems following the manufacturers' protocol with DNase treatment. Due to the difficulty in quantifying cfRNA reliably, we used fixed volumes in reactions [4,62]. For gene fragments, we designed custom Taqman probes using the Custom TaqMan ® Assay Design Tool from Applied Biosystems (Foster City, CA, USA), using the sequences corresponding to the respective probe sets (Supplementary Table S2 and Figure S1). The reference gene used for mRNA analysis was 18S rRNA, as previously described [63].
For miRNA detection, we used the Megaplex RT Primers Human Pool A v2.1 for RT, and specific Taqman probes as described in the text. As the problems of defining suitable reference genes for miRNA detection in plasma are well documented [4], we measured levels of three previously described reference miRNAs (miR-24, miR-16 and miR-191) [64,65], in a cohort of 63 samples (control/stage 0 (n = 22), stage I/II (n = 26), stage III/IV (n = 15)). Using the NormFinder algorithm we identified miR-24 and miR-191 as the most stable combination of reference genes in our samples (Supplementary Figure S3).
Custom primers and Taqman probes for Y3P1, Y4P1 and Y4P25 were designed using the Custom Taqman ® Small RNA Assay Design Tool from Applied Biosystems. We performed ddPCR using QX200TM Droplet DigitalTM PCR system (Bio Rad, Hercules CA, USA), following the manufacturers' protocol. Data analysis was performed by QuantaSoft analysis software from Bio-Rad. Expression levels were compared using the Kruskal-Wallis multiple comparison test, and the Mann-Whitney independent t-test to carry out a pairwise comparison between individual groups (Graphpad Prism v. 5.0, La Jolla, CA, USA). ROC analysis and comparisons were carried out using the method of DeLong et al., as implemented in MedCalc v. 14.8 software [66].

Conclusions
We have carried out a comprehensive, non-biased elucidation of the circulating transcriptome of melanoma patients and identified a number of promising candidate biomarker RNA species, not only miRNAs. These candidates were validated in independent cohorts by ourselves, however it is clear that further studies should be carried out by independent research groups in order to strengthen our findings and facilitate the translation of this knowledge into the clinic.
What is obvious is that many virtually unexplored classes of the circulating transcriptome are yet to be fully assessed for their ability to serve as useful cancer biomarkers. As a consequence, while the discovery of circulating miRNAs represented an important event in the history of the liquid biopsy field, it is clear that there is much that we have still to explore.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-6694/11/1/70/s1, Figure S1: Mapping positions (GRCh37) of selected gene fragment probe sets used to design Taqman probes. Figure S2: Variation of Ct values in levels of miRNA reference genes in selected cohort measured by qRT-PCR. Figure S3: ROC analysis of mRNA expression in different melanoma stage cohorts. Table S1: Summary of NGS results from pooled RNA samples. Table S2: Characteristics of selected differentially protein-encoding gene fragment probe sets. Table S3: YRNA reads from NGS. Table S4: Clinical details of patients used for NGS cohort. Table S5: Clinical details of patients used in mRNA validation cohort. Table S6: Clinical details of patients used in miRNA validation cohort. Table S7: Clinical details of patients used in yRNA validation cohort.