Different inflammatory, fibrotic, and immunological signatures between pre-fibrotic and overt primary myelofibrosis

Abstract

Primary myelofibrosis (PMF) is a myeloid proliferative neoplasm (MPN) characterized by bone marrow fibrosis. Pre-fibrotic PMF (pre-PMF) progresses to overt PMF. Megakaryocytes play a primary role in PMF; however, the functions of megakaryocyte subsets and those of other hematopoietic cells during PMF progression remain unclear. We, therefore, analyzed bone marrow aspirates in cases of pre-PMF, overt PMF, and other MPN using single-cell RNA sequencing. We identified 14 cell types with subsets, including hematopoietic stem and progenitor cells (HSPC) and megakaryocytes. HSPC in overt PMF were megakaryocyte-biased and inflammation/fibrosis-enriched. Among megakaryocytes, the epithelial-mesenchymal transition (EMT)-enriched subset was abruptly increased in overt PMF. Megakaryocytes in non-fibrotic/non-PMF MPN were megakaryocyte differentiation-enriched, whereas those in fibrotic/non-PMF MPN were inflammation/fibrosis-enriched. Overall, the inflammation/fibrosis signatures of the HSPC, megakaryocyte, and CD14+ monocyte subsets increased from pre-PMF to overt PMF. Cytotoxic and dysfunctional scores also increased in T and NK cells. Clinically, megakaryocyte and HSPC subsets with high inflammation/fibrosis signatures were frequent in the patients with peripheral blood blasts ≥1%. Single-cell RNA-sequencing predicted higher cellular communication of megakaryocyte differentiation, inflammation/fibrosis, immunological effector/dysfunction, and tumor-associated signaling in overt PMF than in pre-PMF. However, no decisive subset emerged during PMF progression. Our study demonstrated that HSPC, monocytes, and lymphoid cells contribute to the progression of PMF, and subset specificity existed regarding inflammation/fibrosis and immunological dysfunction. PMF progression may depend on alterations of multiple cell types, and EMT-enriched megakaryocytes may be potential targets for diagnosing and treating the progression.

Introduction

Philadelphia-negative myeloproliferative neoplasms (MPN) are myeloid hematopoietic stem cell (HSC) neoplasms characterized by the overproduction of myeloid, erythroid, and megakaryocytic cells, resulting in polycythemia vera (PV), essential thrombocythemia (ET), and primary myelofibrosis (PMF).¹ PV, ET, and PMF share the clinical and molecular features that may result in disease progression to acute myeloid leukemia and harbor specific driver mutations (JAK2, CALR, or MPL) that activate JAK2 signaling.² The major criteria for diagnosing PMF include megakaryocyte proliferation with bone marrow (BM) fibrosis, the absence of diagnostic criteria for ET or PV, and the presence of driver mutations or additional high-molecular-risk mutations.^3-5 JAK inhibitors can relieve symptoms in PMF but do not entirely resolve BM fibrosis,⁶ suggesting that the pathogenesis of PMF is more complex than that of other MPN.⁷

PMF arises from a single HSC with driver mutations that endow it with a selective advantage, thereby promoting myeloid cell proliferation and BM fibrosis.⁸ Atypical megakaryocytes are the histological hallmark of PMF and play a vital role in its development.⁴ PMF-megakaryocytes are characterized by enrichment of inflammatory and immunoregulatory signals that alter cross-talk between BM cells, thereby exacerbating genetic instability in PMF-HSC^9,10 and promoting fibrosis.^5,10 However, the specific megakaryocyte subsets relevant to PMF remain undetermined. Furthermore, the contribution of other myeloid cells, besides megakaryocytes, to PMF is largely uncharacterized. Immune evasion and dysregulation of the immune system contribute to the clonal evolution of PMF-HSC and BM fibrosis.^11-13 These immunological alterations have been identified in the peripheral blood (PB) cells of patients with PMF; however, they do not reflect the complete immunological landscape of PMF-BM.

The 2016 World Health Organization classification categorizes PMF into pre-fibrotic PMF (pre-PMF, grade 0 or 1 fibrosis) and overt PMF (overt PMF, grade 2 or 3 fibrosis).³ As constitutional symptoms and hematologic pathology worsen, pre-PMF progresses to overt PMF, resulting in poorer survival.^3,14-16 Subtle genomic differences in the hematopoietic clones (unfavorable karyotypes and high-molecular-risk mutations) and bulk transcriptomic differences in inflammatory signatures have been noted between pre-PMF and overt PMF.^14-16 However, the precise molecular mechanisms underlying the progression from pre-PMF to overt PMF remain unexplored.

JAK-STAT activation has been identified as a driver mechanism for MPN; nevertheless, PMF remains the most heterogeneous disease among MPN and is further complicated by co-existing inflammation.¹⁷ After acquiring the driver mutations, PMF passes decades of latency before the development of disease manifestations.¹⁸ Studies suggest the existence of multiple disease-modifying factors that could result in diverse cell states or subpopulations in the BM cells of patients with PMF.

Single-cell RNA sequencing (scRNA-seq) can define the individual transcriptomes of admixed cells in tissues, further dissecting subpopulations and enabling precise inference of the functions and interactions of the cells that cannot be distinguished by traditional bulk analyses.¹⁹ To the best of our knowledge, only a handful of scRNA-seq datasets for PMF-BM are currently available. A patient with overt PMF was assessed using a BM biopsy, with primary analyses of non-hematopoietic cells.¹⁰ Furthermore, cases of PMF were studied by simultaneous mutation and scRNA-seq analyses; however, the proof-of-principle method utilized in that study could not intersect multiple phenotypic readouts across a range of cell subtypes,²⁰ warranting further investigations. In this study, we hypothesized that specific cell subtypes, besides megakaryocytes, might exist in overt PMF, and that their distinctive molecular signatures might contribute to PMF progression. To elucidate these, we analyzed BM aspirates from pre-PMF and overt PMF using scRNA-seq and identified cellular subpopulations with different inflammatory, fibrotic, and immunological signatures between pre-PMF and overt PMF.

Methods

Bone marrow samples

BM aspirates were collected from 33 patients with MPN (5 with ET, 1 with PV, 5 with pre-PMF, 12 with overt PMF, 6 with post-ET myelofibrosis, and 4 with post-PV myelofibrosis). Among the 12 patients with overt PMF, six were treated with the JAK inhibitor ruxolitinib, while the remaining six were not treated with any JAK inhibitor. To minimize the risk of PB dilution in BM aspirate samples, we repeated the aspiration at a different site if a dry tap was encountered during the procedure. We also performed microscopic examination of the BM aspirates to assess their quality and composition. The patients’ BM fibrosis was confirmed using BM biopsies, indicating that our aspirates represented the remaining fluidic areas surrounding the fibrotic PMF-BM. This study was approved by the institutional review board of the Catholic University of Korea (KC20TISI0206).

Single-cell RNA-sequencing library preparation

Individual cells in BM aspirates were isolated using density gradient centrifugation with a Ficoll-Paque Plus medium (GE Healthcare). After removing red blood cells using a red blood cell lysis buffer (Miltenyi Biotec), the cells were counted and stored at -80°C until their utilization. The scRNA-seq library was prepared using the Chromium instrument system with a Single Cell 3’ v3 Reagent kit (10x Genomics), according to the manufacturer’s protocol. scRNA-seq libraries were sequenced on an Illumina NovaSeq platform. Raw sequencing data have been deposited in the Sequence Read Archive under accession number PRJNA1070224. Additional details of the scRNA-seq library preparation are described in the Online Supplementary Methods.

Single-cell RNA-sequencing data analysis

The sequenced data were processed into expression matrices using CellRanger (10x Genomics). Sequencing reads were mapped to the GRCh38 reference genome. Bioinformatics processing of the scRNA-seq data was performed using Seurat.²¹ After removing low-quality cells, the data were log-normalized, and highly variable features were identified based on a variance stabilizing transformation method. All individual datasets were integrated using Harmony.²² Principal component analysis and graph-based clustering were performed on the integrated datasets, and the clustering data were applied to the uniform manifold approximation and projection. Each cell cluster was annotated for its cell type using SingleR and well-known cell-type-specific markers.

Gene expression analysis was performed to identify significantly differentially expressed genes within each cluster using the ‘FindMarkers’ function in Seurat. Gene signature scores were calculated using UCell. Gene set enrichment analysis was performed using fgsea with MSigDB. Single-cell reference mapping compared the megakaryocyte subset abundance between PMF and other MPN. Receptor-ligand interactions were analyzed using CellChat to examine cellto-cell communication between different cell types in PMF-BM.²³ Additional details of the scRNA-seq data analyses are provided in the Online Supplementary Methods.

Statistical analysis

Fisher exact test was used for categorical variables. The Mann-Whitney U test was used for continuous variables. Statistical analyses were performed using SPSS (IBM). GraphPad Prism software was used to create graphs. Statistical significance was set at a P value <0.05 for all the analyses.

Results

Bone marrow from patients with overt primary myelofibrosis is enriched with megakaryocytes and hematopoietic stem and progenitor cells

We performed scRNA-seq on whole BM mononuclear cells isolated from 17 patients with PMF (5 with pre-PMF and 12 with overt PMF) (Table 1). After quality control, we obtained 98,677 cells from pre-PMF (19,452 cells) and overt PMF (79,215 cells) clustered in 14 cell populations: hematopoietic stem and progenitor cells (HSPC), megakaryocytes, erythroid cells, five myeloid lineages (CD14⁺ monocytes, CD16⁺ monocytes, myeloid dendritic cells, plasmacytoid dendritic cells, and neutrophils), five lymphoid lineages (T cells, NK cells, pre-B cells, B cells, and plasma cells), and mesenchymal stromal cells (Figure 1A).

We observed differences in the abundance of cells between pre-PMF and overt PMF. In overt PMF, we found increased HSPC (9.3% of BM-mononuclear cells) and megakaryocytes (11.7%), compared to those in pre-PMF (HSPC: 3.4%, megakaryocytes: 5.9%) (Figure 1B-D). However, these differences were not statistically significant. On the other hand, the populations of myeloid dendritic cells, neutrophils, and plasma cells decreased significantly as the diseases progressed (Figure 1D). Regarding clinical variables, megakaryocytes and HSPC were enriched in the BM of patients with PMF with PB blasts ≥1% compared to those with PB blasts <1% (megakaryocytes: 16.2% vs. 4.5%, P=0.015; HSPC: 12.4% vs. 3.3%; P=0.027). Megakaryocytes were also enriched in the BM of patients at high risk according to the Mutation-enhanced International Prognostic Score System (MIPSS) compared to those at low or intermediate risk (11.5% vs. 4.2%; P=0.021). In line with a previous finding,²⁴ it was observed that lymphocyte levels decreased in patients treated with ruxolitinib compared to the levels in those not treated. We found no significant differences in cell type abundance according to PMF driver mutations, age, leukocyte count, hemoglobin concentration, platelet count, or prognostic groups of the IPSS and Dynamic IPSS-plus (DIPSS-plus). Taken together, BM cellular abundance of HSPC and megakaryocytes was noted during PMF progression and was associated with an increase in PB blasts.

Hematopoietic stem and progenitor cells in overt primary myelofibrosis are megakaryocyte-biased with inflammatory and fibrotic activation

The sub-clustering of HSPC revealed nine subsets (Figure 2A; Online Supplementary Table S1): HSC (HSC1 and HSC2 expressing AVP and HLF), megakaryocyte-erythroid-mast progenitors (MEMP expressing CLC and FCER1A), megakaryocyte-erythroid progenitors (MEP1 and MEP2 expressing GATA1 and KLF1), early erythroid (expressing ALAS2 and GYPA), granulocyte-monocyte progenitors (GMP expressing MPO and AZU1), lymphoid progenitors (expressing CD247 and THEMIS), and proliferating cells (expressing TOP2A and MKI67).

In HSC1, gene set enrichment analysis identified the enrichment of inflammatory (tumor necrosis factor a signaling and interferon-g response) and fibrotic signaling (transforming growth factor b [TGF-b] and coagulation) (Online Supplementary Figure S1A). Pseudo-bulk differentially expressed gene analysis identified eight overt PMF-specific genes, of which BACH2, ANXA2, and ANO2 were associated with megakaryocyte differentiation (Online Supplementary Figure S1B; Online Supplementary Table S2).²⁵ The abundance of HSC1 was not different between pre-PMF and overt PMF (Figure 2B), suggesting that HSC1 might be related to PMF development.

The MEP1 subset was more abundant in overt PMF than in pre-PMF (11.0% vs. 4.2%; P=0.029) (Figure 2B). Additionally, MEP1 was more abundant in patients with PMF with PB blasts ≥1% than in those with PB blasts <1% (13.6% vs. 4.1%; P=0.001). Gene set enrichment analysis of MEP1 revealed enrichment in megakaryocyte-lineage differentiation (Online Supplementary Figure S1C).²⁶ Eight overt PMF-specific differentially expressed genes were detected in MEP1 (Online Supplementary Figure S1D; Online Supplementary Table S2). There were no differences in HSC1 and MEP1 between the ruxolitinib-exposed and unexposed patients, in terms of quantitative or qualitative (inflammatory and fibrotic signaling) measures (Online Supplementary Figure S1E). Previously, megakaryocyte-biased HSPC expressing PF4, MPIG6B, VWF, SELP, and GP9 were observed in the PB of patients with PMF.⁵ In our study, the proportions of these HSPC were greater in overt PMF than in pre-PMF (27.7% vs. 8.3%, P=3.7x10-⁹⁹) (Figure 2C). Notably, the megakaryocyte signatures were not limited to a specific HSPC subset, suggesting that the megakaryocyte bias may be widespread among HSPC (Figure 2C; Online Supplementary Figure S2).

Table 1.Clinicopathological features of the patients with primary myelofibrosis included in the study.

MPIG6B and PF4 expression was observed in uncommitted (HSC1 and HSC2) and lineage-committed subsets (MEP1 and MEP2) (Figure 2D), consistent with a previous report.⁵ The HSPC expressing MPIG6B (59.6% vs. 28.9%, P=1.7x10^-81) and PF4 (22.0% vs. 10.7%, P=1.5x10^-16) were more prevalent in overt PMF than in pre-PMF (Online Supplementary Figure S2). Collectively, HSC1 and MEP1 subsets in overt PMF showed a greater bias toward megakaryocyte production compared to the same subsets in pre-PMF. These subsets may have a preference for expanding or interacting with other BM cells to promote megakaryocyte differentiation, and BM inflammation and fibrosis.

Figure 1.Single-cell profiling of bone marrow aspirates from patients with primary myelofibrosis. (A) The single-cell RNA sequencing data of bone marrow aspirates from 17 patients with primary myelofibrosis (PMF) were integrated. Two-dimensional uniform manifold approximation and projection (UMAP) visualization of 98,677 bone marrow mononuclear cells identified 14 cell types after unsupervised clustering. Each point represents a single cell and is colored based on cell types. (B) UMAP plot colored by the clinical groups. (C) Cell type composition for each sample. Each box’s color is consistent with the cell type. (D) Box plots representing the proportion of each cell type between pre-PMF (N=5) and overt PMF (N=12). The mean and 95% confidence interval are represented with black lines. HSPC: hematopoietic stem and progenitor cell; PC: plasma cell; MSC: mesenchymal stromal cell; pDC: plasmacytoid dendritic cell; mDC: myeloid dendritic cell; MK: megakaryocyte; NK: natural killer; pre-PMF: pre-fibrotic PMF.

Figure 2.Hematopoietic stem and progenitor cell subsets and their gene signatures. (A) Uniform manifold approximation and projection (UMAP) plot colored by hematopoietic stem and progenitor cell (HSPC) subsets. (B) Box plots representing the proportion of each HSPC subset. The mean and 95% confidence interval are represented with black lines. (C) Megakaryocyte (MK) signature scores are shown in the UMAP plot by clinical groups (upper panels) and the violin plot by HSPC subsets (lower panels). (D) The expression levels of MPIG6B and PF4 are shown by clinical groups. MEP: megakaryocyte-erythroid progenitor; MEMP: megakaryocyte-erythroid-mast progenitor; HSC: hematopoietic stem cell; GMP: granulocyte-monocyte progenitor; PMF: primary myelofibrosis; pre-PMF: pre-fibrotic PMF.

Increased signatures of inflammation and fibrosis in megakaryocytes of overt primary myelofibrosis

We identified five megakaryocyte subsets (MK1-5) (Figure 3A; Online Supplementary Table S1). The proportion of MK5 was higher in overt PMF than in pre-PMF (Figure 3B). We compared the module scores to identify the quality difference among the megakaryocyte subsets and found that megakaryocyte differentiation, fibrosis, TGF-β, and cytokine scores increased from pre-PMF to overt PMF (Figure 3C). These signatures increased in overt PMF of all megakaryocyte subsets, except MK2. Notably, MK5 was most specific to overt PMF and was scarce in pre-PMF (1.3% vs. 0.2%; P=0.044) (Figure 3B). Furthermore, the MK5 subset was enriched in patients with PMF with PB blasts ≥1% compared to those with PB blasts <1% (2.0% vs. 0.12%; P=8.2x10-⁵). MK5 specifically expressed epithelial-mesenchymal transition (EMT)-related genes, including TTK, ITGA6, and ILK (Online Supplementary Figure S3A).²⁷ There were no significant differences in megakaryocyte subsets between patients exposed to ruxolitinib and those not exposed to this drug (Online Supplementary Figure S3B, C).

Figure 3.Megakaryocyte subsets and their gene signatures. (A) Uniform manifold approximation and projection (UMAP) plot colored by megakaryocyte (MK) subsets. (B) Box plots representing the proportion of each MK subset. The mean and 95% confidence interval are represented with black lines. (C) Signature scores are shown in the violin plot by MK subsets. (D) The heatmap shows the 19 differentially expressed genes in MK between pre-PMF and overt PMF. Red and blue colors indicate upregulated and down-regulated genes, respectively. PMF: primary myelofibrosis; pre-PMF: pre-fibrotic PMF; TGFβ: transforming growth factor beta.

We performed a pseudo-bulk gene expression analysis on all megakaryocytes and identified 19 overt PMF-specific differentially expressed genes (Figure 3D; Online Supplementary Table S2). Upregulated genes, including COL24A1, CXCL2, and LTBP1, were related to fibrosis and enriched in the MK3 subset. For example, CXCL2, which is known to be associated with pulmonary fibrosis,²⁸ was nearly absent in pre-PMF but was highly expressed in MK3 of overt PMF (0.9% vs. 10.9% of megakaryocytes; P=2.3x10^-12). LTBP1 is an extracellular matrix protein associated with fibrillin-microfibrils²⁹ and can induce TGF-β activation,³⁰ which is essential for the development of fibrosis in PMF.³¹ LTBP1 expression increased with the progression from pre-PMF to overt PMF (15.7% vs. 38.9%; P=7.8x10^-28).

Furthermore, we projected the megakaryocytes from 16 patients with non-PMF MPN (5 with ET, 1 with PV, 6 with post-ET myelofibrosis, and 4 with PMF) to compare mega-karyocyte subset differences among MPN. Most megakaryo-cytes (96.6%) in the ET/PV cases were assigned to the MK2 subset. In contrast, the megakaryocyte subset distribution of post-ET/PV myelofibrosis was widespread and similar to that of PMF (Online Supplementary Figure S3D). Notably, MK5, the most overt PMF-specific subset, showed a similarly high distribution between overt PMF and postET/PV myelofibrosis (Online Supplementary Figure S3E). Collectively, scRNA-seq identified MK5 as the most overt PMF-specific subset with an increased fibrosis signature.

Overt primary myelofibrosis immune cells show immune dysfunction and suppression signatures

We identified 15 subsets of T and NK cells in the BM aspirates (Figure 4A, B; Online Supplementary Table S1): naïve T (CD8⁺, CD4⁺, and CD3^-CD4⁺), helper T, regulatory T (Treg), mucosal-associated invariant T (MAIT), gd T, cytotoxic memory (Tmem), cytotoxic terminal effector (Teff), CD56^bright NK, NK (NK1, NK2, NK3, and NK4), and proliferating cell subsets. The subset distribution of T and NK cells was not significantly different between pre-PMF and overt PMF, nor between ruxolitinib-exposed and unexposed patients, except for NK1: the NK1 subset was significantly decreased in the ruxolitinib-exposed patients (P=0.009) (Online Supplementary Figure S4A). However, functional module analyses showed a significant increase in cytotoxic scores in most T-cell subsets with the progression from pre-PMF to overt PMF (Figure 4C). Similarly, the dysfunctional module score of T-cell activation, measured using PDCD1, LAG3, TIGIT, CD244, and CTLA4 expression analyses, increased in the gd T and Teff subsets with progression (Online Supplementary Figure S4B). For example, overt PMF exhibited a higher prevalence of LAG3-expressed Teff cells than pre-PMF (P=0.004) (Figure 4D). These were not significantly different according to ruxolitinib exposure.

We identified eight CD14⁺ monocyte subsets among the myeloid lineages (Figure 5A, B; Online Supplementary Table S1). Of these subsets, mono3 increased from pre-PMF to overt PMF (24.5% vs. 31.7% of CD14⁺ monocytes) (Figure 5C). Mono3 may be a variant of the monocytic myeloid-derived suppressor cell (M-MDSC)³² that has immunosuppressive functions. Mono3 expressed M-MDSC markers³² with higher MHC-II expression than that of conventional M-MDSC (Figure 5B). Furthermore, we discovered that the interferon signature in overall CD14⁺ monocytes was significantly higher in overt PMF than in pre-PMF (P=5.9x10^-131) (Figure 5D). The scRNA-seq analysis revealed altered immune and inflammatory signaling in overt PMF compared to pre-PMF, potentially leading to reduced immune activity.

Prediction of cell-cell communication

Ligand-receptor interaction analysis using scRNA-seq revealed a higher number and greater strength of ligand-receptor interactions in overt PMF than in pre-PMF (Online Supplementary Figure S5A). Among the 90 putative interactions identified, 38 were significantly enriched in overt PMF, with ten exclusive interactions (Figure 6). The difference in interactions between pre-PMF and overt PMF was primarily due to increased interactions among HSPC, megakaryocytes, T cells, and monocytes (Online Supplementary Figure S5B), which were identified as overt-PMF-driving cells using scRNA-seq. Overt-PMF-enriched interactions were largely categorized into four signatures (Online Supplementary Figure S6): megakaryocyte differentiation (GP1BA and VWF),⁵ pro-inflammatory/fibrotic signaling (CD34, CD40, EPHB, and TGF-b),^19,31,33,34 immunological effector/dysfunction signaling (PROS, SN, THBS, TIGIT, LCK, CCL, and PAR),^35-37 and tumor-associated signaling (ESAM, JAM, and HSPG).^38,39 Ephrins and Eph receptors are known mediators of fibrosis.³⁴ MK3 was identified as the most prominent source of the ephrin ligand (EFNB2) in overt PMF, affecting T cells (Online Supplementary Figure S6). This interaction was not observed in pre-PMF. CD34 and CD40 signals, which are known to be involved in inflammatory disease development,^19,33 were identified exclusively in overt PMF (Figure 7A, C). MK2 interacted with HSPC, M-MDSC, and myeloid dendritic cells through the TGFB1-(TGFBR1+TGFBR2) axis in pre-PMF (Figure 7B). MK5 strongly expressed TGFB1, and HSPC strongly expressed ACVR1, another receptor of TGF-b signaling, in overt PMF (Figure 7C). This prediction suggests that the TGFB1-(TGFBR1+TGFBR2) axis is activated in prePMF, with MK2 as a hub, whereas the TGFB1-(ACVR1+TG-FBR1) axis is activated in overt PMF, with MK5 as a hub (Online Supplementary Figure S7). TIGIT, a marker of T-cell exhaustion,⁴⁰ inhibits immune cell responses at multiple steps of the tumor‐immunity cycle. TIGIT signaling in pre-PMF was predicted to originate from Tmem, whereas that in overt PMF was predicted to originate from Tmem and Treg (Online Supplementary Figure S6). T-cell communication partners in TIGIT signaling were much more diverse in overt PMF than those in pre-PMF, primarily targeting HSPC and MK3 (Figure 7C).

Discussion

Previous investigations of the molecular pathogenesis of PMF primarily relied on megakaryocyte alterations. This informational gap led us to analyze pre-PMF, overt PMF, and other MPN (ET/PV) using scRNA-seq. Our results indicated that the differences between pre-PMF and overt PMF were attributable to the inflammation/fibrosis and immunological alterations of multiple cellular subsets, rather than one or two. First, specific subpopulations of HSPC, megakaryocytes, monocytes, and lymphoid cells increased during the progression from pre-PMF to overt PMF. Second, pro-inflammatory/fibrotic and immunological dysfunction signatures increased during the progression. Third, no single decisive subpopulation emerged during the progression of pre-PMF to overt PMF. These gradual alterations support the idea that pre-PMF and overt PMF are in a disease continuum, with many disease-progressing factors involved in the pathogenesis.

The DIPSS-plus system uses eight prognostic survival factors of patients with PMF, including age >65 years, constitutional symptoms, and PB blasts ≥1%.⁴¹ PB blast increase is directly associated with leukemic transformation, plausibly indicating poor survival; however, other possibilities have remained uncertain. Approximately 20% of patients with PMF die from leukemic transformation, with most of them succumbing to BM failure and other complications.⁴² We found that the PB blast ≥1% predictor was related to increases in HSPC and megakaryocytes in PMF-BM, particularly in MEP1 and MK5 subsets. MEP1 was significantly enriched in overt PMF, with a corresponding increase in inflammatory and fibrotic functions. MK5 highly expressed EMT-related genes according to our data. EMT and inflammation cooperate in the progression of organ fibrosis.⁴³

Figure 4.T and natural killer subsets. (A) Uniform manifold approximation and projection (UMAP) plot colored by T and natural killer (NK) subsets. (B) Dot plot of subset-specific marker genes per T or NK subset. Dot intensity (from blue to red) represents the average expression value of all cells per T or NK subset, whereas dot size represents the proportion of cells expressing the genes. (C) The expression level of the cytotoxic signature score is shown in the violin plot by T and NK subsets. (D) Box plots represent the proportion of dysfunctional gene-expressing cytotoxic terminal effector cells by clinical group. The mean and 95% confidence interval are represented with black lines. Treg: regulatory T cells; MAIT: mucosal-associated invariant T cells; Tmem: cytotoxic memory T cells; Teff: cytotoxic terminal effector T cells; PMF: primary myelofibrosis; pre-PMF: pre-fibrotic PMF.

Figure 5.CD14⁺ monocyte subsets. (A) Uniform manifold approximation and projection (UMAP) plot colored by CD14⁺ monocyte subsets. (B) Dot plot of subset-specific marker genes per monocyte subset. Dot intensity (from blue to red) represents the average expression value of all cells per monocyte subset, whereas dot size represents the proportion of cells expressing the genes. Monocytic myeloid-derived suppressor cell (M-MDSC) markers are highlighted in a red box. (C) Box plots represent the proportion of the Mono3 (M-MDSC) subset by clinical group. The mean and 95% confidence interval are represented with black lines. (D) The expression level of the interferon signature score is shown in the violin plot by monocyte subsets. PMF: primary myelofibrosis; pre-PMF: pre-fibrotic PMF.

EMT is identified as a hallmark signature of murine PMF,⁴⁴ indicating that the PB blast ≥1% predictor may be associated with inflammation and fibrosis in PMF. These data suggest that the blast increase may result from HSPC and megakaryocyte proliferation in PMF-BM, which is related to non-leukemic hematologic complications rather than the leukemic transition itself.

Somatically mutated megakaryocytes induce or alter the development and progression of MPN, where megakaryocyte-derived TGF-|3 plays a primary role in HSC proliferation and BM fibrosis.³¹ Consistent with this, we observed increased TGF-β signaling in overt PMF compared with pre-PMF, particularly because of high TGFB1 expression in MK5. Megakaryocytes are essential in the pathogenesis of MPN; however, how megakaryocytes differ between yet-to-be-fibrotic ET/PV and PMF and between PMF and post-ET/ PV-MF remains unclear. We found that the proportion of MK5 was much greater in overt PMF than in pre-PMF, and that there was a higher expression of the TGF-β signature in overt PMF. Moreover, MK5 abundance in non-fibrotic ET/ PV was lower than in PMF, whereas post-ET/PV myelofibrosis showed a similar MK5 abundance to that in overt PMF. These findings suggest that MK5 may be one of the key determinants in the development and progression of BM fibrosis, irrespective of MPN subtypes.

Figure 6.Cell-cell communication analyses. A bar plot represents signaling pathways, ranked according to the differences in overall information flow within the inferred networks between pre-fibrotic primary myelofibrosis (pre-PMF) and overt PMF. The overall information flow of a signaling network is calculated by summarizing all the communication probabilities within the network. Signaling pathways that are more enriched in overt PMF are marked in pink, whereas those more enriched in pre-PMF are marked in blue.

Dysregulation of the immune system contributes to the expansion and survival of the neoplastic myeloid clone in MPN.¹³ For example, differentiation of monocytes to dendritic cells is reduced in PB, whereas MDSC increase in PMF.^11,45 Furthermore, reductions in CD56^bright NK cells and CD3⁺ T cells and an increase in Treg cells have been reported in MPN.^12,46 However, our study did not identify altered NK- or T-cell distribution in PMF-BM samples. Instead, we observed increased cytotoxicity and dysfunction scores in the T-cell subpopulations of overt PMF. Ligand-receptor pair analysis showed that TIGIT (TIGIT-NECTIN2 axis) and PAR (mainly GZMA-F2R and GZMA-PARD3 axes) harboring NK- and T-cell suppressive functions were highly expressed in the cytotoxic cells of PMF interacting with HSPC. Furthermore, we observed an increase of M-MDSC (immune-suppressive) in myeloid cells in overt PMF. These results indicate that the PMF-BM microenvironment alters the immune and inflammatory responses of T cells and monocytes, which may lead to reduced immune activities in the BM of subjects with PMF.

JAK inhibitors targeting mutant hematopoietic clones have improved the symptoms, splenomegaly, and survival of patients with PMF. However, these drugs are not capable of curing progressed PMF, particularly when used as monotherapy.⁴⁷ Inflammation and fibrosis are drivers of PMF pathogenesis; therefore, novel drugs targeting these factors, along with immune modulation, are currently being investigated in clinical trials.¹⁷ In this sense, our molecular data may precisely identify target subpopulations and molecules, enabling the combination of targeted therapy with JAK inhibitors.

Simultaneous BM aspirate and biopsy analysis is essential for a precise delineation of the interactions between stromal cells and hematopoietic cells in the BM. However, our study focused on BM aspirates, primarily because of the challenges of obtaining simultaneous BM biopsies suitable for scRNA-seq. Consequently, the analysis of non-hematopoietic BM cells involved in PMF, such as fibroblasts and myofibroblasts, was precluded. It is likely that the BM aspirates were hemodiluted with contributions from the PB, probably due to extramedullary hematopoiesis. Current consensus suggests that extramedullary hematopoiesis in PMF results from the sequestration, accumulation, and proliferation of circulating progenitor cells.^48,49 Thus, even with hemodilution, the aspirates would still reflect the neoplastic nature of the BM in PMF patients. However, the impact of extramedullary hematopoiesis on the BM microenvironment in PMF requires further investigation. Our study was limited by the complexity of treatments received by the patients, influenced by the rarity of the disease, its broad spectrum, and the variability of therapeutic options. Many patients with overt PMF in our study were treated with a JAK inhibitor (ruxolitinib), which would be expected to alter inflammatory signaling in these patients. However, our findings in overt PMF were predominantly due to PMF progression rather than to ruxolitinib treatment. The small sample size may have led to an underestimation of the impact of ruxolitinib treatment. Therefore, further studies using extended serial BM sampling in a larger cohort of uniformly treated patients are necessary to identify preventive and therapeutic targets for PMF. Additionally, the scRNAseq methodology we employed did not allow genotyping of driver mutations; therefore, the observed differences between pre-PMF and overt PMF could actually be due to the size of the mutated clones. Finally, there were no BM aspirate controls from healthy individuals.

Figure 7.Representative interaction pathways. (A, B) Circle plots represent the inferred interaction pathways of (A) CD34 and CD40, which were exclusively identified in overt primary myelofibrosis (PMF), and those of (B) transforming growth factor-β, identified in both pre-PMF and overt PMF. The edge width represents the communication probability (strength of the interactions) between cell populations. Edge colors are consistent with the signaling source. (C) The expression of ligands and receptors for EPHB, CD34, CD40, TGF-β, and TIGIT in each cell subset from cases of pre-PMF (blue) and overt PMF (red) is shown. Genes corresponding to ligands in each signaling pathway are indicated in bold. Major sources and targets of each signaling pathway are highlighted with boxes. The MK5 subset of megakaryocytes from pre-PMF samples was excluded because of the low number of cells. Details are provided in the Online Supplementary Methods. MK: megakaryocyte; GMP: granulocyte-monocyte progenitor; MEP: megakaryocyte-erythroid progenitor; MEMP: megakaryocyte-erythroid-mast progenitor; HSC: hematopoietic stem cell; MSC: mesenchymal stromal cell; NP: neutrophil; mDC: myeloid dendritic cell; CD16⁺: CD16⁺ monocyte; PC: plasma cell; ERY: erythroid; M-MDSC: monocytic myeloid-derived suppressor cell; CD14⁺: CD14⁺ monocyte; CD56^bright: CD56^bright natural killer cell; NK; NK: natural killer cell; Teff: cytotoxic terminal effector T cell; T: Tmem: cytotoxic memory T cell; MAIT: mucosal-associated invariant T cell; Treg: regulatory T cell.

In summary, our study revealed the single-cell transcriptome signatures and cellular subsets of megakaryocytes, HSPC, and immune cells, characterized according to PMF progression. No overt PMF-specific cell subset emerged during the progression. PMF progression may rely on multiple cell type alterations. Megakaryocytes, HSPC, monocytes, and lymphoid cells contributed to the progression, and there was subset specificity regarding inflammation/ fibrosis and immunological dysfunction. Our results may aid in defining the molecular diagnosis for PMF progression and discovering potential target subsets in PMF, such as the EMT-enriched MK5 subset.

Footnotes

• Received April 12, 2024
• Accepted September 27, 2024

Correspondence

*S-HJ and S-EL contributed equally as first authors

S.H. Lee suhulee@catholic.ac.kr

Y-J. Chung yejun@catholic.ac.kr

S-H. Jung hyun@catholic.ac.kr

Disclosures

No conflicts of interest to disclose.

Funding

This study was supported by grants from the National Research Foundation of Korea (2017R1A2B2002314, 2019R1A5A2027588, 2017R1E1A1A01074913, RS-2022-00165497, RS-2024-00450891, and 2019R1C1C1004909). We also appreciate support from the Basic Medical Science Facilitation Program through the Catholic Medical Center of the Catholic University of Korea funded by the Catholic Education Foundation, and KREONET (Korea Research Environment Open NETwork), which is managed and operated by KISTI (Korea Institute of Science and Technology Information).