Analysis of genomic aberrations and gene expression profiling identifies novel lesions and pathways in myeloproliferative neoplasms.

Polycythemia vera (PV), essential thrombocythemia and primary myelofibrosis, are myeloproliferative neoplasms (MPNs) with distinct clinical features and are associated with the JAK2V617F mutation. To identify genomic anomalies involved in the pathogenesis of these disorders, we profiled 87 MPN patients using Affymetrix 250K single-nucleotide polymorphism (SNP) arrays. Aberrations affecting chr9 were the most frequently observed and included 9pLOH (n=16), trisomy 9 (n=6) and amplifications of 9p13.3-23.3 (n=1), 9q33.1-34.13 (n=1) and 9q34.13 (n=6). Patients with trisomy 9 were associated with elevated JAK2V617F mutant allele burden, suggesting that gain of chr9 represents an alternative mechanism for increasing JAK2V617F dosage. Gene expression profiling of patients with and without chr9 abnormalities (+9, 9pLOH), identified genes potentially involved in disease pathogenesis including JAK2, STAT5B and MAPK14. We also observed recurrent gains of 1p36.31-36.33 (n=6), 17q21.2-q21.31 (n=5) and 17q25.1-25.3 (n=5) and deletions affecting 18p11.31-11.32 (n=8). Combined SNP and gene expression analysis identified aberrations affecting components of a non-canonical PRC2 complex (EZH1, SUZ12 and JARID2) and genes comprising a 'HSC signature' (MLLT3, SMARCA2 and PBX1). We show that NFIB, which is amplified in 7/87 MPN patients and upregulated in PV CD34+ cells, protects cells from apoptosis induced by cytokine withdrawal.

Introduction

The Ph-negative myeloproliferative neoplasms (MPNs) including polycythemia vera (PV), essential thrombocythemia (ET) and primary myelofibrosis (PMF) are clonal hematopoietic disorders characterized by the overproduction of blood cells belonging to the myeloid, megakaryocytic and/or erythroid lineages.[1] Although clinically distinct, these diseases are characterized by a recurrent mutation, JAK2V617F, which is present in ∼95% of patients with PV, ∼65% with PMF and ∼55% with ET.[2] JAK2 is one member of the Janus family of non-receptor tyrosine kinases that transduces signals downstream of type I and II cytokine receptors via signal transducer and activators of transcription (STATs), and is an important regulator of cell survival, proliferation, differentiation and apoptosis.[3, 4] The JAK2V617F mutation leads to autophosphorylation of JAK2 and constitutive activation of downstream pathways including the Janus kinase (JAK)/STAT, phosphatidylinositol 3 kinase/AKT (v-akt murine thymoma viral oncogene homolog 1) and mitogen-activated protein kinase (MAPK)/extracellular signal-regulated kinase pathways.[5, 6, 7, 8, 9] Although overexpression of JAK2V617F in a murine bone marrow transplant models leads to a PV-like disease and megakaryocytic hyperplasia, such mice do not develop thrombocytosis and questions regarding the effect of allele dosage on disease phenotype remain unanswered. The finding of Theocharides et al.[10] that 70% of JAK2V617F-positive MPN patients went on to develop JAK2V617F-negative acute myeloid leukemia (AML), suggests that additional predisposing factors may precede the JAK2 mutation. Recently, several acquired inactivating mutations in TET2,[11] ASXL1(ref. 12) and IDH1/2(ref. in MPN progenitors were identified, and it is likely that additional genes are involved in the initiation of a pre-leukemic clone or as disease modifiers, in co-operation with JAK2V617F.

Previous studies using standard karyotyping revealed abnormalities in 50% of PMF patients, 15% of PV and <5% of ET patients.[14, 15, 16] Although such conventional chromosome analysis has proven useful for diagnosis and prognostic information in PMF patients and other myeloid malignancies, the sensitivity of these approaches is limited. Indeed, studies by Tefferi using comparative genomic hybridization (CGH) oligonucleotide arrays to identify genomic aberrations in a panel of MPN patients, detected copy number aberrations (CNAs) in 29% and 18% of patients with PV and ET, respectively, while cytogenetics for these patients appeared normal.[17] Recently, high-resolution single-nucleotide polymorphism (SNP) arrays have emerged as a powerful, high-throughput technology to allow detection of micro-deletions and duplications and inference of loss of heterozygosity (LOH) from genotype calls, in cases where conventional cytogenetic technologies and array CGH are unsuccessful. Since then, SNP arrays have been successfully used to identify cryptic oncogenic deletions and duplications affecting putative oncogenes and tumor suppressors in various hematological malignancies. For example, the identification of abnormalities on chromosome 4q24 affecting the TET2 gene in patients with MDS, MPN and AML was facilitated by use of SNP array analysis.[11, 18] Similarly, PAX5 was identified as a target for mutation in B-cell acute lymphocytic leukemia using SNP arrays.[19] Recent studies using high-resolution SNP arrays identified novel, recurrent aberrations in a large cohort of MPN patients and revealed associations between specific chromosomal anomalies and disease progression.[20]

Here, we report an integrative analysis of chromosomal abnormalities and gene expression changes in MPN patients. We profiled CNAs in peripheral blood granulocytes of 87 patients with MPN and identify new regions affected by recurrent chromosomal aberrations. In a subset of the profiled patients, genomic copy number changes could be correlated with gene expression in granulocytes, leading to the identification of novel genes deregulated in MPN. We show that patients affected by chromosome 9 abnormalities, which represent the most frequent CNA, are characterized by deregulation of set of genes including JAK2, STAT5B and MAPK14, suggesting that these genes may be involved in disease pathogenesis. We also performed meta-analysis of five independent MPN expression data sets to identify genes commonly deregulated in MPN. Comparison of genes identified in these meta-analysis with those affected by CNAs, revealed genes that may have important roles in disease pathogenesis including EZH1, MLLT3, SOCS3 and SMARCA2. Finally, we show that the transcription factor NFIB, which is located adjacent to JAK2 and upregulated in our MPN meta-analysis, leads to increased net cell growth and decreased apoptosis in the murine hematopoietic cell line, Ba/F3-EPOR.

Patients and methods

Patient samples

A total of 87 MPN patients (32 PV, 24 PMF and 31 ET) were included in this study (Supplementary Table S1). Samples were obtained with informed consent at the Dana–Farber Cancer Institute as described previously.[8]

SNP arrays

DNA was extracted from peripheral blood granulocytes using the QIAamp DNA Blood Maxi Kit (Qiagen, Valencia, CA, USA) and analyzed with GeneChip Mapping 250K arrays according to the manufacturer's protocol (Affymetrix, Santa Clara, CA, USA).[11] Briefly, 250 ng of genomic DNA was digested by Nsp restriction enzymes and ligated with the appropriate adaptor for amplification with Platinum Pfx DNA polymerase (Invitrogen, Carlsbad, CA, USA). Three 100 μl PCR reactions were carried out for each adaptor-ligated DNA sample. PCR products were pooled and purified. In all, 40 μg of PCR product was fragmented with DNase I and visualized on a 4% TBE agarose gel. Fragmented PCR products were subsequently end-labeled with biotin. Hybridization and washing was performed using the Affymetrix Fluidics Station 450 and Gene chip and chips were scanned on an Affymetrix GeneChip Scanner G7 (Affymetrix).

JAK2V617F allele burden analysis

The JAK2V617F allele status was assessed using a quantitative reverse transcriptase-PCR assay as previously described.[21] In brief, a standard curve for the proportion of JAK2V617F/total JAK2 against the dCT (Ct-JAK2V617F-CtJAK2WT) was generated using DNA from a healthy control and a PV erythroid colony homozygous for the JAK2V617F mutation. This method was then used to determine the JAK2V617F allele burden (%) in DNA extracted from individual patient granulocytes.

Data analysis

CNA analysis

These data were processed and analyzed using a custom developed procedure (Supplementary Figure S1). Genotype calls were generated by GTYPE (Affymetrix) from CEL files and 250K library files for Nsp arrays. Raw data were normalized using invariant set normalization method. The intensity value for each SNP was extracted using a model-based perfect match/mismatch algorithm in dChip (http://biosun1.harvard.edu/complab/dchip/). The log2 ratio of the intensity for each SNP for each sample relative to data from 60 unrelated CEU subjects of HapMap project (http://www.hapmap.org) was calculated. The circular binary segment algorithm implemented in BioConductor DNAcopy package (http://bioconductor.org) was applied to the log2 ratio data to segment the chromosome into contiguous regions of equal copy number based on a maximal t-statistics with the use of a permutation reference distribution (using coordinates of the SNP based on Genome Build HG-17).[22] Segmental value was calculated by averaging copy number change values for the SNPs residing in that segment. The baseline value of genome-wide segments in the sample with a bias in regional changes (in case of the samples with more regional gains than losses, or vice versa or with only regional gains or losses) may be shifted because of the preceding steps of global normalization and median centering. To correct the baseline, the segmental values were adjusted by subtracting the median of segmental change values for all the SNPs. Considering a large segment with low-level copy number change as important as a small segment with high-level change, the adjusted segmental change value was further weighted by multiplying log10 of the length of that segment. A segment in the MPN sample was considered significant if its segmental value is >3 s.d. of segmental values in a range of Q1−3*(Q2−Q1) to Q2+3*(Q2−Q1) where Q1 or Q2 are the 25% or 75% percentile of segmental values in all the samples, respectively. The significant segments having an overlap of >75% with CNA regions in normal subjects reported in the Database of Genomic Variations (http://projects.tcag.ca/variation/) were filtered out for further analysis. The frequencies of copy number changes in either direction for significant segments among MPN samples were summarized for each disease using the boundary information of significant segments. To assess the significance of genomic aberrations in a group of samples, the permutation was performed to assess the frequency of genomic aberrations occurring in the genome by a random chance. The genome and aberrant regions were fragmented into 50 kb blocks and the fragmented blocks for each sample were randomly assigned to any block in the genome. The frequency of aberrations occurring in every block in the genome at random was calculated and compared with the actual frequency of aberrations. These steps were repeated 1000 times. If the random chance of obtaining a proportion higher than actual frequency of genomic aberrations for a given region is <5% (equivalently to 50 times from 1000 permutations), the recurrent aberration for a region was determined to be significant. The Fisher's exact test was performed to examine whether a particular chromosomal aberration was associated with JAK2 mutation status or disease type. To determine if the genomic instability profiles were segregated by disease type, hierarchical clustering of genome-wide copy number change profiles was performed in PV, ET and PMF patients using the distance of Pearson correlation co-efficiencies.

LOH analysis

LOH analysis was performed comparing MPN genotypes with 60 CEU unrelated normal controls derived from the HapMap project (http://www.hapmap.org). A Hidden Markov-based algorithm was used to infer LOH from unpaired samples implemented dChip.[23] The regions with a minimum of 80% probability of LOH call were identified for each sample and LOH regions for MPN and HapMap normal samples were then summarized and compared.

RNA extraction and microarray

Total mRNA was extracted from peripheral blood granulocytes of 38 selected MPN patients and 11 normal control patients using Trizol (Invitrogen). In all, 10 μg was hybridized to HGU133A V2 arrays (Affymetrix). Data were analyzed using Genespring 11.0.2 (Agilent Technologies Inc., Santa Clara, CA, USA) following RMA normalization and genes with a fold change of ⩾1.5 (P<0.05) were considered statistically significant. For MPN meta-analysis, a one-sided t-test (P<0.05) was performed between disease and control samples for data sets 1, 2, 3 and 5 (PV data set 1: CD34+ cells from 10 PV and 9 normal controls,[24] PMF data set 2: CD34+ cells from 3 PMF and 3 normal controls,[25] ET data set 3: platelets from 6 ET and 5 normal controls,[26] ET data set 4: CD34+ derived megakaryocyte cells from 2 ET and 2 normal controls,[27] MPN data set 5: granulocytes from 93 MPN patients and 11 normal controls, unpublished data, includes 38 MPN patients described above). For data set 4, an average of log ratios of two samples (⩾2-fold) was used due to limited sample size. Genes deregulated in the same direction (average of ⩾1.5-fold and present in at least 3/5 data sets) were considered significant. For analysis of gene expression between patients with chromosome 9 aberrations and normal chromosome 9 cytogenetics, a one-sided t-test was performed and genes with a fold change ⩾1.5 (P<0.05) were considered significant.

Bioinformatic analysis

Gene biofunctions and network analysis was performed using Ingenuity Pathway Analysis (Ingenuity Systems, Redwood City, CA, USA).

Cell culture

293T cells were grown in Dulbecco's modified Eagle's medium containing 10% (v/v) fetal bovine serum (FBS), 50 U/ml penicillin and 50 μg/ml streptomycin. Ba/F3 cells expressing the murine erythropoietin receptor (Ba/F3-EPOR) were previously described[28] and cultured in RPMI containing 10% (v/v) FBS supplemented with recombinant erythropoietin hEPO (1 U/ml) (Cell Sciences, Canton, MA, USA). HEL cells were grown in RPMI containing 10% (v/v) FBS. TF-1 cells were cultured in RPMI containing 10% (v/v) FBS supplemented with human granulocyte-macrophage colony-stimulating factor (2 ng/ml) (Peprotech, Rocky Hill, NJ, USA). UKE-1 cells (Gift of Gabi Vohwinkel, University Hospital Eppendorf, Hamburg, Germany) were grown in Iscove's modification of Dulbecco's medium containing 10% (v/v) FBS, 10% (v/v) horse serum and hydrocortisone (10−6 ) (Sigma-Aldrich, St Louis, MO, USA). For JAK2 inhibition studies, HEL and UKE-1 cells were treated with either JAK inhibitor I (Calbiochem, Gibbstown, NJ, USA) or NVP-BSK805 Inhibitor (1–2 μ) (Novartis, Schaumberg, IL, USA), respectively, for 16 h, using dimethylsulphoxide as a control.

Viral constructs and transduction

Full-length human NFIB complementary DNA (cDNA) was obtained from Open Biosystems (Huntsville, AL, USA) and cloned into the two promoter EF1α-cytomegalovirus-green fluorescent protein (GFP) lentiviral vector (Gift of Dr Linzhao Cheng, John Hopkins). Full-length human JAK2WT and JAK2V617F cDNA was subcloned from the MSCV-Neo/Puro retroviral vector[8] into the two promoter A2-GFP lentiviral vector. Viral particles were produced by transient transfection of 293T cells with pRSV-REV (Rous sarcoma virus-REV), pMDLg/pRRE (gag/pol elements) and pMD.G (envelope elements) plasmids with Fugene6:DNA (Roche, Mannheim, Germany). After 24 h, transfection reagent was removed and replaced with RPMI and incubated at 32 °C for 48 h. Supernatants were pooled and filtered through a 0.45μm syringe filter for immediate use. For transduction of Ba/F3-EPOR or TF-1 cells, 2 × 106 cells were pelleted and resuspended in 2 ml of lentivirus with polybrene (4 μg/ml) and spinoculated at 1010 × g for 2 h at 25 °C. Infected cells were then cultured for 48 h and the percentage of GFP was determined by flow cytometry.

cDNA preparation and quantitative reverse transcriptase- PCR

cDNA was synthesized from 1 μg of total RNA using Iscript cDNA Synthesis Kit (Bio-Rad, Hercules, CA, USA). Quantitative reverse transcriptase-PCR reactions were performed using Quantitect SYBR Green PCR Kit (Qiagen) using the Mx3000P Real-Time PCR System (Stratagene, Agilent Technologies, La Jolla, CA, USA). Quantitect primers (Qiagen) used for quantitative quantitative reverse transcriptase-PCR are located in Supplementary Table S2.

Cytokine withdrawal

Sorted GFP+ Ba/F3-EPOR cells were cultured in RPMI plus erythropoeitin (EPO) (1 U/ml) and cumulative cell growth counts were performed over 10 days by Trypan blue exclusion. For EPO withdrawal studies, cells were washed three times in phosphate-buffered saline and seeded at 1 × 105 cells/ml with varying concentrations of EPO (1, 0.1 and 0.01 U/ml) for 2 days. The percentage of dead cells was calculated by Trypan blue exclusion.

Results

Identification of novel MPN-associated disease genes using a combination of SNP and gene expression profiling

To identify genes recurrently affected by genomic aberrations in MPN, we hybridized DNA from peripheral blood granulocytes of 87 MPN patients (32 PV, 24 PMF and 31 ET) to Affymetrix 250K SNP arrays. Analysis of genomic profiles for CNAs and copy neutral LOH/uniparental disomy confirmed previously reported chromosomal abnormalities in MPN, including recurrent gain of 1q (n=4), trisomy 8 (n=2), 9pLOH (n=16), trisomy 9 (n=6) and deletions in 20q11.21–13.13 (n=4)[20, 29] (Figure 1a; Supplementary Table S1). We also found that several genes mutated in MPN were also affected by genomic aberrations in our cohort. For example, analysis of the TET2 gene, which is mutated in ∼12% of patients with MPN,[11] revealed small deletions at chr4q24 affecting TET2 in two PV patients (Supplementary Figure S1). We also observed loss of 20q11.21–q13.12 in one PV (s328) and one PMF patient (s212), a region including the ASXL1 gene, which was recently shown to be mutated in ∼11% of myelodysplastic syndromes (MDS) and in ∼7% of patients with MPN ET and PMF,[30] gain of ch11q23.3 encompassing CBL in one ET patient (s403) and loss of chr7q35 harboring EZH2 in one PMF patient (s212) (Supplementary Table S1).

Figure 1

CNAs in MPN (a). Ideogram representing regions of chromosomal gains (red) and loss (green) identified using Affymetrix GeneChip Mapping 250K arrays in 87 patients with MPN. Individual lines represent a single patient (b). Top: genomic profile of chromosome 17 in one patient (s531) showing a gain at 17q21.1–21.31 encompassing 144 genes including EZH1 and ITGA2B, which are also transcriptionally upregulated in MPN patients (Supplementary Table S6). Bottom: genomic profile of chromosome 6 in one patient (s242) showing a gain at 6p22.3 encompassing JAR1D2. Red lines represent the mean log2 ratio of the intensity of the samples relative to 60 unrelated HapMap normal controls.

Overall, our analysis identified CNAs in 88% of MPN patients, (77/87), however, the majority (∼90%) of recurrent aberrations occurred at a frequency of <5%. The most frequent abnormalities involved chromosome 9, including 9pLOH (n=16) as well as gains of whole chromosome 9 (n=6) and partial amplifications including 9p13.3–23.3 (n=1), 9q33.1–34.13 (n=1) and 9q34.13 (n=6). Recurrent gains of 1p36.31–36.33 (n=6), 17q21.2–q21.31 (n=5), 17q25.1–25.3 (n=5) and deletions of 18p11.31–11.32 (n=8) and 15q15.3 (n=7) were also observed (Table 1). In addition to these larger aberrations, a subset of CNAs identified in our screen involved smaller regions of recurrent gain or loss, typically containing <10 genes. We hypothesized that these areas may contain putative oncogenes or tumor suppressors, as in the case of the TET2 tumor suppressor. Such candidates include FLI1 and ETS1, located in a 754–841-kb region of loss (3/87), TCBA1 located in a 600-bp region of loss (3/87) and CNTNAP2 located in a 35-kb region of loss (2/87) (Supplementary Table S1).

Table 1

SNP arrays identify regions of recurrent gain or loss in MPN patients

Chr	Ab	Start	End	Dist (kb)	Freq	COSMIC database	MPN vs control (>1.5-fold, P<0.05)	MPN meta-analysis (>1.5-fold, P<0.05)
1	Gain	1p36.33	1p36.31	4338	6/87	PRDM16	LRRC47	—
1	Gain	1q32.1	1q44	45 208	4/87	ELK4, MDM4, SLC45A3	ADSS, ANGEL2, ARF1, ARID4B, BTG2, CNIH4, FAIM3, FBXO28, HEATR1, LGALS8, NUAK2, PFKFB2, RAB3GAP2, RCOR3, TOMM20, TRAF3IP3, TSNAX, ZC3H11A	AKT3, GALNT2, LGALS8, MTR, PSEN2, RAB3GAP2, RPS6KC1, SMYD2
8	Gain	8p23.1	8p23.1	428	4/87	—	—	—
9	Gain	9p24.3	9p23	10 489	6/87*	JAK2	JAK2, SMARCA2	JAK2, SMARCA2
9	Gain	9p23	9p13.3	16 366	7/87*	CDKN2A, MLLT3, NFIB	—	ADFP, MLLT3, NFIB
9	Gain	9p13.3	9q33.1	83 625	6/87*	FANCC, FANCG, GNAQ, NR4A3, OMD, PAX5, PTCH, SYK, TAL2 XPA	AGTPBP1, ANXA1, CTNNAL1, DCTN3, DDX58 KLF4, NDUFB6, NOL8 NR4A3, SMU1, SPTLC1 T DRD7, TLN1	ABCA1, CDC14B, DAPK1 DDX58, KIAA0367, NTRK2, TMOD1, TPM2, TXNDC4
9	Gain	9q33.1	9q34.13	10 771	7/87*	CEP1	CDK5RAP2, FBXW2, GAPVD1, GOLGA1, PBX3 POLE3, PRPF4, PTGS1, RABGAP1, STOM	PTGS1, RABGAP1
9	Gain	9q34.13	9q34.13	1097	10/87*	—	—	ANGPTL2
9	Gain	9q34.13		11 137	12/87*	ABL1, BRD3, FNBP1, NOTCH1, NUP214, SETTSC1	ASB6, CDK9, GPR107, LCN2, RALGDS SET, SPTAN1, SSNA1	—
12	Gain	12p13.33	12p13.33	1008	4/87	—	—	—
15	Loss	15q15.3	15q15.3	172	7/87	—	—	—
17	Gain	17q21.2	17q21.31	3984	5/87	BRCA1, ETV4	ATP6V0A1, EZH1, FMNL1 HEXIM1, ITGA2B, MAP3K14, PSME3	EZH1, ITGA2B
17	Gain	17q25.1	17q25.3	9770	5/87	ASPSCR1, CANT1	GGA3, MAFG, SEC14L1, SGSH, SLC16A3 ST6GALNAC2,WBP2	BAIAP2, EXOC7, SEC14L1 SOCS3
18	Loss	18p11.32	18p11.31	813	8/87	—	—	—
20	Loss	20q11.23	20q13.12	9917	4/87	MAFB, TOP1	BLCAP, BPI, CTNNBL1, MMP9, SFRS6	C20orf67, EYA2, HNF4A L3MBTL, SRC, TGM2, UBE2C
20	Gain	20q13.33	20q13.33	3353	4/87	SS18L1	LIME1, PCMTD2, RPS21	DATF1

Abbreviations: CNA, copy number aberration; COSMIC, Catalogue of Somatic Mutations in Cancer; ET, essential thrombocythemia; MPN, myeloproliferative neoplasm; PMF, primary myelofibrosis; PV, polycythemia vera; SNP, single-nucleotide polymorphism.

Frequently occurring CNAs identified using 250K SNP Array Analysis in MPN patients (n=87). Asterisks indicate that these frequencies include six cases of whole chromosome 9 gain. Genes affected by CNA, somatic mutations in cancer (COSMIC) and gene expression analysis (MPN vs normal; genes deregulated in peripheral blood granulocytes from 37/87 MPN patients (10 PV, 17 ET and 10 PMF) compared with 11 normal controls (⩾1.5-fold, P<0.05; Supplementary Table S4) and meta-analysis of five MPN array data sets (one-sided t-test, P<0.05, average of 1.5-fold in 3/5 data sets; Supplementary Table S6) are shown in the far right columns.

In order to further investigate the disease potential of these genomic gains and losses, we compared genes affected by CNAs in our MPN cohort (⩾2 patients; 4448 genes) (Supplementary Table S3) with genes deregulated in peripheral blood granulocytes from a subset of these patients (n=37, 10 PV, 17 ET and 10 PMF, Supplementary Table S1) compared with normal controls (n=11) (⩾1.5-fold, P<0.05; 1447 genes, Supplementary Table S4) and cancer-associated genes obtained from the Catalogue of Somatic Mutations in Cancer (COSMIC) (426 genes) (Supplementary Figure S2). Of the 426 COSMIC genes, 75 were affected by CNAs and 42 were deregulated in expression (Supplementary Table S5). Nine genes fulfilled all criteria and these included JAK2 and Translocated Promoter Region, which were both affected by copy number gains and overexpressed in MPN patients compared with normal controls. Using this tri-level analysis, we identified aberrations affecting components of a non-canonical PRC2 complex including upregulation of SUZ12 (1.8-fold, P<0.05), amplification and upregulation of EZH1 (5/87; 1.6-fold, P<0.05) and amplification of JARID2 (2/87) (Figure 1b). Several genes constituting part of a HSC ‘gene signature'[31] were also identified in our analysis, including MLLT3 and SMARCA2, amplified as a result of trisomy 9 in six patients and PBX1 amplified in three patients on chr1q22–1q25.3, strengthening the hypothesis that aberrant self-renewal may underlie a common MPN progenitor (Table 1).

Finally, we compared genes affected by CNAs with those deregulated in a comprehensive MPN meta-analysis. Analysis of five individual MPN data sets comparing patient and normal control samples identified 574 upregulated genes and 636 downregulated genes. Comparison of this data set with genes affected by genomic aberrations in our MPN patient cohort (CNAs present in ⩾2 patients), identified 81 upregulated genes in regions of chromosome gain and 20 downregulated genes in regions of loss (Supplementary Table S6). Ingenuity pathway analysis of these genes revealed an enrichment of genes belonging to various signaling pathways including the JAK/STAT, interleukin-8 and erythropoietin pathways, reinforcing the role of dysregulated JAK/STAT signaling in MPN (Table 2). Furthermore, genes previously implicated in hematological disease such as PBX1, PDE4C, PTGS1 and SRC as well as in cell growth and proliferation and apoptosis were also enriched (Table 3). Taken together, these data indicate that genomic aberrations present in MPN account for a proportion of gene expression changes, which may have a role in disease pathogenesis.

Table 2

MPN-associated pathways

Pathway	P-value	Genes	Aberration	Fold change
JAK/STAT signaling	3.4 × 10−4	AKT3	Gain (3/87)	1.8
		JAK2	Gain (6/87)	1.9
		SOCS3	Gain (6/87)	1.9
		JAK3	Gain (2/87)	1.5
IL-8 signaling	4.4 × 10−4	AKT3	Gain (3/87)	1.8
		ANGPT1	Gain (2/87)	2.8
		AZU1	Gain (3/87)	1.5
		IKBKB	Gain (2/87)	1.6
		PTK2B	Gain (2/87)	1.8
		SRC	Loss (4/87)	−1.6
Erythropoietin signaling	5.5 × 10−4	AKT3	Gain (3/87)	1.8
		JAK2	Gain (6/87)	1.9
		SOCS3	Gain (6/87)	1.9
		SRC	Loss (4/87)	−1.6

Abbreviations: IL, interleukin; JAK, Janus kinase; MPN, myeloproliferative neoplasm; SNP, single-nucleotide polymorphism; SOCS, suppressor of cytokine signaling; STAT, signal transducer and activators of transcription.

Ingenuity pathway analysis was performed using genes commonly identified by SNP array (n=87) and MPN meta-analysis (n=114) (89 genes).

Table 3

MPN-associated bio-functions

Functions	P-value	Genes	Aberration	Fold change
Hematological disease (4.7 × 10−5–1.2 × 10−2)
Erythrocytosis	7.6 × 10−4	JAK2	Gain (6/87)	1.9
		SOCS3	Gain (6/87)	1.9
		SRC	Loss (4/87)	−1.6
Leukemia	1.7 × 10−3	AKT3	Gain (3/87)	1.8
		JAK2	Gain (6/87)	1.9
		JAK3	Gain (2/87)	1.5
		MTR	Gain (3/87)	1.9
		PBX1	Gain (3/87)	1.7
		PDE4C	Gain (2/87)	2.7
		PTGS1	Gain (7/87)	2.1
		SRC	Loss (4/87)	−1.6
Cell growth and proliferation (1.5 × 10−5–1.2 × 10−2)
Colony formation	9.2 × 10−4	IFI16	Gain (3/87)	1.7
		UBE2C	Loss (3/87)	−1.9
		CLU	Gain (2/87)	2.0
		ELF3	Gain (3/87)	1.6
		JAK2	Gain (6/87)	1.9
		LAPTM4B	Gain (2/87)	2.0
		PDLIM2	Gain (2/87)	2.8
		SRC	Gain (4/87)	−1.6
Proliferation of cells	1.3 × 10−3	AKT3	Gain (3/87)	1.8
		ANGPT1	Gain (2/87)	2.8
		BLZF1	Gain (3/87)	1.6
		CD244	Gain (3/87)	2.2
		HRPT2	Gain (3/87)	2.2
		CLU	Gain (2/87)	2.0
		CTSB	Gain (2.87)	2.7
		ELF3	Gain (3/87)	1.6
		ENPP2	Gain (2/87)	1.5
		HNF4A	Loss (3/87)	−1.8
		HPGD	Gain (2/87)	2
		IFI16	Gain (3/87)	1.7
		IKBKB	Gain (2/87)	1.6
		ITGA2B	Gain (5/87)	2.3
		JAK2	Gain (6/87)	1.9
		JAK3	Gain (2/87)	1.5
		MAPRE	Loss (2/87)	−1.5
		NFIB	Gain (7/87)	2.3
		PBX1	Gain (3/87)	1.7
		SMARCA2	Gain (6/87)	2
		SOCS3	Gain (6/87)	1.9
		SRC	Loss (4/87)	−1.6
		TGM2	Loss (4/87)	−3.3
		UBE2C	Loss (3/87)	−1.9
Cell death (2.5 × 10−6–1.3 × 10−2)
Apoptosis	3.2 × 10−4	ANGPT1	Gain (2/87)	2.8
		CLU	Gain (2/87)	2.0
		DAPK1	Gain (6/87)	1.8
		DNAJB9	Loss (2/87)	–1.6
		IKBKB	Gain (2/87)	1.6
		PTK2B	Gain (2/87)	1.8
		SOCS3	Gain (5/87)	1.9

Abbreviations: IL, interleukin; JAK, Janus kinase; MPN, myeloproliferative neoplasm; SNP, single-nucleotide polymorphism; SOCS, suppressor of cytokine signaling.

Ingenuity pathway analysis was performed using genes commonly identified by SNP array analysis (n=87) and meta-analysis (n=114) (89 genes).

Aberrations involving chromosome 9 are associated with dysregulated gene expression

Previous studies show that wild-type JAK2 can exert a dominant-negative effect over mutant JAK2V617F, suggesting that there may be a selective advantage for elevating the ratio of JAK2V617F/JAK2WT.[32] We observed that the majority of patients with a high JAK2V617F allele burden (>85% JAK2V617F) exhibit chr9pLOH (86% 12/14), confirming previous studies that JAK2V617F mutant allele status can be increased through LOH because of mitotic recombination[8] (Supplementary Table S1). We also found that the majority of patients with trisomy 9 had an elevated JAK2V617F allele burden (average 67%), consistent with duplication of the mutant allele.

As trisomy 9 and 9pLOH represent the predominant CNAs in MPN, we used gene expression profiling to investigate the differences between patients with normal (n=25) versus abnormal chromosome 9 cytogenetics (n=13). We identified a set of 493 genes, including 210 upregulated and 283 downregulated genes that were differentially expressed between these two groups (⩾1.5-fold, P<0.05) (Supplementary Table S7). Among those genes upregulated in patients with chromosome 9 aberrations, were JAK2, STAT5B, MAPK14 and CD177 (Figure 2).

Figure 2

Chromosome 9 aberrations are associated with elevated JAK2V617F allele burden and gene expression signature. Gene expression differences in patients with normal (n=25) versus abnormal (+9, 9pLOH) chromosome 9 cytogenetics (n=13) (Supplementary Table S7). A subset of genes upregulated (>2.7-fold, P<0.05; 26 genes) or downregulated (>3-fold, P<0.05; 26 genes) in patients with chromosome 9 anomalies compared with normal chromosome 9 cytogenetics are displayed in a heatmap.

PMF and ET patients display higher frequencies of chromosomal aberrations compared with PV

Although PV, ET and PMF patients may all harbor the JAK2V617F allele, their clinical manifestations differ. This may be at least in part explained by the gene dosage hypothesis, which is based on the observation that homozygous JAK2V617F mutant erythroid colonies are found in the majority of patients with PV but rarely in ET patients.[33] This suggests that higher kinase activity (JAK2V617F homozygosity) gives rise to an erythroid/PV phenotype, whereas lower kinase activity (JAK2V617F heterozygosity) leads to an ET phenotype. In support of this model, we observed that JAK2V617F homozygosity was associated with PV (64% 9/14) and PMF (36% 5/14) but not ET patients (Supplementary Table S1).

In addition to the role of JAK2V617F allele dosage, we hypothesized that specific genomic aberrations may have a role in disease phenotype. To test this, we performed unsupervised hierarchical clustering of SNP profiles for MPN patients heterozygous (allele burden range 1–85%) for the JAK2V617F allele (n=63) (Supplementary Figure S3). Although genomic instability profiles of these MPN samples did not segregate according to disease subtype, we found that PMF patients were more likely to have large (>10 Mb) chromosomal aberrations compared with ET or PV patients (P=0.009) (Figure 3). In addition, ET and PMF patients had more moderate sized (1–10 Mb) chromosomal aberrations than PV patients (P=0.01).

Figure 3

PMF and ET patients display higher frequencies of chromosomal aberrations compared with PV. Frequency plots of chromosomal gains (red) or loss (green) in individuals with PV (n=23), PMF (n=15) and ET (n=25). Y axis represents the frequency of copy number change in each patient subgroup. T-test was used to identify significant differences in the frequency of chromosomal aberrations between these three patient subtypes.

NFIB overexpression is associated with chromosome 9 amplification and protects hematopoietic cells from cytokine withdrawal-induced apoptosis

Accumulating evidence suggests that additional genetic lesions cooperate with JAK2V617F to induce MPN. We identified NFIB as a gene upregulated in our MPN meta-analysis, which includes a recently published data set of CD34+ cells of PV patients compared with normal CD34+ specimens.[24] In our SNP array study, we show that NFIB, which lies 9.1MB centromeric to the JAK2 locus, is amplified as a result of trisomy 9 in 6/87 MPN patients and gain of 9p23–13.3 in 1/87 MPN patients. NFIB expression was not significantly changed following JAK2 inhibition in HEL cells and modestly upregulated following JAK2 inhibition in UKE-1 cells, which are both homozygous for the JAK2V617F mutant allele. NFIB expression was not altered by overexpression of either WT or mutant JAK2 in TF-1 cells, which express wild-type JAK2 (Supplementary Figure S4). In contrast, PIM1 was modulated in response to JAK2 inhibition or overexpression in keeping with its positive regulation by JAK/STAT signaling. Taken together, these data suggest that aberrant JAK2 signaling does not account for the upregulation of NFIB in MPN.

To explore the putative oncogenic potential of NFIB, we utilized the erythropoietin-dependent Ba/F3-EPOR murine hematopoietic cell line. Stable overexpression of NFIB in Ba/F3-EPOR cells protected cells apoptosis induced by EPO withdrawal (Figure 4a) Furthermore, NFIB-infected Ba/F3-EPOR cells displayed a modest, but consistent increase in cell number compared with vector (P<0.02, Figure 4b). Taken together, these data suggest that ectopic NFIB expression may cooperate with aberrant JAK2 signaling to dysregulate normal hematopoiesis.

Figure 4

NFIB protects cells from cytokine withdrawal-induced cell death and promotes cell growth (a). Baf3/EPOR cells transduced with control (EF) or EF-NFIB lentivirus for 48 h and GFP+ cells were replated in RPMI with 1 U/ml, 0.1 U/ml and 0.01 U/ml EPO for 48 h. Cell viability was determined by Trypan blue exclusion (b). Cumulative cell counts were performed over 10 days. Data are expressed as the average of three independent experiments±s.e.m. (**P<0.05).

Discussion

The recent use of genome-wide methods to interrogate copy number alterations and gene expression changes in hematologic malignancies has been critical to our understanding of the mechanisms underlying these heterogeneous diseases. In this study, we employed a combination of high-density SNP arrays and gene expression arrays to identify novel oncogenic lesions and dysregulated pathways in a cohort of 87 MPN patients. In addition to identifying copy number variations affecting genes TET2, ASXL1, EZH2 and CBL, which have previously been reported in MPN,[34] we identified new MPN-associated aberrations, which may have a role in disease initiation and/or progression.

The most frequent aberrations observed in our patient cohort involved gains of whole chromosome 9 or 9pLOH, suggesting that genes on this chromosome may provide a selective clonal advantage. Previous studies have shown that wild-type JAK2 is capable of inhibiting STAT5 activation induced by mutant JAK2V617F.[6] This inhibition is presumably due to competition between the two JAK2 forms for cytokine receptor binding because wild-type JAK2 does not affect autophosphorylation of JAK2V617F.[8] Since JAK2 is present on chromosome 9, we hypothesized that gain of chromosome 9 may provide an alternative mechanism for cells to acquire an extra copy of mutant JAK2. In support of this, we observed that a majority of patients trisomy 9 were associated with high JAK2V617F allele burden, suggesting that JAK2 was mutated in this extra chromosome 9. Although the presence of other disease-associated genes on chromosome 9 may also contribute to the pathogenesis of MPN, these data support the hypothesis that acquisition of an additional JAK2V617F mutant allele, either by 9pLOH or trisomy 9, confers a selective, competitive edge.

The finding that JAK2 and STAT5B are significantly upregulated in patients with trisomy 9 and 9pLOH, which are comprised largely of PV patients (57%), reinforces the role of aberrant JAK/STAT signaling in this disease subtype. MAPK14 (also known as p38α), which is also upregulated in patients with chromosome 9 abnormalities, is a kinase that lies downstream of JAK2, and is activated in response to pro-inflammatory stimuli and hematopoietic growth factors. Recent studies reveal a role for MAPK14 in red blood cell production, since mice lacking MAPK14 are anemic due to failed definitive erythropoiesis.[35] Given the association of high-level JAK2V617F expression with a PV-like phenotype in murine models of MPN,[36] it is tempting to speculate that aberrant expression of MAPK14 in patients with two copies of JAK2V617F may activate an erythroid differentiation pathway, contributing to the development of PV.

In addition to gaining extra copies of JAK2V617F, modulation of JAK2 kinase activity may also have a role in disease progression. The suppressor of cytokine signaling (SOCS) proteins are negative regulators of cytokine signaling and are expressed after cytokine stimulation, directly interacting with JAK/STAT pathway members to modulate the responsiveness of cells to further signaling.[37] Although it is postulated that SOCS proteins are tumor suppressors, there is also evidence to suggest that overexpression of SOCS may favor tumor progression. Persistent expression of SOCS3 is observed in AML and coincides with constitutive activation of JAK-STAT pathways.[38] Interestingly, our study identified SOCS3 within a region of chromosomal gain on chromosome 17 in 6/87 patients and SOCS3 was also upregulated ∼2-fold in our meta-analysis. Although SOCS3 is a potent negative regulator of erythropoietin signaling through the EPO receptor and wild-type JAK2, a study by Hookham et al.,[39] revealed that SOCS3 is unable to negatively regulate mutant JAK2V617F and in contrast, enhanced the proliferation of Ba/F3 cells in the presence of EPOR. Furthermore, SOCS3 degradation was inhibited by JAK2V617F and this correlated with SOCS3 tyrosine phosphorylation, a feature also observed in the PBMC of patients homozygous for JAK2V617F. Indeed, we found that all patients harboring a gain of chromosome 17q leading to amplification of SOCS3, had intermediate JAK2V617F allele burden (average 55.4%), supporting a model in which overexpression of SOCS3 may cooperate with JAK2V617F to enhance myeloproliferation.

Although aberrant JAK2 signaling is clearly a recurrent theme in MPN, it is likely that other JAK2V617F-independent events contribute to disease pathogenesis. Recent studies have identified several disease-associated mutations such as TET2 and ASXL1 that may initiate a premalignant clone, since they appear to precede the JAK2V617F mutation.[2, 11] In this study, we show that NFIB, which we previously reported to be overexpressed in CD34+ cells from nine PV patients compared with normal controls, is amplified in six patients via gains of whole chromosome 9 and in one patient by gain of region 9p23–13.3. Previous studies reporting overexpression of NFIB in a small subset of PV patients (two with 9pLOH and two without), hypothesized that this was a consequence of 9pLOH.[40] It was later shown, however, that NFIB was overexpressed in eight MPN patients positive for endogenous erythroid colony formation, however, none of these patients had 9pLOH.[41] Although the ectopic NFIB expression may be a consequence of trisomy 9, it is also possible that NFIB may be upregulated by other factors at the level of gene transcription. Indeed, we observed that only two genes present in the 9pLOH region, JAK2 and C9orf46, which encodes a transmembrane protein with unknown function, were upregulated in these patients compared with patients with normal chromosome 9 cytogenetics. NFIB encodes a CAAT box-binding transcription factor and previous studies showed that overexpression in 32D cells increases resistance to TGF-β1 inhibition.[40] NFIB expression was not positively regulated by JAK2 in cell lines, however, overexpression in Ba/F3-EPOR cells conferred a modest growth advantage, and impaired apoptosis induced by EPO withdrawal. The role of the related NFIA gene in erythropoiesis is well established, and a recent study revealed structural alterations in NFIA in 2% of chronic myeloid diseases, theoretically leading to inactivation of the protein.[42] The functional role of NFIB in normal and malignant hematopoiesis remains unknown. Preliminary experiments suggest that NFIB is expressed at low levels in human CD34+ cells, however, expression is upregulated during normal megakaryocte differentiation (data not shown). Future studies will seek to understand the effect of NFIB overexpression on normal hematopoiesis and in the context of the JAK2V617F overexpression.

In addition to NFIB, overlap of genes identified by our SNP array and MPN meta-analysis identified several other potentially interesting disease candidates. For example, MLLT3 and SMARCA2, both of which are located on chromosome 9 and upregulated in MPN patients, were identified as part of a gene set expressed by ‘normal' quiescent HSCs.[31] MLLT3 is a positive regulator of erythroid and megakaryocyte differentiation during normal hematopoiesis[43] and is also implicated in leukemogenesis, as a translocation fusion partner of MLL in AML.[44] EZH1, which is located in a region of chromosome gain (17q21.1–21.31) in 5/87 patients and is upregulated in the MPN meta-analysis, is an H3K27 methyltransferase that forms part of a non-canonical PRC2 complex and has a role in preserving H3K27me3, preventing derepression of PRC2 targets involved in ESC pluripotency.[45] Recent studies identified loss of function mutations in EZH2 in MDS and MPN, suggesting a tumor-suppressor activity. Although the functional consequences of EZH1 overexpression in MPN remains to be elucidated, these studies highlight an emerging theme in myeloid malignancies involving the dysregulation of genes involved in epigenetic regulation, such as ASXL1, UTX, SETD2, JARID1C and TET2.

Recent studies have demonstrated that JAK2V617F may affect genomic stability by displacing HP1α from chromatin[46] and by deregulating DNA double-strand break repair mechanisms.[47] Plo et al.,[47] showed that overexpression of JAK2V617F and to a lesser extent wild-type JAK2 in Ba/F3-EPOR cells, led to an increase in spontaneous homologous recombination and RAD51 nuclear foci formation. JAK2V617F expression was associated with genetic instability, as assessed by the higher occurrence of centrosome abnormalities and sister chromatid exchange. Furthermore, increased homologous recombination and RAD51 formation was observed in patient CD34+ cells from JAK2V617F-positive PV and PMF patients.[47] In support of these studies, we observed significantly more genomic aberrations in patients with intermediate JAK2V617F allele burden (1–85%) compared with patients with wild-type JAK2 (<1% JAK2V617F) (data not shown). Unexpectedly, however, patients with high JAK2V617F allele burden (>85%) displayed fewer aberrations. The ability of mutant JAK2 to induce hyper-recombination is thought to facilitate mitotic recombination that may lead to LOH and account for disease heterogeneity.[47] Whether genomic instability induced by JAK2V617F facilitates leukemic transformation, however, remains contentious. For example, in a study of 17 patients with JAK2V617F-positive MPN that later transformed to AML, 9 patients went on to develop JAK2V617F-negative AML.[10] Furthermore, a recent study by Thoenssin et al.,[48] found no significant association between JAK2V617F mutation status and time to leukemic transformation in patients with MPN-blast phase, suggesting the existence of a pre-JAK2V617F clone.

Although unsupervised clustering of PV, ET and PMF patients did not reveal distinct aberration patterns associated with these MPN subtypes, we noted that PMF patients displayed more chromosomal aberrations that PV or ET patients. These observations are consistent with other studies showing that patients with PMF exhibit increased copy number changes and abnormal cytogenetics compared with PV or ET.[17] The consequences of this increased genomic instability remains unclear, however, a recent study investigating the chromosomal abnormalities of patients with leukemic transformation of MPN, revealed that time from diagnosis to leukemic transformation was shorter for PMF than ET or PV patients.[48] Further analysis on a larger cohort may help us determine the aberrations associated with a specific MPN subtype.

In conclusion, we have used a combination of high-resolution SNP arrays and gene expression arrays to identify novel potential disease candidates that may have a role in MPN. Future studies will seek to address the specific contribution of aberrated genes and signaling pathways identified in this study to the development of MPN.