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.

Methods

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 log2 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 log2 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 Entropy<1.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<−0.5 with FDR<0.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 tests<=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 Magia2. Magia2 (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 real-time-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.

Figure 1

Experimental design and analysis procedures flow chart. miRNA and gene expression data in CD34+ cells of PMF patients and healthy controls were considered; CTR BM and CTR PB controls were considered separately in parallel analyses, and results were merged at the final steps. MAGIA integrated analysis outputted the subset of predicted miRNA-target relations supported by expression data; these were used to enrich KEGG pathway-derived miRNA-gene networks based on pathways annotations and on validated miRNA-target relations. Topological pathway analyses by Micrographite identified most modulated paths, showing significant gene/miRNA expression variations and changes in relations strength, then a non-redundant comprehensive meta-pathway was derived; iterative analysis of the meta-pathway identified the most modulated network for each comparison. Magia2 analysis identified mixed TF-miRNA-gene circuits. Results of parallel comparisons were merged in an integrated network.

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 MAGIA2. 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. The resulting networks of the two analyses were merged achieving a set of 181 non-redundant genes and miRNAs and their interactions that represent the compendium of the deregulated pathway-guided interactions in PMF. The network contains a selection of 401 gene-gene consolidated interactions from KEGG pathways, 45 miRNA-gene connections already validated with gene reporter assays and only one in silico predicted anti-correlated miRNA-gene connection (Figure 2a). See http://compgen.bio.unipd.it/pmf-net/ for an interactive, searchable and zoomable version of the PMF network model.

Figure 2

PMF network model. (a) PMF network integrating pathway-derived miRNA-gene networks deregulated in PMF and mixed TF-miRNA-gene circuits discovered by reverse engineering of expression data. The network (see Supplementary Figure 3 for a larger version) gives a non-redundant and comprehensive picture of most modulated paths in the two PMF vs CTR comparisons, of the impact of miRNAs on pathway genes, and of connected TF-miRNA-gene mixed circuits discovered in the study. Genes are reported as round rectangles, transcriptional factors as diamonds and miRNAs as triangles. Node colors represent the fold-change (FC) of the gene expressions in the PMF vs CTR BM (node inner color) and PMF vs CTR PB (node border color). The type of edges depends on the type of interaction: arrow for activation, T arrow in case of inhibition, and arrow line for miRNA-target predicted/supported interactions. The light blue shade indicates the part of the network resulting from miRNA and gene topological pathway analysis. The yellow shade indicates mixed TF-miRNA-gene circuits, inferred by Magia2 analysis, connected to the path-derived network. (b) Cluster analysis of expression profiles of miRNAs and genes included in the final PMF network do not show clustering of PMF patients by mutation. Samples are colored according to the carried mutation as shown in the legend (3N indicates triple negative). Sample clustering was obtained according to Euclidean distance and complete clustering. See Supplementary Figure 4 for the corresponding heatmap. (c) Summary of most deregulated pathways represented in the miRNA-gene network, and connections thereof. See http://compgen.bio.unipd.it/pmf-net/ for an interactive, searchable and zoomable version of the PMF network model.

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. Magia2 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 Magia2 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

miRNAs impacting in pathway-derived miRNA-gene networks deregulated in PMF and participating to mixed TF-miRNA-gene circuits

Main pathways and genes	miRNAs
Ca(2+) signaling	miR-125b-5p
	miR-133a
CDC42	miR-185-5p
	miR-29b-3p
	miR-29c-3p
HIF-1 alpha hypoxia pathway	miR-1244
NF-kappaB signaling	miR-124-3p
	miR-146a-5p
	miR-146b-5p
PTEN signaling	miR-106b-5p
	miR-19a-3p
	miR-20a-5p
	miR-21-5p
	miR-221-3p
	miR-222-3p
	miR-26a-5p
	miR-494
	miR-519d
	miR-93-5p
PTEN-GEF crosstalk	miR-17-5p
PTEN-CDC42 crosstalk	miR-29a-3p
RAS-PTEN crosstalk	miR-217
RAS/RAS signaling	miR-195-5p
	miR-7-5p
	let-7a-5p
	miR-122-5p
	miR-126-3p
	miR-145-5p
RAS/RAF signaling, GEF activity	let-7d-5p
RHOA	miR-31-5p
RHOA, MAPK and NF-kappaB signaling	miR-155-5p
Transcription factors	miRNAs
MYC	miR-19b-3p
ATF2, CEBPA, MYC	let-7d-5p
FOXJ2, MYC, IRF1	miR-106a-5p
FOXJ2, MYC	miR-20b-5p
REL, CEBPA, FOXJ2, MYC, IRF1	miR-17-5p

Abbreviations: miRNA, microRNA; PMF, primary myelofibrosis; TF, transcription factor.

The upper part of the table indicates the pathways or genes targeted by miRNAs in PMF; the lower part of the table indicates the transcription factors for which the interplay with miRNAs has been in silico inferred.

Table 1 summarizes miRNAs and TFs involvement. Magia2 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, and in Supplementary Figure 4, cluster analysis of samples based on expression profiles of miRNAs and genes included in the PMF network of Figure 2a did not highlight any correlation with the presence of a specific mutation. Therefore, results of cluster analysis suggest that our results provide a comprehensive picture of deregulated elements and relations thereof, that are common to all PMF patients, independent of the underlying driver mutation. Accordingly, fold-changes against control cells have been recalculated stratifying patients by driver mutations; these data do not show mutation-specific contributions in the here presented network (Supplementary Figure 5).

Confirmation of miRNA and gene differential expression in PMF

Genes and miRNAs are included in the PMF network only if they have shown coordinated expressions across patients and opposite behaviors if compared to healthy donors. Specifically, modulated genes and miRNA are included in the PMF network if they are either linked according to validated relations (pathway-derived part of the network and TF-gene or TF-miRNA relations in the part containing MAGIA2 circuits) or linked by predicted miRNA-target relations that are anyway significantly supported by expression data analysis.

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-value<0.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.

Table 2

miRNAs and genes included in the network for which differential expression in PMF patients compared to healthy controls has been confirmed

PMF vs CTR	CD34+			GR and Plasma
	Arrays	RNA-seq	RT-PCR	RT-PCR (PMF vs CTR; RQ mean±dev)
miR-106a-5p	UP			UP GR (1.9±1.2 vs 1.08±1.0)
miR-145-5p	UP			UP GR (1.7±0.9 vs 1.0±0.5)
miR-146b-5p	UP		UP
miR-155-5p	UP		UP
miR-17-5p	UP			UP GR (1.5±0.7 vs 1.1±0.9)
miR-185-5p	UP			UP GR (3.0±3.1 vs 1.5±1.2)
miR-195-5p	UP		UP	UP GR (15.3±10.5 vs 1.8±2.1) UP Plasma (12.4±17.5 vs 1.37±0.95)
miR-19a-3p	UP		UP	UP GR (15.3±10.5 vs 1.8±2.1) UP Plasma (12.4±17.5 vs 1.37±0.95)
miR-19b-3p	UP	UP	UP
miR-21-5p	UP		UP
miR-29a-3p	UP	UP	UP	UP Plasma (3.0±3.1 vs 1.5±1.2)
miR-29c-3p	UP		UP
miR-494-3p	UP	UP	UP
ANGPT1	UP		UP
ARHGEF7	DOWN		DOWN
CDC42	DOWN		DOWN
MEF2D	DOWN		DOWN
SMAD7	DOWN		DOWN

Abbreviations: miRNA, microRNA; PMF, primary myelofibrosis.

The second column indicates if miRNAs and genes resulted up- or down-regulated according to array data used for network reconstruction. The following columns indicate for which elements a significant differential expression has been previously confirmed using RNA-seq in CD34+ cells,[22] and using RT-PCR in CD34+ cells,[20] in granulocytes (GR) and in plasma samples.

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.

Figure 3

Details of the PMF network model showing deregulated miRNAs and genes that participates to specific connected pathways, linked in turn to biological processes and functions germane to the disease. (a) Akt signaling; (b) Rho GTPases, CDC42, PLD2 and PTEN; (c) HIF-1a pathway; (d) Calcium signaling.

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 Ca2+-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 PMF compared to controls and not present in the network (PMF equal to CTR PB controls and only slightly upregulated in PMF in comparison to CTR BM controls) evidencing that genetic alterations of PMF patients do not affect significantly the transcriptional behavior of the CALR gene. However, many CALR interactors and genes participating in CALR-involving pathways are represented in the PMF network: ATF4 and ATF6 participate to unfolded proteins response (Protein processing in endoplasmic reticulum), and transduce to the nucleus ER stress, eventually inducing apoptosis; CD209 and RAC1 participate to the Phagosome; several genes, including FOS/JUN and NFKB are in the Chagas disease path; ATF2/CREB, NFAT, EP300 and RAS genes share with CALR the HTLV-1 infection path.

The Ca2+-related part of the net is connected to the TFs ATF2 and CREB. Several circuits show cross-talk of ATF2 and let-7d-5p, both downregulated in PMF, to 12 correlated genes, mostly downregulated in PMF (RNF7, HIC1, HIC2, PRKACG, PRKX, PRKACA, PRKACB, MEF2D, ARMC8, KLF9, ZNF516). Among them HIC1/2 are tumor suppressor genes,[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.