Bone marrow samples were collected from seven patients with MPN as per guidelines and approval by Institutional Ethics Committee at St. John's Medical College and Hospital, Bengaluru (approval no. 103/2016). The median age of patients was 39 years (32–47 years). The samples were collected between June 2016 to July 2018. Patients >18 years old and classified according to the World Health Organization criteria for MPNs (20) were included in the study. CD34 cells were isolated from bone marrow samples using CD34 Microbead kit from Miltenyi Biotec, Inc., as per the manufacturer's instructions.

RNA was isolated from CD34 cells using TRIzol^® (Invitrogen; Thermo Fisher Scientific, Inc.). Whole transcriptome analysis of the seven samples was performed using the Whole Transcript (WT) PLUS Reagent kit (Applied Biosystems; Thermo Fisher Scientific, Inc.). The assay workflow in the order included first-stand cDNA synthesis, second-strand cDNA synthesis, cRNA amplification, cRNA purification and quantification, 2nd cycle ss-cDNA synthesis, template RNA removal, ss-cDNA purification and quantification, fragmentation, terminal labelling and hybridization to WT. All the steps were performed as per manufacturer's instructions. In total, 100 ng RNA was used as the input.

The quality check of CEL file was performed using Transcriptome Analysis Console software version 4.0.1 (Thermo Fisher Scientific, Inc.) and the files were normalized using the Oligo-R-based Bioconductor package (22). The processing was done using RMA method for R (23).

The normalized RMA files were checked for any discrepancies using the top ten housekeeping gene and their expression values across the sample. Principal Component Analysis was performed to cluster the samples. The samples that were clustered together based on the housekeeping genes were used for further analysis.

The normalized text file was then subjected to the AltAnalyze (http://www.altanalyze.org/) automated tool for analysis. Differentially expressed genes in the JAK2V617F negative vs. JAK2V617F positive samples were identified using ±1.5 as the fold-change cut-off with a statistical significance of <0.05. Prune ontology was calculated using the z score (initial filtering 1.96 cut-off). The database used for analysis was Ensembl (https://asia.ensembl.org/). The pathway network was generated automatically using the AltAnalyze tool.

Differentially expressed genes were calculated using ±1.5 fold-change and P≤0.05. Database for Annotation, Visualization and Integrated Discovery (DAVID) (https://david.ncifcrf.gov/) tool version 6.8 was used for Gene Ontology (GO), and functional and pathway enrichment analyses (24,25). The protein-protein interaction (PPI) network was constructed using STRING database version 11 (https://string-db.org/) using minimum confidence 0.4 (26). The genes which were enriched for DNA modification and stemness properties based on GO analysis were fed into StemChecker tool (http://stemchecker.sysbiolab.eu/) for analyzing stemness markers and transcription factors for enriched genes (27). Hierarchical clustering and heatmaps were created using the ClustVis online tool (28).

RNA was isolated from cells using TRIzol^® (Invitrogen; Thermo Fisher Scientific, Inc.). In total, 1 µg RNA was converted into cDNA using the MMLV enzyme (Invitrogen; Thermo Fisher Scientific, Inc.) as per the manufacturer's instructions. Gene expression was evaluated by using the TB Green Premix Ex Taq (Takara Bio, Inc.) as per manufacturer's instructions and qPCR setup was run on a 7500 Fast Real-Time PCR system (Applied Biosystems; Thermo Fisher Scientific, Inc.). GUSB was used as housekeeping gene. Primer sequences for selected genes are presented in Table I. The relative quantification was calculated using the 2^−ΔΔCq method (29). MPN mutation panel, including CALR type I and II and myeloproliferative leukemia virus (MPL) mutations (W515L, W515K, W515A and S505N), were analyzed using the TRUPCR MPN panel kit (3B BlackBio Biotech India Ltd.) as per manufacturer's instructions.

The bone marrow cells post red blood cell lysis (150 mM ammonium chloride, 10 mM potassium bicarbonate, 0.1 mM ethylenediaminetetraacetic acid) were labelled with CD34 FITC (clone 581 cat. no. 555821; BD Biosciences) at room temperature for 30 min in dark and analyzed using a FACS Calibur instrument (BD Biosciences). The results were analyzed using FCS4 Express software version 6.06.0025 (De Novo Software).

The statistical significance for CD34(+) cells was calculated using an unpaired student's t-test using SPSS v16 (IBM). P<0.05 was considered to indicate a statistically significant difference. All Q-PCR experiments were performed in triplicates (n=3). The differential expression of genes was calculated using one-way ANOVA and Benjamini and Hochberg false discovery rate (FDR).

The role of JAK2V617F mutation is well established in driving a subset of the myeloproliferative neoplasms (6–9); however, the etiology of the JAK2V617F negative MPNs is unknown. To understand this, JAK2V617F mutant positive and negative MPNs were compared using microarray-based transcriptional analysis. The hematopoietic cancer stem cells are determinant of the cellular hierarchy of tumor progression in hematological malignancies; however, the underlying mechanism is poorly understood in most cases, except CML wherein the role of translocation and resultant chimeric protein is a driving force in hematopoietic stem cells (HSCs) (30). To understand the changes that may be driving the JAK2V617F negative MPNs, bone marrow derived CD34(+) cells were targeted in the present study. CD34(+) fractions were separated from patient samples and subjected to microarray analysis. Tables II and III present the baseline characteristics and mutation profile of the patients involved in the present study. Patients who were JAK2V617F negative were also negative for CALR type 1 and 2 mutations and MPL mutations (W515L, W515K, W515A and S505N) (Table III).

The microarray analyses of CD34(+) cells isolated from JAK2V617F negative and positive samples resulted in 21,448 gene signatures. The analysis resulted in 472 upregulated and 202 downregulated genes. The maximally downregulated genes, with >3-fold difference were 3-oxo-5-α-steroid 4-dehydrogenase 1 (SRD5A1), matrix metalloproteinase (MMP7), Wilms tumor protein (WT1) and Ankyrin repeat and death domain-containing protein 1B (ANKDD1B), while the maximally upregulated genes with >3-fold difference were X-ray repair cross-complementing protein (XRCC6), TUBB, histone-lysine N-methyltransferase SUV39H1 and pantetheinase VNN1 as shown in Table IV. The expression of SUV39H1 and MYB was validated in five samples used for microarray analysis as shown in Fig. S1 where in the expression of both the genes was higher in JAK2V617F negative samples compared with control samples.

Hierarchical clustering of significantly upregulated and downregulated genes was performed to understand the correlation between samples (Fig. 1). Hierarchical clustering separated seven samples into two different groups consistent with the experimental design. Red bars indicated a high level of expression, blue bars indicated a low level of expression and white bars indicated no significant difference (Fig. 1). There was a distinctive difference in the expression of genes in JAK2V617F positive and negative samples (Fig. 1). The expression pattern of genes was homogeneous within the mutation negative samples; however, there was heterogeneity within mutation positive samples.

To understand the important biological processes that determine the phenotypes in these two subsets of MPNs, GO was performed. The enriched genes and common pathways were analyzed using DAVID (Fig. 2). The distribution of the transcripts in the GO terms revealed a distinct set of biological processes (BP) (Fig. 2A), including ‘mitotic nuclear division’, ‘DNA replication’ and genes involved in ‘cell division, migration and proliferation’. Notably, all the BPs that were identified were indicative of a DNA-related processes, such as histone exchange and G₂/M transition. The cellular components category also indicated an enrichment of genes associated with the nuclear compartment, such as ‘nucleus’, ‘nucleoplasm’ and ‘nuclear matrix’ (Fig. 2B), suggesting alterations in events occurring in the nuclear compartment. The most significant GO terms describing molecular functions for the enriched genes were ‘protein binding’, ‘chromatin binding’ and ‘ATP binding’ (Fig. 2C).

Following the aforementioned observations, the interactions between the enriched proteins were analyzed using the STRING database to understand the important molecular processes that may be activated in the JAK2V617F (−) MPNs. The PPI network constructed for the upregulated genes is shown in Fig. 3A. Notably, the STRING analysis revealed the existence of a distinct network of genes represented by tight network clusters. A detailed analysis of the network indicated the involvement of genes in two main biological processes: Cell cycle (red) and nucleobase modification (blue) (Fig. 3A). The central tight network cluster included genes involved in the cell cycle and its regulation, such as CDK1, TUBB, G2 and S phase-expressed protein (GTSE1), borealin (CDCA8), SUV39H1 and centrin 2 (CETN2) (Fig. 3A). Nucleobase modifying proteins, such as nuclear pore complex protein NUP85, NUP62, MYB, SS18, nucleolar and coiled-body phosphoprotein NOLC1, SWI/SNF complex subunit SMARCC2 and SMARCE1, were also linked in the network (Fig. 3A). SUV39H1 encodes a histone methyltransferase that tri-methylates lysine 9 of histone H3, which leads to transcriptional gene silencing (31). SMARCC1 and SMARCE1 are part of the SWI/SNF chromatin remodeling complex involved in transcriptional activation of target genes in an ATP-dependent manner (32). The majority of the proteins in the network were involved in acetylation and phosphorylation (Fig. 3A). Proteins such as hypoxia up-regulated protein HYOU1, Isoleucine-tRNA ligase (IARS), tyrosine t-RNA synthetase YARS, dual specificity protein kinase TTK, acyl-coenzyme A synthetase ACSM3, CDK1, XRCC and carbamoyl-phosphate synthetase CAD marked for nucleotide binding and acetylation indicating towards the role of DNA modification in this subset of MPNs. However, unlike other hematological malignancies, such as acute myeloid leukemia (AML) (33), little information is available on epigenetic modulations and their role in regulating hematopoiesis and fibrosis in MPNs. A heatmap was constructed to visualize the differences in gene expression between samples for genes comprising the aforementioned tight cluster (Fig. 3B). There was a categorical difference between samples with the red bars indicating a high expression of genes in mutant negative samples and blue bars indicating decreased expression. A bar graph representing the fold-change between mutant positive and negative samples shows the significant relative difference in expression of selected few genes involved in cell cycle and nucleotide modification, such as MYB, SUV39H1, XRCC, SMARCC2, SMARCCE1 and CDK1 (Fig. 3C).

A network analysis was also performed for downregulated genes, which highlighted the genes involved in developmental process, response to stimuli, multicellular organismal process and regulation of apoptosis (Fig. S2). Key genes involved in this network included muscarinic acetylcholine receptor (CHRM3), IL6, Kalirin (KALRN), MMP7 and ankyrin repeat domain-containing protein ANKRD1.

Analyses using the DAVID and STRING databases based on GO and network enrichment for the upregulated genes suggested an association between the significantly upregulated genes and stem-like properties. Enriched genes which clustered into DNA modification and stemness properties (based on DAVID GO analysis) were further analyzed using StemChecker for evaluating the stemness signatures and transcription factors associated with these genes. Seven genes (12%) were linked to the hematopoietic signature (Fig. 4A) while forty-seven (76%) genes were related to the embryonic signature (Fig. 4A). The abilities of self-renewal, maintenance and differentiation of stem cells serve as a core reservoir of cancer initiation and development and tumor growth (34). Elevated numbers of embryonic or lineage-specific stem cells can lead to deregulation and altered differentiation of the progeny (34). The expansion of one specific lineage in mutant negative MPNs might be attributed to an altered hematopoietic stem cell. Notably, 76% of genes matched with embryonic stem cells (ESC) (Fig. 4A). ESC belongs to the primitive class of stem cells and are capable of differentiating into various cell types (35). Owing to their diverse nature, even a minute perturbation in these cells can lead to notable changes in cellular physiology (35). The transcription factors (TFs) associated with stem cells were also analyzed (Fig. 4B). Overexpression of stem-cell specific TFs may contribute to the pathological self-renewal characteristics of cancer stem cells (36). The majority of the genes mapped to E2F4 and NANOG TF studies by Boyer et al (37) (Figs. 4B and S3) which are involved in cell cycle and pluripotency respectively (38,39). To further evaluate the correlation between upregulated genes and stemness and chromatin modification, a network was constructed using STRING (Fig. 5A). The central tight cluster mainly comprised of genes implicated in the cell cycle, nucleobase modification and chromatin modification (Fig. 5A). These genes included MYB, tubulin γ, probable ATP-dependent RNA helicase DDX41, TUBB, IARS, SS18, CDCA8 and TTK (Fig. 5A). The heat map representing the relative expression of genes involved in the cluster indicated a higher expression of these genes in mutant negative samples (Fig. 5B). There was heterogeneity in mutant positive samples, but mutant negative samples had a homogenous profile. The bar graph representing the fold-changes of a select few key chromatin modifiers from the PPI network shows the significant relative difference in expression between mutant positive and negative samples such as chromatin modification protein (YEATS2), MYST-associated factor (MEAF6), methyl CpG binding protein (MBD1), nuclear receptor corepressor (NCOR2), PHD finger protein (PHF10) and SS18 (Fig. 5C). SUV39H1 is a well-known epigenetic marker and has an established role in MPN (40); however, the role of MYB is unknown in this context. Therefore, the expression of MYB was evaluated in JAK2V617F negative samples. There was a ~10-fold increase in MYB expression in the patient samples compared with controls (Fig. 6).

On the basis of the aforementioned observations, the stem cells in patient samples were further evaluated. CD34 was selected as the marker for interest as it is the most characterized and reported HSC marker shown to be involved in various hematological malignancies (41). Using flow cytometry, it was demonstrated that there was a significant increase in CD34 counts in JAK2V617F negative samples compared with positive samples (P<0.05; Fig. 7C), which concurred with the aforementioned bioinformatic analysis.