Twenty-two NB primary tumour samples, equally distributed between MYCN-amplified (MNA) and MYCN-non amplified (MNnA) samples, were obtained at diagnosis before any treatment from children admitted to the Department of Oncology at Alder Hey Children's NHS Trust, Liverpool. Histopathological diagnosis was confirmed using immunohistochemistry. Details of the patients are described in Table I. Institutional written informed consent was obtained from the patients parents or legal guardians.

The study underwent ethical review and approval according to the local institutional guidelines (Alder Hey Children's NHS Foundation Trust Ethics Committee, approval number 09/H1002/88).

Human DRG total RNA (Clontech Laboratories, Inc., Mountain View, CA, USA), pooled from normal tissues of 21 male/female Caucasians who had died suddenly, was used as a normal control. DRG were utilised as a normal control because of their neural crest-derived histogenesis (23) and since NB can occasionally arise from these organs (24).

Tumour samples were frozen in liquid nitrogen immediately after surgical removal and stored at −80°C. Total RNA was isolated by tissues ground under liquid nitrogen using 1 ml of TRIzol LS reagent (Life Technologies, Carlsbad, CA, USA) per 50–100 mg of sample according to the manufacturers protocol. RNA concentration and purity were measured by NanoDrop 2000 (Thermo Fisher Scientific, Inc., Waltham, MA, USA), having 260/280 ratio values in the range of 1.9–2.1.

Total RNA samples (5 µg) were enriched for small RNAs up to 200 bp by size selection using PureLink miRNA Isolation kit (Life Technologies) and RNA quality was assessed, before and after RNA enrichment, by the Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) using a Agilent RNA 6000 Nano and a Agilent Small RNA kit, respectively. Enriched RNA samples were processed using the Small RNA Expression kit according to the manufacturers protocol (Small RNA Expression kit, rev. B; Applied Biosystems, Foster City, CA, USA). Briefly, 20 ng of RNA was first hybridized and ligated with the adapter mix ‘A’, subsequently reverse transcribed and treated with RNAse H. The cDNA libraries were then PCR amplified, purified and size-selected by PAGE, resulting in libraries containing inserted small RNA sequences of 20–40 bp length. Size, integrity and purity of the libraries were verified by the Agilent 2100 Bioanalyzer, Agilent DNA 1000 kit. The cDNA libraries were barcoded using the SOLiD RNA barcoding kit and amplified onto beads using emulsion PCR. Templated beads were deposited on slides and analysed using the Applied Biosystems SOLiD 4 sequencer.

The qualified clean reads were mapped and analysed with the ‘small RNA’ bioinformatics pipeline from the Thermo Fisher LifeScope version 2.5.1 software (https://www.thermofisher.com/it/en/home/technical-resources/software-downloads/lifescope-genomic-analysis-software.html), using as a target the human genome GRCh38/hg38 and the dataset of mature and precursor miRNA sequences (miRBase, release 21.0) (25). Any sequence match against repetitive elements of the genome (SINEs and LINEs), and against non-miRNA small RNAs (snoRNAs, piRNAs, tRNAs and rRNA fragments) were filtered out from the results.

Sequence counts were extracted and reformatted with Genomnia Perl scripts from the pipeline output. Differential expression analysis was performed with the edgeR Bioconductor statistical library version 3.1.0 on R (26) incorporated in a proprietary Genomnia pipeline. TMM-normalized sequence counts in the libraries were transformed in counts per millions (CPM) according to the formula: CPM = (transcript normalized counts/total miRNA matches) × 1,000,000. After having estimated the tagwise dispersion, genewise exact test (26) as implemented in edgeR was used to measure the significance of differential expression, using the miRNA ‘Pseudo-counts’. Sequences were filtered out to remove poorly expressed microRNAs (i.e. miRNAs with <10 CPM in more than half of the samples). miRNAs were considered significantly differentially expressed if the false discovery rate (FDR; P-value corrected for multiple comparison with the Benjamini-Hochberg procedure) obtained with this method was <0.05. The MNA vs. MNnA comparison was also performed with the edgeR generalized linear model (GLM) mode optimised for the analysis of data with sequence spikes (27), in order to confirm the obtained results with a second robust statistical approach well suited to handle outlier data.

For each differentially expressed miRNA, a list of experimentally validated targets was extracted from miRTarBase release 6 (http://mirtarbase.mbc.nctu.edu.tw/) (28). In addition, de novo target prediction was performed using the miRanda software (29) on selected downregulated and upregulated microRNAs from the comparison MNA vs. MNnA. The 3UTR sequences of all human transcripts present in Ensembl (release 78) were used as a target set for this analysis, and only predictions with total score ≥155 and total energy ≤-20 were considered as putative miRNA targets. Annotation and enrichment of functional pathways and ontology terms associated with the selected target genes were evaluated using DAVID Bioinformatics Resources 6.7 (30). The same target gene list was used as the starting dataset for the generation of a Functional Interaction network analysis and related Gene Ontology enrichment analysis with the Cytoscape Reactome FI plugin (31,32). Promoter regions of the differentially expressed miRNAs in the comparison MNA vs. MNnA were scanned for motifs by FIMO 4.10.2 software (http://meme-suite.org/tools/fimo) (33). A sequence of 5 Kb was considered as putative upstream regulatory region for each miRNA or miRNA cluster based on the histone acetylation track (H3K27Ac) reported in UCSC. Finally, the MYCN binding motifs available in JASPAR CORE database (http://jaspar.genereg.net/) were used in the promoter analysis.

The alignment files in BAM format corresponding to the same biological group (MNA, MNnA and Control) were merged and converted to sequence alignment/map (SAM) format with samtools (34). The files were then processed and analysed with the miRDeep2 software for miRNA prediction (35). The differential expression analysis for isomiRs was carried out with edgeR from Bioconductor, with the same analytical strategy previously described. A set of potentially novel microRNAs was also identified, by selecting among miRDeep2 predictions the hairpins with a probability >60% to be a true positive. Their expression was also evaluated in all samples. miRNA extraction, library preparation, sequencing and bioinformatics analysis were performed by Genomnia s.r.l.

Quantitative real-time PCR (qPCR) was performed to validate some results obtained from the high-throughput sequencing. In brief, step-loop reverse transcription and qPCR were carried out using TaqMan MicroRNA assays (Life Technologies) specific for human miR-18a-5p, miR-217, miR-491-5p and miR-628-3p, according to the manufacturers protocol. Ct values for miRNAs were normalized against U6 small nuclear RNA (RNU6) levels, used as an internal control. The relative expression of each miRNA was calculated by the comparative Ct method and expressed as fold-change. All samples were run in triplicate reactions in the StepOnePlus Real-Time system (Life Technologies) instrument. Results are expressed as means ± standard deviation (SD), and a P<0.05 was considered statistically significant.

The MYCN-amplified KCNR cell line was grown in RPMI-1640 medium supplemented with 10% fetal bovine serum (FBS), 1% v/v non-essential amino acid solution, 2 mM glutamine, 100 U/ml penicillin and 100 µg/ml streptomycin. HEK (human embryonic kidney) 293T cells were cultured in high-glucose DMEM, supplemented with 10% FBS, 2 mM L-glutamine, 100 U/ml penicillin and 100 µg/ml streptomycin. All cell lines were maintained in a humidified incubator at 37°C and 5% CO₂, and appropriately sub-cultured before confluence was reached.

The 3UTR region of the human MYCN gene (3HSAR055440 derived from http://utrdb.ba.itb.cnr.it/) was amplified from genomic DNA using a proof reading Phusion High-Fidelity PCR Master Mix (Finnzymes, Espoo, Finland) with the following primers: MYCN 3UTR-F 5-GCTCTAGAACGCTTCTCAAAACTG GACAGTCAC-3 and MYCN 3UTR-R 5-GCTCTAGAAACATGAGGTATTTCAAAGTGCTATAAGATGC-3. The 880-bp amplified product was inserted into the XbaI site (underlined primer sequences) of the pRLTK vector (Promega, Madison, WI, USA) immediately downstream of the Renilla luciferase gene. The luciferase reporter construct, termed pMYCN-3UTR, was validated by specific restriction digestions as well as bi-directional sequencing by Bio-Fab Research (Rome, Italy). HEK293T cells were plated at a density of 1.5×10⁵/well in 24-well plates and co-transfected after 24 h with 50 ng of pMYCN-3UTR Renilla luciferase construct, 10 ng of reference firefly luciferase reporter (pGL3-SV40; Promega) and 100 nM of mir-628-3p (mirVana miRNA mimic - MI0003642/MIMAT0003297 by Thermo Fisher Scientific) or miR-Ctr (miRIDIAN microRNA Mimic Negative Control #1; Dharmacon, Lafayette, CO, USA) molecules using 2 µl of Lipofectamine 2000 (Invitrogen). Cells were lysed 72 h after transfection and luciferase signals were detected by the Dual-luciferase reporter assay (Promega) using the GloMax Discover System (Promega), following the manufacturers instructions. Renilla luciferase activity was normalized to the firefly luciferase activity for each reaction. Two independent experiments were carried out, each in triplicate. All the luciferase data are expressed as the mean ± SD normalized to the negative control miRNA.

KCNR cells were seeded at 6.0×10⁵/well in 12-well plates and transfected the day after with miR-628-3p or miR-Ctr (50 mM final concentration) using Lipofectamine 2000, according to the manufacturer's instructions. After 72 h, protein expression was determined by western blotting as previously described (36) with a primary antibodies against MYCN (1:1,000; Santa Cruz Biotechnology, Santa Cruz, CA, USA). Tubulin (Sigma-Aldrich, St. Louis, MO, USA) was used as a normalization control for equal loading. Chemiluminescent signals were revealed by ChemiDoc XRS+ System (Bio-Rad Laboratories, Hercules, CA, USA) and densitometric analysis was performed by Image Lab 5.1 software (Bio-Rad Laboratories).

Comprehensive miRNA, isomiR and novel miRNA identification, annotation and analysis were developed in the framework of this project and are summarized in Fig. 1. The miRNA transcriptome was profiled by sequencing of mature miRNA molecules in a panel of primary NB tumours, 4 with MYCN-amplification (MNA) and 4 without MYCN-amplification (MNnA) and from normal dorsal root ganglia (DRG). The clinical characteristics of the patients are reported in Table I. Small RNA libraries were prepared and deep sequenced by using a SOLiD 4 sequencer platform. Approximately 71 million reads mapped with high quality against the human reference genome (GRCh38/hg38, repeat masked), with an average of 8 million reads per library. The reads corresponding to annotated miRNAs were identified by mapping all reads against the human precursor and mature sequences included in miRBase v21.0. Fig. 2A shows that a total of 532 distinct known miRNAs had on average at least 10 read counts in the NB tumour samples, while 512 miRNAs were detected in the DRG sample; 382 of these were expressed in both conditions, with a 72% overlap for NB miRNAs and 75% overlap for DRG miRNAs. Thus, about a quarter of the identified known microRNAs were uniquely expressed in NB or control samples. NB-specific miRNAs might therefore play an important functional role in tumour development and progression, through the regulation of downstream targets and the biological pathways affected by these genes.

Based on the high-throughput sequencing of small RNAs, we performed differential expression analysis of miRNAs comparing the two libraries from NB tumours (MNA and MNnA collectively considered) and DRG, after removing poorly expressed miRNAs. Of a total of 198 differentially expressed miRNAs in the dataset, 128 (64.6%) displayed statistically significant differential levels in NB tumours in comparison with DRG, after FDR correction. Notably, out of these 128, 47 miRNAs (36.7%) were expressed at higher levels in NBs (Table IIA), while 81 (63.3%) were expressed at lower levels (Table IIB). Fig. 2B shows the distribution of the different mature miRNA species in abundance classes comparing NB tumours (MNA and MNnA collectively considered) against DRG, allowing a survey of the entire miRNA population. Significantly dysregulated miRNAs are indicated in red. Among the differentially expressed miRNAs, some have been previously associated with NB development, such as upregulated miR-876-5p, miR-19a, as well as downregulated miR-183-5p, miR-423-3p and miR-92b-3p (18,36–39). Other miRNAs have been implicated in the pathogenesis of many types of human cancer but not in NB malignancy to date. In particular, of main interest are upregulated miR-127-5p, miR-369-3p and miR-655-3p encoded in the large 14q32 miRNA cluster and associated with oligometastases through suppression of shared target genes involved in cytoskeletal organization, cell motility and TGF-β signalling pathways (40). Also the marked downregulation of miR-328-3p in NB samples in comparison to DRG suggests its possible role as a tumour suppressor also in this malignancy as already observed in other types of cancers, such as acute myeloid leukaemia, chronic myelogenous leukaemia and glioblastoma (41,42). Therefore, future experiments will address the role of this miRNA in the proliferation, invasion and metastasis of NB tumour cells. Another new attractive result is the significant downregulation of miR-1247-5p, whose suppression has been demonstrated to play an important role in pancreatic cancer (43,44). Notably, the molecular mechanisms underlying the role of miR-1247-5p have been correlated with its capacity to directly target and silence the expression of neuropilins (44), which are transmembrane glycoproteins recently implicated in tumour functions and found as promising drug target candidates in leukaemia and lymphoma (45). Furthermore, miR-1247-5p also belongs to the 14q32 miRNA cluster (but in reverse strand in respect to the above mentioned miR-127-5p, miR-369-3p and miR-655-3p) and its aberrant expression in colorectal cancer has been linked to hypermethylation, confirming the importance of epigenetic modifications in miRNA-mediated cancer development (46). Since alterations in DNA methylation are frequent in NB and can also modulate sensitivity to cytotoxic agents (47), the analysis of the miR-1247-5p promoter region in NB samples could provide new insight into the mechanisms that correlate epigenetic processes and NB development.

To investigate the biological implications of the dysregulated miRNAs, we assessed miRNA target genes within the regulatory network for enrichment in Gene Ontology (GO), one of the most useful methods for functional annotation and classification of gene products. GO categories of the targets of the differentially upregulated and downregulated miRNAs (NB vs. DRG) were organized into three different functional groups, i.e. Biological process, Cellular component and Molecular function (Fig. 3), among which a significant over-representation of transcription factor activity/binding, protein kinase activity, chromatin modifications, metabolic processes, cell communication, plasma membrane, cytoplasm/intracellular organelle part, cell proliferation and programmed cell death regulation were found. To further explore the biological pathways potentially affected in NB development, targets of the differentially upregulated and downregulated miRNAs (NB vs. DRG) were subjected to KEGG pathway analysis. Several cancer-related networks, including TGF-β, neurotrophin, MAPK, ErbB, Toll-like receptor and p53 signalling pathways, were significantly enriched (Tables III–IV). Indeed, the expression levels of TGFBR2, SMAD2 and SMAD4, three important components along the TGF-β cascade, were found to be downregulated in NB cells (48). VEGF, Wnt, JAK-STAT and GnRH pathways were found among the targets of the underexpressed miRNAs in accordance with their growth-inhibitory potential on NB (49–51). Moreover, some biological processes linked to cell cycle, adherents junction and focal adhesion could have also been altered. KEGG database interrogation also suggested that many of the miRNA targets contribute to the development of other human tumours, including colorectal, pancreatic cancer, melanoma and gliomas.

The miRNA expression profiles of MNA and MNnA tumours largely overlapped as shown by a Venn diagram (Fig. 4A), suggesting that only a small number of miRNAs are specific for each subgroup. After removing poorly expressed miRNAs, 70 molecules were found to be differentially expressed between MNA and MNnA tumours. Of these, 23 had a P<0.05 after FDR correction: 17 miRNAs were significantly upregulated in MNA samples, while 6 miRNAs were significantly downregulated (Table V). The most homogeneously expressed molecules among different tumours inside each subgroup were miR-16-2-3p, miR-18a-5p, miR-217, miR-338-5p, miR-489-3p, miR-491-5p, miR-628-3p, miR-628-5p and miR-3065-3p (data not shown). Among these, miR-18a-5p, miR-217, miR-491-5p and miR-628-3p were selected as candidates to validate the miRNA expression data by qPCR experiments in the 8 deep-sequencing-analysed NB tumours and in 14 additional NB samples as well as in 4 human NB cell lines (MNA LAN5 and KCNR and MNnA SH-SY5Y and SK-N-SH). qPCR demonstrated that the expression levels of the selected dysregulated miRNAs showed strong overlap with the deep-sequencing findings, thus confirming the overexpression of miR-18a-5p and miR-217 as well as the downregulation of miR-491-5p and miR-628-3p in the MNA (Fig. 4B) in comparison with the MNnA tumours. The same trend was found in MNA vs. MNnA NB cell lines (data not shown). The differential expression of a substantial number of miRNAs between favourable and unfavourable NBs was consistent with previous reports (Table V), reinforcing their fundamental role in NB development (18,48,52,53,56,62,63,74,76). In particular, the present study confirms how the majority of the differentially expressed miRNAs positively correlate with MYCN amplification, suggesting that miRNA induction might be a general mechanism caused by MYCN oncogene in high-risk NB tumours. Additionally, bioinformatic analysis of the promoter regions of the overexpressed miRNAs identified MYCN-binding sites (38,54) in 5 out of 17 miRNAs (Table V), and this is consistent with a direct transcriptional regulation of these molecules by MYCN. Overexpression of miR-17-5p, miR-18a-5p, miR-20a-5p and miR-92a-3p, belonging to the polycistronic miR-17-92 cluster at 13q31.3, was consistently higher in MNA cell lines compared to the MNnA tumours as previously reported (55,56) and in line with the evidence that MYCN directly activates miR-17-92 expression in NB cell lines (55). Indeed, the miR-17-92 cluster has a unique promoter control region and encodes for a precursor transcript that contains six tandem-loop hairpin structures, which in turn produce the six mature molecules miR-17, miR-18a, miR-19a, miR-20a, miR-19b-1 and miR-92a-1 (57,58). These miRNAs may act independently or synergistically to target multiple mRNAs encoding for proteins implicated in cancer-related pathways (57,58). Given the miR-17-92 oncogenic role, inhibition of the clustered oncomiRs by delivering single-stranded antisense oligonucleotides (59) may represent a novel and attractive treatment approach for NB therapy. Regarding this aspect, Fontana et al (60) demonstrated that in vitro and in vivo treatment of MNA NB cells with a specific antagomiR against miR-17-5p was able to dramatically reduce proliferation and to activate massive apoptosis through p21 and BIM over-expression, respectively. Concerning some downregulated miRNAs in MNA tumours, our data are in agreement with recently published studies. In particular, miR-491-5p was found to be significantly downregulated at NB metastatic sites compared with primary tumours by Guo et al (61), and miR-628 was confirmed to be moderately expressed in favourable NBs and nearly absent in unfavourable cases, thus, supporting its tumour-suppressive function and its potential as a good candidate in discriminating between these two NB subsets (38,62,63). Indeed, among the downregulated miRNAs, miR-628-3p has one pairing region at position 718–724 of the MYCN 3 untranslated region (3UTR; Fig. 5A), as predicted by Miranda software (best score ≥120). To verify whether miR-628-3p was able to directly target MYCN mRNA, HEK293T cells were co-transfected either with a reporter construct containing 880 base pairs of the human MYCN 3UTR downstream of the Renilla luciferase open reading frame (Fig. 5A) with either miR-628-3p mimic or negative control miRNA (miR-Ctr). Overexpression of miR-628-3p upon transfection was confirmed by qPCR assays (data not shown). Luciferase activity was specifically reduced to ~35% (P<0.01) in cells transfected with miR-628-3p compared to the control group (Fig. 5B), thus indicating that miR-628-3p functionally interacts with the MYCN 3′UTR and decreases the expression of the corresponding protein. Indeed, miR-628-3p mimics transiently introduced into KCNR cells, an in vitro model of MNA NB, were able to slightly reduce endogenous MYCN protein levels as compared with miR-Ctr positive samples (Fig. 5C), suggesting a functional role of miR-628-3p in MYCN gene expression. Finally, our analysis identified a set of differentially expressed miRNAs, such as miR-130b-5p, miR-217, miR-3065-3p and miR-3182, known as oncomiRs in several other tumours (64–67), but whose role in NB has not yet been investigated (Table V). Notably, miR-130b and miR-301b belong to the same genomic cluster located on chromosome 22 and are significantly upregulated in triple-negative (lacking the expression of the oestrogen receptor, the progesterone receptor and the human epidermal growth factor receptor 2) breast cancer (64), having a direct role in cyclin G2 regulation. miR-338 and miR-3065 belong to the brain-specific genomic cluster located on chromosome 17 in an intronic region within the apoptosis-associated tyrosine kinase (AATK) gene, thus having a tumour-suppressor function (68). Accordingly, it was described that miR-338-3p strand has the potential to inhibit cell proliferation, cell migration and invasion in human NB cells by affecting the PTEN/Akt pathway (69). miR-3065 was recently mapped to the antisense strand of miR-338 (70) and regarded as a novel promising candidate contributing to clear cell renal cell development (71). Furthermore, our miRNA promoter in silico analysis showed that downregulated miR-338-5p, miR-489-3p, miR-628-3p and miR-3065 harbour putative MYCN binding motifs in their regulatory regions (Table V), suggesting an inverse correlation and possible feedback regulatory circuits between MYCN amplification and the reduced levels of these miRNAs. To this concern, recent studies have supported the dual role of MYCN in controlling gene expression, being not only a transcriptional activator but also a mediator of transcriptional silencing by facilitating DNA methylation (54,72,73). Thus, it will be interesting to characterise the molecular mechanisms that regulate miRNA expression.

Using the miRTarBase web tool, we found a relevant number of cancer-related genes, such as PTEN, KRAS and EZH2 among the validated target genes for the differentially expressed miRNAs. Indeed, since the availability of experimentally validated targets in miRTarBase was variable and possibly biased towards the more extensively investigated miRNAs, we decided to extend our analysis and perform a de novo target prediction using the miRanda software. In particular, we considered miRNAs with a homogeneous expression level between tumours inside each group: upregulated miR-16-2-3p, miR-18a-5p, miR-217 and downregulated miR-338-5p, miR-489-3p, miR-491-5p, miR-628-3p, miR-628-5p and miR-3065-3p. The 3UTR sequences of all human transcripts present in Ensembl were screened for the identification of putative microRNA binding sites: a total of 867 putative target genes were identified for the upregulated miRNAs, while 2148 target genes were identified as potentially modulated by the downregulated miRNAs. The predicted targets were grouped into functional modules on the basis of their functional interactions, and for each module an enrichment analysis was performed with the Cytoscape Reactome FI plugin. Fig. 6 shows the main functional module putatively regulated by miR-16-2-3p, miR-18a-5p and miR-217, which are overexpressed in MNA tumours. The different interaction networks are enriched in gene products implicated in cancer-related biological processes, such as DNA-repair and apoptosis, as well as axon guidance, adenylate cyclase activity and the FGFR signalling pathway. Target genes associated to each specific overexpressed miRNA are reported in Table VI. The GO categories associated with the predicted miRNA targets showed significant overlap with biological processes that are over-represented in the MYCN target genes (38), confirming the existence in NB of a miRNA:MYCN related network able to broadly perturb gene expression (38,74,75). Fig. 7 shows the main functional module putatively regulated by miR-338-5p, miR-489-3p, miR-491-5p, miR-628-3p, miR-628-5p and miR-3065-3p which are underexpressed in MNA samples. The enriched interactions consist of proteins involved in synaptic transmission as well as TLR and EGFR signalling pathways. Target genes associated to each specific downregulated miRNA are reported in Table VII. The complete list of network modules and GO BP - pathway enrichment is included in the Table VIII. Our results suggest that the selected differentially expressed miRNAs are relevant in NB biology and could be of potential importance to stratify high-risk and low-risk tumours. Further investigations will be necessary to validate the specific target genes and to assess the downstream effects.

Recent findings identified a large number of sequences with different 5 or 3 flanking regions, called isomiRs (76). From a biological point of view, the most interesting isomiRs are the 5′-isomiRs (having different 5-ends compared to their corresponding canonical miRNAs) as they are potentially able to recognize different targets with respect to the main isoform. A list of isomiRs (length and sequence variants) was generated from the miRDeep2 output with proprietary scripts, and their expression was evaluated in all samples. As expected, compared to the total number of isomiRs detected in our dataset, 5-isomiRs were less frequent: out of 18,809 isomiRs with at least 10 read counts in at least one library, 3,769 were 5-isomiRs (20%), as described in Fig. 8A. The presence of miRNA length variants suggests that these isomiRs might have different regulatory functions or efficiencies, this allowing a finely tuned control of global gene expression. Furthermore, in silico data analysis identified 18 putative novel miRNAs expressed in the tumour samples. Two of these were differentially expressed in MNA vs. MNnA tumours. In particular, a hairpin with genomic coordinates chr4:180531388-180531451 was detected as significantly overexpressed in the MNA tumours (P=0.001; FDR=0.02) using the ‘GLM’ analysis, and it demonstrated a uniform expression within each biological subgroup. The secondary stem-loop structure of this putative novel miRNA is shown in Fig. 8B. According to the miRanda predictions, this candidate miRNA has ~900 putative targets, which are enriched in transcriptional regulation and protein phosphorylation ontology terms. Interestingly, a significant number of putative target genes, highlighted in red in Fig. 8C, are correlated with the mTOR signalling cascade, outlining the evidence for the aberrant activation of this pathway in NB (77–79). For this reason, the putative new miRNA may be an interesting candidate for future functional studies in order to determine its biological relevance in mTOR signalling, which in turn might lead to the identification of further molecular targets in NB tumour.

In conclusion, in the present study, NGS technology and bioinformatics analysis allowed the identification of a specific set of miRNAs that are differentially expressed in NB tumours compared to DRG, or in NB tumours with or without MYCN amplification. A relevant percentage of these microRNAs has been already found to be involved in NB, another percentage has been shown to be involved in the pathogenesis/progression of other tumours, while other miRNAs are new. Target prediction and functional annotation highlighted the possible roles of the identified microRNAs, and in vitro assays provide initial evidence of the contribution of miR-628-3p to MYCN regulation. New potential miRNA candidates not previously reported were also found, but further studies are needed to clarify their role in NB development.

Altogether, these findings may provide further insight concerning NB development and may also contribute to the identification of novel potential therapeutic targets in this highly lethal childhood tumour.