Prognostic value of natural killer T cell related genes in acute myeloid leukemia

Abstract

Background

Acute myeloid leukemia (AML) is a hematological malignancy characterized by complex immune microenvironment. This study aims to identify immune-related prognostic biomarkers in AML.

Methods

Multiple public sequencing datasets were utilized to analyze differentially expressed genes (DEGs) in AML. Single-sample gene set enrichment analysis (ssGSEA) and weighted gene co-expression network analysis (WGCNA) were also performed. Immune cell infiltration was assessed at the single-cell level. NKT cell marker genes were intersected with the most AML-relevant module genes to identify key genes. Prognostic genes were screened using the Cox Lasso regression model, and their prognostic value was evaluated with Cox random forest and Kaplan-Meier survival analyses. Gene expression was validated using RT-qPCR and Western blot, and immune cell levels were analyzed by flow cytometry.

Results

A total of 1,919 common DEGs were obtained between AML and controls. WGCNA revealed that the brown module was most strongly associated with AML. Single-cell analysis showed that NKT cell infiltration was significantly reduced in AML patients, consistent with ssGSEA results. Forty intersecting genes were identified between NKT cell marker genes and brown module genes. Cox Lasso regression identified 10 prognostic genes (FGFBP2, GZMB, GZMH, IKZF3, IL2RB, KLRB1, KLRC2, RHOF, RUNX3, and STAT4). A risk score model based on these genes stratified AML patients into high-risk and low-risk groups, with significant differences in survival prognosis between the two groups. RT-qPCR and Western blot analyses showed that these genes were significantly downregulated in AML patients. Flow cytometry results revealed significantly lower levels of NKT and CD8 + T cells in AML patients compared to controls.

Conclusion

This study identified key prognostic genes in AML and highlighted the critical role of NKT cells in AML pathogenesis. The study provides new insights and potential biomarkers for understanding AML biology, prognosis, and therapeutic targets.

Introduction

Acute myeloid leukemia (AML) is an aggressive and highly heterogeneous hematologic malignancy characterized by the abnormal proliferation of myeloid precursor cells in the bone marrow, which suppresses normal hematopoiesis [1]. AML is the most common form of acute leukemia in adults. Despite significant improvements in survival rates for some patients due to targeted therapies and hematopoietic stem cell transplantation, the overall 5-year survival rate remains below 30% [2, 3]. The effectiveness of AML treatment is influenced by various factors, including the genetic heterogeneity of the tumor, complex molecular mechanisms, and significant alterations in the immune microenvironment [4]. Immune evasion and a suppressive tumor microenvironment are considered major obstacles to AML treatment, although their specific molecular mechanisms are not yet fully understood [5].

In AML patients, the function of effector immune cells is significantly impaired, creating an immunosuppressive state that facilitates the survival and proliferation of AML cells [6, 7]. Moreover, the high degree of individual heterogeneity among AML patients poses challenges for prognosis. Traditional prognostic models based on clinical features and single biomarkers fail to adequately capture the molecular and immunological complexities of AML, limiting their ability to guide precision treatment effectively [8]. With the widespread application of machine learning in biomedical research, integrating big data analysis and multi-omics data has made it possible to construct precise prognostic models for AML [9], offering new avenues for identifying prognostic factors and therapeutic targets.

Based on this, the present study aims to systematically analyze the molecular mechanisms, immune microenvironment changes, and prognostic features of AML patients using multi-level data integration and experimental validation approaches. First, by integrating public transcriptome datasets, we identified differentially expressed genes (DEGs) in AML and used weighted gene co-expression network analysis (WGCNA) to construct key gene modules associated with AML. Second, we explored changes in the immune microenvironment of AML patients through single-sample gene set enrichment analysis (ssGSEA) and single-cell RNA sequencing. Third, we applied Cox Lasso regression and Cox random forest algorithms to screen key prognostic genes and construct a risk score model to evaluate their clinical significance. Finally, RT-qPCR, Western blot, and flow cytometry were used to experimentally validate changes in key genes and immune cell populations. By combining multi-level and multi-technical approaches, this study aims to elucidate the core pathological mechanisms of AML and lay the groundwork for future individualized therapies and biomarker research.

Materials and methods

Data collection and differential expression analysis

The GSE114868 dataset [10] includes gene expression profiles of bone marrow mononuclear cells from 194 AML patients and 20 healthy controls. The GSE37642 dataset [11] contains gene expression data and clinical information for bone marrow mononuclear cells from 136 AML patients. Gene expression data from these datasets were normalized using the oligo and affy packages. The TCGA dataset includes gene expression data, clinical information, mutation information and for 179 AML patients, while the gene expression profiles of 337 normal whole blood samples were obtained from the GTEx cohort [12]. Additionally, single-cell RNA sequencing (scRNA-seq) data were derived from the GSE235857 dataset [13], which contains peripheral blood mononuclear cell (PBMC) data from 6 AML patients and 6 controls.

Differentially expressed genes (DEGs) between AML samples and normal controls were identified using the limma package to compare gene expression differences. The following criteria were applied to define DEGs:|log2FC| > 1 and p < 0.05. To improve the efficiency of DEG identification, DEG analyses were performed separately for different datasets (TCGA, GTEx, and GSE114868), and overlapping genes across these datasets were determined using Venn diagram analysis. The overlapping genes were subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses using the clusterProfiler package.

Weighted gene Co expression network analysis

To uncover the potential molecular mechanisms and immune-related modules in AML, weighted gene co-expression network analysis (WGCNA) was performed to construct a gene co-expression network. Analysis was conducted using the WGCNA package in R, with an appropriate soft threshold (power) selected to ensure the scale-free topology of the network. Genes were grouped into multiple modules based on their expression patterns. The correlation between each module and AML clinical traits was calculated to identify the module most strongly associated with AML.

Single sample gene set enrichment analysis

To further investigate the immune microenvironment in AML, single-sample gene set enrichment analysis (ssGSEA) was employed to evaluate the infiltration levels of various immune cells. Known immune cell marker gene sets from the ImmPort database were used to represent different immune cell types. ssGSEA calculates the enrichment score (ES) of specific immune cell-related gene sets in each sample to estimate the relative abundance of different immune cell populations.

Analysis of single-cell RNA sequencing data

Immune cell infiltration in AML patients was assessed at the single-cell level using the GSE235857 dataset. The Seurat package was employed for quality control of single-cell data, retaining cells with 500 to 6,000 detected genes and a mitochondrial gene expression proportion of less than 20%. Log-normalization was applied to normalize the gene expression matrix of each cell. Principal component analysis (PCA) was performed to reduce dimensionality, and the top 20 principal components were selected for downstream analysis. Further dimensionality reduction was conducted using t-SNE, and clustering analysis was performed using the FindClusters function in Seurat. Immune cell subsets were annotated by analyzing the gene expression characteristics of each cluster in combination with known immune cell marker genes. The proportion of each immune cell population in AML patients and healthy controls was calculated to evaluate the infiltration patterns of various immune cells in AML.

To identify intersection genes, NKT marker genes were compared with genes from the AML-associated module identified by WGCNA. Genes expressed in both datasets were selected as intersection genes.

Cox Lasso regression model and prognostic gene screening

Cox Lasso regression was used to identify prognostic genes in AML. Lasso regression applies L1 regularization to select genes associated with prognosis, enabling the construction of a risk score model based on prognostic genes. The model was initially trained using differentially expressed candidate genes and subsequently validated using the Cox proportional hazards regression model. Further validation was performed using Kaplan-Meier survival analysis, stratifying AML patients into high-risk and low-risk groups to evaluate survival differences. The diagnostic value of the gene-based risk score for AML prognosis was assessed using ROC curve analysis. The prognostic impact of the identified genes on AML patients was further evaluated using Cox random forest analysis and Kaplan-Meier (K-M) survival analysis.

Clinical sample collection

Peripheral blood PBMC samples were randomly collected from 10 AML patients and 10 age- and sex-matched healthy individuals as controls. All sample collection and processing complied with the Declaration of Helsinki and ethical guidelines and was approved by the Ethics Committee of Rizhao People’s Hospital (IEC-Form-016-6.1). Written informed consent was obtained from all participants.

Inclusion Criteria: Diagnosis consistent with the World Health Organization (WHO) criteria for acute myeloid leukemia (2016 revision); Cytological confirmation of AML through bone marrow aspiration; Age ≥ 18 years, with no restrictions on race or sex; No history of other malignancies or significant organ dysfunction; No prior chemotherapy or other AML treatments before sample collection.

Exclusion Criteria: Presence of major infections or other immune system diseases; Pregnant or breastfeeding women.

RT-qPCR

Total RNA was extracted from PBMCs using TRIzol reagent (Invitrogen), and RNA concentration was measured with a NanoDrop spectrophotometer. For complementary DNA (cDNA) synthesis, 1 µg of total RNA from each sample was used with the PrimeScript™ RT Reagent Kit (TaKaRa) according to the manufacturer’s instructions. Specific primers for target genes and the housekeeping gene GAPDH were designed using the Primer-BLAST tool (Table S1). PCR reactions were performed using TB Green™ Premix Ex Taq™ II (TaKaRa) on an ABI StepOnePlus™ Real-Time PCR system. Ct values for all samples and genes were recorded, and relative gene expression levels were calculated using the 2^−ΔΔCt method. Each sample was analyzed in at least three technical replicates.

Western blot

Proteins were extracted from PBMCs using RIPA lysis buffer supplemented with protease and phosphatase inhibitors. Protein concentration was measured using a BCA protein assay kit (Beyotime, China). A total of 30 µg of protein was loaded onto an SDS-PAGE gel for electrophoresis and subsequently transferred onto a PVDF membrane. The membrane was blocked with 5% skim milk for 2 h and then incubated overnight at 4 °C with specific primary antibodies (Abcam, CA, UK). The membrane was further incubated with HRP-conjugated secondary antibodies (Abcam) for 2 h. ECL detection reagents were applied evenly to the membrane, and chemiluminescent signals were captured using a ChemiDoc imaging system. β-actin was used as the internal reference protein, and the grayscale values of target protein bands were measured using ImageJ software (NIH, MD, USA). Each sample was analyzed in at least three technical replicates.

Flow cytometry analysis

To further validate the changes in immune cell levels in AML patients, flow cytometry was used to analyze immune cells in peripheral blood. Fluorescently labeled antibodies included CD45 KO, CD3 APC, and CD8 FITC for CD8 + T cells, and CD45 KO, CD3 APC, and CD56 PerCP for NKT cells. A total of 1 × 10⁶ PBMCs were resuspended in 100 µl PBS in flow tubes, followed by the addition of antibody mixtures and incubation at 4 °C in the dark for 30 min. Cells were washed with 1 ml PBS, centrifuged (300 g, 5 min), and resuspended in PBS. Immune cell frequencies and numbers were measured using a BD FACSCanto™ II flow cytometer. All flow cytometry results were analyzed using FlowJo software (Tree Star) and compared with normal controls. Each sample was analyzed in at least three technical replicates.

Statistical analysis

All statistical analyses were performed using R (v4.1.0) and GraphPad Prism (v9.3.0, GraphPad Software, USA). Student’s t-test or one-way ANOVA was used to compare results between groups. A p-value less than 0.05 was considered statistically significant.

Results

Identification of differential genes

In the transcriptomic data of AML and normal control groups, 10,401 differentially expressed genes (DEGs) were identified from the TCGA and GTEx datasets (Fig. 1A), and 2,870 DEGs were identified from the GSE114868 dataset (Fig. 1B). Intersection analysis identified 1,919 common DEGs shared across these datasets (Fig. 1C and D). Gene Ontology (GO) enrichment analysis revealed that these common DEGs were significantly enriched in processes such as leukocyte activation and immune response (Fig. 2A). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis showed that these DEGs were significantly involved in pathways such as cytokine-cytokine receptor interaction, chemokine signaling pathway, and NOD-like receptor signaling pathway (Fig. 2B).

Fig. 1

Differential analysis between AML and controls. (A) Volcano plot of differentially expressed genes in TCGA and GTEx. (B) Volcano plot of differentially expressed genes in GSE114868. (C) The intersection of two sets of differentially expressed genes. (D) Heatmap of expression of common DEGs

Fig. 2

Enrichment analysis of common DEGs. (A) The main GO enrichment terms enriched by common DEGs. (B) The main KEGG enrichment terms enriched by common DEGs

Co-expression network and immune cell infiltration

Using the common DEGs, a co-expression network was constructed through WGCNA. By selecting a soft-thresholding power of β = 14 (R² = 0.85), four gene modules were identified (Fig. 3A and B), each consisting of genes with similar expression patterns. Further analysis of module-clinical trait correlations revealed that the brown module had the strongest association with AML (Fig. 3C). The ssGSEA was used to assess differences in immune cell infiltration levels between AML patients and healthy controls. The results showed that multiple immune cells had significantly lower enrichment scores in AML patients compared to healthy controls (Fig. 4A and B).

Fig. 3

Weighted gene co-expression network analysis (WGCNA) reveals key modules associated with AML. (A) Soft-thresholding powers analysis. (B) Cluster diagram of gene modules. (C) Module-trait relationship heatmap showing the correlation between gene modules and AML clinical traits

Fig. 4

Single-sample gene set enrichment analysis (ssGSEA) reveals altered immune infiltration in AML. Comparison of immune cell infiltration scores between AML patients and healthy controls in TCGA (A) and GSE114868 (B)

Screening important immune cells and genes in AML

Single-cell RNA sequencing data were used to analyze the immune microenvironment of AML patients and healthy controls, revealing significant alterations in immune cell composition in AML patients. Through t-SNE dimensionality reduction and cell clustering analysis, 28 clusters were identified (Fig. 5A). Based on the expression levels of cell markers, 15 immune cell types were annotated, including MAIT, NK, CD8 + Trm, CD8 + Tem, Bn, CD8 + Teff, CD4 + Tn, CD8 + Tn, natural killer T cells (NKT), platelet, plasma, CD4 + Tcm, DC, FCGR3A + monocyte, and CD14 + monocyte (Fig. 5B and C).

Fig. 5

Single-cell RNA sequencing analysis reveals immune microenvironment changes in AML. (A) t-SNE plot of cell clusters. (B) Violin plots of expression of marker genes in cell clusters. (C) t-SNE plot of immune cells. (D) t-SNE plot of immune cells in AML and control groups. (E) Percentage of immune cells in AML and control groups. (F) Differences in immune cells between AML and control

By comparing immune cell expression levels between AML patients and controls, we found that CD4 + Tcm and CD8 + Trm were significantly increased in AML, whereas NKT, MAIT, NK, and CD8 + Tn cells were significantly reduced (Fig. 5D, E and F). Further analysis of intercellular communication among different cell subsets revealed enhanced interactions between CD8 + Tn and CD8 + Teff with NKT cells in AML patients (Figure S1). Notably, the significant reduction of NKT cells in AML is consistent with the ssGSEA results.

Prognostic gene screening and risk scoring model construction

NKT marker genes were intersected with the genes from the brown module identified in the WGCNA analysis, resulting in 40 overlapping genes (Fig. 6A). Using the Cox Lasso regression model, 10 genes closely associated with AML prognosis were selected with the optimal λ value determined at 0.04 during cross-validation: FGFBP2, GZMB, GZMH, IKZF3, IL2RB, KLRB1, KLRC2, RHOF, RUNX3, and STAT4 (Fig. 6B and C).

Fig. 6

Cox Lasso regression identifies prognostic genes. (A) The intersection of NKT marker genes and brown module genes in WGCNA analysis. (B) Cross-validation for the selection of the optimal λ value. (C) Gene coefficient profiles as λ increases in Lasso regression

Based on the expression levels of these 10 prognostic genes, risk score models were constructed in both the TCGA (Fig. 7A, B and C) and GSE37642 (Fig. 7D, E and F) datasets. AML patients were stratified into high-risk and low-risk groups. Kaplan-Meier survival analysis showed that the overall survival of the high-risk group was significantly shorter than that of the low-risk group. Additionally, time-dependent ROC curve analysis demonstrated that the risk score model achieved high predictive accuracy for survival at 1 year, 3 years, and 5 years. These findings indicate that the model is an effective prognostic tool for predicting AML patient survival. To further explore the relationship between risk score models and molecular genetic characteristics of AML patients, we analyzed the distribution of gene mutations in AML patients in different risk groups (high-risk group and low-risk group) in the TCGA database (Fig. 7G). The results showed that the mutation frequency of genes such as DNMT3A, and NPM1 was higher in the high-risk group than in the low-risk group. This indicates that the risk scoring model constructed in this study can further differentiate patient subgroups at the molecular genetic level.

Fig. 7

Construction and validation of risk score model for AML patients. (A) Distribution of risk scores and survival status of AML patients in TCGA. (B) Kaplan-Meier survival curves of AML patients stratified into high-risk and low-risk groups based on the risk score in TCGA. (C) Time-dependent ROC curves for the risk score model at 1 year, 3 years, and 5 years in TCGA. (D) Distribution of risk scores and survival status of AML patients in GSE37642. (E) Kaplan-Meier survival curves of AML patients stratified into high-risk and low-risk groups based on the risk score in GSE37642. (F) Time-dependent ROC curves for the risk score model at 1 year, 3 years, and 5 years in GSE37642. (G) The distribution of gene mutations between high-risk and low-risk groups

The contribution of the prognostic genes to patient survival was further evaluated using the Cox random forest algorithm. The results showed that GZMB and FGFBP2 significantly influenced patient prognosis in both the TCGA (Fig. 8A) and GSE37642 (Fig. 8B) datasets. Kaplan-Meier survival analysis revealed that FGFBP2, GZMB, GZMH, IL2RB, KLRB1, KLRC2, RHOF, RUNX3, and STAT4 significantly affected prognosis in the TCGA dataset (Figure S2A), while FGFBP2, GZMB, GZMH, IKZF3, IL2RB, and KLRB1 significantly impacted prognosis in the GSE37642 dataset (Figure S2B).

Fig. 8

Evaluation of prognostic genes. Cox random forest of prognostic genes in TCGA (A) and GSE37642 (B)

Experimental verification

RT-qPCR and Western blot analyses were performed to validate the expression levels of prognostic genes in AML patients and healthy controls. RT-qPCR results showed that the mRNA levels of FGFBP2, GZMB, GZMH, IKZF3, IL2RB, KLRB1, KLRC2, RHOF, RUNX3, and STAT4 were significantly lower in AML patient samples compared to controls (Fig. 9A). Western blot analysis demonstrated that protein expression levels were consistent with RT-qPCR results, with significantly reduced expression of the prognostic genes in AML patients compared to healthy controls (Fig. 9B).

Fig. 9

Detection of prognostic genes. (A) The mRNA levels of marker genes in AML and controls detected by qRT-PCR. (B) Protein levels of marker genes in AML and controls detected by western blot. *P < 0.05, **P < 0.01, ***P < 0.001

Flow cytometry further confirmed changes in immune cell levels, revealing that the levels of NKT cells and CD8 + T cells were significantly lower in AML patients compared to controls (Fig. 10).

Fig. 10

The level of immune cells in AML and controls detected by flow cytometry

Discussion

This study systematically analyzed the molecular mechanisms, immune microenvironment changes, and prognostic factors of AML. By integrating bioinformatics, machine learning, and experimental validation, the study uncovered the roles of key genes and immune cells in AML, providing new insights for the diagnosis, prognosis evaluation, and treatment of AML.

DEGs identified in the transcriptomic data of AML patients and healthy controls were significantly enriched in immune-related pathways, including leukocyte activation, cytokine-cytokine receptor interaction, chemokine signaling pathway, and NOD-like receptor signaling pathway. Suppression of the immune microenvironment enables malignant cells to evade immune surveillance, and numerous new therapeutic approaches targeting immune evasion aim to restore anti-leukemic immune activity [14]. Altered cytokine signaling in the immune microenvironment contributes to disease progression and treatment resistance [15]. Cytokines and chemokines within the bone marrow microenvironment play critical roles in supporting AML cell survival, promoting resistance to conventional chemotherapy and targeted therapies, and ultimately leading to disease relapse [16]. Activation of the NOD-like receptor signaling pathway induces chronic inflammation, which promotes AML development and progression [17]. These functional abnormalities suggest that the immune system in AML patients may be impaired, further facilitating tumor cell growth and proliferation.

At the immune microenvironment level, ssGSEA analysis and single-cell RNA sequencing revealed significantly reduced infiltration of several key immune cell subsets in AML patients, particularly NKT cells and CD8 + T cells. This finding indicates that immune evasion in AML may be closely associated with the depletion of effector immune cells. Consistent with our results, the levels of CD8 + effector subsets and NKT cells in AML patients were significantly lower than in healthy controls [5]. Suppression of T-cell proliferation and NKT cell numbers in AML has been linked to an immunosuppressive microenvironment and worse prognosis [18]. Decreased intracellular signaling and reduced cytotoxicity in NK and NKT cells in AML patients are associated with an increased risk of infection and mortality [19]. These findings demonstrate that the immune microenvironment in AML is significantly dysregulated, with immune evasion emerging as a critical mechanism driving AML progression.

Furthermore, WGCNA identified the brown module as the most strongly associated with AML, suggesting its potential central role in AML pathogenesis. Although hub genes identified by WGCNA are valuable, we specifically chose genes from the intersection with NKT cell markers and further refined these through Cox Lasso regression to prioritize genes closely linked to AML prognosis and immune dysregulation. By intersecting NKT cell marker genes with genes in the brown module, 40 overlapping genes were identified. Cox Lasso regression further selected 10 genes closely associated with AML prognosis. Among these, FGFBP2 [20] and KLRC2 [21] were found to be significantly downregulated in AML. FGFBP2, GZMB, and GZMH are cytotoxic factors highly expressed in NKT cells and play essential roles in NK cell-mediated immune functions [22, 23]. High expression of GZMB has been linked to poorer clinical outcomes in AML [24]. IL2RB is associated with immune evasion and suppression [25, 26]. Increasing evidence suggests that KLRB1 plays a critical role in tumor immunity, with its high expression in most cancers correlating with favorable prognosis [27]. In addition, the in vitro expansion of NKG2C + NK cells had been proposed as a simple strategy to enhance NK cell antitumor cytotoxicity for immunotherapy [28]. Contrary to our findings, a study by Wen et al. [29] reported that RHOF is overexpressed in AML and associated with poor prognosis. Similarly, RUNX3 is highly expressed in AML cells and linked to adverse outcomes in AML patients [30]. STAT4, which is upregulated in AML, has also been associated with unfavorable prognosis [31].

These genes are not only closely associated with critical biological processes such as immune regulation and cytotoxicity but also exhibit significant prognostic value. The constructed risk score model effectively stratified AML patients into high-risk and low-risk groups and demonstrated high accuracy in predicting survival at 1, 3, and 5 years. Additionally, the Cox random forest model further confirmed the critical prognostic contributions of genes such as FGFBP2 and GZMB, suggesting that these genes may serve not only as prognostic markers but also as potential therapeutic targets.

Despite the significant progress made in this study, several limitations remain. First, the heterogeneity of public datasets may influence the analytical results, necessitating validation with more multicenter datasets. Second, although the expression patterns of key genes were experimentally validated in clinical samples, further in vitro functional experiments, such as gene knockdown or overexpression assays, cell proliferation, apoptosis assays, and functional immune assays, are required to elucidate the precise biological roles and mechanisms of these genes in AML pathogenesis. Such experiments will be conducted in our future studies. Additionally, due to the limited sample size and the relatively narrow scope of clinical information collected in this study, we were unable to systematically correlate gene expression and immune cell distribution with detailed clinical characteristics. Future studies will incorporate larger patient cohorts and comprehensive clinical data to further explore these correlations and validate the clinical applicability of our findings. Finally, while the risk score model demonstrated high predictive performance, its feasibility and robustness in clinical applications need to be validated in larger cohorts of AML patients, particularly including patients undergoing targeted therapies. Further prospective studies focusing specifically on targeted therapy-treated AML populations are warranted to confirm whether this prognostic score maintains validity and practical utility in guiding treatment decisions in this subgroup.

Conclusion

This study systematically revealed key changes in the immune microenvironment of AML and constructed an efficient prognostic model based on NKT cell-associated genes. The inactivation of NKT cells and their related genes may represent a critical mechanism of immune evasion in AML. Future immunotherapeutic strategies aimed at restoring NKT cell function hold promise as novel treatment approaches for AML patients. This research not only provides important theoretical support for the study of AML pathogenesis but also offers new targets and tools for the development of diagnostic and therapeutic strategies in the context of precision medicine for AML.