Identification of potential therapeutic target genes and miRNAs for primary myelofibrosis with microarray analysis.

The aim of the present study was to identify potential therapeutic target genes and miRNAs for primary myelofibrosis (PMF). The dataset GSE53482 was downloaded from the Gene Expression Omnibus database. The differentially expressed genes (DEGs) and differentially expressed miRNAs (DEMs) of peripheral blood (PB) cluster of differentiation (CD)34+ cells from PMF patients (PB-PMF group) and peripheral blood CD34+ cells from healthy individuals (PB-control group) were analyzed using the Linear Models for Microarray Data package in R. The Kyoto Encyclopedia of Genes and Genomes was used for pathway enrichment analysis. MiRNA-gene joint enrichment analysis was performed by ENViz and a miRNAs-gene regulatory network was constructed. A total of 1,182 DEGs (773 upregulated and 109 downregulated) and 48 DEMs (28 upregulated and 20 downregulated) were identified. According to the pathway enrichment analysis, a number of DEGs were enriched in metabolic pathways, including IDH1 and DNMT1. Other DEGs were enriched in the citrate cycle (tricarboxylic acid cycle; IDH1 and IDH3A) and certain DEGs were enriched in pyrimidine metabolism, including CARD8. For downregulated genes, certain DEGs were enriched in the spliceosome, including SF3B1 and CDC40. Furthermore, hsa-miR-127-3p, hsa-miR-140-3p and hsa-miR345 were associated with cell cycle-related biological processes, signal transduction and cell surface receptor signaling pathway. The DEM-DEG regulatory network indicated that hsa-miR-543 regulated 113 genes, including CARD8 and TIFA. The present study identified a number of genes, including IDH1, DNMT1, SF3B1 and CARD8, and miRNAs, including hsa-miR-127-3p and hsa-miR-140-3p, which may be therapeutic targets in the treatment of PMF.

Introduction

Primary myelofibrosis (PMF) is a clonal stem cell disorder currently classified as a myeloproliferative neoplasm (MPN), and is typically characterized by megakaryocyte proliferation and a typical morphology, extensive bone marrow fibrosis and splenomegaly secondary to extra-medullary hematopoiesis (1). Ineffective erythropoiesis and extra-medullary hematopoiesis are the primary causes of anemia and organomegaly, respectively (2). PMF is a rare hematological malignancy, with an estimated incidence of 0.3–1.5 cases per 100,000 individuals annually (3). Although hematopoietic stem cell transplant and splenectomy are the conventional treatment options for PMF patients, the 5-year survival rates in the transplant cohorts are 69% for low-risk, 52% for intermediate (int)-1, 50% for int-2, and 32% for high-risk patients (4). The mortality rate and median survival for all patients after splenectomy is 14% and 18.5 months (range 0.1–71.9 months), respectively (5). Thus, there is a need for novel, effective treatment strategies for PMF.

Previously, the genetic underpinnings of MPN were identified and a number of studies have indicated that the Philadelphia chromosome, which has been demonstrated to harbor an oncogenic BCR-ABL1 fusion transcript, is the disease-causing mutation in chronic myelogenous leukemia (6–8). In addition, Janus kinase 2 (JAK2) gain of function (JAK2V617F) and myeloproliferative leukemia virus (MPL) gene mutations have been described in ~50% and 5–10% of PMF patients, respectively (9). Norfo et al (10) identified biomarkers of PMF, including FGR, LCN2, and OLFM4, and using miRNA-gene expression integrative analysis, revealed that miR-155-5p affects the in vitro expansion of the megakaryocyte lineage through the modulation of JARID2 expression. This suggests that the miR-155-5p/JARID2 axis may contribute to the abnormal megakaryopoiesis typical in PMF. Nevertheless, the disease genes and its related miRNA molecular signature of PFM remain to be completely identified and the cause of disease is not clearly understood. Therefore, it is essential to comprehensively describe the alterations in gene expression profiles.

In order to reduce the variables, the present study only investigated the expression microarray data of peripheral blood. The microarray dataset GSE53482 deposited by Norfo (10) was used to identify differentially expressed genes (DEGs) and differentially expressed miRNAs (DEMs). Then miRNAs-target gene network analysis combined with functional enrichment analysis was performed. The aim of the present study was to investigate PMF-associated miRNAs and key genes, which may serve key roles in the treatment of PMF.

Materials and methods

Microarray data

The expression profile dataset of GSE53482 was downloaded from Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) (10). A total of 73 chips were available, including 42 peripheral blood (PB) cluster of differentiation (CD) 34+ cells from PMF patients (PB-PMF group), 16 peripheral blood CD34+ cells from healthy individuals (PB-control group), 15 bone marrow (BM) CD34+ cells of healthy individuals (BM-control group). In order to reduce the variables, only the expression microarray data from PB tissue were analyzed. Raw data from mRNA and miRNA profiles were collected using GPL13667 (HG-U219) Affymetrix Human Genome U219 Array and GPL14613 (miRNA-2_0) Affymetrix Multispecies miRNA-2_0 Array (Affymetrix Inc., Santa Clara, CA, USA), respectively.

Data preprocessing

The downloaded raw gene expression data of mRNA and miRNA were preprocessed. Background correction, quantile normalization, probe summarization and log2-transformation were performed using a Robust microarray analysis package (version 3.3.2; Bioconductor, Buffalo, NY, USA) (11). If one probe corresponded to multiple genes, the expression value of this probe was removed. However, if multiple probes corresponded to one given gene, the mean expression value was defined as the final expression value of the gene.

DEGs and DEMs screening

DEGs and DEMs of the PB-PMF and PB-control groups were analyzed by the Linear Models for Microarray Data (LIMMA) package (R) (version: 3.30.3; Bioconductor) (12). In order to circumvent the multiple testing problems, which may induce false-positive results, the Benjamini-Hochberg procedure was used to control false discovery rate (FDR) by adjusting the raw P-values (13). log2FC (fold change)≥1 and FDR<0.05 were set as the screening criteria for DEGs. However, log2FC≥0.5 and FDR<0.05 were set as the screening criteria for DEM (12).

Functional enrichment analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was conducted with R package clusterProfiler (version: 3.2.14; Bioconductor) (14). KEGG is a database resource for understanding high-level functions and utilities of the biological system, including the cell, organism and ecosystem, from molecular level information (15). The P-values were adjusted using the Benjamini-Hochberg procedure. P<0.05 was considered to indicate a statistically significant difference.

MiRNA-gene joint enrichment analysis

In order to investigate the biological function of DEMs, miRNA-gene joint enrichment analysis was performed by Enrichment Analysis and Visualization (ENViz) in Cytoscape (version 3.4.0; Cytoscape Consortium, San Diego, CA, USA) which was driven by three input matrices: the primary data matrix, the annotation matrix and the pivot data matrix (16). The multiple types of data were integrated in the same sample and the graphical results were presented (16). Primary data (gene expression profile), pivot data (miRNA expression profile) and annotation data were entered into the procedure. The enrichment analysis for the primary element, which was most associated with the pivot element, was performed through ENViz to annotate the biological function of the pivot element.

Construction of miRNAs-target gene regulatory network

MiRNAs are single-stranded RNA molecules of 21–23 nucleotides in length that regulate gene expression generally by inducing cleavage of mRNAs or inhibiting translation in a sequence-specific manner (17). miRanda (www.microrna.org) (18), MirTarget (http://cbit.snu.ac.kr/~miTarget/) (19), PicTar (http://pictar.bio.nyu.edu/) (20), TargetScan (www.targetscan.org) (21), miRecords (http://mirecords.biolead.org/) (22) and miRwalk (www.umm.uni-heidelberg.de/apps/zmf/mirwalk) (23) were used for the miRNA target gene prediction. MiRNA-genes present in at least three of the aforementioned databases were used to construct the miRNA-gene regulatory network.

Statistical analysis

The expression levels of genes are presented as mean ± standard deviation. Statistical analysis was performed using SPSS 20.0 (IBM SPSS, Armonk, NY, USA). Student's t-test was used to analyze differences between PB-PMF and PB-control groups. P<0.05 was considered to indicate a statistically significant difference.

Results

Screening of DEGs and DEMs

A total of 1,182 DEGs including 773 upregulated and 109 downregulated genes were selected from the PMF samples (Table I). Following DEMs screening, 48 DEMs were identified, including 28 upregulated and 20 downregulated DEMs (Table II).

Table I.

Number of differentially expressed genes and microRNAs in the PB-PMF group compared with PB-control group.

Category	Genes	MicroRNAs
Upregulated	773	28
Downregulated	109	20
Total	1,182	48

PMF, primary myelofibrosis; CD, cluster of differentiation; PB-PMF group, peripheral blood CD34+ cells from PMF patients; PB-control group, peripheral blood CD34+ cells from healthy individuals.

Table II.

Differentially expressed microRNAs in the peripheral blood cluster of differentiation (CD)34+ cells from primary myelofibrosis patients.

MicroRNAs	logFC	FDR
hsa-miR-127-3p	3.387003333	1.59E-11
hsa-miR-487b	2.898958717	4.77E-09
hsa-miR-409-3p	2.38774578	1.90E-06
hsa-miR-1246	2.117570289	0.001769255
hsa-miR-486-5p	1.987400714	0.000180305
hsa-miR-1308	1.918602857	1.74E-09
hsa-miR-494	1.770673137	7.10E-06
hsa-miR-134	1.503354598	3.53E-05
hsa-miR-382	1.503337152	0.000110142
hsa-miR-503	1.425249491	8.13E-05
hsa-miR-486-3p	1.392454196	0.047462229
hsa-miR-18a	1.270815722	4.58E-08
hsa-miR-1979	1.267982054	0.012984921
hsa-miR-92a-1	1.184540893	0.003833481
hsa-miR-152	1.141451113	0.011362506
hsa-miR-379	1.042055815	9.74E-07
hsa-miR-432	0.957774043	9.27E-05
hsa-miR-19b	0.940775327	1.48E-05
hsa-miR-195	0.940193085	1.37E-05
hsa-miR-376c	0.849398149	0.007550278
hsa-miR-485-5p	0.811446693	0.000110142
hsa-miR-543	0.754953378	0.003583896
hsa-miR-34a	0.754130661	0.000852361
hsa-miR-21	0.710161631	0.023095868
hsa-miR-18b	0.633363635	0.000111584
hsa-miR-370	0.604671741	0.016171501
hsa-miR-20b	0.510119307	0.002615483
hsa-miR-146b-5p	0.500727976	0.023667101
hsa-miR-744	−0.53118132	1.74E-05
hsa-miR-1915	−0.55023808	0.000708478
hsa-miR-423-5p	−0.55960845	0.000997801
hsa-miR-3196	−0.59688598	0.028223221
hsa-miR-140-3p	−0.62122708	0.001454615
hsa-miR-23a	−0.63403583	0.003583896
hsa-miR-345	−0.64862991	0.01242366
hsa-miR-1275	−0.66625937	0.003833481
hsa-miR-638	−0.74490649	0.005061941
hsa-miR-4281	−0.76486024	0.007550278
hsa-miR-2861	−0.78904119	0.010372911
hsa-miR-324-3p	−0.81814792	0.009448033
hsa-miR-27a	−0.85423333	0.007550278
hsa-miR-197	−0.9544247	0.000228802
hsa-miR-762	−0.95937914	0.003304652
hsa-miR-766	−1.04226524	0.000725484
hsa-miR-1469	−1.04372292	9.21E-05
hsa-miR-181a-2	−1.23940304	0.003054099
hsa-miR-3185	−1.31621408	0.007373289
hsa-miR-663	−1.55496702	3.62E-07

FC, fold-change; FDR, false discovery rate.

Functional enrichment analysis results

According to KEGG pathway enrichment analysis, the upregulated DEGs were significantly enriched in metabolic pathways (P=3.21E-04), cell cycle (P=1.31E-05), ribosome biogenesis in eukaryotes (P=2.83E-03) and DNA replication (P=3.21E-04). The downregulated DEGs were significantly enriched in the Toll-like receptor signaling pathway (P=0.025), osteoclast differentiation (P=0.002), protein processing in the endoplasmic reticulum (P=0.016) and the transforming growth factor-β signaling pathway (P=0.016) (Fig. 1). As presented in Table III, for upregulated genes a number of DEGs were significantly enriched in metabolic pathways (IDH1 and DNMT1; P=3.21E-04), the citrate cycle (TCA cycle; IDH1 and IDH3A; P=1.47E-02) and within the pyrimidine metabolism (CARD8; P=2.92E-02). However, for downregulated genes, a number of DEGs were significantly enriched in the spliceosome (SF3B1 and CDC40; P=4.54E-02).

Figure 1.

Pathway enrichment analysis for up- and downregulated genes. The abscissa represents P-value following log2-transformation, the ordinate represents the enriched biological pathways with log2FC≥1 and P<0.05, and the x-axis represents the number of genes. TGF-β, transforming growth factor-β; TCA, tricarboxylic acid cycle (citric acid cycle).

Table III.

Functional annotation for differentially expressed genes in primary myelofibrosis.

Category	ID	Description	Count	Gene symbol	P-value
Upregulated genes	hsa04110	Cell cycle	21	BUB1, BUB1B, CCNA2, CCNB1, CCNB2,CDC6, CDK1, CHEK1, CUL1, E2F5	1.31E-05
	hsa01100	Metabolic pathways	81	ACADM, ACO1, B3GALNT1, B3GALT6, BCAT2, IDH1, CYP2R1, DNMT1, EXTL2, GALK2	3.21E-04
	hsa03030	DNA replication	9	FEN1, MCM2, MCM3, MCM4, MCM6, PCNA, POLE2, RNASEH2A, RPA2	3.21E-04
	hsa03008	Ribosome biogenesis in eukaryotes	12	EFTUD1, GNL3L, IMP3, MPHOSPH10, POP5, POP7, PWP2, REXO2, RPP40, UTP14C	2.83E-03
	hsa00020	Citrate cycle (TCA cycle)	6	ACO1, DLAT, IDH1, IDH3A, PDHB, SUCLA2	1.47E-02
	hsa00280	Valine, leucine and isoleucine degradation	7	ACADM, ALDH3A2, BCAT2, HIBADH, HIBCH, HMGCL, MUT	2.12E-02
	hsa00240	Pyrimidine metabolism	11	CARD8, PNP, PNPT1, POLE2, POLR1B, POLR2H, POLR3B, RRM2, TK1, TYMS, UMPS	2.92E-02
	hsa04146	Peroxisome	9	ACSL3, ACSL5, FAR2, HMGCL, HSD17B4, IDH1, PEX3, PRDX1, SLC27A2	4.46E-02
Downregulated genes	hsa04380	Osteoclast differentiation	10	CTSK, CYLD, FOSL2, IFNGR1, IL1A, JUND, NFKBIA, PIK3CA, RELB, SQSTM1	2.31E-03
	hsa03018	RNA degradation	6	BTG1, BTG3, CNOT8, DCP1A, DDX6, XRN1	1.57E-02
	hsa04141	Protein processing in endoplasmic reticulum	9	DDIT3, EIF2AK3, ERO1L, ERO1LB, SEC24A, SEC24B, UBE2G1, UBE2J1, YOD1	1.57E-02
	hsa04120	Ubiquitin-mediated proteolysis	8	BIRC3, ITCH, RNF7, SIAH1, SMURF1, UBE2G1, UBE2H, UBE2J1	1.57E-02
	hsa04350	TGF-β signaling pathway	6	EP300, ID1, ID2, ID3, SMAD7, SMURF1	1.57E-02
	hsa04722	Neurotrophin signaling pathway	7	CDC42, FOXO3, IRS2, MAPKAPK2, NFKBIA, NFKBIE, PIK3CA	2.48E-02
	hsa04620	Toll-like receptor signaling pathway	6	CCL3, CTSK, IL8, NFKBIA, PIK3CA, TLR1	2.48E-02
	hsa05142	Chagas disease (American trypanosomiasis)	6	CCL3, GNA15, IFNGR1, IL8, NFKBIA, PIK3CA	2.48E-02
	hsa04144	Endocytosis	9	ADRBK1, CDC42, CHMP2B, CHMP4B, ITCH, SMAD7, SMURF1, STAM, VPS37B	2.48E-02
	hsa04062	Chemokine signaling pathway	8	ADRBK1, CCL20, CCL3, CDC42, FOXO3, IL8, NFKBIA, PIK3CA	4.00E-02
	hsa03040	Spliceosome	6	CDC40, RBMXL1, SF3B1, SLU7, SYF2, WBP11	4.54E-02

TGF-β, transforming growth factor-β; TCA cycle, tricarboxylic acid cycle (citric acid cycle).

MiRNA gene joint enrichment analysis

In order to study the effect of DEMs and DEGs in different biological processes (BP), DEMs and DEGs joint enrichment analysis was performed. DEMs were primarily enriched in the biological processes associated with the cell cycle, immune response and response to stimulus. It was also demonstrated that the expression of hsa-miR-127-3p, which exhibited a significantly higher expression, was negatively associated with four cell cycle-related BP terms, including cell cycle, mitotic cell cycle, mitotic cell cycle process and cell cycle process. Notably, hsa-miR-140-3p was negatively associated with the cell cycle and positively associated with different biological pathways, including signal transduction and cell surface receptor signaling pathway (Fig. 2).

Figure 2.

miRNA-GO enrichment network. Cumulative enrichment score for each GO term is the sum of enrichment scores for all edges connected to the GO term node. Unconnected parents are colored gray. Pivot nodes (miRNAs) are colored gray. Edges are color-coded by direction of gene-pivot association. Red edges indicate that the genes related to the enriched GO term are positively regulated by the pivot miRNA, and blue edges indicate that the genes related to the GO term are negatively regulated by the pivot miRNA. The colour of boxes (yellow to red) indicates the enrichment score from 8.09 to 152.09. Black arrows indicate the important GO terms. GO, gene ontology; miRNA, microRNA.

DEM-DEG regulatory network analysis

Based on the aforementioned 7 miRNA target gene prediction methods, the regulation associated between DEMs and DEGs was identified (Fig. 3). The regulatory network included 27 miRNAs, 427 genes and 820 edges. The top 5 miRNAs with the highest connective degree in the DEM-DEG regulatory network included hsa-miR-543, hsa-miR-494, hsa-miR-20b, hsa-miR-23a and hsa-miR-152 (Table IV). For example, hsa-miR-543 regulated 113 genes, including CARD8, PELI1 and TIFA (Fig. 3).

Figure 3.

Differentially expressed microRNAs-differentially expressed gene regulatory network. The triangle nodes represent miRNAs and the circle nodes represent genes. Red squares and circles represent the key DEGs and miRNAs, respectively.

Table IV.

Top 5 microRNAs with the highest connective degree in the differentially expressed microRNAs-differentially expressed genes regulatory network.

MicroRNA	Connective degree
hsa-miR-543	113
hsa-miRt-494	87
hsa-miR-20b	73
hsa-miR-23a	62
hsa-miR-195	54

There was a regulatory association between miRNA and transcription factors (TF) in the DEM-DEG regulatory network (Fig. 4), which included 13 miRNAs and 9 TFs. Of these, MEF2C and BCL11A were regulated by >2 DEMs; for example MEF2C connected with hsa-miR-21, hsa-miR-23a, hsa-miR-18a and hsa-miR-18b.

Figure 4.

Differentially expressed microRNA-transcription factor regulatory network. The triangle nodes represent microRNAs, and the circle nodes represent transcription factors. Red squares indicate the transcription factors that are regulated by more than two differentially expressed microRNAs.

Discussion

PMF is a Philadelphia-negative myeloproliferative neoplasm with a heterogeneous clinical presentation (24,25). In the present study, a total of 1,182 DEGs and 48 DEMs were identified from PB CD34+ cells of patients with PMF. Among them, a number of upregulated DEGs were significantly enriched in metabolic pathways, including IDH1 and DNMT1. Certain DEGs were markedly enriched in pyrimidine metabolism, including CARD8. A number of downregulated DEGs were enriched in the spliceosome, including SF3B1.

IDH1 is located on chromosome 2q33.3 and encodes isocitrate dehydrogenase 1, which catalyzes the oxidative decarboxylation of isocitrate to α-ketoglutarate (26). IDH1 serves an important role in cytoplasmic micotinamide adenine dinucleotide phosphate production. IDH mutations are frequent in blast-phase myeloproliferative neoplasms and therefore, may contribute to leukemic transformation. Previous findings have reported higher frequencies of IDH1/IDH2 and LNK mutations in blast-phase PMF, suggesting a pathogenetic contribution to disease progression (27). The mutant IDH1 was associated with a worse prognosis in cytogenetically normal acute myeloid leukemia with a NPM1+/FLT3− molecular profile (28,29). The phenotypic and prognostic effects of the IDH1 mutation among patients with PMF have been previously shown (27). Clinical observations underscore the potential relevance of investigating other mutations or epigenetic abnormalities that functionally mimic the mutant IDH, in JAK2-mutated PMF (30).

DNMT1 encodes DNA (cytosine-5-)-methyltransferase 1, which belongs to a family of DNA methyltransferases (DNMT) including DNMT1, DNMT3A and DNMT3B, and catalyzes the addition of methyl groups to cytosine residues of CpG nucleotides (31). DNA methyltransferase is among the most common occurring mutation identified in myeloid malignancies (32). DNA methylation of DNMT1 consists of an enzymatic addition of a methyl group at the carbon 5 position of cytosine in the context of the sequence 5-cytosine-guanosine (CpG), mediated by DNMTs (33). Although the mechanisms leading to aberrant DNA hypermethylation remain to be fully elucidated, increased levels of DNMT1 are observed in malignant myeloid blasts compared with normal bone marrow mononuclear cells, suggesting that DNMT overexpression contributes to gene promoter hypermethylation and in turn, to leukemogenesis (34). Furthermore, DNMT family genes were analyzed in PMF, chronic myelomonocytic leukemia and advanced phases of myeloproliferative neoplasms (35).

SF3B1 is a critical component of the splicing machinery, which catalyzes the removal of introns from precursor mRNA and it was observed to be enriched in spliceosome in the present study. One of the most notable findings arising from the recent next-generation sequencing data has been the identification of mutated splicing factor 3b1 (SF3B1) as a putative candidate driver gene of the common lymphoid malignancy, chronic lymphocytic leukemia (36). A clinical, histopathological and genetic association study of 155 patients indicated that SF3B1 mutations in PMF may encode a key component of the mRNA splicing machinery (37). SF3B1 mutations were initially identified in myelodysplastic syndromes and other solid tumors suggesting it is serving an important role in cancer biology (38). SF3B1 mutations are associated with disease subtype (38), progression (39,40), chemotherapy resistance (41) and longer overall patient survival (42).

CARD8 is the role of caspase activating and recruitment domain 8 and encodes a protein whose exact function remains unclear (43). In the current study, it was observed that CARD8 was upregulated and enriched in the pyrimidine metabolism. CARD proteins are generally involved in signal transduction pathways, important for apoptosis, inflammation, immune function and nuclear factor-κB activation (1). Furthermore, a previous study demonstrated that the presence of somatic mutations and expression analysis of CARD family gene may be potentially relevant to PMF pathogenesis, at least in rare cases (1). In the present study, the involvement of CARD8 was predicted by the bioinformatics methods mentioned. This, combined with previous studies, indicates that it may be a therapeutic target for PMF.

In the present study, hsa-miR-127-3p was negatively associated with the four cell cycle-related BP terms, which may suggest that the high expression of hsa-miR-127-3p inhibits the proliferation of stem cells in peripheral blood (44). However, hsa-miR-140-3p was identified to be the most active miRNA in the current study. While it was negatively associated with cell cycle pathways, it was positively associated with signal transduction and cell surface receptor signaling pathway. A previous study demonstrated that miR-140-3p is differentially expressed in hematopoietic lineages, and it may be associated with platelet production and activation (45). Hsa-miR345, hsa-miR20b and hsa-miR486-5p were positively associated with cell cycle-related terms such as mitotic cell cycle process and mitotic cell cycle. In particular, hsa-miR345 was positively associated with DNA replication and metabolic pathways (46). This suggested that these miRNAs play key roles in PMF. A clinical study demonstrated that in patients with PMF, the number of peripheral blood CD34+ cells was significantly higher than normal (47). Therefore, the abnormal expression of these miRNAs in patients with PMF affected the CD34+ cell cycle in peripheral blood, thus affecting the occurrence and development of PMF. Although the target genes and miRNAs for PMF were identified, validation of the results is required in the future.

In conclusion, the DEGs and DEMs were obtained in PMF. A number of significant genes, including IDH1, DNMT1, SF3B1 and CARD8 and miRNAs, including hsa-miR-127-3p and hsa-miR-140-3p, may serve pivotal roles in the development of PMF. These genes and miRNAs may be used as specific therapeutic targets in the treatment of PMF, which may be validated by experiments in future studies.