Bioinformatics Analysis of Key Genes and Pathways Associated with Thrombosis in Essential Thrombocythemia.

BACKGROUND Essential thrombocythemia (ET) is a form of chronic myeloproliferative neoplasm (MPN), and thrombosis is an important complication. This study aimed to use bioinformatics analysis to identify differentially expressed genes (DEGs) in ET associated thrombosis. MATERIAL AND METHODS Two datasets were identified from the Gene Expression Omnibus (GEO) database to investigate the expression profiles in ET. The GSE103176 dataset included 24 patients with ET and 15 healthy individuals with samples from CD34+ bone marrow cells. The GSE54644 dataset included 47 patients with ET and 11 healthy individuals with samples from peripheral neutrophils. GEO2R was used to screen DEGs, followed by over-representation analysis. Protein-protein interaction (PPI) network analysis and module analysis were performed using the STRING database and Cytoscape software. Hub genes were identified using the cytoHubba Cytoscape plugin, and maximal clique centrality (MCC) was identified. The MCODE Cytoscape plugin was used to identify network clusters, or highly interconnected regions. RESULTS There were 586 and 392 DEGs identified from the GSE103176 and GSE54644 datasets, respectively. The upregulated DEGs for CD34+ cells were predominantly enriched for granulocyte activation or related pathways for biological process (BP), and secretory vesicle for the cellular component (CC). The top hub genes within CD34+ cells included CXCL1, CAMP, HP, MMP8, PTX3, ORM1, LYZ, LTF, PGLYRP1, and OLFM4. CONCLUSIONS Bioinformatics analysis identified DEGs and hub genes that interacted with CD34+ cells and neutrophils that may predict an increased risk of thrombosis in patients with ET. These preliminary findings should be validated using next-generation sequencing (NGS) and clinical studies.

Background

Essential thrombocythemia (ET) is a form of chronic myeloproliferative neoplasm (MPN) that is associated with arterial and venous thrombosis. According to the revised World Health Organization (WHO) classification, BCR-ABL-negative classical MPN comprises polycythemia vera (PV), ET, and primary myelofibrosis (PMF) [1]. The incidence of ET is reported to be 0.84 per 100,000 per year [2,3]. Thrombosis, transformation to leukemia and myelofibrosis, remain the most common complications for patients with ET [4]. ET is characterized by an increased risk of arterial and venous thrombosis [5,6], with a prevalence of thrombotic events being 20.7% [7].

Prognostic models were established to predict the risk of thrombosis in patients with ET. The International Prognostic Score for Essential Thrombocythemia (IPSET) model was developed to predict survival and occurrence of thrombosis based on a cohort study, including 867 patients, assigning patients with ET into low, intermediate and high-risk groups [8]. The prognostic assessment is based on advanced age (≥60 years), leukocyte count ≥11×109/L and a history of thrombosis, according to IPSET. In 2012, the International Prognostic Score of Thrombosis in Essential Thrombocythemia (IPSET-thrombosis) was developed, in which age >60 years was scored as one point, a history of thrombosis was scored as two points, mutation in JAK2 V617F was scored as two points, and cardiovascular risk factors were scored as one point, to assess the risk of thrombosis [9]. Patients with ET are stratified into low, intermediate, and high-risk groups identified by an accumulated score of <2, 2, and >2 points, respectively. The revision of prognostic models has been indicated for patients with ET with the mutation status of JAK2 the risk of thrombosis is comparable regardless of the presence of CV risk factors [10]. The revised IPSET-thrombosis guidelines stratify patients into the following risk groups: very low risk (age ≤60 years, no prior history of thrombosis, and no JAK2 mutation); low risk (age ≤60 years, no prior history of thrombosis, and JAK2 mutation); intermediate risk (age >60 years, no prior history of thrombosis, and no JAK2 mutation); and high risk (prior history of thrombosis and/or age >60 years with JAK2 mutation).

Although leukocytosis was not included in the IPSET-thrombosis model as a risk factor, it is reported to be significantly associated with increased arterial thrombotic events [11,12]. According to a meta-analysis including 41 studies [13], for patients with ET, the relative risk of thrombosis was 1.65 (95% CI, 1.13–1.86) in the presence of leukocytosis. While the adverse effect of leukocytosis was significant in arterial thrombosis (RR 1.45; 95% CI, 1.13–1.86) instead of venous thrombosis (RR 1.14; 95% CI, 0.65–1.98), leukocytosis was shown to be associated with thrombotic events [14], including granulocyte activation. Therefore, bioinformatics analysis of the neutrophil transcriptome may identify key genes and pathways involved in thrombosis in patients with ET.

Given that gene abnormalities are associated with hematopoietic stem cells in the pathogenesis of ET, investigation of the expression profiles of CD34+ bone marrow cells and neutrophils may provide information on the genes involved in thrombosis in ET. Therefore, this study aimed to use bioinformatics analysis to identify differentially expressed genes (DEGs) in thrombosis associated with ET. Bioinformatics analysis of datasets from bone marrow CD34+ cells and peripheral neutrophils were used to study the expression signature of differentially expressed genes (DEG) in patients with ET compared with healthy individuals.

Material and Methods

Data source

After the review of all gene expression datasets of essential thrombocythemia (ET) from the GEO database (), GSE103176 and GSE54644 were selected for analysis. The GSE103176 dataset was obtained from the GPL13667 platform (HG-U219) Affymetrix Human Genome U219 Array, and the GSE54644 was obtained from the GPL4685 platform (U133AAofAv2) Affymetrix GeneChip HT-HG_U133A Early Access Array. GEO103176 included expression profiles of bone marrow CD34+ cells from 24 ET and 15 healthy individuals. GSE54644 consisted of peripheral neutrophils samples from 47 patients with ET and 11 healthy individuals.

Identification of differentially expressed genes (DEGs) in essential thrombocythemia (ET)

The GEO2R tool () was used to distinguish DEG in ET by comparison with healthy controls. An adjusted p-value <0.01 and |log2FC| >1 were defined as the criteria for DEG.

Over-representation analysis of DEGs

The human ConsensusPathDB (CPDB) () integrated network, including genetic, signaling, gene regulation, as well as complex protein–protein interactions (PPIs) was used [8]. Gene Ontology (GO) enrichment, The Reactome () bioinformatics analysis database, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed using the human ConsensusPathDB (CPDB) () integrated network, with a cutoff p-value of 0.01.

The protein–protein interaction (PPI) network analysis

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database () was used to predict the protein–protein interaction (PPI) network information for DEGs. After mapping the DEGs to te STRING database, minimal criteria for extracting PPI pairs was 0.4. The top ten nodes were calculated and ranked by the maximal clique centrality (MCC) method and were analyzed using Cytoscape software and the cytoHubba plugin. Module analysis was performed by using the MCODE plugin (molecular complex detection) in Cytoscape. The parameters was set as degree cutoff=2, node score cutoff=0.2, k-core=2, and maximum depth=100.

Results

According to previously published data [16], the GSE103176 dataset included 11 male and 13 female patients with ET, the median age at onset was 61 years (range, 37–72 years), the median hemoglobulin is 14.6 g/dL (range, 11.9–18.1 g/dL), the median platelet (PLT) count was 677×109/L (range, 481–1168×109/L), the median white blood cell (WBC) count was 8.7×109/L (range, 5.9–18.2×109/L), 17 patients had JAK2 V617F mutations, and seven patients carried the CALR mutation. However, data on clinical features and molecular status were not available for GSE54644. A total of 586 DEGs were identified from GSE103176, including 318 upregulated DEGs and 268 down-regulated DEGs. For GSE54644, in the dataset for peripheral neutrophils, 360 upregulated DEGs and 32 down-regulated DEGs were identified. All DEGs were screened by comparing samples from patients with ET, as CD34+ cells or peripheral neutrophils, with healthy donors.

Venn diagram analysis (Funrich 3.1.3) was performed to evaluate the intersection of DEGs involved across cell subpopulation. There were 10 common upregulated DEGs and three common down-regulated DEGs for both subpopulations. Since the portion of overlapping genes was small, this demonstrated the spectrum of DEG validation across the cell groups.

Over-representation analysis for DEGs associated with CD34+ cells

The Gene Ontology (GO), Reactome, and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were performed using the human ConsensusPathDB (CPDB) online tool (). The results for upregulated DEGs that originated from CD34+ cells in ET are shown in Figure 1. GO analysis showed that upregulated DEGs originating from CD34+ cells of ET were significantly enriched in the components of the biological process (BP) that included neutrophil activation and granulocyte activation. Cellular component (CC) analysis showed that upregulated DEGs were significantly enriched in the secretory granule, secretory vesicle, and filocin-1-rich granule. Molecular function (MF) analysis showed that upregulated DEGs were enriched in immunoglobulin binding, haptoglobin binding, oxygen carrier activity, CCR1 chemokine receptor binding. In the KEGG pathway analysis, upregulated DEGs were significantly enriched in the phagosome, and NF-kappa B signaling pathways.

Figure 1

Analysis of upregulated differentially expressed genes (DEGs) in CD34+ cells in essential thrombocythemia (ET). The distribution of DEGs for the corresponding Gene Ontology (GO) term is shown in the upper plot. The Y-axis represents the gene count. The left inferior plot shows the enriched biological process (BP) terms by GO analysis. The right inferior plot shows the enriched pathways by the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. The colored dots represent the p-value, and the size of the dots correlates with the involved gene counts. The X-axis represents the gene enrichment factor.

Reactome analysis showed that the upregulated DEGs were significantly enriched in neutrophil degranulation, the innate immune system, erythrocyte oxygenation, and carbon dioxide release Down-regulated DEGs were enriched in biological process (BP), including B cell activation, lymphocyte differentiation, lymphocyte activation and regulation of cellular biosynthesis. Molecular function (MF) analysis showed that down-regulated DEGs were enriched in DNA binding, RNA binding, MHC class II protein complex binding. Cellular component (CC) analysis showed that down-regulated DEGs were enriched in the cell nucleus, nuclear lumen, nucleoplasm, and the beta-catenin/TCF complex.

Over-representation analysis for DEG originating from peripheral neutrophils

The results of enrichment analysis for upregulated DEGs of peripheral neutrophils showed a distinct signature when compared with CD34+ hematopoietic cells (Figure 2). The significantly enriched BPs included organelle organization, cell activation, protein modification process, vesicle-mediated transport, indicating that upregulated DEGs of neutrophils were not significantly enriched in granulocyte-activating pathways. Enriched CC included the cytosol, nucleus, and cytoplasm. Enriched MF included chromatin binding, ATP binding, adenyl nucleotide binding, transferase activity. KEGG analysis showed enriched DEGs in starch and sucrose metabolism, the phosphatidylinositol signaling system, lipid metabolism, and ubiquitin-mediated proteolysis. However, down-regulated DEGs of neutrophils were enriched in the cellular response to oxygen-containing compound, vesicle-mediated transport, antigen receptor-mediated signaling pathway, and the cytokine-mediated signaling pathway.

Figure 2

Analysis of upregulated differentially expressed genes (DEGs) in neutrophils in essential thrombocythemia (ET). The distribution of DEGs for the corresponding Gene Ontology (GO) term is shown in the upper plot. The Y-axis represents the gene count. The left inferior plot shows the enriched biological process (BP) terms by GO analysis. The right inferior plot shows the enriched pathways by the Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. The colored dots represent the p-value, and the size of the dots correlates with the involved gene counts. The X-axis represents the gene enrichment factor.

Protein–protein interaction (PPI) network and hub genes in CD34+ cells in ET

After mapping the DEGs to the STRING database, the PPI network of upregulated DEGs was established, containing 284 nodes and 546 edges (Figure 3). The top hub genes ranked by the maximal clique centrality (MCC) method were selected using the cytoHubba Cytoscape plugin (Table 1). After a careful review of the top hub genes in ET, the DEGs that were selected as hub genes included CXCL1, CAMP, HP, MMP8, PTX3, ORM1, LYZ, LTF, PGLYRP1, and OLFM4.

Figure 3

The protein–protein interaction (PPI) network of upregulated differentially expressed genes (DEGs). The top 20 nodes with the highest degree of connectivity are colored in red.

Table 1

The top 20 genes in the network ranked by maximal clique centrality (MCC).

Rank	Name	Score
1	CXCL1	2.44E+18
2	CAMP	2.44E+18
3	HP	2.44E+18
4	MMP8	2.44E+18
5	PTX3	2.44E+18
6	ORM1	2.44E+18
7	LYZ	2.44E+18
8	LTF	2.44E+18
9	PGLYRP1	2.44E+18
10	OLFM4	2.44E+18
11	TCN1	2.44E+18
12	CRISP3	2.44E+18
13	QPCT	2.44E+18
14	SLPI	2.44E+18
15	ARG1	2.44E+18
16	RETN	2.44E+18
17	LCN2	2.44E+18
18	CHI3L1	2.44E+18
19	BPI	2.43E+18
20	ORM2	2.43E+18

MCC – maximal clique centrality.

Module analysis

The MCODE Cytoscape plugin was used to identify clusters, or highly interconnected regions, in the network for module analysis of upregulated DEGs 12 significant clusters were established, and the gene enrichment analysis was conducted. The results of KEGG analysis are shown in Table 2. For Cluster 1/2/4/6, the genes were significantly enriched in granulocyte-activating or related pathways. Cluster 8 were significantly enriched in association with negative regulation of blood coagulation.

Table 2

Enrichment analysis for the upregulated differentially expressed gene (DEG) clusters.

ID	l included	GO analysis for BP	GO analysis for CC
Cluster 1	MMP9, TCN1, ARG1, RETN, OLFM4, DEFA4, CKAP4, BPI, QPCT, CTSS, CHI3L1, CRISP3, CAMP, CDA, CYSTM1, LCN2, SLPI, ORM2, LYZ, CD36, ORM1, LTF, PTX3, HBB, FCER1G, C3AR1, PGLYRP1, PTAFR, MGAM, PTPRJ, CXCL1, MMP8, CTSH, HP, TMEM30A	Neutrophil degranulation; neutrophil mediated immunity; granulocyte activation; myeloid cell activation involved in immune response	Secretory vesicle; tertiary granule lumen; cytoplasmic vesicle lumen; cytoplasmic vesicle
Cluster 2	STOM, SLC2A3, C5AR1, MNDA, CD14, S100A12, NMU, FCGR3B, FCGR3A, OLR1, MCEMP1, CLEC4D, CCL4L1, HCK, PTGS2, CD68, MS4A6A, S100A9, S100A8, FPR2, FGR, HCAR3, CYBB, HCAR2, CEACAM1, SELP, CD177, AQP9, CXCR2, SNAP23, VEGFA, FCAR, HK3, PF4, TYROBP, CEACAM8	Regulated exocytosis; neutrophil degranulation; neutrophil mediated immunity; granulocyte activation	Secretory vesicle; tertiary granule lumen; cytoplasmic vesicle lumen; cytoplasmic vesicle
Cluster 3	AHSP, HBA2, HBD, HBA1, HBG2, HBM, HBG1, ALAS2, HBQ1	Oxygen transport; gas transport; hydrogen peroxide catabolic process; hydrogen peroxide metabolic process	Cytosolic part; endocytic vesicle lumen; cytosol; cytosolic small ribosomal subunit
Cluster 4	GNS, GCA, CXCL2, CCL3, VNN2, NCF2, LILRB3, PRF1, GZMB, RNASE3, IGSF6, ALOX5AP, FGL2, ALOX5, TLR8, CXCL8, FCN1, CPPED1, CX3CR1, CD1D, CCL5, CCL4, LILRB2, ACTR2, SDCBP, CD86	Myeloid leukocyte activation; granulocyte activation; leukocyte activation involved in immune response; myeloid cell activation involved in immune response	Secretory vesicle; azurophil granule lumen; cytoplasmic vesicle; secretory granule lumen
Cluster 5	CD8A, GZMK, EPB42, ENTPD1, FYB, GNLY, STAT4, GZMA, IL2RB, GYPB, CA1, GZMH, KLRB1, GIMAP4	Interleukin-12-mediated signaling pathway; response to interleukin-12	–
Cluster 6	CPNE3, ITGAM, FCGR2A, CLEC4A, KLRK1, LAMP2, CEACAM6, CD83, CASP3, FPR1	Myeloid leukocyte activation; neutrophil degranulation; neutrophil mediated immunity; granulocyte activation	Azurophil granule membrane; cytoplasmic vesicle membrane; secretory vesicle; lysosomal membrane
Cluster 7	HIF1A, CDKN1A, CLEC4E, CLEC7A, CDH1, LGALS3, ITGB1	Columnar/cuboidal epithelial cell maturation; regulation of protein localization; epithelial cell maturation; lymphocyte activation;	Perinuclear region of cytoplasm; lamellipodium; membrane raft;
Cluster 8	TGM2, PROS1, VCAN, SERPINA1, THBS1	Negative regulation of blood coagulation; negative regulation of coagulation; negative regulation of hydrolase activity; negative regulation of response to wounding	Secretory granule lumen; secretory vesicles
Cluster 9	DEFA3, DEFA1B, DEFA1	–	–
Cluster 10	APOBEC3B, SAMHD1, APOBEC3A	–	–
Cluster 11	SGPP1, UGCG, GALC	–	–
Cluster 12	FBXO9, PJA2, RNF138	–	–

uDEG – differentially upregulated genes; GO – gene ontology; BP – biological process; CC – cell component.

Discussion

Essential thrombocythemia (ET) is a form of chronic myeloproliferative neoplasm (MPN) that is associated with arterial and venous thrombosis. Previously published data from a large-scale cohort study in Sweden indicated that the incidence of thrombosis in MPN has not significantly reduced during the past decades [5]. With the development of microarray technology, studies have been conducted to investigate the expression profiles of ET and the pathways associated with the pathogenesis of ET [9,10], and myeloid transformation [11]. However, the mechanisms associated with the increased risk of thrombosis in ET remain to be determined.

Leukocytosis is a significant risk factor for thrombosis in ET and other forms of chronic MPN [12-15]. In addition to the increase of cell counts, granulocytes show functional activation in patients with ET and are characterized by increased leukocyte alkaline phosphatase (LAP) and CD11b proteins on the cell membrane [16]. CD11b protein is encoded by the ITGAM gene and is implicated in cell adhesion of granulocytes. Coucelo et al. reported that increased expression of CD11b proteins correlated with thrombosis in patients with ET, as shown by cytofluorometric analysis [17]. The possible mechanism is the promotion of platelet attachment to the endothelium [18]. CD11b has an important role in the formation of leukocyte-platelet aggregates [19-21]. CD11b was also reported to be associated with thrombosis in other pathological conditions [22-24].

The profiles of CD11b/LAP mRNA expression in CD34+ bone marrow cells and neutrophils from this study are shown in Table 3. CD11b mRNA was not significantly upregulated in peripheral neutrophils, compared to healthy controls. Therefore, the increase of CD11b proteins of neutrophils in ET was not due to intrinsic upregulated mRNA expression but obtained from other sources. Also, the expression of CD11b mRNA in CD34+ cells was significantly upregulated and enriched in secretory vesicles (by CC analysis). Considering the increased number of CD34+ cells in peripheral blood of patients with ET compared with healthy controls [25], this may facilitate the interaction between peripheral CD34+ cells and neutrophils. It was assumed that CD11b were overexpressed in CD34+ cells and that the gene products were transported to neutrophils in a secretory vesicle-dependent manner, leading to the extrinsic gain of CD11b protein in neutrophils. Activated neutrophils promoted the formation of thrombosis in ET. However, further studies are required to validate the upregulated expression of CD11b mRNA by low-throughput methods, and the existence of secretory vesicles containing CD11b genes products originated from CD34+ cells. The secretory vesicles containing CD11b protein may potentially predict the risk of thrombosis and provide novel therapeutic targets.

Table 3

The profiles of leukocyte alkaline phosphatase (LAP) and CD11b expression in CD34+ cells and peripheral neutrophils.

Gene	Cell origin	Dataset	logFC	Adjusted p-value	p-value
CD11b	CD34+ cells	GSE103176	2.715432	4.18E-09	1.02E-11
	peripheral neutrophils	GSE54644	−0.00447	9.87E-01	9.83E-01
LAP	CD34+ cells	GSE103176	0.25435	1.50E-03	1.90E-04
	peripheral neutrophils	GSE54644	−0.4378	4.09E-01	3.14E-01

LAP – leukocyte alkaline phosphatase; FC – fold change.

The generation of neutrophil extracellular trap (NET) by activated neutrophils has previously been reported as an important mechanism involved in vascular thrombosis [26-28]. Also, increased expression of NET has previously been shown to promote deep vein thrombosis in mice with conditional knock-in of JAK2 V617F, and this effect was reduced in the presence of a JAK2 inhibitor [29]. In the present study, among the identified hub genes, the protein encoded by the CAMP gene binds to bacterial lipopolysaccharides (LPS), resulting in antibacterial activity [30]. The CAMP gene-encoded peptide is cleaved into two chains, antibacterial peptide FALL-19 and antibacterial peptide LL-37 [30,31]. As a specific ligand of integrin α(M)β2 (Mac-1, CD11b), LL-37 induces induced a potent α(M)β2-dependent cell migratory response and results in the activation of α(M)β2 on neutrophils [32,33]. LL-37 was shown to participate in facilitating the formation of neutrophil extracellular traps (NETs) and protection from bacterial nuclease [34,35]. Therefore, CAMP may play a synergistic role with CD11b in promoting thrombosis in ET.

According to the UniProt database, other hub genes are involved in the inflammatory response (LTF, PTX3, HP, ORM1, and LYZ), neutrophil activation and chemotaxis (CXCL1), innate immunity (PGLYRP1), proteolysis of the extracellular matrix and leukocyte migration (MMP8/9), and neutrophil degranulation (OLFM4). The association of the expression of these genes with the risk of thrombosis in ET has not been previously investigated. In the present study, the role of the identified hub genes requires required further validation. The present study provided insights for potential markers of thrombosis in patients with ET. However, both clinical factors and genetic background have an impact on thrombotic events in patients with ET.

The CALR mutation does not have a significant role in the prediction of thrombosis [42], and gene carriers have a lower incidence than patients with ET and JAK2 mutations, partly due to the tendency to cluster with lower risk factors. These genetic aberrations seem to be independent with leukocytosis demonstrated by multivariate analysis, suggesting that leukocytosis and the related expression changes may explain only part of the mechanism of thrombosis in ET. Therefore, it is hoped that the combination of the gene expression signature, genetic variables, and clinical factors will be better predictors of the risk of thrombosis for patients with ET in the future.

It has previously been reported that an increased white blood cell (WBC) count was associated with the risk of bleeding in ET [36], which appears to be a paradoxical phenomenon. In the present study, module analysis showed that Cluster 8 enriched genes were involved in the negative regulation of blood coagulation and were located in secretory vesicles, suggesting that similar vesicle-based transportation may exist and leads to both granulocyte activation and a coagulation disorder. Because the occurrence of bleeding is associated with leukocytosis [43], the abnormal expression of genes in Cluster 8 may explain this phenomenon. However, in a previously published study that reported GSE103176, no information on bleeding events was available to determine whether genes of Cluster 8 were abnormally expressed in cases of bleeding cases. In the future, comparison between patients with ET without bleeding or thrombosis, ET with thrombosis, and ET with bleeding may determine the key pathways involved and provide possible therapeutic targets for both thrombosis and bleeding in ET.

Because thrombotic events may occur after between five or ten years following the initial diagnosis of ET, data on clinical outcome is difficult to collect, s shown by previous publications on this topic. Therefore, studies are needed to directly compare the expression profiles between patients with ET and thrombosis, ET without thrombosis, and healthy controls, as data from these studies are not available. Also, following the findings from this preliminary bioinformatics analysis, future studies are required to investigate expression profiles and clinical outcome in patients with ET using long-term follow-up.

Conclusions

This study aimed to use bioinformatics analysis to identify differentially expressed genes (DEGs) in thrombosis associated with essential thrombocythemia (ET). In the present study, bioinformatics analysis identified DEGs and hub genes that interacted with CD34+ cells and neutrophils that may predict an increased risk of thrombosis in patients with ET. These preliminary findings should be validated using next-generation sequencing (NGS) and clinical studies.