MicroRNA signature and integrative omics analyses define prognostic clusters and key pathways driving prognosis in patients with neuroendocrine neoplasms

Neuroendocrine neoplasms (NENs) are mutationally quiet (low number of mutations/Mb), and epigenetic mechanisms drive their development and progression. We aimed at comprehensively characterising the microRNA (miRNA) profile of NENs, and exploring downstream targets and their epigenetic modulation. In total, 84 cancer‐related miRNAs were analysed in 85 NEN samples from lung and gastroenteropancreatic (GEP) origin, and their prognostic value was evaluated by univariate and multivariate models. Transcriptomics (N = 63) and methylomics (N = 30) were performed to predict miRNA target genes, signalling pathways and regulatory CpG sites. Findings were validated in The Cancer Genome Atlas cohorts and in NEN cell lines. We identified a signature of eight miRNAs that stratified patients in three prognostic groups (5‐year survival of 80%, 66% and 36%). Expression of the eight‐miRNA gene signature correlated with 71 target genes involved in PI3K–Akt and TNFα–NF‐kB signalling. Of these, 28 were associated with survival and validated in silico and in vitro. Finally, we identified five CpG sites involved in the epigenetic regulation of these eight miRNAs. In brief, we identified an 8‐miRNA signature able to predict survival of patients with GEP and lung NENs, and identified genes and regulatory mechanisms driving prognosis in NEN patients.

The low incidence, anatomical dispersion, and molecular and clinical heterogeneity of NENs have hindered an in-depth characterisation of the molecular landscape of these tumours.These neoplasia, and particularly NETs, are mutationally quiet, harbouring relatively few somatic mutations as compared to other malignancies [7,8].Studies in sporadic NENs have reported few recurrently mutated genes: MEN1 (20%) and ARID1A (8%) in lung NETs [8], MEN1 (40%) and ATRX/DAXX (35%) in pancreatic NETs [9] and CDKN1B (10%) in small intestine NETs [10].This would suggest that there are other alternative mechanisms implicated in the development and progression of these tumours.Consistent with this, several recent reports suggested that epigenetic changes may play a key role in the pathogenesis of NENs [10][11][12].Nevertheless, a better characterisation of these nonmutational alterations is required to improve our understanding of the molecular mechanisms involved in the biological and clinical behaviour of NENs.Alterations in epigenetic regulators such as microRNAs (miRNAs) or DNA methylation could be relevant in this context as they can modulate gene expression and, consequently, specific oncogenic or tumour suppressive pathways.Moreover, as these changes are reversible, epigenetic modulation could potentially be exploited for therapeutic purposes.
In this context, the aim of this study was to comprehensively characterise the miRNA profile of NENs and to explore, through integrative omics analyses, downstream targets and epigenetic modulation to identify key pathways involved in NEN biological and clinical behaviour.

Patients and samples
The study population included 85 patients with G1-G3 NENs of gastroenteropancreatic (GEP; N = 34) and lung (N = 51) origin diagnosed between 2000 and 2015 (small cell neuroendocrine carcinoma of lung origin was excluded) in the Hospital Universitario Virgen del Rocıó (HUVR) and Hospital Universitario Central de Asturias (HUCA).Archival formalin-fixed paraffin-embedded (FFPE) tumour and nontumour tissue samples from each patient were retrospectively collected and assessed by expert NEN pathologists prior to molecular analysis.
The study was carried out in accordance with the Declaration of Helsinki and approved, together with the protocols, by the Hospital 12 de Octubre ethics committee (Committee Register Number: 18/449).All patients provided written informed consent prior to the study entry.

Nucleic acid extraction
Total DNA and RNA were extracted using standard methods (QIAamp DNA FFPE Tissue Kit (QIAGEN, Valencia, CA, USA) and Recover All TM Total Nucleic Acid Isolation Ambition Kit (Applied Biosystems, Foster City, CA, USA)) according to the manufacturers' instructions.Total RNA was extracted from cell lines using the miRNeasy Mini Kit (QIAGEN) according to the manufacturer's instructions.DNA and RNA integrity, purity and concentration were assessed with the Qubit ® Fluorimeter (Life Technologies, Carlsbad, CA, USA) or Nanodrop 1000 (NanoDrop, Wilmington, DE, USA) instruments, respectively.

Analysis of miRNA expression profiles
Quantitative Real-Time PCR (qRT-PCR) was performed to determine the expression levels of mature human miRNAs.Complementary DNA (cDNA) was obtained by reverse transcription (RT) of extracted RNA from tumour and nontumour patient samples using the cDNA miScript II Reverse Transcription kit (QIAGEN).Expression levels of miRNAs were determined using the Human Cancer Pathway Finder miScript miRNA PCR Array (MIHS-102ZE; QIAGEN), which included 84 of the most relevant cancer-related human miRNAs.Reaction conditions for amplification were as follows: first step of 95 °C 15 min, then 40 cycles of three-step 94 °C 15 s, 55 °C 30 s and 70 °C 30 s.For miRNA expression analysis, we have used nontumour paired tissues to normalise the miRNA data for each tumour sample.Expression of miRNAs for each NEN patient was assessed by a relative quantification approach, in which the amount of each miRNA is expressed in relation to the geometric average of five reference housekeeping genes (small nucleolar RNAs: SNORD61, SNORD68, SNORD95, SNORD96A, RNU6B) using the 2 ÀΔΔCt method [13]: ΔC t = C t (miRNA) − C t (average of the housekeeping genes); ΔΔC t = ΔC t (tumour) − ΔC t (nontumour); 2 ÀΔΔCt = 'fold change': relative expression of the miRNA in tumour sample with respect to the normal tissue.Expression of eight selected miRNAs (miR-17-5p, miR-18a-5p, miR-19a-3p, miR-20a-5p, miR-20b-5p, miR-92a-3p, miR-203a-3p and miR-210-3p; reference numbers 002308, 002422, 000395, 000580, 001014, 000431, 000507, and 000512, respectively; Thermo Scientific) was also quantified in 40 NENs with available tumour and nontumour samples using TaqMan qPCR as an alternative method (Table S1).Geometric average of SNORD95 and SNORD61 (reference numbers Hs03464473_s1 and Hs03298134_s1, respectively; Thermo Scientific) was used to normalise the data.In all, 10 ng of total RNA was reverse transcribed using the TaqMan miRNA reverse transcription kit in a total volume of 15 μL, according to the manufacturer's protocol.Thereafter, 0.7 μL of cDNA was used for TaqMan MicroRNA Assays.The reactions were incubated at 95 °C for 10 min, followed by 40 cycles of two-step 95 °C 15 s and 60 °C 1 min.All reactions were performed in an Applied Biosystems ViiA 7™ Real-Time PCR instrument (BioRad, Hercules, CA, USA).

Gene expression and DNA methylation profiling
Clariom S Human microarrays (Affimetrix, Santa Clara, CA, USA) representing more than 23 000 unique transcripts were performed in all tumour samples with sufficient remaining RNA (N = 63).After assessing the RNA integrity, labelling and hybridisation were performed according to Affimetrix protocols.Briefly, 50 ng of total RNA was amplified and labelled using the Pico reagent kit (Affimetrix) and then hybridised to Clariom S Human microarrays (Affimetrix).Washing and scanning were performed using the Affymetrix GeneChip System (GeneChip Hybridization Oven 645, GeneChip Fluidics Station 450 and Gene-Chip Scanner 7G).

DNA methylation profiling
Total DNA was extracted from paired FFPE tumour (N = 30) and nontumour (N = 10) tissue samples using the QIAamp DNA FFPE Tissue Kit (QIAGEN) according to the manufacturer's recommendations.DNA purity and concentration were assessed with the Qubit ® Fluorimeter (Life Technologies).Quality control and normalisation were performed with the Infinium HD FFPE QC assay (Illumina, San Diego, CA, USA) kit.Samples that passed the quality control were subjected to bisulfite conversion using 250 ng of DNA and the EZ-96 DNA Methylation Kit (Zymo Research, Irvine, CA, USA).Converted DNA was processed with the Infinium HD FFPE Restore Protocol (Illumina), and the ZR-96 DNA Clean & Concentrator-5 (Deep Well) (Zymo Research) kit to restore FFPE DNA.Methylation analysis was performed using Infinium MethylationEPIC BeadChip (Illumina), which allows extensive coverage of CpG islands and gene enhancers.Fluorescence detection and scanning were performed with the iScan System (Illumina) platform.
The obtained raw data files (.IDAT) were processed and filtered using the ChAMP pipeline available as a R package.The beta-values obtained were transformed into M-values to normalise distribution of the data.Methylation levels of CpG sites (beta-values) were obtained for each sample and were compared between clusters defined by miRNAs expression in tumour and nontumour samples.

qRT-PCR of selected miRNAs and target genes
MiR-17-5p and miR-19a-3p expressions in transiently transfected cell lines were assessed as described in the analysis of miRNA expression profiles.Selected miRNA target genes were analysed in transiently transfected cell lines by RT of 1000 ng of RNA using the High-Capacity cDNA Reverse Transcription Kit (Thermo Scientific) according to the manufacturer's instructions.cDNA was amplified using Hs01093259_m1).ACTB (ID: Hs01060665_g1) was used as a reference gene to normalise cycle thresholds (C t s).Relative miRNA/mRNA expression levels were determined using the 2 ÀΔΔCt method [13].All reactions were performed in an Applied Biosystems ViiA 7™ Real-Time PCR machine (BioRad).

Statistical analysis
Descriptive statistics were used to characterise most relevant clinical parameters.The distribution of quantitative variables among study groups was evaluated by parametric (Student's t-test or ANOVA) and nonparametric (Kruskall-Wallis or Mann-Whitney) tests, as appropriate.OS was calculated as the time elapsed from the date of diagnosis to the date of death from any cause or last follow-up of living patients.The Kaplan-Meier method was used to estimate time-dependent variables (OS), and differences observed among patient subgroups were assessed by the log-rank test.The prognostic value of all miRNAs was analysed by a Cox univariate regression model considering their expression as a continuous variable.P-values were corrected for multiple testing by the false discovery rate (FDR) method, considering FDR ≤ 0.05 to indicate significance.Finally, a Cox multivariate proportional hazards model was used to further assess the prognostic value of selected miRNAs adjusted for other potentially confounding variables, such as gender, age, grade, TNM stage and primary tumour site (GEP or lung).
Spearman correlation analyses between selected miRNAs and in silico predicted target genes were performed.In addition, Cox regression analyses were performed to interrogate the role of significantly correlated target genes in OS (see Section 2.9).Correlations between miRNAs and CpG sites at regulatory regions (−50 kbp upstream from transcription start sites (TSS)) [14] were also assessed using Spearman correlations, with the Benjamini-Hochberg method used to adjust for multiple testing.MiRNA-CpG sites with an FDR < 0.05 were considered significant.These CpG sites were further explored, and differences between tumour and normal samples, and between prognostic clusters were assessed with Student's t-tests and ANOVA tests, respectively.P-values < 0.05 were considered significant.All statistical analyses were performed using SPSS 21.0 (IBM) and R software (v4.1.2).

Bioinformatic analysis
In silico predicted targets for the eight miRNAs significantly associated with patient OS were identified using target genes consistently predicted in all three different databases explored: TargetScan Release 7.2 [15], miRDB [16] and Diana-Tarbase v7.0 [17].
A heatmap showing the Spearman correlation coefficients among the 84 cancer-related miRNAs expressions was created using the MORPHEUS software (https://software.broadinstitute.org/morpheus/).Rows and columns of the heatmap were hierarchically clustered with 1-Spearman rank correlation as distance.On the contrary, a heatmap showing the expression levels of the eight miRNAs and the prognostic groups was obtained by COMPLEXHEATMAP [18] of the R software package using the UPGMA hierarchical clustering algorithm with Euclidean distance over the standardised Log (miRNA expression + 1).

Public database analysis
Normalised gene expression and miRNAs data from The Cancer Genome Atlas (TCGA) cohorts (Ovarian Serous Cystadenocarcinoma (TCGA, Nature 2011); Breast Invasive Carcinoma (TCGA, Nature 2012); Glioblastoma (TCGA, Nature 2008); Kidney Renal Clear Cell Carcinoma (TCGA, Nature 2013); Lung Squamous Cell Carcinoma (TCGA, Nature 2012); Colorectal Adenocarcinoma (TCGA, Nature 2012); and Pediatric Rhabdoid Tumor (TARGET, 2018)) were analysed using cBioPortal for Cancer Genomics (https://www.cbioportal.org/)[23].The expression of selected genes was evaluated using z-scores of either log RNA-seq RPKM (Reads Per Kilobase Million), log RNA-seq V2 RSEM or log microarray.Accordingly, miRNA expression was evaluated using z-scores obtained from microarray expression data.Pearson's correlation analysis was assessed for all miRNA-target gene pairs of interest.Values of P < 0.05 were considered significant.Only samples with both miRNA and mRNA expression data were included in the analysis for each TCGA study.

Characteristics of the study population
The main clinical and pathological features of the study population are summarised in Table S1.With a median follow-up of 131 months (range 26-229 months), 41 patients (48.2%) have died.The median OS of the study population was 91 months (range 0-229 months), with a 5-year OS rate of 65%.

MiRNA and mRNA integrative expression analyses enabled the identification of key genes and pathways involved in NEN's biology
To further investigate which genes and pathways are regulated by the eight prognostic miRNAs, an integrative gene expression analysis of NEN samples (N = 63) was performed (Table S1).First, we performed correlation analysis between the miRNA expression and their in silico predicted target genes (N = 62) (Table S3; see Section 2).Seventy-one of these target genes were found to be significantly correlated (P < 0.05) with their corresponding miRNA following integrative expression analyses performed on tumour samples from our study population; 41 of them were inversely correlated (r < 0), and 30 were directly associated (Table S4).ORA of these 71 target genes revealed involvement of the precited target genes in key signalling pathways including Hallmarks TNFα Sig. via NF-kB and KEGG PI3K-Akt signalling pathways, among others (FDR < 0.25) (Table S5).All of these dysregulated pathways are implicated in biological functions related to almost all hallmarks of cancer [25] (Fig. 1B).
We next interrogated the prognostic impact of the 71 target genes and found that 28 of them were significantly associated with OS (P < 0.05) (Table 1).The genes associated with poor prognosis were positively correlated with miRNA expression (r > 0) and those associated with better prognosis were negatively correlated with their corresponding miRNA (r < 0).Interestingly, a String network of these 28 prognostic genes revealed seven functionally related genes involved in mitotic spindle, G2/M checkpoint, apoptosis and cell cycle and proliferation (Fig. S3).

In vitro and in silico validation of target genes of the eight prognostic miRNAs
To further validate a direct correlation between the eight miRNAs and some of their target genes observed in patient samples, we explored target gene expression upon downregulation of miR-17-5p and miR-19a-3p in BON-1 pancreatic and H727 lung NET cell lines.Basal expression levels of miR-17-5p and miR-19a-3p in NET cell lines were similar to those obtained in the good prognosis cluster of NENs, whereas a nonsignificant trend was observed towards higher miRNA levels in tumour samples of patients of the poor prognosis cluster (Fig. S4).As expected, inhibition of miR-17-5p in the H727 cell line significantly increased the expression of its targets CRY2 and AGFG2.A nonsignificant increase trend was observed for MINK1 (Fig. 2A).In BON-1 cell lines, a significant increase of MINK1 and a nonsignificant trend in CRY2 and AGFG2 was observed upon miR-17-5p inhibition (Fig. 2A).Similarly, inhibition of miR-19a-3p in H727 cells lead to a significant increase of ALG2 and a slight nonsignificant trend of REEP3 levels (Fig. 2A).In BON-1 cell lines, a nonsignificant reduction of miR-19a-3p was observed  [25] in which the eight prognostic miRNAs are involved based on an ORA of the 71 target genes that significantly correlate with these miRNAs.
upon inhibitor transfection which led to a significant ALG2 increase and a nonsignificant increase trend for REEP3 (Fig. 2A).In addition, in silico validation of the 28 prognostic target genes was performed using publicly available data.Similar positive and negative correlations between the miRNAs and their target genes as those documented in our NEN cohort were observed in seven TCGA cohorts, further supporting the miRNA/target gene functional interactions (Fig. 2B).

An 8-miRNA signature stratifies prognosis of NEN patients
Given the biological and independent prognostic relevance of the eight miRNAs, we defined a signature based on their expression profile.Hierarchical clustering of NEN patients based on the expression levels of the eight miRNAs in the signature permitted the identification of three distinct expression clusters grouped by survival, with 5-year OS rates of 80% (good prognosis, N = 35), 66% (intermediate prognosis, N = 27) and 36% (poor prognosis, N = 22; P < 0.001; Fig. 3).Moreover, Cox multivariate analysis showed that prognostic clusters defined by the 8-miRNA signature were significantly associated with survival (P < 0.05), independently of other relevant clinicopathological variables in NENs such as age, gender, primary tumour location, tumour differentiation and stage (Table 2).

Expression of the eight prognostic miRNAs is regulated at the epigenetic level
To test whether dysregulation of the 8-miRNA signature could be mediated by upstream alterations in CpG methylation, we performed DNA methylation analyses on NEN samples (N = 30; Table S1) and analysed the correlation between miRNA expression and the methylation levels of CpG sites within 50 kbp of the TSS of each miRNA.Interestingly, the most significantly correlated (FDR < 0.05) CpG sites were Table 1.List of predicted target genes of the 8-miRNA signature with significant prognostic impact.Cox univariate regression model outcomes using OS and gene expression data of predicted target genes as continuous variables in 63 NENs are shown.Hazard ratios (HR) and P-values are shown as well as the adjusted statistical significance (FDR).Twenty-eight predicted target genes had significant prognostic impact (P < 0.05): 12 were associated with a poorer prognosis and 16 were associated with better outcome.Overall, negatively correlated genes from our study were also negatively correlated in the TCGA samples (in green), and directly correlated genes from our study were also positively correlated in the TCGA (red).OV, ovarian serous cystadenocarcinoma; BRCA, breast invasive carcinoma; GBM, glioblastoma; KIRC, kidney renal clear cell carcinoma; LUSC, lung squamous cell carcinoma; COAD, colorectal adenocarcinoma; MRT, malignant rhabdoid tumor.

Gene
located within the MIR17HG gene (oncomiR-1), which encodes five of the eight miRNAs of the signature (Table S6   Moreover, these methylation differences across prognostic clusters were not observed when we analysed paired healthy tissue (Fig. 4).Accordingly, a robust negative correlation was observed between MIR17HG gene expression, and its promoter methylation levels in 8 TCGA cohorts which further supports the epigenetic regulation of oncomiR-1 (Fig. S5).Similar observations were made for miR-210-3p and two regulatory CpG sites that are involved in its epigenetic regulation (Fig. S6).

Discussion
Relevant advances in the molecular characterisation of NENs have put the spotlight on epigenetic mechanisms as decisive drivers of these tumours [26,27].Although methylation marks and histone modifications are the best known mechanisms of epigenetic regulation, miR-NAs are also key molecules involved in the regulation of gene expression [28].Our study reports a deep characterisation of the expression profile of 84 cancer-related miRNAs in a cohort of 85 patients with NENs.To this end, we identified for the first time a signature based on the expression profile of eight miRNAs (miR-17-5p, -18-5p, -19a-3p, -20a-5p, -20b-5p, -92a-3p, -203a-3p, -210-3p), that is significantly associated with patient prognosis, independent of other well-established prognostic factors.Specifically, the 8-miRNA signature permitted the stratification of three distinct prognostic groups, with 5-year OS rates of 80%, 66% and 36%, respectively (P = 0.046).Additional gene expression and DNA methylation analysis for further down-and upstream integrative characterisation of dysregulated miRNAs identified key genes and crucial regulatory mechanisms potentially involved in the prognostic outcomes of NEN patients.Interestingly, five of the eight miRNAs of the signature (miR-17-5p, -18a-5p, -19a-3p, -20a-5p and -92a-3p) are members of the miR-17-92 cluster (also known as OncomiR-1), localised in chromosome 13, which plays a critical role in cancer [29,30].Members of this cluster and paralogs (miR-20b-5p) have been found to be involved in angiogenesis, proliferation, autophagy, epithelial-mesenchymal transition, invasiveness, migration, cell cycle and apoptotic processes, and are overexpressed in many haematological and solid tumours.They have also been related to biological aggressiveness and chemotherapy resistance [31].The downstream molecular mechanisms of miR-17-92 are nevertheless poorly understood, although some validated target genes are related to the inhibition of EGR (early growth response) [32,33], MYC oncogene or E2F family transcription factors [31], as well as other signalling pathways such as JAK-STAT, RAS/ MAPK, PI3K/Akt/mTOR or TGFβ [31,34].
Knowledge of the role of oncomiR-1 in NENs is limited [31,35].Coordinated overexpression of miR-17, -20 and -92 has been reported in pancreatic acinar and endocrine tumours [36], as well as in angiogenic islets of the RT2 murine model [37].However, the association of miR-17-92 with clinical parameters has not been fully explored.Only miR-19 overexpression has been related with poorer survival in patients with lung NENs [38].Our study thus represents the first report of the role of the miR-17-92 cluster as an oncogene, associated with a poorer prognosis in NEN patients.Similarly, little is known about the role of miR-203a in NENs besides its tumour suppressor role in Merkel cell carcinomas [39] and its association with cell growth, migration and invasiveness in other tumours [40,41].In the case of miR-210, Zimmermann et al. [42] described the overexpression of this miRNA in metastasis vs. primary GEP NENs and in samples with higher proliferative rates.Moreover, its role in angiogenesis is also well-known: miR-210 levels increases in hypoxia and induce HIF1A and VEGF expression by inhibiting VEGF's negative regulators or through the JAK/STAT pathway [43].
We also identified in our study a set of 71 target genes, the expression levels of which were significantly correlated (P < 0.05) with their respective miRNA.Forty-one target genes were inversely correlated (r < 0), suggesting that their corresponding miRNAs regulate their expression through 3 0 UTR mRNA binding and degradation.On the contrary, 30 target genes were directly correlated (r > 0), which indicates alternative regulation mechanisms.MiRNA/target gene associations were also validated in silico using TCGA databases from several solid tumours and NEN cell lines in vitro, where we observed upregulation of targets upon miRNA depletion.ORA revealed that these putative targets and, consequently, the miRNAs that regulate them are involved in classic oncogenic processes described in the hallmarks of cancer and in successfully targeted pathways in NEN patients, namely KEGG PI3K-Akt.Interestingly, epigenetic dysregulation through DNA methylation and miRNA expression modulation could represent a novel mechanism of PI3K-Akt-mTOR pathway alterations, which harbours mutations in about 15% of pancreatic NETs and has proven to be a successful target in a wide spectrum of NETs as demonstrated by the mTOR inhibitor Everolimus [44].This supports the biological relevance of these miRNAs and suggest that miRNA dysregulation could play an important role in driving patients' clinical course.
Finally, we reported here that miR-17-92 cluster expression levels are potentially regulated by DNA CpG methylation, which leads to three clusters of patients with different miRNA expression and prognosis.This is in accordance with recent studies pointing to the relevance of epigenetic alterations in NEN development and progression [10,11,26,27,45].Interestingly, several authors have shown that the miR-17-92 cluster expression can be epigenetically regulated in the context of alveolar differentiation [46][47][48].Moreover, miR-17-92 seems to be progressively repressed during lung differentiation, with several studies reporting that miR-17-92 cluster overexpression inhibits differentiation, which is consistent with what we have observed in NENs.However, further in vitro studies are required to better understand the role of miR-17-92 in dedifferentiation [49,50] and tumour progression, as such efforts could potentially lead to novel therapeutics approaches.
We acknowledge some limitations to our study, including the fact that our cohort was retrospective and heterogeneous in terms of primary tumour site, grade and stage at diagnosis.Nevertheless, the sample size here was not negligible considering the low incidence of the disease, and the patient cohort was wellcharacterised from a clinical and pathological perspective.Furthermore, the results obtained were robust, were validated in in vitro and in silico studies and are biologically plausible as confirmed by integrative transcriptomic and methylomic analyses.Our findings are also supported by preclinical and clinical data published in the scientific literature.Nevertheless, external validation in a prospective independent cohort is warranted to confirm the clinical potential of the described miRNA signature.

Conclusions
In summary, NENs are mutationally quiet and increasing evidence suggests that epigenetic changes play a major role in NEN development and progression.Our work describes, for the first time, a signature based on the expression levels of eight miRNAs (miR-17-5p, -18-5p, -19a-3p, -20a-5p, -20b-5p, -92a-3p, -203a-3p, -210-3p) able to predict survival of patients with GEP and lung NENs, independent of other wellestablished prognostic factors.Integrative analyses of miRNA profiles, transcriptomics and methylomics identified key target genes and dysregulated pathways driving prognosis of NEN patients, along with critical CpG sites involved in the epigenetic regulation of these eight miRNAs.Although our findings require external validation in independent series of patients and further functional analyses to better understand the mechanisms by which the eight miRNAs exert their actions in these malignancies, the miRNA signature identified here constitutes a promising tool to facilitate a more accurate prognostic stratification of NEN patients.Furthermore, we have demonstrated that these miR-NAs regulate oncogenic pathways and key cellular processes that could potentially be explored in the future for therapeutic purposes.prognostic impact.Despite not being statistically significant, miR-19a-3p, miR-20b-5p and miR-92a-3p show a clear trend towards a poorer prognosis with higher miRNA levels.(B) Cox univariate regression model data for the 8-miRNA signature using qPCR expression data as a continuous variable and overall survival.Hazard ratios (HR) and p-values are shown.Seven of the eight miRNAs showed a significant prognostic impact on OS, which confirmed our earlier observations with the PCR array.p < 0.05 was considered significant.Fig. S3.Functional interactions of 28 prognostic target genes.A string functional network is shown based on the 28 target genes that were significantly associated with OS (p < 0.05).Line width increases with higher interaction robustness.Interactions with medium confidence or higher are shown (0.400).Seven genes (E2F3, CCND2, E2F8, KIF23, DUT, PRC1 and RHOB) are functionally related (PPI enrichment p-value: 0.0517) despite being regulated by different miRNAs.These genes are involved in mitotic spindle, G2/M checkpoint, apoptosis and cell cycle and proliferation.Fig. S4.Basal miR-17-5p and miR-19a-3p expression levels in H727 and BON-1 cell lines are similar to good prognosis NEN patients.MiR-17-5p and miR-19a-3p expression levels in basal H727 and BON-1 NET cell lines and in 40 NEN patients from our study which had available RNA were assessed by RT-PCR.No significant differences in miR-17-5p and miR-19a-3p levels were observed between NET cell lines and good prognosis NEN patients.Conversely, poor prognosis patients showed a trend towards higher levels of these two miRNAs than NET cell lines (miR-17-5p: H727 vs NENs, p = 0.063; BON-1 vs NENs, p = 0.066; miR-19a-3p: H727 vs NENs, p = 0.085; BON-1 vs NENs, p = 0.088).This is in accordance with the fact that these two cell lines are derived from well-differentiated NETs, whereas our patient cohort included a wider spectrum of NENs including high grade poorly differentiated NECs.MiRNA expression levels are expressed as 2 −ΔΔCt on the y-axis using H727 cell line as reference.Student's t-tests between the cell lines and the patients from the different prognostic clusters were performed to evaluate differences.The methylation levels (β-value) of the significantly correlated CpG sites across the 8-miRNA prognostic clusters were assessed in both healthy and tumoural samples.ANOVA or Kruskall-Wallis test were performed between clusters as appropriate and pairwise comparisons were performed using Tukey's range or Dunn's Test, respectively.Promoter-like 2 CpG site methylation levels show increasing levels in the three clusters in healthy tissue.Conversely, decreasing levels are observed in the tumoural tissue (p = 0.0094).Similar results were obtained for Promoter-like 3 CpG site, with increasing methylation across clusters in healthy tissue, whereas decreasing methylation is observed in tumoural tissue (p = 0.0049).Table S1.Clinico-pathological features of the study population.The clinic-pathological characteristic of our study population depicted by the different omic approaches.Table S2.Technical validation of the eight prognostic miRNAs using qRT-PCR.qRT-PCR of the eight selected miRNAs was performed on 39 NEN patients in our cohort.A Spearman correlation analysis between PCR array and RT-PCR data (N = 39) for the eight prognostic miRNAs was performed.The r value of the Spearman correlation coefficient and the p-value are shown.With the exception of miR-20b-5p, a strong (r > 0.4) and significant (p < 0.05) correlation was observed between the miRNA array and qRT-PCR miRNA expression levels for all miRNAs of the signature.p < 0.05 was considered significant.Table S3.In silico predicted target genes of each selected prognostic miRNA.Genes regulated by more than one miRNA are in bold font.
Table S4.List of predicted target genes of the eight selected miRNAs with significant expression correlation with their respective miRNA.Spearman correlation coefficients (r) using PCR array miRNA and gene expression data of 62 NENs are shown as well as pvalues and corresponding miRNA.Predicted target genes are ordered by decreasing r.Only predicted target genes with significant correlation are listed in the table (N = 71).Forty-one of the 71 target genes were inversely correlated (r < 0) in accordance with the classical miRNA-target gene interaction, while 30 showed a direct correlation (r > 0), which suggests an alternative regulation.p < 0.05 was considered significant.* 'Genes that correlate with several miRNAs.Table S5.Overrepresentation analysis (ORA) using the 71 target genes of the eight selected miRNAs and the Hallmarks and KEGG gene sets revealed relevant signalling pathways in NENs.Overrepresentation analysis was performed with Enrichr using the 71 target genes that were significantly correlated with their miRNAs, and the gene sets from Hallmarks and KEGG.Gene sets with FDR <0.25 were considered significant.Next, significant pathways, the target genes involved in each pathway and their miRNAs were related (based on their biological functions) to the Hallmarks of cancer which was recently updated by Hanahan and colleagues.Table S6.Significantly correlated CpG sites present differentially methylated levels in NEN tumours.

Fig. 1 .
Fig. 1.Eight miRNAs are associated with poor prognosis and key oncogenic pathways in NENs.(A) Prognostic impact of miRNAs in NENs.P < 0.05 was considered significant.CI, confidence interval.(B) The sector diagram represents the hallmarks of cancer[25] in which the eight prognostic miRNAs are involved based on an ORA of the 71 target genes that significantly correlate with these miRNAs.

Fig. 2 .
Fig.2.Validation of target genes of the eight prognostic miRNAs.(A) MiR-17-5p and miR-19a-3p regulate the expression of their predicted target genes in NENs cell lines.MiR-17-5p inhibitor was transfected for 48 into H727 and BON-1 cell lines using lipofectamine (n = 3).Treated BON-1 and H727 cells showed a reduced expression of miR-17-5p at 48 h upon transfection.Accordingly, expression of candidate target genes CRY2 and AGFG2 was significantly increased in H727 cells and expression of MINK1 was significantly increased in BON-1 cells.A nonsignificant increase of MINK1 and of CRY2 and AGFG2 was observed in H727 and BON-1 cell lines, respectively.Similarly, miR-19a-3p inhibitor was transfected for 48 h into H727 and BON-1 cell lines using lipofectamine (n = 3).Treated BON-1 and H727 cells showed a reduced expression of miR-19a-3p at 48 h upon transfection.Accordingly, the expression of candidate target genes ALG2 and REEP3 was increased in the BON-1 and H727 cell line, although it did not reach statistical significance for REEP3.No differences were observed in RHOB.MiRNA or target gene expression levels are expressed as 2 ÀΔΔCt on the y-axis using the negative control at each time as reference.The dotted line represents the NC (Negative Control) value.Mean AE SEM (Standard Error of the Mean) is shown.Student's t-tests between the negative control and inhibitor-transfected cells at each time point were performed to evaluate differences, with P-values <0.05 considered statistically significant (*P < 0.05; ***P < 0.001; ns, nonsignificant).(B) Correlations between miRNAs and target gene expression were validated in silico in different solid tumours.We performed Pearson's correlation analyses for each miRNA/target gene pair in seven solid tumour TCGA cohorts.Heatmaps represent the miRNA-target gene coefficient of our NEN cohort and of the seven TCGA cohorts.Inversely (left panel) and directly (right panel) correlated miRNA-target genes are shown in the two heatmaps.The seven TCGA databases and our NEN cohort are shown in columns and miRNA-target genes are shown in rows.Individual Pearson's correlation coefficient values are colour-coded, ranging from red (positive correlation; max = 1) to green (negative correlation; min = −1).Overall, negatively correlated genes from our study were also negatively correlated in the TCGA samples (in green), and directly correlated genes from our study were also positively correlated in the TCGA (red).OV, ovarian serous cystadenocarcinoma; BRCA, breast invasive carcinoma; GBM, glioblastoma; KIRC, kidney renal clear cell carcinoma; LUSC, lung squamous cell carcinoma; COAD, colorectal adenocarcinoma; MRT, malignant rhabdoid tumor.

Fig. 3 .
Fig.3.8-miRNA signature defines prognostic groups of NEN patients.(A) We conducted a heatmap using the expression profile of the eight miRNAs in tumour samples of the NEN patient cohort.Red and blue colours represent high and low expression, respectively.Hierarchical clustering was performed for patients (columns) and miRNAs (rows).The upper tree shows three clusters defined by the 8-miRNA signature expression profile and survival analysis: poor (left), intermediate (middle) and good prognosis (right).Clustering based on the expression of miR-NAs is shown in the tree on the right.Patients were annotated with primary tumour site, grade, stage, differentiation and the presence or absence of a functional syndrome.diff, differentiated; G1, grade 1; G2, grade 2; G3, grade 3; (B) Kaplan-Meier curves of NEN patient prognostic clusters according to the 8-miRNA expression profile identified three distinct prognostic groups, with 5-year survival rates of 80% (good prognosis, depicted in blue), 66% (intermediate prognosis, depicted in grey) and 36% (poor prognosis, depicted in red) (P = 0.00014).OS in months is represented on the x-axis and survival probability is shown on the y-axis.Values of P < 0.05 were considered statistically significant.

Fig. 4 .
Fig. 4. 8-miRNA signature is epigenetically regulated by CpG methylation.Methylation levels (β-value) of CpG sites which significantly correlated with miRNA expression are shown across the prognostic clusters (defined by miRNA expression profile) in both healthy (n = 10) and tumoural tissue (n = 30).ANOVA or Kruskall-Wallis test were performed between clusters as appropriate and pairwise comparisons were performed using Tukey's range or Dunn's Test, respectively.(A) The Proximal (Prox.)enhancer-like 1 CpG site showed no methylation differences in healthy tissue, whereas decreasing methylation was observed for the 'Good prognosis' to 'Poor prognosis' clusters in tumoural tissue (P = 0.048).(B) The Proximal enhancer-like 2 CpG site showed stable/slightly increasing levels for 'Good prognosis' to 'Poor prognosis' clusters in paired healthy tissue.Conversely, a significant decrease in methylation was observed from 'Good prognosis' to 'Poor prognosis' cluster in tumoural tissue (P = 0.003).(C) Finally, the Promoter-like 1 CpG site showed no differences in the normal tissue, and a slight decrease, particularly in the poor prognosis cluster, in the tumoural tissue (P = 0.019).P-values <0.05 considered statistically significant (*P < 0.05; **P < 0.01; ns, nonsignificant).

Fig. S6 .
Significantly correlated CpG sites present differentially methylated levels in NEN tumours.(A) Spearman correlation analysis was performed using PCR-array miRNA expression and CpG methylation levels of sites within 50 kbp of miR-210-3p (regulatory region) in 30 NENs.The Spearman correlation coefficients (r) are shown as well as p-values and adjusted statistical significance (FDR).Promoter-like 2 and 3 CpG sites (Cg01325426; Chr11:536330-536332 and Cg22174486; Chr11:536816-536818, respectively) are inversely correlated (both FDR =0.07523) with miR-210-3p in accordance with the classic CpG methylation regulatory mechanism.These sites are in a CpG island in a promoter-like region within the HRAS gene.(B)

Table 2 .
Multivariate analysis of prognosis clusters defined by the 8-miRNA signature.Cox multivariate proportional hazards model using OS and the clusters defined by the 8-miRNA signature as well as other possible confounding variables in 84 NENs is shown.Hazard ratios (HR) and P-values are shown.The clusters defined by the signature were associated with a poorer outcome in NEN patients independently of other relevant variables such as gender, age, primary tumour location, differentiation and stage.CI, confidence interval; Ref, reference.
Mean AE SEM (Standard Error of the Mean) is shown.P-values <0.05 considered statistically significant.Glioblastoma Multiforme, GBM; Kidney Renal Clear Cell Carcinoma, KIRC; Lung Squamous Cell Carcinoma, SLCC; Ovarian.Serous Cystadenocarcinoma, OV).As we observed in our own cohort, MIR17HG promoter methylation was inversely correlated with MIR17HG expression.