Transcriptome profiling reveals target in primary myelofibrosis together with structural biology study on novel natural inhibitors regarding JAK2

Introduction

Primary Myelofibrosis (PMF), one of the myeloproliferative neoplasms (MPNs) which arises from clonal proliferation of hematopoietic stem cells and leads to progressive bone marrow (BM) fibrosis, is currently classified with polycythemia vera (PV) and essential thrombocythemia (ET) under the broad WHO category of MPNs [1]. The main pathological manifestations are marked reactive bone marrow fibrosis, osteosclerosis, angiogenesis, extramedullary hematopoiesis (EMH) and abnormal cytokine expression. Patients with PMF frequently complain about fatigue as well as symptoms due to splenomegaly, such as ascites, fever and night sweats [2]. Besides, constitutional symptoms including severe anemia, marked hepatosplenomegaly, bone pain, infarct, pruritus, cachexia, thrombosis and bleeding were also reported appearing in PMF [3–5].

The overall prognosis of patients with PMF is generally poor, with a short median survival and poor life quality. The estimated median survival time was 15 years for patients younger than 60, and 6 years for patients older than 60 [5, 6]. Poor life quality may due to constitutional symptoms and cachexia. Unfortunately, current chemotherapy of PMF including JAK2 inhibitors do not provide a promising prospect view [7, 8]. Fedratinib (TG101348) is a JAK2-selective inhibitor which demonstrated clinical benefits in patients with MF in early-phase clinical trials, and is approved by the Food and Drug Administration (FDA) as the first drug to treat myelofibrosis [9]. Currently, only this drug has shown a relatively efficient efficacy, the disadvantage of this inhibitor also remains obvious and needs to be improved. Therefore, there is an urgent need to discover more innovative drug candidates regarding PMF in order to improve this situation.

In recent years, with the development of bioinformatics and microarray technology, they were widely employed to analyze malignant neoplasms. It helps us to study initiation, progression and metastasis of carcinoma under molecular level, making it possible to analyze the genetic alteration and molecular mechanisms in the development of PMF. Bioinformatics analysis allows researchers reveal molecular therapeutic target and provide a theory basis through a systematic, effective and accurate manner. Natural compounds and their derivatives always provide unique chemical structures as well as potential biological functions in today’s pharmacologic market, due to their malleable and readily available property, they have made a great contribution to medication screening [10, 11]. Small molecules and natural products screening are an essential aspect if not the first means to tackle an emergent or uncontrollable disease, these molecules and approaches have been used to improve chemotherapy as well as overall cancer treatment [12, 13]. Consequently, bioinformatics combined with structural biology were implemented in this study to accelerate the discovery of PMF drugs.

In the aspect of MPN, it is well known that the accurate differential diagnosis is the key to conducting prognosis as well as therapy [14]. Firstly, this study aimed at identifying hub genes through bioinformatics method in the occurrence of PMF. Totally 3 messenger RNA microarray datasets (GSE26049, GSE53482, GSE61629) involving primary myelofibrosis were downloaded from the Gene Expression Omnibus (GEO) to identify driver genes and key pathways causing the progression of PMF. Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to discover molecular function changes and abnormal signaling pathways by showing their biological process, cellular component and molecular function. Protein-protein interaction network analysis (PPI) was then carried out to visualize the connection between different genes. After identifying the hub genes causing the occurrence of PMF, this study further employed a series of structural biological and chemical methods (such as virtual screening, molecular docking, toxicity prediction, molecular dynamics simulation etc.) to screen and identify novel lead compounds that could have biological effects targeting hub genes in the treatment of PMF. This may contribute to new ideas and resources for drug discovery in the pharmaceutical market. This study provided a list of drug candidates as well as their pharmacologic properties from ZINC15 database, which could offer a solid practical foundation for gene product inhibitors’ research. The whole diagram and framework of this study were shown in Figure 1.

Results

Discussion

Primary Myelofibrosis (PMF) is a myeloproliferative neoplasm characterized by stem cell-derived clonal myeloproliferation, which is frequently but not always accompanied by gene mutation, abnormal cytokine expression, bone marrow fibrosis, anemia etc., [4, 17, 18]. PMF starts insidiously and progresses slowly, often with a long asymptomatic period before diagnosis, and patients could suffer from a series of symptoms [19–22]. Currently, the treatment of PMF primarily emphasizes on palliative treatment, aiming at relieving anemia, splenomegaly, constitutional symptoms, and bone pain. The reason caused this situation may be blamed for lacking effective diagnostic methods and biomarkers at the early stage of the disease [23, 24]. Accordingly, a comprehensive understanding of the molecular mechanism, progression as well as pathogenesis of PMF is imperative to formulate efficiently diagnostic and therapeutic strategies. Bioinformatics analysis combined with structural biology method was applied in this study to fully investigate the PMF from the carcinogenesis to the treatment, and assess the properties of existed drug Fedratinib in the meantime.

In the current study, 54 PMF samples and 72 normal samples were extracted from 3 mRNA microarray series from GEO. Firstly, raw data (“.CEL” format) was downloaded for bioinformatics analysis. Then this study preprocessed the raw data as well as discussing the reason for data preprocessing: GEO do not take responsibility for uploaders’ chips quality, so it remains unknown what algorithm the author used in their “series_matrix.txt” file, different “series_matrix.txt” file might be generated from different algorithms, platforms or time. Some matrices’ gene expression levels may fluctuate at ten whereas others approach tens of thousands, which invisibly increases the difficulty to compare more than one dataset for study. The purpose of data standardization was to eliminate high co-relevance generated by systematic bias while preserving high correlation of gene expression level caused by the real biological reason. Through RMA algorithm, standard gene expression profiles data were obtained by background correction and standardized pretreatment. After that, RNA degradation image of each GSE series as well as NUSE box plot was plotted to verify the chips’ quality, and 2 samples were removed due to their poor chip quality. Next, “batch effects” were eliminated for the gene expression profiles, which was due to sample preparation or array variation (charge, type, and/or platform), it could mask or confound the real difference among biology [15, 16, 25]. Non-biological experimental variation or “batch effects” was commonly observed across multiple batches of microarray experiments, which rendered the task of combing data from different batches difficult. In order to check batch effects among these 3 series, principal component analysis (PCA) was performed, and PCA image was plotted, results illustrated that these 3 datasets were clustered in different places, which meant they had a different pattern of gene expression. PCA analysis was conducted again after batch effects were removed, and results displayed that the points of 3 datasets interleaved with each other, with almost no outliers. PCA analysis provided solid evidence for either the existence or elimination of batch effects. Consequently, eliminating batch effects was essential for researchers to analyze the real difference in the gene expression profiles.

Next, bioinformatics analysis was conducted to find out the hub genes of PMF. A total of 810 DEGs, 337 DEGs and 1742 DEGs were identified from 3 datasets, respectively. Altogether, 95 mutual DEGs were obtained among those 3 datasets, as shown in Venn plot. After GO analysis of abnormal expression genes, we detected that those up-regulated genes were mainly associated with erythrocyte differentiation, killing of cells of other organism, ATP binding, protein binding and hemoglobin complex, which explained why the fast multiplication of tumor cells as well as the destroy of normal cells. Down-regulated genes were primarily involved in humoral immune response, transcription factor binding and plasm membrane. Our results suggested that the abnormal changes including erythrocyte differentiation, neutrophil degranulation may occur in hemoglobin complex, which agreed with the previous study that obstacle to myeloproliferative differentiation may cause aplastic anemia in the progression of PMF [26]. Some studies have also demonstrated that cell membrane ion channels, such as iron ion binding, participate in cell signal transduction, proliferation, apoptosis as well as regulation of gene expression in tumor levels [27, 28]. Furthermore, analyses of KEGG and GSEA revealed that the mutual up-regulated DEGs were chiefly enriched in the cell cycle, JAK-STAT signaling pathway and TNF signaling pathway. Previous studies have reported that tumor necrosis factor (TNF) interacts with tumor cells to trigger cytolysis or cell death. TNF could also promote inflammatory responses, which acts via the TNF receptor (TNFR) and is part of the extrinsic pathway for triggering apoptosis [29–31]. Advanced studies have further reported that dysregulation of the cell cycle accelerated the carcinogenesis of tumor [32–34], and drugs acting on the cell cycle may benefit patients [35], which was subscribed to our findings that cell cycle was aberrantly activated in PMF. JAK-STAT signaling pathway would be discussed in detail later.

With the aim of screening hub genes among DEGs identified in our former work, the 95 mutual DEGs were analyzed with construction of PPI network based on the STRING database. 11 different algorithms were applied in Cytoscape software (“cytohubba” plug-in), the final hub genes were selected by intersecting the top 25 genes in each algorithm, which could provide a solid reliability of the selected hub genes, including LCN2, JAK2, MMP8, CAMP, DEFA4, LTF, MPO, HBD, STAT4, EBF1. Particularly JAK2 and LCN2 ranked prominent in these 11 algorithms, suggesting that they were the most essential genes in the occurrence of PMF.

LCN2, located on chromosome 9q34, is a secreted protein that belongs to the lipocalins, a group of transporters of small lipophilic molecules such as iron, steroids, fatty acids and lipopolysaccharides in circulation [36]. The altered expression of LCN2 could trigger disease in several pathologic organs, including liver injury, steatosis, kidney injury, brain injury, cardiomyopathy, musculoskeletal disease and cancer of several organs [37]. Previous research cloned the full length of human LCN2 cDNA and results demonstrated that this lipocalin was mainly expressed in myeloid cells [38]. Existed studies accessed LCN2 protein levels in the serum between PMF samples and healthy donors [39], results indicated that the level of LCN2 secreted protein was significantly higher in PMF patients than in healthy donors, suggesting that LCN2 could be considered as a biomarker. In the course of last decades, LCN2 had been well studied as a potential biomarker whether for the kidney injury or to estimate the outcome of different diseases. In one multicenter, prospective cohort study, a significant correlation between urinary injury biomarkers (LCN2, KIM-1, IL-8, L-FABP and albuminuria) and the disease outcome indicated that LCN2 only distinguished progression alone from non-progression, elucidating that LCN2 could be potentially served to recognize patients who were at high risk of progression [40]. Those studies provided solid evidence for the utility of LCN2 as a biomarker, and LCN2 could also be a promising diagnostic and therapeutic target of PMF in our study.

Until now, four JAKs have been recognized in mammals: JAK1, JAK2, JAK3 and TYK2, of which Janus Kinase 2 (JAK2) is a cytoplasmic tyrosine kinase engaged with numerous signaling pathways, referring interleukin-3, granulocyte colony-stimulating factor, granulocyte-macrophage colony-stimulating factor, receptors for erythropoietin and thrombopoietin [41]. The pathogenetic contribution of JAK2 is currently thought to involve the upregulation of JAK-STAT signaling pathway [42]. JAK-STAT signaling pathway is a signal transduction pathway stimulated by cytokines, which participated in many necessary biological processes, such as cell proliferation, differentiation, apoptosis, and immune regulation. The mechanism of JAK signal transmission is that binding of cytokines to the corresponding receptors leads to the dimerization of receptor molecules, making the receptor-coupled JAK kinases close to each other and activated by Interactive tyrosine phosphorylation. After activation of JAK, the tyrosine residues on the receptor are phosphorylated, and then these phosphorylated tyrosine sites form a “berthing site” (docking site) with the surrounding amino acid sequence, which also contains the SH2 domain recruited to this docking site. Finally the STAT protein which bound to the receptor is phosphorylated by kinase JAK, and the activated STAT protein enters in nucleus in the form of dimer to bind to the target gene and regulates gene transcription [43–45]. None of references pointed that JAK2 had any pivotal cellular activity to maintain the functions of the body, which elucidated the specificity of JAK2. Therefore, JAK2 is an indispensable therapeutic target without participating in any other pivotal cellular activity, which could be inhibited to prevent the JAK-STAT signaling pathway and thus prevent the occurrence of PMF.

At present, there is no effective treatment for PMF, namely there is no radical way to cure this disease [6]. Based on the anomalous activation of JAK-STAT signaling pathway in PMF, several studies have reported their discovery of JAK2 inhibitors [46, 47]. Unfortunately, current drug therapy in PMF including JAK2 inhibitor lacks disease-modifying activity. Despite the fact that great progress with inhibitors has been made regarding JAK2 in medication design and development, merely Fedratinib, which was selected as the reference drug in this study, has displayed a relatively mature research until now [9].

Though Fedratinib, which is currently being tested as a potential anti-cancer compound, has therapeutic potential in anti-neoplastic fields, it still had great therapeutic limitations. Relevant research showed that common adverse events with Fedratinib treatment were anemia, gastrointestinal symptoms, increased levels of serum creatinine, liver transaminases as well as pancreatic enzymes. Besides, Wernicke encephalopathy of unknown reason were found when treatment with Fedratinib [9]. Consequently, discovering more compounds targeting JAK2 is of profound significance for pharmacology as well as clinical application. Based on the previous outcomes in this study, structural biology was performed to screen natural inhibitors targeting JAK2 for further investigation.

In structural biology study, 17799 biogenic-for sale-named compounds were obtained from ZINC15 database for virtual screening, followed by ADME, TOPKAT, CDOCKER and Molecular Dynamics Simulation. Libdock module was firstly conducted in order to screen appropriate compounds which could dock with JAK2 through a fast docking analysis from the tremendous ligands. Libdock score represents the energy optimization and stability of the conformation. Compounds with a higher libdock score indicates a better energy optimization and a more stable conformation. Altogether, 10102 compounds were found eligible to bind with JAK2 stably after libdock module. Among those, 667 compounds were calculated to have higher libdock score than the reference ligand Fedratinib (libdock score: 129.056), suggesting that these 667 compounds may form a more stable conformation together with a better energy optimization compared to Fedratinib. Based on libdock score, the top 20 natural compounds were selected and pooled for following study.

Pharmacologic properties including ADME and toxicity were predicted to access these selected compounds. Taking all the properties as well as safeties into account, two compounds, ZINC000013513540 and ZINC000004099068, were identified as ideal lead compounds. Their properties of being soluble in water and easily absorbed by the intestine could alleviate the gastrointestinal symptoms compared to Fedratinib. Additionally, Fedratinib was predicted with hepatotoxicity while these two compounds were non-hepatotoxicity, which could reduce the damage to the liver. As for CYP2D6, cytochrome P450 (CYP450) is the principal enzyme involved in drug metabolism. Drugs behaved as inhibitors of CYP450 could weaken the activity of drug enzymes and thus slow down the drug metabolism. CYP2D6 is one of the enzymes in CYP450, and these compounds were predicted with non-CYP2D6 inhibition, suggesting that these two compounds could make the drug enzyme function greatest in drug metabolism. Furthermore, they had weak plasma protein binding property, the free part of ligand could behave the best pesticide effect, which demonstrated the good selectivity of the compounds regarding JAK2: compounds won’t be bound by other substances in the body such as plasma protein and thus reduced the combination with JAK2. However, only this molecule was analyzed, which was inadequacies in this study, it still remained unknown whether the drugs could be combined by other substances in the body, since different selectivity of drugs resulted in different dosage in future vitro/vivo experiment. Further study ought to be focused on the selectivity of the inhibitor. Moreover, these two compounds were accessed with non-Ames mutagenicity, low rodent carcinogenicity and less developmental toxicity potential compared to other compounds, which strongly suggested their perspective application in drug development. Although the remaining compounds on the list (Tables 5, 6) possessed relatively negative effects and toxicity, they still had potential application for the drug design and refinement by adding or removing specific groups and atoms to obtain better pharmacologic properties. Thus, ZINC000013513540 and ZINC000004099068 were finally selected as ideal lead compounds and were applied in further analysis.

Docking mechanism as well as chemical bond interaction of candidate compounds with JAK2 was then analyzed and visualized. CDOCKER module computation results demonstrated that these two compounds had the same or even higher binding affinity with JAK2 compared to Fedratinib. Subsequently, the chemical structures of the two compounds as well as Fedratinib with JAK2 were conducted with intuitive visualization (Figures 6, 7). The two compounds combined with JAK2 had more chemical bonds than Fedratinib did, which could maintain a more solid structure and thus prevent the function of JAK2 by contributing a competitive inhibition. In addition, results of pesticide effect (LD50 and LOAEL) of these two compounds indicated that they were relatively safe drugs as Fedratinib did, that they won’t have rapid adverse reactions due to high dosage of these drugs. Meanwhile, their properties of non-inhibition with liver drug enzyme CYP2D6 also proved the safety of these compounds. Consequently, the computational predictions could provide a reference for future experimental verification.

Lastly, molecular dynamics simulation was performed to verify their stability under natural circumstances. RMSD and energy values were determined as evaluation indexes. Results revealed that the RMSD trajectories of complexes reached equilibrium after 16, 18 ps, respectively. Energy curves got stabilized with time, hydrogen bond heat map illustrated that the hydrogen bonds, which contributed largely to the stability of complex, could exist stably with the progression of molecular dynamics. In conclusion, these two compounds, ZINC000013513540 and ZINC000004099068, could interact with JAK2 stably, and their complexes could exist under natural environment steadily, as well as inhibiting the function of JAK2 just the same as Fedratinib did.

Based on all the results above, drug development, including modification and refinement, could be prospectively conducted to make ligands and receptors bind more stably. It is also noteworthy that agonists and inhibitors always share a similar skeleton in the chemical structure [10, 11], adding or removing different groups or atoms may cause the opposite effects. Combining the advantages of pharmacologic properties, highly binding affinity as well as stabilization with JAK2 of these two compounds, they could provide a valuable resource for JAK2-related medication development.

Overall, in this study, from the identification of DEGs by bioinformatics to the study of inhibitors by structural biology, each step had been adequately explained. This study provided a guideline for preprocessing the GSE series’ data, as well as screening lead compounds which could have potential effects regarding receptor. The JAK-STAT signaling pathway and other relevant signaling pathway and the mechanism of the function of inhibitors identified in this study were illustrated in Figure 8. Firstly, we fully discussed the selection, downloading and preprocessing of GSE series’ data. Through a series of analytical methods, ten genes (LCN2, JAK2, MMP8, CAMP, DEFA4, LTF, MPO, HBD, STAT4, EBF1) were identified as hub genes responsible for the carcinogenesis of PMF. Of which, two genes, JAK2 and LCN2 were the most essential genes and fully discussed. Next, structural biology was performed to screen natural compound inhibitors regarding JAK2 based on the meager research status. Molecular conformation, pharmacologic properties, binding affinity as well as stability of each candidate compounds were fully investigated to verify their superiority compared to Fedratinib. ZINC000013513540 and ZINC000004099068, were finally selected as the ideal lead compounds, which may have potential effect in current chemotherapy of PMF. It is also worth noting that there is no single drug that could be marketed directly, unless through thousands of refinements, improvements, and finally clinical tests. Based on the two compounds selected in this study, further research could focus directly on the improvement and modification of them.

Although this study was conducted by elaborate design and precise methods have been carried out, it still had some limitations. Further experiments including animal test and cell experiment, were required to verify our results more firmly. Besides, more substances, like the plasma protein analyzed in this study, should be assessed in further studies in order to test whether these compounds were effectively selective inhibitor of JAK2.

Conclusions

This study combined two parts of bioinformatics and structural biology. First, after analysis of bioinformatics, this study found out two significant hub genes: JAK2 and LCN2, which were remarkably highly expressed and were responsible for the progression of PMF. Next, based on the JAK2, structural biology method was conducted, after a series of computer-aided structural and chemistry techniques (e.g., virtual screening, molecular docking, pharmacologic properties prediction, molecular dynamics simulation), two compounds, ZINC000013513540 and ZINC000004099068, were selected as safe and effective drug candidates, they were promising drugs in the treatment of patients with PMF, which could also contribute greatly to JAK2 inhibitors’ development and refinement.

Materials and Methods

Supplementary Materials

Supplementary Figures