Single-Cell Analyses Reveal Megakaryocyte-Biased Hematopoiesis in Myelofibrosis and Identify Mutant Clone-Specific Targets

1Advances in single-cell technologies have recently provided new insights into the cellular and molecular diversity and pathological mechanisms underlying many diseases, including cancers, pre-malignant and non-malignant conditions. Parallel interrogation of mutation status and the transcriptome at a single-cell level provides an unprecedented opportunity to identify cancer-cell-specific targets. Single-cell resolution also uniquely enables the identification of rare cell types and analysis of combinatorial patterns of gene expression, both of which are necessary to reconstruct differentiation trajectories and to accurately define cellular heterogeneity between populations, such as normal and malignant tissues, as well as to identify the mediators of interactions between different cell types. For example, pathological fibrosis underlies many prevalent diseases, including cancer, where fibrosis is well recognized to be important for disease progression and metastasis. It is broadly proposed that pro-fibrotic mediators secreted by cancer cells and infiltrating immune cells activate non-malignant stromal cells, such as myofibroblasts, to deposit collagen fibrosis. However, an understanding of the specific cellular populations that mediate fibrosis in a given disease model, their molecular features, and the cellular pathways through which they are generated is necessary for these cells to be therapeutically targeted.

2Myelofibrosis is the most severe of the “myeloproliferative neoplasms” (MPNs), a group of heterogeneous disorders that result from somatic mutations in hematopoietic stem and progenitor cells (HSPCs) affecting Janus Kinase (JAK) signaling. The most common driver mutation is JAK2V617F, occurring in ∼60% of myelofibrosis patients, with mutations affecting calreticulin (mutCALR) found in the majority of other patients. Myelofibrosis can occur as a primary disorder (PMF) or develop secondary to the other MPNs polycythemia vera (post-polycythemia vera myelofibrosis [PPV-MF]) or essential thrombocythemia (post-essential thrombocythaemia myelofibrosis [PET-MF]). Myelofibrosis is characterized by a progressive bone marrow fibrosis that destroys the hematopoietic microenvironment, resulting in the cardinal disease features of cytopenias, mobilization of HSPCs to peripheral blood, extramedullary hematopoiesis, and a high propensity for leukemia. Survival is typically 5–10 years from diagnosis and is not substantially improved by currently available drug therapies. Megakaryocytes, the platelet-producing cells in the bone marrow, are dramatically increased in number in myelofibrosis and are one of the key cellular drivers of the destructive bone marrow remodelling by releasing excess pro-fibrotic cytokines and growth factors. In normal hematopoiesis, megakaryocyte progenitors (MkPs) have a low proliferation rate, typically undergoing less than 8 cell divisions before mitotic arrest and the onset of polyploidization, and megakaryoyctes are relatively rare cells in healthy bone marrow.

3The cellular and molecular pathways that give rise to the dramatically increased megakaryocyte numbers and megakaryocyte dysfunction leading to tissue fibrosis are unclear. In traditional models of normal hematopoiesis, megakaryocytes are said to arise from a bipotent progenitor shared with the erythroid (red cell) lineage—the megakaryocyte-erythroid progenitor (MEP). Recent advances in single-cell technologies, including single-cell transplantation and lineage tracing studies of unperturbed hematopoiesis, have revealed that hematopoiesis occurs over a continuum rather than by distinct, oligopotent intermediate steps and also that a proportion of hematopoietic stem cells (HSCs), at least in the murine system, are megakaryocyte-biased but retain the capacity for multilineage reconstitution. Lineage-committed megakaryocytes arising directly from HSCs, sometimes without cell division, have also been reported.

4Targeting megakaryocytes in myelofibrosis has been shown to ameliorate the disease in mouse models and early-phase human studies, but technical challenges have precluded the extensive study of megakaryopoiesis in myelofibrosis patients. These challenges include the rarity of megakaryocytes in healthy bone marrow, gaps in our knowledge of the cellular pathways of megakaryopoiesis, and their extreme cell size and fragility. In addition, the severe fibrosis typically prevents bone marrow aspiration (“dry tap” aspirate). However, bone marrow HSPCs are mobilized to the peripheral blood in myelofibrosis. In this study, we used this phenomenon to capture peripheral blood HSPCs and perform the first in-depth single-cell analysis of abnormal megakaryocyte differentiation and function in patients with myelofibrosis, suggesting key cellular and molecular targets. Using multiparameter immunophenotyping, functional studies, high-throughput single-cell RNA sequencing (scRNA-seq), targeted single-cell mutational analysis with simultaneous scRNA-seq (TARGET-seq), and single-cell proteomics, we identify potential targets for the inhibition of pathological megakaryocyte differentiation and megakaryocyte-induced fibrosis and validate G6B as a cell surface marker that may enable specific ablation of myelofibrosis cells using immunotherapy. This study illustrates the power of single-cell “multi-omics” in the characterization of cellular heterogeneity in cancers associated with aberrant fibrosis, including the identification of potential therapeutic pathways and cancer-cell-specific targets.

5Multi-parameter flow cytometric analysis of the CD34+ lineage (lin)− HSPC compartment in peripheral blood samples from healthy mobilized apheresis donors and patients with myelofibrosis (Table S1) was performed to compare frequencies of the classically defined HSPC subsets (Figure 1A). This demonstrated reduced lymphoid-primed multipotent progenitors (LMPPs) and increased multi-potent progenitors (MPPs; Figure 1A) in myelofibrosis patients. The cell-surface antigen CD41 has previously been reported to identify cells primed for megakaryocyte differentiation. A 5-fold increase in the percentage of CD41+ cells was detected within both CD38-negative, early stem/progenitor (HSC and MPP) and CD38-positive, down-stream progenitor (MEP and common myeloid progenitor [CMP]) cell fractions (Figures 1A and 1B), suggesting a bias toward megakaryocyte differentiation originating during the earliest phases of HSC lineage commitment. Morphological analysis of CD38−CD41+ and CD38+CD41+ cells from the CD34+lin−CD45RA− compartment showed undifferentiated blast cell morphology and not mature megakaryocytes (Figure S1A).

6The CD41+ fraction of human CD38-positive CD34+ lin− CD45RA− HSPCs contains megakaryocyte-biased progenitors with significant erythroid differentiation potential as well as unipotent MkP. However, the phenotype of CD41+ cells within the CD38-negative HSC/MPP compartment has not previously been defined. We, therefore, sought to determine whether the CD41+ HSCs and MPP cells isolated from healthy donors retained a capacity for multi-lineage differentiation or were lineage-committed MkP. CD34+ Lin– CD38– CD45RA– CD90+ CD41– (CD41– HSC), CD34+ Lin– CD38– CD45RA– CD90– CD41– (CD41–MPP), and CD34+Lin– CD38– CD45RA–CD41+ (CD41+HSC/MPP) cells were isolated by fluorescence-activated cell sorting (FACS) for liquid culture differentiation assays. When stimulated with thrombopoietic cytokines, CD41+ HSC/MPP cells showed accelerated megakaryocyte differentiation with a substantially higher proportion of cells expressing the mature megakaryocyte surface antigen CD42, a large cell size, and proplatelet extensions at early time points as compared to CD41– HSCs and MPPs (Figure 1C). In parallel megakaryocyte, erythroid, and myeloid differentiation assays, CD41+ HSC/MPP showed a similar potential for CD11b/CD14+ myeloid differentiation and a superior potential for CD71+/glycophorin A erythroid differentiation than CD41– fractions (Figure 1D).

7In comparison to those from healthy donors, CD41– HSC/MPP cells from myelofibrosis patients showed a megakaryocyte versus erythroid differentiation bias (Figure 1E), in keeping with the clinical phenotype of myelofibrosis patients in which excessive megakaryocyte numbers occur in parallel with anemia. In single-cell clonogenic assays supportive of myeloid and erythroid (but not megakaryocytic) colony formation (methocult), CD41+ and CD41– fractions of HSCs and MPPs gave rise to expected colony frequencies with myelofibrosis CD41- HSC/MPP showing a bias toward myeloid versus erythroid colonies (Figure S1B). Together, these results support that in myelofibrosis, HSPCs are biased toward megakaryocyte-lineage differentiation from the earliest stem cell compartment, even before expression of canonical megakaryocytic markers.

8To identify the cellular and molecular basis for megakaryocyte-biased hematopoiesis in myelofibrosis without bias from pre-selected cell surface antigens, high-throughput scRNA-seq was performed on 135,929 individual CD34+ lin– HSPCs from patients with JAK2V617F+ or mutCALR+ myelofibrosis (93,157 cells, n = 15) according to World Health Organization (WHO) criteria and age-matched healthy donors (42,772 cells, n = 6) by using the 10x Genomics Chromium platform (Table S2). Filtering, quality control and doublet exclusion was performed (Table S3). Healthy donor control and myelofibrosis cells were aggregated and individual donor effect was regressed out, following which no batch effect remained (Figure S1C). A contaminating population of plasmacytic dendritic cells was identified and removed from all down-stream analysis (Figure S1D). Following these steps, 120,196 cells (82,255 myelofibrosis and 37,941 control cells) were used for down-stream analyses (Table S3).

9Dimensionality reduction and unsupervised clustering were performed using a uniform manifold approximation and projection (UMAP) method combined with the Louvain community-detection clustering method to enable identification of distinct cell populations while preserving inter-cluster relationships (Figure 2A). A total of 8 clusters were identified and manually annotated by correlation of differentially expressed genes for each cluster with reference marker genes for each lineage (Figure 2A; Figure S1E; Table S4). “Lineage signature” gene sets were established by an analysis of published datasets to identify genes selectively expressed in erythroid, myeloid, lymphoid, and megakaryocyte lineage progenitors and uncommitted HSPCs and were superimposed on the UMAPs (Figure 2B; Figures S1F and S2; Table S5).

10Most HSPCs expressing megakaryocyte signature genes were derived from myelofibrosis patients with very few from healthy donors (Figures 2B–2D; Figure S1F). The average proportion of HSPCs with a megakaryocyte gene signature was ∼5% for individual myelofibrosis patients, 11-fold higher than for healthy donors (4.98% [range, 0.52%–19.81%] for myelofibrosis versus 0.44% [range, 0%–1.51%] for controls, p < 0.05). Down-sampling the myelofibrosis cells to create a dataset with equal numbers of control and myelofibrosis HSPCs confirmed that over 93% of HSPCs in the MkP cluster originated from myelofibrosis donors, whereas less than 10% lymphoid progenitor cells were from myelofibrosis patients. The fractions of myeloid and erythroid progenitor cells were not significantly different (Figures 2D and S1F), supporting a strong megakaryocyte bias and reduction in lymphoid differentiation in myelofibrosis.

11To study differentiation trajectories, cells were ordered in gene expression space by using force-directed graphs (FDG), and the lineage signature gene scores were superimposed on the graphs (Figures 3A–3D). The HSPC signature highlighted cells at the origin of the trajectory, and erythroid, megakaryocyte, myeloid, and lymphoid trajectories formed distinct paths, with MPP cells primarily located in intermediate positions in the trajectories (gray cells, Figures 3A–3D, left plots; Figure S3). Expression of megakaryocyte signature genes (purple) was observed along a prominent distinct trajectory originating from the apex of the HSPC cluster in the FDG trajectory in myelofibrosis patients (Figures 3B and 3D). In contrast, in healthy donors, very few MkP cells were observed in a distinct differentiation trajectory (Figures 3C and 3D). Notably, most of these cells derived from two of the six healthy donors (donors ID06 and ID09; Figure S3C). Together with functional data (Figure 1), these data support a model in which a direct route for MkP production from HSPCs is aberrantly expanded in JAK2V617F and mutCALR-driven myelofibrosis. Our observations were consistent across all clinical and molecular patient subgroups (Figures S3A–S3C), and trajectory analyses using diffusion maps created with Scanpy, an alternative toolkit, confirmed findings with our in-house analysis pipeline (Figure S3D).

12We next sought to identify molecular regulators that might specifically drive aberrant megakaryocyte differentiation in myelofibrosis and may potentially be targeted without major toxicity to the other blood cell lineages. Trajectory analysis of all myelofibrosis HSPCs performed using the Scanpy toolkit demonstrated a distinct trajectory for megakaryocyte differentiation through “pseudotime” from an HSC origin (HSC → HSPC2 → Mega; Figure 4A). Expression patterns of 1,639 human transcription factors were analyzed, and transcription factors showing progressive changes over the megakaryocyte and erythroid differentiation trajectories, either increased or decreased expression, were identified. Expected differential expression patterns of transcription factors known to be involved in megakaryocyte versus erythroid lineage specification were observed, e.g., progressive increase in GATA1 and GATA2 and antagonistic expression of two key regulators of megakaryocyte-erythroid cell fate decision, namely FLI1 and KLF1 (Figures 4B and 4C). Additional genes not previously implicated as regulators of megakaryocyte versus erythroid differentiation showed striking differential expression between the erythroid and megakaryocyte trajectories, including YBX1, HMGA1, PLEK, SOX4, and MYC (Figures 4B and 4C), suggesting additional targets for strategies to specifically inhibit pathological megakaryopoiesis while preserving erythropoiesis in myelofibrosis patients.

13To evaluate the pathological role of the expanded MkPs in driving bone marrow fibrosis, we next examined potential mediators of fibrosis among HSPCs. Fibrosis regulators were identified from previously published datasets studying lung and liver fibrosis as well as bone marrow fibrosis. Genes detected at expression levels over 1 (using log-transformed unique molecular identifier [UMI]) in our HSPC dataset were selected for a “fibrosis signature” gene score (Table S5). Superimposition of this score on the UMAP for all healthy donor and myelofibrosis HSPCs clearly highlighted the MkP cluster cells as being the key regulators of fibrosis among all HSPCs (Figure 5A).

14Healthy donor and myelofibrosis MkPs were then extracted for further analyses (Figure S4A). Gene set enrichment analysis (GSEA) of all expressed genes pre-ranked according to their differential expression between myelofibrosis and control MkP showed significant enrichment of metabolic (e.g., oxidative phosphorylation and fatty acid metabolism) and inflammatory pathways (e.g., transforming growth factor α [TNF-α] signaling and interferon gamma response) in myelofibrosis MkPs (Figure 5B). A subset of myelofibrosis MkPs (cluster cluster 6) showed a similar gene expression profile to control MkPs (Figures 5C and 5D). These cells had a high expression of key inflammatory mediators previously implicated in myelofibrosis, including PDGFA, CCL5, and CXCL5 (Figures 5C and 5D). However, the majority of myelofibrosis MkPs had a unique gene expression profile with overall upregulation of genes normally expressed at low levels in healthy donor MkPs (Figures 5C and S4B). This finding suggests that megakaryocyte-induced fibrosis in myelofibrosis is due to both expansion of a population of megakaryocytes analogous to those in normal bone marrow as well as the generation of an aberrant population, an observation that would not have been possible without analysis at the single-cell resolution.

15To further dissect cellular and molecular heterogeneity among myelofibrosis MkP, unsupervised clustering using Louvain community detection was performed on myelofibrosis MkP. Nine sub-clusters were identified with distinct expression of fibrosis and proliferation-associated genes (Figures 5D, 5E, and S5A). Genes encoding key mediators of fibrosis (TGFB1, TNF, and LTBP1, which encodes a protein that targets the latent form of transforming growth factor beta [TGF-β] and contributes to its activation; ), were most highly expressed in MkP clusters 6–8, whereas MkP clusters 2 and 8 showed the highest expression of the proliferation marker MKI67 and a G2M gene signature (Figures 5E and S5B). AURKA was selectively expressed in two clusters, with particularly high expression in the minor cluster 8 (Figure 5E). This is of interest as AURKA is the target for alisertib (MLN8237), which was recently demonstrated to promote megakaryocyte polyploidization and ameliorate the myelofibrosis phenotype in mouse models, with some efficacy also in patients with myelofibrosis.

16Normal megakaryocytes have a low proliferation index, and healthy donor MkPs showed low expression of the proliferation marker MKI67. By contrast, MKI67 was strongly expressed in the majority of myelofibrosis MkPs and the MkP cluster showed the highest expression of MKI67 among all myelofibrosis HSPC clusters (Figure S5C) as well as an enrichment of a G2M checkpoint gene signature (Figure S5B; Table S5), suggesting that increased proliferation of MkPs may contribute to the pathological accumulation of megakaryocytes in myelofibrosis, in addition to megakaryocyte-biased hematopoiesis.

17Increased expression of megakaryocyte genes in the myelofibrosis aggregate was noted not just within the MkP cluster but also within clusters of uncommitted HSPCs and other lineage-affiliated clusters (Figure 6A). This included intracellular proteins (PF4 and VWF) and also cell surface antigens (ITGA2B [CD41] and MPIG6B [G6B]). Increased expression of MPIG6B, encoding the G6B protein, was particularly striking (Figure 6A). G6B is an immunoreceptor tyrosine-based inhibition motif (ITIM)-containing inhibitory receptor, considered to be exclusively expressed on mature megakaryocytes in normal hematopoiesis. As the majority of healthy donor CD34+ lin– HSPCs did not express megakaryocyte genes and because mature megakaryocytes normally lose expression of CD34 during differentiation, we hypothesized that aberrant co-expression of stem and progenitor and megakaryocyte surface antigens may enable the selective identification of myelofibrosis clone-derived HSPCs.

18Patients with myelofibrosis have distinct genetic subclones of HSPCs, including residual wild-type clones (non-mutated) as well as clones with co-mutations in addition to driver mutations (JAK2V617F or mutCALR). To determine whether the increase in the expression of megakaryocyte-associated genes was specific to mutant clone-derived HSPCs or due to cell-extrinsic signals affecting both mutated and un-mutated HSPCs, data generated by combined high-sensitivity mutational analysis and parallel transcriptome profiling (TARGET-seq; ) were analyzed. A total of 2,734 cells were examined—678 healthy donor cells plus 2,056 myelofibrosis cells (388 JAK2 wild type and 1,668 JAK2V617F mutated). Expression of megakaryocyte genes, in particular G6B, was significantly higher in JAK2V617F-mutated HSPCs than in either wild-type cells from the same patients or in wild-type cells from healthy donors (32% versus 22.9% versus 14.7%, respectively; p < 0.001; Figure 6B). Wild-type cells from myelofibrosis patients also showed increased frequency of G6B expression, albeit to a lower degree than JAK2V617F-positive cells, in keeping with cell-extrinsic signals also contributing to this aberrant megakaryocyte differentiation (Figure 6B).

19]The high-throughput TARGET-seq and 10x Chromium datasets included all CD34+ lin– cells. Expression levels of G6B were also increased specifically in individual JAK2 mutant CD38– early stem and progenitor cells (HSC/MPP; Figure 6C) identified from the index sorting data of the TARGET-seq cells, and in 100-cell “mini-bulk” preparations of FACS-isolated immunophenotypic CD34+ lin– CD38–CD45RA– CD90+ HSCs, CD34+ lin–CD38– CD45RA–CD90– MPPs, and CD41+ HSC/MPPs (Figure 6D).

20To examine where cells with specific genotypes fell on the HSPC trajectory, the datasets of myelofibrosis HSPCs analyzed by high-throughput 10x Genomics and TARGET-seq were integrated using Harmony (Figure S6A). FDG trajectory analysis showed that both wild-type and mutant progenitors fell in all 3 of the lineage progenitor trajectories (myeloid, erythroid, and megakaryocyte; Figure S6A). However, a higher proportion of the cells in the megakaryocyte and myeloid trajectories were mutant versus wild type than the erythroid trajectory (Figure S6A, right plot). In two patients with 3+ co-mutations in addition to the driver JAK2V617F mutation, the increase in G6B was observed in all genetic sub-clones detected (Figure S6B).

21High-throughput, single-cell proteomics by mass cytometry time of flight (CyTOF) was performed to simultaneously measure 20 surface proteins in multiple samples in parallel by using barcode multiplexing (Figure 7A; Table S7). G6B was consistently detected at substantially higher levels in patients with primary and secondary myelofibrosis and with JAK2V617F and mutCALR driver mutations than in healthy donors (Figures 7A and 7B). In addition, high cell surface G6B expression was also detected exclusively on JAK2V617F-mutated MPN cell lines (HEL and SET2) and not on the other leukemia cell lines, namely K562, HL60, JURKAT, and MARIMO, and HEK (human embryonic kidney) cells (Figure S7A). G6B expression was noted in both the CD41-positive and -negative cell fractions in myelofibrosis by FACS (Figure 7B, right plots).

22To examine G6B expression in bone marrow megakaryocytes in situ, immunohistochemical staining was performed on trephine biopsy sections from healthy donors and patients with mutCALR and JAK2V617F+ myelofibrosis, confirming the expected expression on control megakaryocytes but with a dramatic increase in G6B+ cells in myelofibrosis (Figure 7C).

23To confirm G6B as a potential marker of HSPCs derived from the malignant clone, G6B-positive and -negative cells were FACS-isolated from healthy donor and myelofibrosis patient MNCs and the expression of mutant versus wild-type JAK2 was determined by quantitative real-time PCR. Expression of mutant JAK2V617F relative to wild-type JAK2 was dramatically higher in G6B-positive cells (Figure 7D).

24To validate that aberrant co-expression of G6B by mutant clone-derived HSPCs may be a potential strategy for immunotherapy, a tool bi-specific antibody was generated to target a stem and progenitor surface antigen (CD34) and G6B and labeled with a pH-sensitive cyanine CypHer5E dye that maximally fluoresces at an acidic pH following cellular internalization. Rapid internalization of the bi-specific by SET-2 cells (a JAK2V617F+ MPN cell line expressing both CD34 and G6B), as indicated by clear intracellular fluorescence, was observed (Figure 7E). Internalization was significantly reduced by two inhibitors of endocytosis, Dynasore (a GTPase inhibitor of dynamin) and Pitstop 2 (an inhibitor of the clathrin terminal domain). Internalization was also observed in CD34-negative, G6B+ HEL cells (Figure 7B) and not in G6B-negative cells nor with an isotype-control antibody, confirming that internalization was mediated by the G6B receptor. This strategy may have the potential to deliver therapeutic antibody conjugates directly to G6B+ mutant clone-derived cells. Together, these data identify G6B as a promising cell surface antigen worthy of further validation as enabling selective targeting of aberrant megakaryocytic differentiation and mutant clone-derived HSPCs in myelofibrosis.

25Bone marrow transplant is currently the only potentially curative treatment for myelofibrosis but is associated with significant risk and most patients are ineligible due to age and comorbidities. The introduction of JAK inhibitors has led to significant improvement in symptomatic management, but most patients continue to experience substantial morbidity and a significant reduction in life expectancy. New approaches to treatment are urgently required. Megakaryocytes are well recognized as the key cellular drivers of disease pathogenesis; however, only one megakaryocyte-targeting therapy—alisertib, a specific inhibitor of aurora kinase—has been developed to date. A major obstacle to identification of novel targets has been the inability to isolate megakaryocytes from patients for detailed study. In the present study, we reasoned that aberrant megakaryopoiesis in myelofibrosis is very likely to be caused by aberrant differentiation of HSPCs, rather than proliferation of mature megakaryocytes alone, and that this process might be amenable to therapeutic targeting to “turn off the supply.” We, therefore, set out to characterize the distinct cellular and molecular features of megakaryocyte differentiation pathways in myelofibrosis by using a combination of single-cell approaches. We demonstrate a dramatic expansion of megakaryocyte differentiation from uncommitted stem and progenitor cells in JAK2V617F-driven hematopoiesis. Furthermore, we identify a number of molecular targets that may inhibit the abnormal megakaryocyte differentiation and potentially ablate mutant clone-derived HSPCs and MkPs.

26Importantly, several key observations were only possible due to the single-cell-level resolution of study, highlighting the power of single-cell technologies in understanding disease pathology and in novel therapeutic target discovery. First, our data indicate that aberrant megakaryopoiesis in myelofibrosis is due to both a dramatic expansion of MkP with a similar transcriptional signature to healthy donor MkPs, as well as generation of a unique, aberrant MkP population. Second, by simultaneously interrogating the mutational status and the transcriptome of individual cells, we demonstrated that certain megakaryocyte surface antigens, in particular G6B, are markedly over-expressed in mutant clone-derived HSPCs compared with wild-type HSPCs from myelofibrosis patients or healthy donor HSPCs. This validates combinatorial targeting of stem cell (e.g., CD34) and megakaryocyte (e.g., G6B) surface antigens, e.g., with bi-specific antibody therapies as a potential strategy worthy of further investigation for selective ablation of the myelofibrosis clone. As none of the currently available treatments for MPNs reliably induce clonal remissions or substantially reduce fibrosis, this work sets the stage for immunotherapeutic targeting of aberrant hematopoiesis in myelofibrosis. Furthermore, the approach we have adopted and the resulting insights are highly relevant to other studies seeking to identify cancer-cell-specific drug targets and cancer-associated fibrosis in other malignancies, as well as non-malignant disorders of tissue fibrosis.

27Further information and requests for resources or materials will be fulfilled by bethan.psaila@ndcls.ox.ac.uk or adam.mead@imm.ox.ac.uk

28HEL, JURKAT, K562, HEK, HL60 and MARIMO cells were obtained from the American Type Culture Collection (ATCC). SET-2 cells were kindly provided by Dr. Jacqueline Boultwood and Dr. Andrea Pellagatti (Radcliffe Department of Medicine, University of Oxford). All cells were maintained in culture in RPMI-1630 supplemented with 10% fetal calf serum (FCS) and 1% penicillin-streptomycin. SET-2 cells were supplemented with 20% FCS.

29A summary of demographic and clinical details of myelofibrosis patients and normal donors used for analysis can be found in Tables S1 and S2.

30Patients and normal donors provided written informed consent in accordance with the Declaration of Helsinki for sample collection, tissue banking and use in research under either the INForMed Study, University of Oxford (IRAS: 199833; REC 16/LO/1376) or Imperial College London (approval reference: R13077; HTA license 12275; REC 12/WA/0196). Cryopreserved peripheral blood mononuclear cells stored in FCS with 10% DMSO were thawed and processed by warming briefly at 37°C, gradual dilution into RPMI-1630 supplemented with 10% FCS and 0.1mg/mL DNase I, centrifuged at 500G for 5 minutes and washed in FACS buffer (PBS + 2mM EDTA + 10% FCS).

31FACS-sorting was performed using a Becton Dickinson Aria III or Fusion 2 and cells isolated into 1.5ml Eppendorf tubes or 96-well plates depending on the experiment. Single color stained controls and fluorescence minus one (FMO) controls were used for all experiments. HSPCs were stained with the following antibody cocktail (see Key Resources Table) for 20 minutes at 4°C and passed through a 70 μm mesh cell strainer if necessary prior to sorting: CD34-APC-efluor780; Linege-BV510; CD38-PE-TxRed; CD123-PeCy7; CD45RA-PE; CD71-AF700; CD41-APC; CD90-BV421. The following antibody cocktail was used to analyze cell differentiation: CD34-APC-efluor780, CD71-AF700, CD36-FITC, CD41 PeCy7, CD42 PE, CD11b-APC, CD14-APC. 7AAD was used for live/dead cell exclusion. For G6B immunostaining, cells were stained with anti-human G6B (17-4) for 30 minutes at 4°C (1:100), washed and stained with goat anti-mouse Alexa Fluor 488 secondary IgG antibody (2:200 ThermoFisher Cat#A10680) for 20 minutes in the fridge and washed prior to staining with fluorescence-conjugated commercial antibodies.

32Cells were isolated by FACS into 1.5 μL eppendorfs, centrifuged at 500G for 5 minutes, resuspended in 100ul culture medium and plated in flat-bottom 96-well plates (Corning). Media used was Stemspan SFEM (StemCell Technologies #09650) + 1% Pen/Strep supplemented with the following cytokines (see also Key Resources)

33Cells were analyzed by FACS on days 6 and 14 (50 μL removed from wells for analysis and replaced with fresh media).

34Cells were FACS-isolated into 1.5ml Eppendorf tubes, centrifuged and resuspended into 200 μL PBS and cytospun at 500RPM for 5 minutes onto Superfrost glass slides. May Grunewald Giemsa stain was prepared as per manufacturers protocol, filtered and slides stained in May-Grunewald for 7 minutes followed by 20 minutes in Giemsa then washed in distilled water, air-dried and coverslip applied.

35Single cells were FACS-isolated into flat bottomed 96-well plates containing 100 μL of MethoCult™ H4435 Enriched (StemCell Technologies Cat#04435). Colonies were visually inspected and classified 11-14 days after plating. Lineage assignment was made by morphological assessment with verification of ambiguous colonies by plucking and FACS analysis.

36Cells were thawed, stained with FACS antibodies and sorted on a BD Aria III or Fusion 2 as described above and as per recommendations in the 10x Genomics Single Cell Protocols – Cell Preparation Guide. 15,000 CD34+ lineage negative cells were sorted into 2 μL PBS/0.05% BSA (non-acetylated) and then the cell number/volume adjusted to the target for loading onto the 10x Chromium Controller. Samples were processed according to the 10x protocol using the Chromium Single Cell 3′ library and Gel Bead Kits v2 (batch 1) or v3 (batch 2) (10x Genomics). Cells and reagents were prepared and loaded onto the chip and into the Chromium Controller for droplet generation. RT was conducted in the droplets and cDNA recovered through demulsification and bead purification. Pre-amplified cDNA was used for library preparation, multiplexed and sequenced on a HiSeq 2500 (batch 1) or a Novaseq S4 (batch 2) aiming to obtain > 50,000 reads per cell. For some samples, a preliminary, low-depth run was done on a MiSeq using MiSeq Nano Reagent Kit V2 (Illumina Cat#102-2001) to estimate the number of cells and total sequencing required.

37The count matrix for 8 myelofibrosis patients and two healthy donors profiled using 3′TARGET-seq were downloaded from GSE122198, normalized by library size and log2-transformed as previously described. Cells were classified into WT-normal (cells from normal donors), WT-patient (non mutant cells present in patient samples) and mutant (cells from patient samples carrying mutations in the genes targeted).

38100 cells from each population were isolated by FACS into 4 μL of lysis buffer containing oligo-dT primer and dNTP mix in 0.2 mL PCR tubes. Cell lysis, RT and PCR preamplification and purification was performed using the Smart-Seq 2 protocol as previously published. Libraries were pooled and tagmentation performed using the Illumina Nextera XT DNA sample preparation kit (Illumina Cat #FC-131-1024), libraries pooled and sequenced on a HiSeq 2000.

39Antibodies were purchased pre-conjugated when commercially available. Non-available antibodies were conjugated to lanthanide metals using Maxpar X8 antibody labeling kit according to the manufacturer protocol (version 10). The antibody cocktail used is listed in Table S7. For barcoding and staining, cells were washed with Maxpar PBS buffer (Fluidigm #201058) and stained with 0.5 μM Cell-ID Cisplatin Viability Stain (Fluidigm #201064) in 200 μL Maxpar PBS for 5 mintutes at room temperature for dead cell exclusion. The reaction was quenched with Maxpar Cell Staining Buffer (CSB, Fluidigm #201063) and cells fixed, permeabilized and barcoded using the Cell-ID 20-Plex Pd Barcoding Kit (Fluidigm #201060) as per the manufacturers user guide. Barcoded cells were washed, combined and stained with the antibody cocktail as per Table S7 for 30 minutes at room temperature. Cells were washed with Maxpar Cell Staining Buffer (Fluidigm #201068), fixed in 1.6% formaldehyde, washed and resuspended in Fix&Perm Buffer (Fluidigm Cat#201067) with Cell-ID intercalator-Ir (Fluidigm #201103B) and incubated overnight at 4°C. The following day, cells were washed and analyzed on a Helios (Fluidigm). The mass cytometer was tuned and QC was run prior to acquiring samples according to the manufacturers’ recommendations.

40Sections of formalin fixed and paraffin embedded (FFPE) bone marrow trephine biopsies were processed as follows: paraffin was removed, then antigen retrieval was performed using citrate (Roche Cell Conditioning 2 Cat#950-123) pre-treatment for 30 minutes. Sections were washed and incubated with G6B antibody diluted 1:100 in Ventana’s DISCOVERY antibody diluent (Roche Cat#760-108) for 60 minutes at room temperature. Secondary detection was performed using UltraMap DAB anti-Ms HRP detection kit (Roche #760-152) for 16 minutes and slides counterstained with hematoxylin (Roche #760-2021) for 4 minutes and Bluing reagent (Roche #760-2037) for 4 minutes.

41For each experiment, MNCs from myelofibrosis patients and healthy donor controls were thawed and combined 1:1 in FACS buffer prior to antibody staining as described above. 50 G6B+ and G6B- cells were sorted into each well of a 96-well PCR plate (10 replicates per population for each experiment), containing CellsDirect One-Step qRT-PCR kit 2X Reaction Buffer and SuperScript III RT/Platinum Taq Mix (Thermo Fisher Cat#11753100), Ambion SUPERase-In RNase inhibitor (Thermo Fisher Cat#AM2694), TE buffer, JAK2 forward and reverse primers and wild-type and JAK2V617F-specific probe mix (see Key Resources Table). RT and PCR were performed as per manufacturer’s recommendations with 18 pre-amplification cycles and then diluted 5x in TE buffer. Taqman RT-PCR was performed in a 20 μL reaction volume using 4 μL of the diluted cDNA, Taqman Fast Advance Mastermix (Thermo Fisher Cat#4444556) and the primers/probes as detailed in the Key Resources Table. Custom Taqman assays were designed (see Key Resources Table) as previously described using RT-PCR primers flanking the mutant region plus two Taqman PCR probes specific for the normal or mutant sequence. An Applied Biosystems 7500 Fast Real-Type PCR system was used with the default PCR conditions, with each replicate run in duplicate. Intra-assay replicates varying more than 5% were excluded.

42The CD34 x G6B bispecific antibody contains a human IgG1 Fc and was produced using ‘knobs-into-holes’ technology, which involves generating a single amino acid substitution in opposite CH3 domains. Each sequence was inserted into the pcDNA 3.4 expression vector (Thermo Fisher). Following preparation of plasmid DNA, each chain was co-transfected into Chinese Hamster Ovarian cells at a 200 mL scale using the ExpiCHO expression system (Thermo Fisher). The Max Titer protocol was followed. Cells were incubated in a 37ᵒC incubator with a humidified atmosphere of 8% CO2 in air on an orbital shaker. On the day after transfection, ExpiFectamine CHO Enhancer and ExpiCHO Feed was added to the flask at the appropriate volume and the flask was transferred to a 32ᵒC incubator with a humidified atmosphere of 5% CO2 in air on an orbital shaker. On Day 5 post-transfection, the second volume of ExpiCHO Feed was added to the flask. On Day 14, the cells were harvested by centrifuging at 18,000 x g for 30 min. The protein was purified from the supernatant using protein A affinity resin (GE HealthCare Life Sciences). Bound protein was eluted with 20 mM Citrate at pH 2.9 and then immediately neutralized with 10% 1 M Tris. The protein was then further purified using size exclusion chromatography (SEC), and characterized by SDS-PAGE gel and analytical scale SEC.

43The CD34 x G6B bispecific and a non-targeting isotype control (DNP-MB273, Acrobiosystems) were conjugated to CypHer5E (GE Healthcare Life Sciences), a red-excitable, pH-sensitive cyanine dye detected in the APC or Ax647 channel that maximally fluoresces at an acidic pH (i.e., after movement from a receptor on the cell surface to acidic endosomes upon internalization). HEL and SET-2 cells were re-suspended in serum-free, no phenol red RPMI and incubated with DMSO or inhibitors for 30 minutes at 37°C, 5% CO2 prior to the addition of either the CD34 x G6B bispecific or isotype control (5 μg/ml), followed by a further incubation for 30 minutes at 37°C, 5% CO2. For flow cytometry, cells were then washed twice with PBS and re-suspended in PBS for flow analysis on a Beckman Flow CytoFLEX cytometer. For live cell imaging, cells were plated onto slides with flat wells (Hendley-Essex) and imaged on a Zeiss inverted confocal LSM870 with an apochromatic 40X oil immersion objective lens. Two inhibitors were used - Dynasore, a GTPase inhibitor of dynamin (Sigma) at 100 μM and Pitstop 2 (Sigma) at 30 μM concentration, that inhibits the clathrin terminal domain as well as clathrin-independent endocytosis.

44Sequencing data in the binary base call (BCL) format were demultiplexed. Unique molecular identifier (UMI) counts were obtained by aligning FASTQ files to the human reference genome (GRCh38 3.0.0) using Cell Ranger software (version 3.0.1) from 10x Genomics. The CellRanger “count” standard pipeline was used to obtain the expression matrix for each individual library for each donor. Cells meeting the following QC parameters (detailed in Table S3) were included in analyses: UMI counts > 1,000 and ≤ maximum UMIs); number of detected genes > 500 and ≤ maximum number of detected genes); the percentage of mitochondrial gene expression < 15% per cell. Genes expressed in at least 10 cells were included. Following application of these filters, 122,154 cells passed quality control (83,753 cells from myelofibrosis patients and 38,401 cells from healthy donors, see Table S3). We scaled UMI counts by normalizing each library size to 10000. The normalized expression values were then log transformed.

45Sparse expression matrices of cells obtaining from the CellRanger output for individual donor were combined. Highly variable genes were identified by fitting mean expression values and the squared coefficient of variation (CV2) calculated with a gamma generalized linear model using the “glmgam.fit” function in the statmod package in R as described previously. Using a mean expression value of > 0.05 and the dispersion score of > 0.05, 800 genes were identified as highly variable genes. 12 ribosomal, mitochondrial and heat shock protein genes were removed, resulting in 788 highly variable genes being used for down-stream analysis. Donor effect was regressed out using Combat from the sva package by regressing out on the donor IDs. Following this, there was no clear batch effect in the dataset (Figure S1C). Principal Component Analysis (PCA) was then performed on normalized expression values on the first 50 principal components (PCs). An elbow plot was inspected to determine the appropriate number of top PCs capturing the most of variances. Using this approach, we selected the first 20 PCs for further analyses. Uniform Manifold Approximation and Projection (UMAP) analysis was performed using the “uwot” function on the embedded matrix derived from the first 20 PCs and 30 neighbors using cosine as the metric parameter. Cells were clustered using the k-nearest neighbor (KNN) approach, using Euclidean metric as the input parameter. The weighted graph was created with the weight values calculated from the normalized shared number of the nearest neighbors. The function ‘cluster_louvain’ from the igraph package was then applied to identify clusters based on the weighted graph. Identified clusters were superimposed on the two-dimensional UMAP. We identified a distinct cluster (n = 1,958 cells) expressing a gene signature corresponding to plasmacytoid dendritic cells (pDC) (Figure S1D) and we removed this cluster from further analyses as they are not hematopoietic stem/progenitor cells but contaminating mature cells that fall into the CD34+Lin– FACS gate. After removing the pDC cluster, 120,196 cells (82,255 from myelofibrosis donors and 37,941 cells from healthy donors) from all 21 donors (myelofibrosis patients + controls) were included in analyses. Highly variable genes were recalculated, and 817 were identified; 12 of which were removed (heat-shock/mitochondrial/ribosomal genes), resulting in 805 highly variable genes being used for all downstream analyses. Eighteen clusters were initially identified, differentially expressed genes inspected and clusters with similar profiles were merged. This resulted in eight major clusters representing HSPC and lineage progenitor populations, as shown in Figures 2A and 2B.

46Lineage signature gene sets were collated by curating known canonical lineage markers selected from multiple recently published hematopoiesis datasets. Genes selectively expressed in hematopoietic lineage progenitors and uncommitted HSPCs were identified (see Table S5). Expression of each gene was plotted individually on the UMAP plot and genes that were highly expressed and most specifically marked distinct lineage clusters or uncommitted HSPCs in the HSPC UMAP were selected for inclusion in the ‘signature gene sets’ (see Figure S2). These gene sets were used to calculate a ‘lineage gene score’ for each cell, based on the average gene expression of each lineage gene set. These scores were superimposed on the UMAP and FDG plots (Figures 2B and 3).

47Differentially expressed genes for each cluster were identified using the nonparametric Wilcoxon test on the log-transformed, normalized UMIs to compare expression level. Fisher′s exact test was used to compare the cell frequency expressing each gene as previously described. P values generated from both tests were then combined using Fisher’s method and were adjusted using the Benjamini-Hochberg (BH) correction. Genes expressed by each individual cluster were compared to all other clusters and differential expression defined as an absolute log2 fold change of ≥ 0.5 and adjusted P value of < 0.05, with the fraction of expressing cell frequency of > 0.3. Differentially expressed genes were ranked using P values and log2FC to select up to 50 genes per cluster (Tables S4 and S6). Clusters were identified by manual inspection of differentially expressed genes for canonical marker genes of blood cell lineages. All heatmaps show scaled (z-score) expression values.

48To generate a dataset with equal numbers of myelofibrosis and control cells, myelofibrosis cells were ‘down-sampled’ to 37,941 cells. This number is equal to the maximum number of cells that passed QC in the healthy control. These cells were integrated with all the cells from healthy donors, and analyses from normalization to clustering as described above were repeated. 722 highly variable genes were detected and 12 ribosomal, mitochondrial and heat shock protein genes were removed, resulting in 710 highly variable genes used for subsequent analyses. Seventeen clusters were identified initially, clusters expressing erythroid, megakaryocyte, myeloid and lymphoid genes signatures were merged (Figure S1F). The percentage of cells in each of these lineage clusters from each donor type was quantified (Figure 2D).

49Using the same method as described above, cells from the two donor types (healthy donor controls and myelofibrosis patients) were integrated separately. 847 and 680 highly variable genes were detected and 18 and 14 clusters for the MF and control respectively. As shown in Figure S4A, clusters 8, 13, 14, 15, and 16 (n = 8134 cells) were identified as MkP clusters in the myelofibrosis UMAP and clusters 13 and 14 (n = 141 cells) in the healthy donor UMAP. We used the AUCell package to calculate the AUCell score for megakaryocyte signature genes for each individual cell within these clusters. Based on the score distribution, we selected AUCell score > 0.4 (Figure S4A) to define MkP. This resulted in 4,134 and 139 MkP cells being identified from myelofibrosis and healthy control donors respectively. To examine the heterogeneity among myelofibrosis MkP, normalization and clustering were performed as described above. 20 PCs were used for analysis and 1195 highly variable genes were identified after removing 11 ribosomal, mitochondrial, and heatshock protein genes. Clustering analysis revealed 9 sub-clusters of myelofibrosis MkP (Figure 5D).

50After removing plasmacytoid dendritic cells for each individual, a standard pipeline described above was used. 20 PCs and 30 neighbors were used for UMAP and clustering analyses. 20 PCs and 5 neighbors were used for the force-directed graph (FDG) analysis (Figure S3).

51ForceAtlas2 software was used to visualize differentiation trajectories over the force-directed graph layout. The KNN weighted graph was used as the input for the software, and analysis performed with the following parameters: iterations = 1,000, edgeWeightInfluence = 1, barnesHutTheta = 2, scalingRatio = 1, gravity = 0.05 and jitterTolerance = 1. The output layout was plotted using ggplot2 and superimposed with normalized lineage signature gene scores (average expression values) for myeloid, erythroid, lymphoid and megakaryocyte lineages (Table S5). Transparent factors for selected colors were calculated from normalized expression values across cell types. ggplot2 was used, adding the dynamic alpha parameter values to “geom_point” to control the transparency of colors as shown in Figures 3A, 3B, and 3C.

52To identify the differentiation paths, the PAGA function in the Scanpy toolkit was used. Clusters were identified using canonical lineage maker genes as hematopoietic stem cell (HSC), hematopoietic stem/progenitor cell (HSPC1, HSPC2, HSPC3), erythroid (Ery), megakaryocyte (Mega) and myeloid (Mye) (Figure 4A), with HSC clearly positioned at the root of the trajectory. As shown in Figures 4A and 4B, the predominant path for megakaryocyte differentiation was HSC → HSPC2 → Mk.

53Scanpy (version 1.4.5) was used for UMAP, FDG, and diffusion maps analyses (Figure S3D) to compare the results of our in-house analysis pipeline with the existing analysis tool. The ‘sc.pp.combat’ function in scanpy was used to correct donor effect by regression out donor IDs. The ‘sc.pp.highly_variable_genes’ function was used by setting ‘n_top_genes’ to the same number of genes used in the in-house pipeline. 20 PCs and 30 neighbors were used to run UMAP, diffusion map, and FDG functions. Megakaryocyte signature genes (Table S5) were superimposed on the UMAP, diffusion map and FDG.

54To examine expression of 1,639 transcriptional factors along the megakaryocyte and erythroid differentiation trajectories, the 82,255 myelofibrosis HSPCs were clustered by the Louvain community-detection method (resolution = 0.5). Seven major clusters were identified and projected on the FDG layout using Scanpy. Differentially expressed genes for each cluster were identified as described above. 16 transcription factor genes showing progressive changes, either increased or decreased expression, along the two trajectories were selected (Figures 4B and 4C).

55The heatmap was generated by calculating the running average of the gene expression for cells along the path using the function ‘scanpy.pl.paga_path’ in Scanpy to investigate dynamic changes in gene expression. The n.avg parameter was set at 5000 for the number of data points to include in the computation of the running average, to ensure a smooth transition.

56GSEA was performed using GSEA software (https://www.gsea-msigdb.org/gsea/index.jsp) with ‘Run GSEAPreranked’ and default parameters. The HALLMARK gene set used for the analysis was downloaded from MSigDB (https://www.broadinstitute.org/gsea/msigdb/collections.jsp). 9,313 expressed genes were used for the analysis. Genes were pre-ranked in order of their differential expression between myelofibrosis and healthy donor megakaryocyte progenitors. The Pathways with FDR q-value < 0.25 were selected showing in Figure 5B.

5737,941 cells (the down-sampled dataset) from 15 myelofibrosis patients processed using the Chromium platform (10x Genomics) and 2,071 cells from 8 myelofibrosis patients (obtained from; GSE122198) were projected into a shared embedding with Harmony, using the top 20 PCA dimensions. The effect of the platform (10X Genomics/TARGET-seq) and donor were simultaneously accounted for and introduced as covariates. Then, dimensionality reduction and force directed graph analysis were performed as described above (Dimensionality reduction, removal of individual donor effect and cell clustering). To quantify the proportion of mutant and wild-type cells in each lineage progenitor cluster, clusters were identified by inspection of differentially expressed genes and super-imposition of the lineage gene set scores on the UMAP and FDG graphs as described above, and mutant and wild-type cells enumerated for myeloid, eythoid and megakaryocyte progenitor clusters. Proportions of wild-type and mutant cells for each of cluster pair were compared using a Chi-square test (p value) with Yates’ continuity correction.

58Flow cytometry data was analyzed using FlowJo software (v10.5.3). Summary charts and associated statistical analyses were performed using GraphPad Prism (v8.1.0). Helios CyTOF Software (v6.7) was used for processing of FCS 3.0 files, normalization to EQ Beads, concatenation of multiple files and debarcoding. Data was then analyzed and histograms and viSNE plots generated using CytoBank.org.

59Statistical tests used, numbers of replicates and definitions of statistical significance are described in the relevant figure legends. All bar charts show mean ± standard error of the mean and were generated using GraphPad Prism (v.8.1.0).

60To compare the proportions of wild-type versus mutant cells in lineage progenitor clusters (Figure S6A), Fisher’s exact test was used. To compare the expression of key genes in different sub-clones within the same patients (Figures 6B, 6C, and S6B), Fisher’s test was used to compare the expression frequencies between the groups and Wilcoxon rank sum test to compare the expression levels. P values were combined using Fisher’s method, and the combined p value for each pairwise comparison are reported in Figure S6.

6110x Genomics single cell RNA-sequencing data has been submitted to GEO database (Accession Number GEO: GSE144568). TARGET-seq single cell RNA-sequencing data is available at Accession Number GEO: GSE122198. The Shiny application for visualization of the data from patients and healthy donors in this study is available at https://github.com/supatt-lab/SingCellaR-myelofibrosis. R scripts used for the analysis are available upon request.