Characterization of the cellular and transcriptomic landscapes of natural killer cells in hematopoietic microenvironment of myelodysplastic syndrome

Abstract

BACKGROUND

Myelodysplastic syndrome (MDS) is a common hematological malignancy, which gradually worsens in patients with chronic cytopenias owing to progressive bone marrow (BM) failure and/or transformation into acute myeloid leukemia.

AIM

To verify the systematic and detailed information of natural killer (NK) cell-mediated pathogenesis in the BM microenvironment.

METHODS

Herein, we enriched BM mononuclear cells from 19 MDS patients and 19 healthy donors (HD), and analyzed the immunophenotype and content of the corresponding resident NK cells (rMDS-NKs, rHD-NKs) using flow cytometry analysis and cell counting. Then, we utilized our well-established “3ILs”-based (rhIL-2, rhIL-15, rhIL-18) strategy for ex vivo amplification and dissection of the expanded NK cells (eMDS-NKs, eHD-NKs). Cell apoptosis and cell cycle analyses were conducted for the comparison of cellular viability, and coculture with diverse tumor cell lines was employed for cytotoxicity assessment. Besides cellular analysis, we performed RNA-sequencing and multifaceted bioinformatics analyses for transcriptomic signature evaluation of eMDS-NKs and eHD-NKs, including gene expression profiling and genetic spectrum.

RESULTS

As to resident NK cells in BM mononuclear cells, resident MDS-NKs showed significant decline in the content of NKG2D⁺ activated and CD25⁺ memory-like NK cell subsets, while minimal differences were observed in total NK cells (CD3^-CD56⁺), total activated NK cells (CD3^-CD56⁺CD16⁺) or relative activated NK cell subsets (NKp44⁺, NKp46⁺). As to expanded NK cells after 14-day’s culture, eMDS-NKs revealed impaired amplification and decrease in the contents of the indicated NK cell subsets (CD25⁺, NKp44⁺, NKp46⁺). Meanwhile, the variations in cellular viability and cytotoxicity were identified, particularly the alterations in cell cycle of eMDS-NKs compared to eHD-NKs. Furthermore, the similarities and diversities of gene expression pattern and genetic spectrum of eMDS-NKs and eHD-NKs were also verified.

CONCLUSION

Collectively, our data indicated the multidimensional alterations and the concomitant conservations of MDS-NKs compared to HD-NKs both at the cellular and transcriptomic levels. These findings would provide novel references for further illuminating the NK cell-mediated pathogenesis of MDS and facilitating the development of diagnostic and therapeutic options for MDS in future.

Key Words: Myelodysplastic syndrome; Natural killer cells; Bone marrow blood; Cell viability; Transcriptomic signatures

Core Tip: Myelodysplastic syndrome (MDS) are myeloid malignancies of hematopoietic stem cells with limited treatment options. Longitudinal studies have highlighted the tumor microenvironment and dysregulated immune response in MDS, yet the pathogenesis of natural killer cells (NKs) are largely obscure. Therewith, we dissected the biological and transcriptomic signatures of NKs in bone marrow of MDS patients, which would facilitate the further exploration of MDS-NK-meditated pathogenesis and treatment options in future.

INTRODUCTION

Myelodysplastic syndrome (MDS) is recognized as a heterogeneous clonal disorder that primarily affects elderly people but also children and young adults, which is characterized by dysplastic features in hematopoietic precursors, deregulation of apoptosis, peripheral blood cytopenias, together with the increased risk for transforming into acute leukemia[1]. Meanwhile, patients with other malignancies that receive radiation therapy and intensive chemotherapy are reported with a high incidence of MDS including therapy-related MDS[2]. Worse still, despite the continuous progress and significant strides in fundamental and clinical investigations, none therapeutic options have been demonstrated to prolong survival apart from allogeneic hematopoietic stem cell transplantation (HSCT)[2]. Hence, to fulfill the robust demands for developing comprehensive and novel therapeutic approaches, there’s an urgency to deepen our understanding upon the leukemogenesis mechanisms and prognostic classifications of MDS and to develop valuable biomarker candidate for precise diagnosis in clinical practice[3].

Natural killer (NK) cells are lymphocytes that belong to innate lymphoid cells characterized by the unique immune response to microbiota and cancer cells without priming and sensitization[4,5]. For decades, we and other investigators in the field are dedicated to fulfilling the demands for large-scale preparation of NK cells and dissecting the cellular and omics signatures of NK cells in the setting to physiological or pathological status. On the one hand, longitudinal studies have reported the high-efficient preparation of NK cells and chimeric antigen receptor-transduced NK cells for therapeutic purposes from multifarious sources, including bone marrow blood[6], peripheral blood[7], umbilical cord blood[6,8], placental blood[8], cell lines[9], CD25⁺ memory-like NK cells[10,11] and even stem cells[12,13]. On the other hand, the pathogenic role of NK cell-mediated immunodysfunction and the interaction with hematopoietic stem cells (HSCs) and relative cells (e.g., mesenchymal stem/stromal cells, dendritic cells, myeloid derived suppressor cells) in the bone marrow microenvironment (BME) of patients with hematologic malignancies (e.g., acute myeloid leukemia, MDS) and relative hematological disease (acquired aplastic anemia) have been extensively illuminated[14-19]. These data consistently put forward the remarkable prospects of NK cells in disease diagnosis and treatment.

Nowadays, impaired NK cell surveillance has been recognized as a hallmark of MDS, which commonly results in the evolution of pathophysiology and disease progression[20]. With the aid of fluorescence in situ hybridization, investigators verified the involvement of NK cells in patients with MDS carrying monosomy 7[21]. Very recently, Xu et al[22] identified the increase expression of hypoxia-inducible factor-1α in peripheral blood and bone marrow of MDS mice for the inhibition of NK cell differentiation (e.g., NKG2D, NKp44, DNAM-1) and cytotoxicity (e.g., KIR expression, and abnormalities in interferon-γ secretion) via modulating JAK1/STAT5/SOCS2 pathway. However, the systematic and detailed characterization of the landscape of cytogenetic abnormalities of NK cells in BME for further unveiling the pathogenesis of MDS is largely unavailable, which in return hinders the development of disease diagnosis and treatment strategies for MDS.

In this study, we primarily compared the cellular phenotypes of both resident and expanded NK cells between MDS patients and healthy donors (HDs), and particularly the decline in the indicated NK cell subsets [NKG2D⁺ and CD25⁺ subsets in resident MDS-NKs (rMDS-NKs); CD25⁺, NKp44⁺, NKp46⁺ subsets in expanded MDS-NKs (eMDS-NKs)] as well as the alterations in cellular viability (e.g., decreased ex vivo expansion, the disturbed sub-stage distribution of cell cycle) and cytotoxicity against Naml6 cells. Besides the cellular signatures, we further verified the similarities and diversities between eMDS-NKs and expanded HD-NKs (eHD-NKs) from the aspects of gene expression profiling and genetic variations. Collectively, these findings outlined the key features of MDS-NKs in BME (cellular immunophenotype, cell viability, cytotoxicity and transcriptome), which would collectively benefit the further exploration of NK cell-based diagnosis and intervention for MDS in future.

MATERIALS AND METHODS

Patients and enrollment

To compare the biological and transcriptomic signatures of HD-NKs and MDS-NKs, we collected bone marrow blood from 19 MDS patients and 19 HDs with the consent and the approval of the Ethics Committee of The Fourth People’s Hospital of Jinan Affiliated to Shandong Second Medical University according to the guideline of Helsinki Declaration (LL20240077, LL20240071). All the MDS patients were with initial diagnosis and initial treatment, and identified according to the inclusion and exclusion criteria (World Health Organization classification, Revised International Prognostic Scoring System risk stratification, and mutation status). The detailed information of HDs and MDSs was available in Supplementary Table 1.

Mononuclear cells isolation from bone marrow

As mentioned above, human bone marrow blood was collected from The Fourth People’s Hospital of Jinan Affiliated to Shandong Second Medical University with the consent of 19 HD and 19 primary MDS patients with the approval of the Ethics Committee. For mononuclear cells (MNCs) isolation from 3-5 mL from blood samples in sterile anti-coagulant tubes, we conducted the Ficoll-based density gradient centrifugation according to the manufacturer’s instructions as we previously described[6,7]. The information of HD and MDS was listed in Supplementary Table 1.

NK cell activation and expansion from MNCs

For in vitro NK cell activation and expansion, 5-10 × 10⁵ HD-MNCs or MDS-MNCs were seeded into 96-well plates in 200 μL NK MACS basal Medium (Miltenyi Biotech, Germany) supplemented with 5% human AB serum, 1% penicillin/streptomycin, 1% NK MACS supplement, and the “3ILs” cocktails (rhIL-2, 100 ng/mL; rhIL-15, 10 ng/mL; rhIL-18, 50 ng/mL)[7] at 37 °C and 5% CO₂ for 14 days. The culture medium was half replaced every two days for 7 days, followed by completely replaced every day for the next 7 days[6,23]. Trypan blue staining was conducted for cell counting, and the number of total CD3^-CD56⁺ NK cells and the activated subsets (CD16⁺, NKG2D⁺, CD25⁺, NKp44⁺, NKp46⁺) was calculated according to of total living cells and the percentages of the corresponding cell population by flow cytometry (FCM) analysis at the indicated time points (day 0, 14). The detailed information of the indicated cytokines was shown in Supplementary Table 2.

FCM assay

FCM assay was performed for the resident NKs at day 0 [resident HD-NKs (rHD-NKs), rMDS-NKs] and the expanded NKs at day 14 (eHD-NKs, eMDS-NKs) of NK cell induction. Briefly, the MNCs (day 0) and MNC-derived cells (day 14) were harvested and resuspended in 1 × phosphate-buffered saline (Solarbio, Beijing, China) for preparation. Then, the enriched cells were incubated with the fluorescence-conjugated antibodies in dark for 30 minutes as we previously reported[18]. Finally, the incubated cells were washed for further analysis by FACS Canto II (BD Biosci, CA, United States) and the FlowJo 10.0 software (Tree Star, OR, United States). The information of the fluorescence-conjugated antibodies was available in Supplementary Table 3.

Cell cycle and cell apoptosis assessment

The cells at day 14 were harvested by 300 × g centrifugation, and the fractions of eHD-NKs and eMDS-NKs were gated by co-incubation with anti-CD3-PE and anti-CD56-APC antibodies. For cell cycle analysis, the cells were turned to 70% (v/v) ethanol (Thermo Fisher Scientific, MA, United States) and propidium iodide staining solution (BD Pharmigen, CA, United States) as reported[18]. Finally, the sub-stage distribution of eHD-NKs and eMDS-NKs was verified by BD LSR II (BD Biosci, CA, United States) and the ModFit software (Verity Software House Co. Ltd, ME, United States). For cell apoptosis analysis, 1 × 10⁶ cells at day 14 were harvested and incubated with the Annexin V Apoptosis Detection Kit (Tianjin Sungene Biotech, Tianjin, China) according to the manufactures’ instructions[6,24]. The percentages of apoptotic cells in eHD-NKs and eMDS-NKs were verified by FACS Canto II (BD Biosci, CA, United States) and the FlowJo 10.0 software (Tree Star, OR, United States).

Cytotoxicity detection of NK cells

The cytotoxicity of eHD-NKs and eMDS-NKs against K562 (human myeloid leukemia cell line), Nalm6 (human lymphocytic leukemia cell line) and U937 (human myeloid leukemia cell line) cell lines was performed as recently described[7,24]. In brief, after incubation with the CellTrace Violet reagent (BV421 Laser channel, Invitrogen, MA, United States), the indicated tumor cell lines were solely (negative control) or co-cultured with eHD-NKs or eMDS-NKs at an effector-to-target ratio (E:T = 1:5) for 6 hours in a 96-well concave plate at 37 °C in a 5% CO₂ incubator. Then, the cells were incubated with the fluorescence-conjugated antibodies (anti-CD3-FITC, ant-CD56-APC, anti-CD-107-PE, and 7-AAD). After mixing with Precision Count Beads (BioLegend, CA, United States), the cells were analyzed by FACS Canto II (BD Biosci, CA, United States) and FlowJo 7.0 software (Tree Star, OR, United States). Cytotoxicity of eHD-NKs or eMDS-NKs was calculated according to the formula: Cytotoxicity = (1 - N₁/N₀) × 100%. N₁ and N₀ represent the total number of living K562 cells in the experimental group and the negative control group, respectively[6,7].

RNA-sequencing and bioinformatics analysis

The eHD-NKs and eMDS-NKs at day 14 were harvested and lysed for total mRNA preparation with the TRIZol reagent (ThermoFisher, United States) according to the manufactures’ instructions. After quality test, mRNAs were sent to BGI Genomics (Shenzhen, China) for RNA-sequencing (RNA-SEQ), and the multifaceted bioinformatics analyses were accomplished with the online platforms (https://biosys.bgi.com/#/report/Login) and KEGG website for signaling pathway analysis (https://www.kegg.jp/)[25,26]. The detailed information of the RNA-SEQ data was available in Supplementary Table 4.

Statistical analysis

Statistical analysis was performed using Graph Pad Prism 7.0 software. Independent sample t test was used for the comparison between the two groups, and one-way analysis of variance was used for comparisons of three groups. The date were shown as mean ± SD or mean ± SEM (n ≥ 3 independent repeated experiments). All values meet the normality and variance equality, and multiple-comparison correction was conducted for multiple NK-cell subsets, time points, viability measures, cell-cycle stages, and cytotoxicity outcomes. For each sample, at least n = 3 repeated wells were performed, and the average value was used for further statistical analysis between independent patient samples. Statistically significant difference was defined only when P < 0.05 (^aP < 0.05, ^bP < 0.01, ^cP < 0.001, NS, not significant).

RESULTS

The rMDS-NKs revealed decline in NKG2D⁺ and CD25⁺ NK cell subsets

To verify the potential variations of immunophenotyping between rDMS-NKs and rHD-NKs, we turned to bone marrow blood-enriched MNCs (BM-MNCs) and detected the content of total NKs and the concomitant subsets by FCM analysis. Accordingly, no significant differences were observed in the content of the CD3^-CD56⁺ total NK cells in the enriched BM-MNCs of rDMS-NKs and rHD-NKs (P > 0.05) (Figure 1A and B). As shown in Figure 1C and D, the proportions of the NKG2D⁺ activated subset and the CD25⁺ memory-like subset were consistently decreased in rDMS-NKs compared to those in rHD-NKs, whereas no differences were verified as to the CD3^-CD56⁺CD16⁺ total activated NKs or the relative activated subsets (CD25⁺, NKp44⁺, NKp46⁺). Additionally, excepted for the increased proportion of CD3⁺CD56⁺ NKT cells in rMDS-NKs, minimal differences were observed in the content of relative immune cells in MNCs such as CD3⁺ total T lymphocytes and the subsets (CD4⁺ T cells, CD8⁺ T cells, CD4⁺Foxp3⁺ Treg cells) (Figure 1E). Taken together, these data indicated the similarities in the content of total NK cells and the indicated subsets between rDMS-NKs and rHD-NKs in BM-MNCs except for the moderate decrease of the NKG2D⁺ activated subset and the CD25⁺ memory-like subset in rDMS-NKs.

Figure 1 Comparison of the immunophenotype and components between resident myelodysplastic syndrome natural killers and resident healthy donor natural killers. A and B: Representative flow cytometry diagrams (A) and statistical analysis (B) of CD3^-CD56⁺ total natural killers (NK) cells in the indicated groups [resident myelodysplastic syndrome NKs (rMDS-NKs), resident healthy donor NKs (rHD-NKs)]; C and D: Representative flow cytometry diagrams (C) and statistical analyses (D) of total activated CD3^-CD56⁺CD16⁺ NK cells and the indicated subsets (NKG2D⁺, CD25⁺, NKp44⁺, and NKp46⁺) in the indicated groups (rMDS-NKs, rHD-NKs); E: Statistical analyses of CD3⁺ total T lymphocytes, CD3⁺CD56⁺ NKT cells, and the T cell subsets (CD3⁺CD4⁺, CD3⁺CD8⁺, and CD4⁺FoxP3⁺ T cells) in the indicated groups (rMDS-NKs, rHD-NKs). Data were shown as mean ± SEM (n = 3 repetitions for each sample). ^aP < 0.05, ^bP < 0.01, NS: Not significant. rMDS-NKs: Resident myelodysplastic syndrome natural killers; rHD-NKs: Resident healthy donor natural killers; NK: Natural killer.

The eMDS-NKs showed impaired ex vivo expansion and relative activation of NK cell subsets

Aiming to further dissect the similarities and differences in content and phenotypes between eMDS-NKs and eHD-NKs, we took advantage of our well-established “3ILs”-based protocol for ex vivo NK cell expansion and activation from BM-MNCs of MDS patients and HDs[7,8]. As shown in Figure 2A, we observed the decrease in the generation of NK cell aggregates in the MDS-NKs group compared to that in the HD-NKs group after the 14-day’s ex vivo expansion. In consistence, the fold change of MDS-NKs at day 14 (eMDS-NKs) in total cell number was significant decreased compared to the HD-NKs (Figure 2B). According to the FCM analysis, we further verified the sharply decline in the percentage of total CD3^-CD56⁺ NK cells in the eMDS-NKs group compared to that in the eHD-NKs group (Figure 2C and D).

Figure 2 Comparison of the immunophenotype and components between expanded myelodysplastic syndrome natural killers and expanded healthy donor natural killers. A: Representative phase contrast images of myelodysplastic syndrome natural killers (MDS-NKs) and healthy donor natural killers (HD-NKs) at the indicated time points (day 0, 7, 14) during the 14-day’s ex vivo expansion and activation. Scale bar = 100 μm; B: Fold change of total natural killers (NK) cells in the MDS-NKs and HD-NKs groups; C: Representative flow cytometry diagrams of CD3^-CD56⁺ total NK cells in the MDS-NKs and HD-NKs groups; D: Statistical analyses of CD3^-CD56⁺ total NK cells between the indicated groups (resident MDS-NKs vs resident HD-NKs, expanded MDS-NKs vs expanded HD-NKs); E and F: Representative flow cytometry diagrams (E) and statistical analyses (F) of total activated CD3^-CD56⁺CD16⁺ NK cells and the indicated subsets (NKG2D⁺, CD25⁺, NKp44⁺, and NKp46⁺) in the indicated groups. Data were shown as mean ± SEM (n = 3 repetitions for each sample). ^aP < 0.05, ^bP < 0.01, ^cP < 0.001, ^dP < 0.0001, NS: Not significant. MDS-NKs: Myelodysplastic syndrome natural killers; HD-NKs: Healthy donor natural killers; NK: Natural killer; eMDS-NKs: Expanded myelodysplastic syndrome natural killers; eHD-NKs: Expanded healthy donor natural killers; rMDS-NKs: Resident myelodysplastic syndrome natural killers; rHD-NKs: Resident healthy donor natural killers.

Subsequently, we tried to compare the phenotypes of NK cell subsets between the indicated groups. As shown by the FCM diagrams and statistical analysis, there were minimal differences in the proportions of CD3^-CD56⁺CD16⁺ total activated NK cells between eHD-NKs and eMDS-NKs (Figure 2E and F). As to the relative NK cell subsets, moderate decreases were observed in the proportion of CD25⁺ memory-like NK cells and the activated NKp44⁺ subset, while there were no statistical differences in the relative activated NKG2D⁺ or NKp46⁺ subsets (Figure 2E and F). Collectively, besides the variations in resident NK cells between rHD-NKs and rMDS-NKs, we further confirmed the impaired ex vivo expansion and the decline in the content of CD25⁺ and NKp44⁺ NK cell subsets in eHD-NKs and eMDS-NKs.

The eMDS-NKs revealed certain variations in cellular viability and conservation in cytotoxicity upon tumor cell lines

Having verified the multifaceted variations of both resident and expanded MDS-NKs, we initially compared the cellular viability between eHD-NKs and eMDS-NKs. With the aid of FCM analysis, we found there were no significant differences in the percentages of Annexin V⁺/7-AAD⁺ apoptotic cells between eMDS-NKs and eHD-NKs (P > 0.05) (Figure 3A and B). According to cell cycle analysis, eMDS-NKs revealed moderate increase in the subset distributed in the G0/G1 stage, whereas with decrease in the G2/M stage compared to the eHD-NKs group (Figure 3C and D). Therefore, besides the decline in ex vivo expansion, eMDS-NKs exhibited alterations in cellular viability like the distribution of the cell cycle sub-stages.

Figure 3 Cellular viability and cytotoxicity of expanded myelodysplastic syndrome natural killers and expanded healthy donor natural killers upon tumor cell lines. A and B: Representative flow cytometry diagrams (A) and statistical analyses (B) of the apoptotic subsets in expanded myelodysplastic syndrome natural killers (eMDS-NKs) and expanded healthy donor natural killers (eHD-NKs); C and D: Representative flow cytometry diagrams (C) and statistical analyses (D) of cell cycle-related sub-stage distribution (G0/G1, S, G2/M) eMDS-NKs and eHD-NKs; E and F: Representative diagrams (E) and statistical analyses (F) of CD107a⁺ natural killer cells in the indicated groups (Ctr, K562, Nalm6, U937) at an effector-to-target ratio (1:5); G and H: The ex vivo cytotoxicity of eMDS-NKs and eHD-NKs according to the BV421-based fluorescence of the remaining surviving tumor cells. Data were shown as mean ± SEM (n = 3 repetitions for each sample). ^bP < 0.01, NS, not significant. eMDS-NKs: Expanded myelodysplastic syndrome natural killers; eHD-NKs: Expanded healthy donor natural killers.

To further verify the cytotoxicity of eHD-NKs and eMDS-NKs, we conducted the classical coculture model of the indicated NK cells with diverse tumor cell lines at the effector-target ratio of 1:5 as previously reported. Compared to the eHD-NKs group, eMDS-NKs showed a moderate lower percentage of CD107⁺ subset against Naml6 cells (a human B-cell precursor leukemia cell line), whereas no statistical differences were observed against the other tumor cell lines, including K562 cells (a human chronic myeloid leukemia cell line) and U937 cells (a malignant tumor cell line of human lymphoid tissues) (Figure 3E and F). According to the BV421-based fluorescence of the remaining surviving tumor cells, we found the cytotoxicity of eHD-NKs and eMDS-NKs revealed minimal differences against the indicated tumor cell lines (P > 0.05) (Figure 3G and H). Taken together, these data indicated the variations in cellular viability (e.g., proliferation, cell cycle) and the relative conservation in ex vivo cytotoxicity upon diverse tumor cell lines. In detail, these data supported relatively preserved tumor-cell killing under these assay conditions rather than a broad impairment of cytotoxic function according to the actual cytotoxicity measurements between eHD-NKs and eMDS-NKs.

The gene expression profiling of eHD-NKs and eMDS-NKs showed multidimensional similarity and diversity

Having clarified the phenotypical and biofunctional characteristics of eHD-NKs and eMDS-NKs at cellular level, we conducted RNA-SEQ and multidimensional bioinformatics analyses upon three representative samples from each group for further dissect their transcriptomic signatures. For the purpose, we collected eHD-NKs and eMDS-NKs at day 14 of culture from 3 MDS patients (denoted as MDS-NK-1, MDS-NK-2, MDS-NK-3) and 3 HDs (denoted as HD-NK-1, HD-NK-2, HD-NK-3) for RNA-SEQ detection, and a total number of 16776 genes were enriched (Supplementary Table 4). According to the box plot of gene expression and bar charts of RPKM (short for Reads per Kilobase per Million mapped reads) distribution, both eHD-NKs and eMDS-NKs revealed similarities and conservations in overall gene expression pattern (Figure 4A and B). As shown by principal components analysis, we intuitively observed the typical genetic relationship of individual eHD-NKs and eMDS-NKs, which was further confirmed by the Pearson correlation analysis (Figure 4C and D).

Figure 4 The gene expression pattern of expanded myelodysplastic syndrome natural killers and expanded healthy donor natural killers. A and B: Box plot (A) and RPKM (short for Reads per Kilobase per Million mapped reads) diagram (B) showed the distribution of genes expression pattern in expanded myelodysplastic syndrome natural killers [myelodysplastic syndrome natural killers-1 (MDS-NKs-1), MDS-NK-2, MDS-NK-3] and expanded healthy donor natural killers [healthy donor natural killers-1 (HD-NK-1), HD-NK-2, HD-NK-3]; C and D: Principal components analysis (short for principal component analysis) and Pearson correlation analysis (D) of the indicated NK cells; E and F: Volcano plot diagram (E) and bar chart (F) showed the down-regulated differentially expressed genes (DEGs) (down), the up-regulated DEGs (up), and the non DEGs (no-DEG) between expanded myelodysplastic syndrome natural killers and expanded HD-NKs; G: HeatMap diagram revealed the hierarchical cluster pattern of DEGs between expanded myelodysplastic syndrome natural killers and expanded HD-NKs. MDS-NKs: Myelodysplastic syndrome natural killers; HD-NKs: Healthy donor natural killers; DEG: Differentially expressed genes.

To further verify the transcriptomic signatures between eHD-NKs and eMDS-NKs, we turned to volcano plot analysis and observed the distribution of 95 differentially expressed genes (DEGs) and the 16681 non-DEGs (Figure 4E and Supplementary Table 4). Of the 95 DEGs, a total number of 35 genes were upregulated and 60 ones were downregulated (|log₂ (fold change)| > 1, false discovery rate < 0.05) in eMDS-NKs compared to the eHD-NKs group (Figure 4F and Supplementary Table 5). Furthermore, with the aid of hierarchical cluster analysis based on FPKM values of the indicated DEGs, we verified that eHD-NKs (MDS-NK-1, MDS-NK-2, MDS-NK-3) and eMDS-NKs (MDS-NK-1, MDS-NK-2, MDS-NK-3) exhibited preferable intergroup relatedness rather than the relatedness between groups (Figure 4G and Supplementary Table 5). Taken together, these findings reflected the conservations in overall gene expression pattern and the concomitant multidimensional alterations in gene expression profiling.

The alterations of biofunctional landscape and genetic spectrum of eMDS-NKs

For the purpose of further illuminating the potential variations of eMDS-NKs in biofunctional landscape, we conducted GOBP (short for gene ontology biological process) analysis of the DEGs, and found that these genes were mainly involved in subcellular structure such as plasma membrane (e.g., integral component of plasma membrane, external side of plasma membrane, plasma membrane raft) and microtubule (e.g., astral microtubule, spindle midzone, mitotic spindle astral microtubule) (Figure 5A). As shown by KEGG analysis of the signaling pathways, the DEGs between eHD-NKs and eMDS-NKs were mainly involved in antigen processing and presentation and diverse metabolism (Figure 5B). From GSEA (short for gene set enrichment analysis) diagrams of eHD-NKs and eMDS-NKs, we further identified the enrichment of gene sets related with certain bioprocesses such as ribosome biogenesis, DNA repair and cell cycle (e.g., cell cycle G2/M transition, G2M checkpoint) (Figure 5C).

Figure 5 differentially expressed gene-based bioprocesses and somatic variations between expanded myelodysplastic syndrome natural killers and expanded healthy donor natural killers. A and B: GOBP (short for gene ontology biological process) analysis (A) and KEGG (short for Kyoto encyclopedia of genes and genomes) signaling pathway analysis (B) of the differentially expressed genes between expanded myelodysplastic syndrome natural killers and expanded healthy donor natural killers according to the website (https://www.kegg.jp/); C: GSEA (short for gene set enrichment analysis) showed the enriched gene sets involved in the indicated bioprocesses; D and E: The distribution of genes with the indicated somatic variations, including variable shear events (D) and differentially variable shear events (E); F: Circos diagrams showed the distribution of genes with somatic variations in chromosomes of expanded myelodysplastic syndrome natural killers and expanded healthy donor natural killers. MDS-NKs: Myelodysplastic syndrome natural killers; HD-NKs: Healthy donor natural killers; GOBP: Short for gene ontology biological process; KEGG: Short for Kyoto encyclopedia of genes and genomes; NES: Normalized enrichment score; FDR: False discovery rate.

Subsequently, we were trying to verify the overall alterations of genetic spectrum eMDS-NKs, and observed the similarities and variations in the content of the major variable shear events, including as_a3ss, as_a5ss, as_mxe, as_ri and as_se (Figure 5D). Among the five subtypes of differentially variable shear events, as_se serves as the dominant subtype between eHD-NKs and eMDS-NKs (Figure 5E). Additionally, the Cicros diagrams intuitively displayed the loci regional distribution of the indicated somatic variations (single nucleotide polymorphisms, insertions and deletions, and gene fusion events), and no obvious trends were observed in the chromosomes of eHD-NKs and eMDS-NKs (Figure 5F). Collectively, these findings suggested the diverse variations in DEG-related bioprocesses and gene sets, together with the conservations in genetic spectrum between eHD-NKs and eMDS-NKs.

DISCUSSION

As mentioned above, MDS, including the lower- and higher-risk counterparts, are a spectrum of heterogenous clonal myeloid malignancies of HSCs and bone marrow failure disorders with limited treatment options including the transplant (e.g., HSCT) and non-transplant (e.g., hypomethylating agents) treatment[27-29]. Longitudinal studies have demonstrated the prominent contribution of the immunosuppressive tumor microenvironment (e.g., immune cells, inflammatory factors) and the resultant dysregulated immune response in inducing MDS blast tolerance and accelerating MDS progression, yet the status and key features of NK cells in the hematopoietic microenvironment are largely unknowable[30]. Therewith, in this study, we systematically compared the cellular signatures of resident and expanded MDS-NKs and HD-NKs from BME, and displayed the unique landscapes of transcriptomic variations between eMDS-NKs and eHD-NKs. Taken together, these findings uncovered the cytogenetic abnormalities of MDS-NKs in bone marrow niche and would contribute to the further innovation of NK cell-related pathogenesis and targeted therapy of MDS.

Currently, increasing evidence has pointed out the critical role of cellular and genetic alterations in hematopoietic stem precursor cells and immune cell-mediated immunodysfunction in BME[1,31]. For instance, diverse genes with cellular and molecular perturbations (e.g., GATA2 and/or TET2 mutations, DDX41 gene variants, RUNX1 deficiency) during malignant transformation are associated with comorbidities or complications during HSCT for hematopoietic reconstitution in MDS patients[19,32-35]. With the aid of single cell RNA-SEQ of entire bone marrow cells, Guo et al[36] profiled overall alteration as well as T-cell hyperactivation and abnormal myelogenesis during deficient hematopoiesis of five MDS patients. In recent years, we and Liu et al[37] reported the multifaceted cellular and molecular alterations of mesenchymal stem/stromal cells in the BME of hematological diseases including MDS, together with the regulatory effects upon the dysfunction of NK cells and T cells[17,38,39]. These data consistently put forward the pivotal role of genetic variations and cell number alteration in the pathogenesis of MDS and the accompanied immune response in the setting of dysimmunity. In this study, we further uncovered the characteristics of resident and expanded MDS-NK cells from the BME, which thus provided useful clues and candidate subjects for further exploring the impact of cytogenetic abnormalities and the resultant deterioration of inflammatory BME upon the clonal architecture of HSCs and the development of MDS[30].

In recent years, a variety of treatment options have been developed and considerable efficacy has been achieved for MDS administration, including supportive care and aggressive treatment (e.g., chemotherapy, HSCT)[40,41]. For example, azacytidine and decitabine are adequate to incorporate into DNA, and result in DNA methyltransferase depletion and the reactivation of key genes and signaling pathways with tumor-suppressive activities[42]. Similarly, with the aid of RNA-SEQ and bioinformatics analyses, we also observed a certain number of DEGs, signaling pathways and the genetic spectrum (e.g., somatic variations) involved in subcellular structure (e.g., plasma membrane, microtubule) and diverse bioprocesses (e.g., DNA impair, cell cycle) between MDS-NKs and HD-NKs. However, considering the substantial biological heterogeneity of MDS, this sample size of RNA-SEQ cohort in this study is limited and insufficient to support broad claims about an MDS-specific transcriptomic landscape. That is, the RNA-seq component including the DEGs in this study is exploratory, which needs further validation in the larger cohort using qPCR or another independent assay. Additionally, due to the limitations of rNKs in BM-MNCs, we turned to eNKs after a 14-day’s ex vivo expansion for RNA-SEQ and bioinformatics analyses. Nevertheless, the eMDS-NKs after “3ILs” stimulation might also reveal culture-induced characteristics rather than faithfully recapitulate the pathological conditions in bone marrow of MDS patients. Therefore, it would be much better for dissecting the similarities and alterations of transcriptomic profiling of freshly isolated resident NK cells (rHD-NKs, rMDS-NKs) in future. In consistence, a variety of these DEGs (e.g., PRF1, GZMB, IFNG, KLRK1) are associated with MDS-NK cell dysfunctions such as killing, proliferation, and survival. For instance, the increased release of IFNG usually triggers the enhanced apoptosis of myeloid precursor cells and the impaired function of immune cells including NK cells[43]. According to the genetic landscape, PRF1 and SEC23B revealed pathogenic germline variants in familial MDS[44]. These findings offered novel insights into the pathogenesis of myeloid malignancies and also supply a framework to prioritize genetic variants for routine diagnostics and MDS management. Nevertheless, the application of the computational and transcriptomic-validation would further help enhance the robust representation, reproducibility, and interpretation of high-dimensional biological data, which would support the broader methodological principle that complex data-driven findings require rigorous validation and careful assessment of generalizability.

Meanwhile, it’s noteworthy that the current management for MDS remains complex and intractable owing to the limitations in our understanding of pathogenesis as well as the deficiency of novel therapeutic approaches[45,46]. Beside the aforementioned NK cell-related pathogenesis, state-of-the-art literatures have also prompted the encouraging outcomes and robust prospects of NK cell-based cytotherapy and comprehensive therapy for patients with MDS or the transformed acute myeloid leukemia[10,47]. Additionally, considering the potential variations and highly heterogeneous features among MDS patients, NK-cell phenotype may differ according to subtype, Revised International Prognostic Scoring System/Molecular International Prognostic Scoring System risk, cytogenetics, mutation profile, treatment history, blast percentage, age, and disease stage. Therefore, it’s of great importance to further summarize these characteristics and explore associations with major NK-cell alterations in future. Collectively, our findings revealed the cellular and transcriptomic characteristics of NK cells in the BME of MDS, and would be conducive to the subsequent illumination of the underlying pathogenesis and the innovation of targeted therapy in future.

CONCLUSION

Overall, we conducted systematic comparisons of the cellular phenotypes and transcriptomic characteristics of both resident (rMDS-NKs, rHD-NKs) and expanded (eMDS-NKs, eHD-NKs) NK cells from the BME between patients with MDS and HDs. The findings in the similarities and the concomitant alterations of MDS-NKs both at the cellular and transcriptomic levels provide new references for identifying exploratory cellular and transcriptomic differences between HD-NKs and MDS-NKs. Collectively, this study provides preliminary biological evidence that warrants validation in larger, molecularly characterized MDS cohorts and functional mechanistic studies.

ACKNOWLEDGEMENTS

The co-authors would like to thank the members in The Fourth People’s Hospital of Jinan Affiliated to Shandong Second Medical University and Qingdao University. We also thank the support from Shandong Provincial Key Medical and Health Laboratory of Blood Ecology and Biointelligence, Jinan Key Laboratory of Medical Cell Bioengineering, and Shandong Health Youth Science and Technology Innovation Team for their assistance.