A data-driven network model of primary myelofibrosis: transcriptional and post-transcriptional alterations in CD34+ cells

microRNAs (miRNAs) are relevant in the pathogenesis of primary myelofibrosis (PMF) but our understanding is limited to specific target genes and the overall systemic scenario islacking. By both knowledge-based and ab initio approaches for comparative analysis of CD34+ cells of PMF patients and healthy controls, we identified the deregulated pathways involving miRNAs and genes and new transcriptional and post-transcriptional regulatory circuits in PMF cells. These converge in a unique and integrated cellular process, in which the role of specific miRNAs is to wire, co-regulate and allow a fine crosstalk between the involved processes. The PMF pathway includes Akt signaling, linked to Rho GTPases, CDC42, PLD2, PTEN crosstalk with the hypoxia response and Calcium-linked cellular processes connected to cyclic AMP signaling. Nested on the depicted transcriptional scenario, predicted circuits are reported, opening new hypotheses. Links between miRNAs (miR-106a-5p, miR-20b-5p, miR-20a-5p, miR-17-5p, miR-19b-3p and let-7d-5p) and key transcription factors (MYCN, ATF, CEBPA, REL, IRF and FOXJ2) and their common target genes tantalizingly suggest new path to approach the disease. The study provides a global overview of transcriptional and post-transcriptional deregulations in PMF, and, unifying consolidated and predicted data, could be helpful to identify new combinatorial therapeutic strategy. Interactive PMF network model: http://compgen.bio.unipd.it/pmf-net/.


INTRODUCTION
Primary myelofibrosis (PMF) is a chronic myeloproliferative neoplasm (MPN) that, with essential thrombocythemia and polycythemia vera, constitutes a heterogeneous group of Philadelphia-negative clonal hematopoietic stem cell (HSC) disorders associated with overproduction of mature myeloid cells. Typical traits of myelofibrosis are an increased proliferation of megakaryocytes, a deposition of fibrosis in the bone marrow, an abnormal stem cell trafficking and an extramedullary hematopoiesis (myeloid metaplasia). Moreover, PMF is associated with an increased risk of thrombosis and/or hemorrhage and a propensity to develop acute myeloid leukemia. 1,2 Important progresses in molecular characterization of MPNs pathogenesis have been done in the last years. Specifically, the discoveries of somatic mutations in JAK2, MPL and CALR genes have improved patients' stratification and molecular characterization highlighting the role of Jak-STAT signaling in MPN pathogenesis. However, several evidences indicate that these mutations are not sufficient for disease initiation and progression and that the phenotypes of the disease are highly heterogeneous, suggesting that other unknown genetic or epigenetic factors might be involved [3][4][5][6][7][8][9][10][11][12][13][14][15] and also that the mutation order could matter: 16 driver mutations can precede or follow additional somatic mutations in other myeloid genes. Interestingly, recently, it has been demonstrated that most of these genes are associated with age-related clonal hematopoiesis in normal elderly subjects, suggesting that pre-malignant clones may be present many years before disease develops and are required, but insufficient, to result in disease. 17,18 microRNAs (miRNAs) have an important role in the regulation of hematopoiesis. [19][20][21][22] Our group demonstrated that miR-16-2 contributes to the expansion of erythroid lineage in polycythemia vera 23 and we showed that miR-155-5p is pathogenically associated to MK hyperplasia in PMF through JARID2 downregulation. 24 Moreover, we recently characterized miRNA and microRNA offset RNAs (moRNA) expression in SET2 cells 25 and in CD34+ stem cells using massive small RNA-seq. In the latter study, we observed specificities in small RNAs expression of PMF cells. 26 Although these findings are supportive of the involvement of miRNAs in PMF pathobiology, our understanding of miRNAs involvement in MPNs is still limited. 27 A deeper characterization of the miRNA-mediated pathogenesis processes 19,28 would be highly desirable in order to identify suitable concurrently targetable pathways amenable to therapeutic intervention. 28,29 Thereby, the aim of this study was to obtain an informative, composite and interactive data-driven picture of pathways and circuits deregulated in PMF. We used a composite pipeline exploiting both 'knowledge-based' and ab initio approaches to discover mixed transcriptional and post-transcriptional deregulated circuits in PMF. This strategy allowed us to describe an unforeseen picture of miRNA and gene regulatory circuits linked to abnormal cellular functions and pathways.

Patient selection and expression data
For this study, we considered miRNA and gene expression data (series GSE41812 and GSE53482) obtained by analyzing CD34+ cells of 42 patients with a diagnosis of PMF and 16 peripheral blood (CTR PB) and 15 bone marrow (CTR BM) samples from normal donors. Primary myelofibrosis patients were in a typical fibrotic stage of the disease according to the WHO, and were molecularly characterized (JAK2V617F, MPLW51 and CALR mutations). See Norfo et al. 24 for additional clinical information about patients.
Gene expression profiling, obtained by Affymetrix HGU219, and annotated by using custom CDF, 30 HGU219_Hs_ENTREZG version 14, was normalized by RMA procedure. A log 2 expression matrix with 18 654 genes and 73 samples was obtained.
MiRNA expression profiling, obtained by Affymetrix miRNA 2.0 microarray, was analyzed by using a modified RMA procedure: the background subtraction was performed considering all species and controls of the chip, while normalization and summarization were calculated only on human probe sets, using a custom-designed CDF. MiRNA names were translated to updated miRBase names filtering out deprecated miRNAs, using the custom annotation package AffyIDs. A log 2 expression matrix of 584 miRNAs and 73 samples was obtained after filtering miRNAs weakly expressed (lowest 25% of average expression level) and/or poorly variable across samples (Shannon Entropyo1.4).

Identification of negatively correlated predicted targets of miRNAs
We identified the predicted regulatory relations supported by our expression data by integrating target predictions and expression profiles. Target predictions were computed with TargetScan-62 algorithm according sequence similarity with conservation (threshold at 0.8).
Pairwise Pearson correlations between miRNA and predicted target gene expression profiles were calculated. The relationships supported by significant anti-correlation (r o−0.5 with FDRo0.1) were considered for the following analyses.
The miRNA/target interactions involving genes annotated in KEGG pathways have been used in the following Micrographite pipeline.

Pathways-derived circuits
Two comparisons were considered for the analysis: PMF vs CTR BM and PMF vs CTR PB. mRNA and miRNA relations involved in PMF were compared to CTR PB or CTR BM and have been identified using the Micrographite pathway analysis. 31 In the following we will report a brief description of the main steps of the method; for details refers to Calura et al. 31 KEGG 32 pathways as stored in graphite Bioconductor package 33 were enriched with validated and predicted miRNA-target interactions. We consider 'validated interactions' those interactions extracted from miRTarBase 34 and miRecords 35 for which reporter assays have been provided in literature, while the 'predicted/supported interactions' are the in silico predicted miRNA-mRNA couples filtered by anticorrelation of expressions as discussed above. In this step, only miRNAs targeting genes already annotated in the cell pathways have been added. Following the Micrographite pipeline, the pathways resulted to be significant at the whole pathway level (q-value for the mean and variance testso =0.1) were used for the second step of the analysis that identifies, within each significant pathway, the subportion of the pathway (hereafter called path) mostly associated to the phenotype. The merge of the top 10 paths is called meta-pathway. The meta-pathway is then re-analyzed in order to find the top scored paths. Finally, the paths that, ranked by scores, belong to upper quartile are considered as the most involved and reported. The Micrographite analysis has been performed independently for each comparison (PMF vs CTR PB; PMF vs CTR BM) and the two final networks have been merged into a single network.
Mixed miRNA-TF-gene circuit identification In parallel, the expression data of genes that are present in the KEGG pathways were used with the miRNAs profile as input of Magia 2 . Magia 2 (http://gencomp.bio.unipd.it/magia2) 36 reconstructs circuits involving positive and negative correlations between miRNAs, transcription factors (TFs) and target genes. It uses experimentally validated TF-miRNA relations obtained by mirGen2.0 and TransmiR, and TF-gene relations obtained by the UCSC human genome annotation.
Target predictions were computed with TargetScan (score threshold at 0.8, top 25% of the predictions distribution) and Pearson's correlation coefficient was used as correlations measure (threshold at 0.5).

Pathways visualization
Network visualization and annotation have been performed using Cytoscape 37 (version 3.1).

Validation of differentially expressed sRNAs
We performed individual miRNAs assay by Taqman quantitative real-time PCR for quantification of abnormally expressed miRNAs in PMF and control granulocytes and plasma. All samples were processed within 4 h after collection. Plasma, collected by peripheral blood centrifugation in tubes containing EDTA (4°C, 3000 revolutions per minute for 15 min, then 6000 revolutions per minute for 5 min), was stored at -80°C in RNase/DNase-free Eppendorf tubes. RNA was isolated and eluted in 20 μl of RNase-free water using a miRNeasy Serum/Plasma Kit (Qiagen). For miRNA expression assays 8 μl eluted RNA was used. 38 cDNA was synthesized from total RNA using miRNA-specific RT primers contained in the TaqMan microRNA Human Assays (Applied Biosystems, Foster City, CA, USA). Briefly, single-stranded cDNA was synthesized from 10 ng total RNA in 15-μl reaction volume with the High-Capacity cDNA Archive Kit (Applied Biosystems) using 1 mM deoxyribonucleoside triphosphates, 50 U Multiscribe reverse transcriptase, 3.8 U RNase Inhibitor and 50 nM of miR-specific RT primers. The reaction was incubated at 16°C for 30 min followed by 30 min at 42°C, and inactivation at 85°C for 5 min. Each generated cDNA was amplified by quantitative realtime-PCR with sequence-specific primers from the TaqMan microRNA Assays on an ABI Prism 7300 real-time PCR system (Applied Biosystems). PCR reactions included 10 μl 2 × Universal PCR Master Mix (No AmpErase UNG), 2 μl each 10 × TaqMan MicroRNA Assay Mix and 1.5 μl reverse-transcribed product; they were incubated in a 96-well plate at 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. Expression variations were calculated using the RQ method and a t-test P-value of 0.05 was used as threshold to identify differentially expressed elements.

RESULTS AND DISCUSSION
A systems biology pipeline identifies deregulated signaling pathways in PMF patients We considered matched expression profiling data of 584 miRNAs and 18 654 genes in 73 CD34+ cells sampled from 42 PMF patients, and 31 sampled from normal subjects (15 from bone marrow (CTR BM) and 16 from peripheral blood (CTR PB). The study grounds on comparative analysis of CD34+ cells of PMF patients compared with control CD34+ from normal donors. In a previous study, specificities of both gene and miRNA expression profiles were observed comparing peripheral blood and bone marrow control samples, 24 confirming known differences between these populations. Thus, peripheral blood and bone marrow control samples were considered separately, for parallel comparisons with PMF samples, to avoid biases due to control samples inhomogeneity, and to obtain more informative results. Then, a union of the comparisons gave a global picture. The flowchart in Figure 1 provides a summary of the procedure.
We obtained 250 interactions in the PMF vs CTR BM (54 miRNAs and 236 gene targets involved) and 442 in the PMF vs CTR PB (43 miRNAs and 428 genes) using MAGIA 2 . Then, the miRNA and gene expression profiles of PMF patients have been independently compared with CTR BM and CTR PB health samples using Micrographite. Supplementary Tables 1-8 contain all the intermediate data of the two analyses whereas a brief comparison of results is summarized in Supplementary Figures 1 and 2 Mixed miRNA-TF circuits connect to altered pathways miRNAs post-transcriptionally regulate gene expression, also by alteration of target mRNAs stability. Moreover, TFs modulate the expression of miRNAs and of target genes at transcriptional level, creating a regulatory network. Magia 2 analysis identified mixed regulatory circuits composed of miRNAs, TFs and target genes, considering in parallel expression profiles in PMF and CTR BM samples and in PMF and CTR PB samples. For each type of control samples used, Supplementary Tables 9-12 report more details on observed correlation values, on the numbers of identified circuits, on predicted interactions and on involved genes and miRNAs. Inferred interplay of miRNAs and TFs in gene/transcripts expression regulation may comprise two types of mixed regulatory circuits: (i) a TF regulates both a miRNA and a gene promoter, and the same miRNA targets the gene mRNA exerting post-transcriptional regulation; (ii) a miRNA post-transcriptionally regulates both a gene encoding a TF and another gene, that on turn is a transcriptional target of the TF. For both mixed circuits topologies, different combinations of positive and negative pairwise correlations should correspond to activatory and inhibitory interactions giving rise to coherent or incoherent feed forward loops. 36,39 Since circuits based on strongest interactions are particularly trustable and may highlight key regulatory relations, we focused on circuits that had absolute correlation at least 0.5 for all the three types of interactions (miRNA-target gene, TF-miRNA, TF-target-gene). PMF vs CTR PB comparison identified 169 interactions, involving 63 nodes, 6 miRNAs and 57 genes. No significant circuits were selected for PMF vs CTR BM comparison.
Pathway annotations are progressively manually updated but, even if enriched with miRNA validated interactions, they are still incomplete. Thus, we enriched the pathway-based networks with nodes and relations coming from selected Magia 2 circuits, that shared with the first network at least one gene or one miRNA. In this way, we obtained the network, consisting of 620 edges and 220 nodes (of which 34 are miRNAs; Table 1); they are reported, with different details and annotation, in Supplementary Figure 3 and are included in the yellow area of the circuit in Figure 2a. Table 1 summarizes miRNAs and TFs involvement. Magia 2 circuits showed that six miRNAs and six TFs seem to play critical roles in the regulation of transcription in PMF ( Figure 2a and Table 1). Almost all of the miRNAs are components of miR-17 family (miR-106a-5p, miR-20b-5p, miR-20a-5p, miR-17-5p and miR-19b-3p). The upregulation of the latter miRNA in PMF CD34+ cells has been previously validated. 24 Only let-7d-5p is a component of another family of miRNAs, the let-7 family, which is present also in Micrographite results. The TFs (MYCN, ATF, CEBPA, REL, IRF and FOXJ2) are important regulators of several target-genes. Relevant contribution in PMF vs CTR PB comparison is done by MYCN gene, in particular the upregulation of the miR-17 family components and the downregulation of let-7d-5p. The bridge elements between the pathway network and predicted circuits are two miRNAs (miR-17-5p and miR-20a-5p) and four genes (ATF, CREB5, NRAS and ARHGEF7). Different genetic lesions characterizing PMF patients (28 JAKV617F, 7 CALR, 4 MPLW51 and five triple negative, 3N, for the former mutations) included in this study are representative of the molecular alterations described in PMF. Given such patient heterogeneity, we checked if some of the genes/miRNAs modulated in PMF were deregulated only, or mainly, in a subgroup of patients. As shown in Figure 2b    Even if the multistep approach we used for the PMF network reconstruction was specifically designed to go beyond the differential expression of single genes or miRNAs, giving a more informative picture of the impact of miRNAs deregulated in PMF to pathways and to mixed connected regulatory circuits it is worth noting that we obtained confirmation of the differential expression in PMF samples compared to controls for 13 miRNA and 5 genes included in the network, giving further support to our results. As detailed in Table 2, in addition to the confirmation of differential expression in CD34+ previously obtained by array and quantitative real-time-PCR 24 and RNA-seq, 26 we also conducted new experiments in granulocytes and plasma (see Methods) to obtain new insights of direct translational relevance. These confirmed a significant upregulation (P-valueo0.05) of miR-195-5p, miR-106a-5p, miR-145-5p, miR-17-5p, miR-185-5p, miR-195-5p in PMF patients granulocytes (GR; considering 50 PMF patients and 10 healthy controls), miR-19b-3p both in PMF granulocytes and plasma (25 PMF patients and 6 controls) and of miR-29b-3p in plasma.
Candidate transcriptional deregulated pathways of PMF Figure 2c summarizes the main alterations of the transcript levels observed in PMF patients, involving key pathways and connected circuits, as discussed below in more detail with the support of Figure 3.
Akt signaling. One of most deregulated pathways is the Akt signaling, a part of the Phosphoinositide-3 kinase/protein kinase-B/ mammalian target of Rapamycin (PI3K/Akt/mTOR) pathway ( Figure 3a). As the Jak/STAT signaling path, the Akt signaling is triggered by receptor tyrosine kinases that increase PI3K levels, activating AKT, phosphorylating several targets including mTOR. Several tumor suppressor genes, as PTEN, interact with this path. PI3K/AKT/mTOR is a central signaling module and several other pathways converge on it regulating cell survival and proliferation. 28 It is constitutively activated in MPNs and also in other cancers. 28 Small molecular inhibitors of this pathway are available in clinics, 23 and evidence of activity of an mTOR inhibitor (Everolimus) in myelofibrosis patients has been provided. 40 The contribution of PI3K/mTOR pathway to MPN pathogenesis has received strong experimental support. 28,41,42 Rho GTPases, CDC42, PLD2 and PTEN. The PMF network shows deregulation of Rho GTPases (Figure 3b), intracellular signal transducers, members of the RAC subfamily of Rho GTPases. Deregulation of GEFs (guanine nucleotide exchange factors), recharging Rho GTPases and determining their spatial-temporal activity is also displayed. Interaction of stem cell factors and adhesion molecules with their receptors (c-Kit, CXCR4 and Integrins) expressed by the HSCs triggers the activation of GEFs. 43 Rho GTPases could affect also cancer progression, 44 both as oncogene and oncosuppressor. Specifically, CDC42 and RAC2 control homeostasis of blood cell production by hematopoietic stem/progenitor (HSC/P) cells balancing HSC/P retention within the bone marrow and migration in the blood. 45 CDC42, highly downregulated in PMF, 20 regulates cell polarity and cell cycle progression. Its deletion affects the number of early myeloid progenitors while suppressing erythroid differentiation. 46 CDC2-deficient mice develop a fatal myeloproliferative disorder manifested by significant leukocytosis with neutrophilia, myeloid hyperproliferation and myeloid cell infiltration into distal organs. 46 Downregulation of RhoA activity, as observed in PMF, resulted in increased HSC engraftment and self-renewal. In the network, CDC42 is connected with the upregulated miR-29a-3p, which we previously demonstrated is able to suppress their expression 24 and which is connected also to the phosphatase PTEN, regulated by RhoGTPases. 47 In turn, PTEN, down-regulated in PMF (and linked to several miRNAs including the upregulated miR-19a-3p, miR-29a-3p and miR-21-5p), modulates intracellular levels of PI3K in cells, functioning as a tumor suppressor. The upregulation in PMF of miRNAs interacting with PTEN was previously established. 24 CDC42 is also connected with miR-29b-3p, miR-29c-3p and miR-185-5p, the latter linked also with the downregulated RHOA. Another element markedly deregulated in the PMF network is Phospholipases D, PLD1 and PLD2 (Figure 3b). Phospholipase D and its product PA (Phosphatidic Acid) have been implicated in many physiological and pathological functions, such as cancer cell invasion and metastasis. [48][49][50][51][52] PLD2 is dual activity enzyme with both lipase and GEF activities for the Rac2 and RhoA GTPases. As demonstrated by Mahankali et al, 53 in the context of leukocyte chemotaxis, both activities are triggered after cell stimulation, but then, the accumulation of RAC2-GTP, due to GEF action, negatively feeds back on PLD2-lipase activity, and the JAK3 tyrosine kinase restores lipase activity while ultimately inhibiting GEF activity. 54 This lipase-GEF duality has important biochemical and cellular implications and new PLD inhibitors have been recently described. 52,54 PA is one of the major lipid second messengers that in turn triggers many signaling pathways such as PI3K/Akt/mTOR, sphingosine kinase 1, Raf-1 protein kinase 53 and small GTPases like Ras and Rac. 55 HIF-1a pathway. In connection with previously mentioned pathways, we found evidence of modulation of the HIF-1a hypoxia pathway ( Figure 3c). As known, both growth factors and oncogenic signaling can trigger an HIF-1a response, with the PI3K/Akt and Ras/Raf/MAPK pathways mediating upregulation of HIF-1a by increasing HIF-1a protein synthesis or stability, resulting in a modified transcriptional activity of HIF-1. In our network, the predicted/supported relation of the still poorly characterized miR-1244 (downregulated in PMF) with NOS3 (nitric oxide synthase 3, endothelial cell) is connected with HIF-1a while the upregulated miR-20a-5p connects HIF-1a to PTEN. NOS3 synthesizes nitric oxide, that is involved in several processes such as muscle relaxation, inhibition of platelet aggregation and leukocyte adhesion, with both vasodilatatory and antithrombotic properties. A recent study showed that NOS3 Glu298Asp polymorphism in homozygosis increased the risk of thrombosis in a cohort of patients with essential thrombocythemia or prefibrotic PMF. 56,57 A link between hypoxia and the activation of the HGF/Met axis, through the oxygen sensor HIF-1a in MPNs, has recently been proposed by Boissinot and colleagues. 58 Oxygen tension is as low as 1% in normal bone marrow 59 and the hypoxic condition may result even more stringent because of the bone marrow hypercellularity in myelofibrosis. In accordance with our results, showing activation of the HIF-1a downstream genes ( Figure 3c) and of NF-κB signaling (Figure 3a), the authors proposed that the activation of the HGF/Met axis, due to bone marrow or/and exposure to inflammation cytokines (NF-κB-mediated activation), leads to enhanced survival and proliferation of mutated myeloid progenitors and increased production of various cytokines responsible for inflammation, neo-angiogenesis and fibrosis. Thus, molecules efficient at blocking NF-κB and HIF-mediated activation might represent an efficient therapy for MPNs by simultaneously blocking the cascade of inflammation cytokines, their signaling and the hematopoietic clone.
Calcium and cyclic AMP signaling. The upper part of the network in Figure 2a includes genes, belonging to Calcium signaling path (Figure 3d). More specifically, this part of the network includes genes encoding several subunits of calcium channel voltage-dependent, and other genes, including ATF/CREB subunits activated by cyclic AMP-dependent protein kinases, tightly related to MAPK and Ras signaling pathways, Chemokine signaling and Apoptosis. ATF/CREB genes participate also to PI3K/Akt signaling pathway and to TNF signaling pathway. The Calcium, and c-AMP networks include, as possible modulators, miR-125b-5p (that might target the downregulated PPP1CA and GRIN2A) and miR-133a (that might target CACNA1C).
CALR mutations are found in PMF and essential thrombocythemia patients. Calreticulin is a multifunctional protein that acts as a major Ca 2+ -binding (storage) protein and chaperone in the lumen of the endoplasmic reticulum and intervene in transcription regulation. 60 CALR has been implicated in the development of different cancers and its effect on tumor formation and progression may depend on cell types and clinical stages. 60 CALR is not differentially expressed in  61 known to be hypermethylated in cancer and in high-risk myelodisplastic syndromes; 62 ZNF516 is a zinc-finger protein that interacts with epigenetic regulation protein such as HDAC1 and KDM1A. The upregulated PLXND1 is one of Semaforin receptors, often highly expressed in cancer, 63 that regulates tumor cell survival by suppressing an apoptotic pathway 64 and seems to controls angiogenesis in cancer. 65 let-7d-5p has also correlated expression with MYC with which its shares KLF13 as putative target. CREB connects to miR-20a 5p and to MYC, and in turn to their common targets.
Predicted circuits open new hypotheses All the deregulated pathways contain links to novel in silico predicted interactions that could help in the candidate selection in research. In support of the results consistency, a subset of de novo predicted interactions is known and validated, that is, the upregulated miR-17-5p, miR-106a-5p and miR-20b-5p with the downregulated CDKN1A. CDKN1A is an important intermediate by which p53 blocks proliferation in response to DNA damage, and it inhibits CDK activity, blocking cell cycle progression. The validated post-transcriptional regulatory relations are also linked to MYC (highly upregulated in PMF) transcriptional activity, which is known to stimulate the expression of onco-miR-17 family. These, in turn suppress IRF1, a tumor suppressor TF that is downregulated in PMF and has important supported targets, as MAP3K2, GLIS3 and LRCH1, all downregulated.
Four upregulated miRNAs (miR-17-5p, miR-106a-5p miR-20a-5p and miR-20b-5p) are negatively correlated with the downregulated TF FOXJ2 (Forkhead box protein J2), forming triangular circuits with four genes including the upregulated PDZ and LIM domain 5 (PDLIM5), recently identified as the main signaling molecule that, in the context of augmented AMPK activity and in connection with Rac1-mediated acting cytoskeleton reorganization, 66 regulates cell migration, a key process in MPNs.

CONCLUSIONS
The network inferred in this work is a unique and integrated signaling pathway of PMF, in which the role of miRNA is to wire, co-regulate and allow a fine crosstalk between all the involved processes. Most deregulated pathways include Akt signaling, linked to Rho GTPases, CDC42, PLD2 and PTEN crosstalk with the hypoxia pathway, which highlights the extreme hypoxic status in PMF bone marrow. Novel hints on cyclic AMP signaling involvement jointly connected with Calcium-linked cellular processes were given. Moreover, mixed circuits involving miRNAs (miR-106a-5p, miR-20b-5p, miR-20a-5p, miR-17-5p, miR-19b-3p and let-7d-5p) connections with key TFs (MYCN, ATF, CEBPA, REL, IRF and FOXJ2) and their common target genes were reported, opening new hypotheses. The comprehensive picture of transcriptional deregulation can be used to improve current understanding of the molecular pathways involved in the disease and delineate new combinatorial therapeutic strategies.

CONFLICT OF INTEREST
The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. The authors declare no conflict of interest.