The clinical samples were bone marrow samples collected from patients with AML and normal donors in the Department of Hematology, The Affiliated Hospital of Guizhou Medical University from 2021 to 2023. The patient’s condition was diagnosed using morphological, cytochemical, and immunotyping. Patients with AML were classified as “newly diagnosed” and “relapse” by diagnosis. Therefore, the clinical samples were divided into three groups: “normal donors” (n = 27), “newly diagnosed” (n = 33) and “relapse” (n = 29). Detailed data are shown in Supplementary Tables S1–S3. Written informed consent was obtained from the individuals for the publication of any potentially identifiable images or data included in this article. The clinical samples in this paper were approved by the Ethics Committee of the Affiliated Hospital of Guizhou Medical University for basic research, and the approval number is 2021 Ethics Approval No. 182.

AML clinical samples were retrieved from TARGET (https://ocg.cancer.gov/programs/target) database as the training set, while data from the TCGA (https://www.cancer.gov/ccg/research/genome-sequencing/tcga) database were utilized as the validation set.

The sangerbox (http://sangerbox.com/) and the packages in R language mentioned below were employed for data visualization.

The “Consensus Cluster Plus” package in R language was applied for unsupervised consensus cluster analysis. Hypoxia signature genes were obtained from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb/) database. We selected ‘REACTOME_CELLULAR_RESPONSE_TO_HYPOXIA.v2022.1.Hs.gmt’ as our target gene set for the study. Hypoxia-related genes associated with AML clinical samples were screened using VLOOKUP function. Cluster numbers were between k = 2-6, and the relative area under the cumulative distribution function (CDF) curve were used to evaluate the clustering stability. Figures 1A–C was plotted using SangerBox.

The differential expression of mRNAs was evaluated using the “Limma” package in R language, with thresholds of p < 0.05 and log2|fold change| >0.585. The expression levels and distributions of DEGs between C1 and C2 were analyzed using the “Pheatmap” package in the R language. The sanger box was utilized to generate Figures 2A, B.

The survival analysis was conducted using the “survival” package in the R language, and the OS of patients belonging to different clusters (C1 and C2) were analyzed and evaluated. Furthermore, a comparison was made between the 1-,3-and 5-year OS of high- and low-risk groups in both the training and validation sets. The Sangerbox was utilized for the visualization of Kaplan-Meier (KM) survival curves in Figures 1D, 3C, 4A.

GO and KEGG were conducted using the “Cluster Profilter” package in R language, false discovery rate (FDR) < 0.05. Figure 2C was generated using sangerbox.

The DEGs was subjected to pathway enrichment analysis using GSEA: (https:// www.broadinstitute.org/gsea/). Normalize enrichment score (NES): The normalized enrichment score after correction was normalized by the data of the gene set; NOM p-val: The p-value obtained by statistical analysis of ES value represents the reliability of the result; FDR q-val: The p-value after multiple hypothesis testing correction represents the probability of false positive results, The, the smaller the p-value, the more significant. |NES|>1, FDR <0.25, p < 0.05 was considered statistically significant. And “Cluster Profilter” in R language was used to draw Figures 2D, E.

Univariate Cox regression analysis was adopted to obtain 34 hypoxia-related genes that exhibited significant associations with OS in AML patients (hazard ratio, HR = 95%, p < 0.05). The “ggforest” package in R language was utilized to construct Figure 3A. LASSO Cox regression analysis was performed on 34 hypoxia-related genes to eliminate any false positive hypoxia-related genes that may be associated with prognosis. The “glmnet” package in R language was utilized to generate Figures 3B, 2 hypoxia-related genes were selected according to the minimum λ value for constructing the risk prognostic model.

The risk scoring formula is as follows: Risk score = ∑ i = 1 n coe f i ∗ x i

The term “Coefi” denotes the coefficient, while “Xi” represents the normalized count of each core gene. Receiver Operating characteristic (ROC) curve was generated using the R language package “time ROC”. The accuracy of the prognostic model in predicting the 1 -, 3 - and 5-year OS of AML patients was assessed by calculating the Area Under Curve (AUC) in both the training and validation datasets. Figures 3D, E was plotted using SangerBox.

The CIBERSORT (https://cibersortx.stanford.edu/) in conjunction with the LM22 feature matrix was applied to analyze the differences in immune infiltration of 22 immune cells among different groups. The Pearson product-moment correlation coefficient was utilized to compute the correlation among immune cells, while the Mantel test was employed to statistically analyze the correlation between the risk score matrix and the immune cell matrix. r = 0-1 represents correlation, with higher values indicating stronger correlation, and p < 0.05 was considered statistically significant. The “CIBERSORT” package was used to conduct Figures 5A–C, while the “ggcor” package in R language was used to generate Figure 5D.

The model was trained and drug sensitivity of the samples was predicted using ridge regression, based on the relationship between gene expression and drug IC₅₀ in the training set. GDSC datebase preformed to analyze drug sensitivity of cancer cells. The sanger box was utilized to draw Figure 6A. The “oncoPredict” packages in R language were used to construct Figure 6B.

5 mL of bone marrow from AML patients or normal donors was collected by bone marrow puncture under routine sterile conditions and preserved with EDTA anticoagulation. Bone marrow was diluted 1:1 in equal volume with saline and was slowly added along the wall to a centrifuge tube pre-loaded with Ficoll (Solarbio Technologies, Beijing, China) separation solution, and Ficoll was 1:1 with diluted bone marrow. After centrifugation at 2,000 rpm for 15 min at room temperature, the intermediate white cell layer was aspirated and transferred to a new centrifuge tube. After centrifugation at 1,500 rpm for 5 min, the supernatant was discarded. The remaining precipitate, namely, bone marrow mononuclear cells, was retained after being washed three times with saline and subsequent discarding of the supernatant.

The authenticity of THP-1 and U937 human leukemia cell lines was confirmed by STR analysis, and they were cultured in a 5% CO₂ incubator at 37°Cusing RPMI1640 medium containing 10% fetal bovine serum. Drug-resistant variants, namely, THP-1R and U937R, were generated by supplementing 1% penicillin (100 units/mL) and streptomycin (100 mg/mL) to the medium along with increasing concentrations of cytarabine (Ara-C). The drug concentration was gradually escalated, repeating this process three to five times at each concentration after the cells have proliferated to a normal shape. The drug induction was maintained for a duration of 6–8 months until the cells achieved a stable state at the final concentration.

The extraction of total RNAs from cells was performed using Trizol reagent (Invitrogen, Carlsbad, CA, United States). The mixture was vigorously shaken for 15 s after the addition of chloroform, followed by incubation at room temperature for 3 min. The samples were centrifuged at 12,000 rpm for 15 min at 4°C, and the supernatant was retained. An equal volume of isopropanol was added, followed by centrifugation at 12,000 rpm for 10 min at 4°C. Subsequently, the supernatant was discarded. The RNA precipitate was washed with 500 µL of 75% ethanol. Subsequently, the samples were subjected to centrifugation at a speed of 7,500 rpm for a duration of 3 min at a temperature of 4°C, followed by removal of the supernatant. The samples were air-dried for 5 min at room temperature to allow ethanol evaporation, followed by addition of DEPC water and measurement of concentration. cDNA was extracted using a reverse transcription kit (MedChemExpressMCE, United States of America). Real-time PCR was performed using SYBR Green PCR Master Mix (MedChemExpressMCE, United States) kit and PRISM 7500 real-time PCR Detection System (Thermo Fisher Scientific, United States). Relative expression of the target genes was calculated using β-actin as the reference through comparative cycle threshold (CT) values (2^−ΔΔCT). The following human primers were used in this paper:

β-actin F 5′-CTA​CCT​CAT​GAA​GAT​CCT​CAC​CGA-3´;

β-actin R 5′-TTC​TCC​TTA​ATG​TCA​CGC​ACG​ATT-3´;

PSMD11 F 5′-AGT​TCC​AGA​GAG​CCC​AGT​CT-3´;

PSMD11 R 5′-TTG​CAC​TGC​CTC​TTC​ATC​GT-3´;

PSMD14 F 5′- GTC​AGT​GTG​GAG​GCA​GTT​GAT​C-3′;

PSMD14 R 5′-CCA​CAC​CAG​AAA​GCC​AAC​AAC​C-3′.

The primary antibodies against PSMD11 and PSMD14 (Affinity Biosciences, United States of America) were diluted at 1:500 and 1:1,000, respectively, the β-actin primary antibody (Wuhan Sanying, China) was diluted at 1:3,000. Protein lysates were extracted from cells by adding 1 mM PMSF to the RIPA lysis buffer (Solarbio Science and Technology). The mixture was vigorously shaken and incubated on ice for 30 min. Subsequently, the supernatant was obtained by centrifugation at 12,000 rpm for 15 min at 4°C. The concentration of protein was determined using the BCA Protein Assay kit (Pierce, Hercules, CA, United States). The proteins were mixed with Loading buffer 1:4 and boiled at 100°C for 10 min 40 μg of proteins were then added to a 10% SDS-PAGE gel and electrophoresed into the separation gel at a constant voltage of 80 V followed by switching to a stable voltage of 120 V. At the end of electrophoresis, the separated proteins were transferred onto PVDF membranes and rotated at 250 mA for 1 h. After shaking with PBS containing 5% skim milk on a shaker for 2 h at room temperature, the membranes were washed. The primary antibodies were then incubated for more than 8 h at 4°C. After washing the membranes, secondary antibodies were incubated for at room temperature for 45 min. All protein bands were visualized using the Enhanced Chemistry kit (7Sea Biotech, Shanghai, China). β-actin was used as the internal reference.

CCK8 assay was used to detect the sensitivity of leukemia cell lines to Ara-C. The cells were inoculated into individual wells of a 96-well plate at a seeding density of 3 ×10⁴ cells/100 μL, with five replicates per experimental group. After subjecting the cells to various concentrations of Ara-c for a duration of 24 h, a volume of 10 μL CCK8 reagent was added into each well, the concentrations of Ara-c in U937 and U937R cell lines were 4, 16, 64, 192, 386, 578, 768 and 1,536 μM, meanwhile the concentrations of Ara-c inTHP-1 and THP-1R cell lines were 0.5, 4, 64, 192, 386, 578, 768 and1536 μM. The absorbance at 450 nm was quantified using a microplate spectrophotometer after co-culthring for 1–2 h. The IC₅₀ value was determined by employing the GraphPad Prism 9.5 software.

The statistical software GraphPad Prism 9.5 was utilized for conducting both the analysis and visualization of data. The Shapiro-Wilk test and the Kolmogorov-Smirnov test were used to test the normal distribution of data. After the data passed the normal distribution test, the unpaired t-test was employed to compare and evaluate the differences between two groups. The experimental data were represented as mean ± standaed deviation (SD). The significance level, denoted by p calue, is interpreted as follows: *p < 0.05, **p < 0.01, ***p < 0.001 ****p < 0.0001; among these values, p < 0.05 is considered statistically significant.

To investigate the association between hypoxia and poor prognosis in AML, AML patient samples were grouped into differentsubtypes. Firstly, a search was conducted in the MSigDB database using the keyword “hypoxia”, resulting in 60 gene sets. After thorough consideration of factors such as the number of genes within each set, correlation with the research field, we ultimately identified a gene set containing 75 hypoxia-related genes as the target gene set. After screening using the VLOOKUP function, 63 of these genes were found to be linked to AML. Subsequently, unsupervised consensus cluster analysis was performed on 187 AML samples downloaded from the TARGET database, which were based on these 63 genes. The clustering stability is assessed by employing the CDF curve and Delta area curve for different cluster numbers (K = 2–6), the selection of K = 2 for clustering was presumed to be the optimal choice (Figure 1A). Accordingly, AML patient samples were classified into two subtypes, namely, C1 (n = 121) and C2 (n = 66) (Figure 1B). The heatmap revealed that the hypoxia-related genes in C1 displayed higher expression compared to those in C2 by analyzing the expression distribution of 63 hypoxia-related genes (Figure 1C). Therefore, C1 was defined as the high-expression hypoxia-related gene group, whereas C2 as the low-expression group. To explore the effects of hypoxia-related genes on the clinical prognosis of AML, survival analysis was performed in the two groups. The KM survival curve indicated a worse OS in the C1 group. (p < 0.05) (Figure 1D).

We screened out the DEGs between the two subtypes, GO, KEGG and GSEA were performed to analyse their biological functions and signalling pathways associated with hypoxia. This analysis allowed us to elucidate the mechanism by which hypoxia affects the clinical prognosis of AML. The “Limma” package in R language was utilized for filtering the DEGs between the two subtypes (|log2 FC|>0.585, p < 0.05). The volcano plot revealed 7,735 DEGs between C1 and C2, with 7,505 genes upregulated and 230 genes downregulated (Figure 2A). The heatmap, further revealed the expression level and distribution of DEGs between C1 and C2, (Figure 2B). GO and KEGG enrichment analyses of DEGs identified the top 15 relevant biological processes (BP) and signaling pathways (Figure 2C). Notably, the top three enriched of GO terms were the cellular protein metabolic process, cellular protein modification process and protein modification process. Parallelly, KEGG-enriched signaling pathways included human T-cell leukemia virus type 1 infection, the T cell receptor signaling pathway, PD-L1 expression and the PD-1 checkpoint pathway in cancer, TH17 cell differentiation and other immune-related signaling pathways.

To compensate for any omissions in the DEGs-based GO and KEGG analysis, especially owing to the threshold setting, and further elucidate the related signaling pathways that may be activated by hypoxia, GSEA was used to compare the gene sets (Figure 2D). GSEA revealed positive correlations between the high expression of hypoxia-related genes and the activation of several signaling pathways, including the cytoplasmic pattern recognition receptor signaling pathway (NES = 2.24), pathways in cancer (NES = 1.55), protein acetylation (NES = 2.21) and mTOR signaling pathways (NES = 1.68). The ridge plot presented the posterior distribution of GSEA-GO and GSEA-KEGG (p < 0.05) (Figure 2E), confirming significant enrichment in the cytoplasmic pattern recognition receptor signaling pathway and protein acetylation pathways, which are consistent with the previously enriched biological processes, such as protein modification processes.

In further investigating the value of hypoxia-related genes in assessing the clinical prognosis of AML, a risk prognostic model of AML was constructed based on hypoxia-related genes. Initially, a univariate Cox regression analysis was performed to select 34 genes that were directly associated with the prognosis of AML from the 63 hypoxia-related genes (HR = 95%, p < 0.05) (Figure 3A). These 34 hypoxia-related genes were further analyzed using LASSO Cox regression to eliminate false positive factors related to prognosis. Using a minimum λ = 0.09, two hypoxia-related hub genes, PSMD11 and PSMD14, were finally selected to construct the risk prognosis model for patients with AML (Figure 3B). The risk scores for patients with AML were calculated by assessing the expression levels of PSMD11 and PSMD14 and applying risk coefficients according to the following formula: Risk Score = PSMD11*0.150603712635304+PSMD14*0.193834319424803.

The training set consisted of 187 AML patient samples downloaded from the TARGET database, which were further categorized into a high-risk group (n = 94) and a low-risk group (n = 93) based on the median risk score. Moreover, survival analysis was performed on the two groups, with the KM survival curve demonstrating a lower OS rate in the high-risk group (p < 0.05) (Figure 3C). Additionally, we analyzed the risk scores, survival status and PSMD11 and PSMD14 expression in the two groups. The high-risk group exhibited poorer clinical outcomes, whereas the expression levels of the two hub genes were significantly upregulated in this group (Figure 3D). Furthermore, the training set model was assessed using ROC-AUC to determine its accuracy. The AUC values for predicting the 1-, 3- and 5-year OS of patients with AML using this model were found to be 0.62, 0.68 and 0.68 respectively (Figure 3E). Thus, collectively, these results indicate that the hypoxia-related genes PSMD11 and PSMD14-based AML risk prognostic model demonstrates good accuracy in predicting prognosis.

To confirm the reliability of the risk prognostic model, external dataset was employed for validation. A total of 151 AML clinical samples downloaded from the TCGA database was used as the validation set and were stratified into high-risk (n = 76) and low-risk (n = 75) groups based on the median risk scores. Similar to the training set, the KM survival demonstrated a lower OS in the high-risk group of the validation set (p < 0.05) (Figure 4A). Analyses were then conducted on the risk scores, survival status and core gene expression levels of AML patient samples in the validation set. The results obtained were consistent with those observed in the training set (Figure 4B). Finally, the AUC values of predicting the prognosis of 1-, 3- and 5-year AML patients were 0.65, 0.56 and 0.62, respectively (Figure 4C). This indicates that the model maintains its prediction ability in the validation set.

The influence of hypoxia on the progression of AML through alterations in the TIME remains unclear. Therefore, the correlation between the risk prognostic model and the TIME of AML was investigated in this study. First, we utilized a heatmap to illustrate the distribution of immune cells between the high- and low-risk groups (Figure 5A). CIBERSORT combined with the LM22 characteristic matrix was adopted to analyze the disparity in immune cell infiltration between the two groups. Stacked bar chart represented the overall calculation of immune cell infiltration for each sample (Figure 5B). Additionally, the box plots demonstrated that the infiltration degree of B cells naive, mast cells resting and T cells CD4 memory resting was higher in the high-risk group, whereas that of eosinophils, macrophage M2 and monocytes was higher in the low-risk group (p < 0.05) (Figure 5C). These findings indicate that the activity of the immune microenvironment was lower in the high-risk group. Finally, the correlation between the risk score matrix and the immune cell matrix was statistically analyzed using Mantel test. The heatmap of the correlation matrix revealed a significant association between the risk score and NK cells as well as mast cells (0.01< Mantel’s p < 0.05, Person’r = 0–1) (Figure 5D).

In order to explore the biological significance of the risk prognosis model, we performed drug predictions in the high- and low-risk groups. The “oncoPredict” package was utilized to identify 22 anti-tumor drugs associated with risk scores based on the GDSC database.The heatmap illustrated the distribution of drug susceptibility between the two groups (Figure 6A). The box plots of 22 drugs demonstrated that the high-risk group was more sensitive to the mTOR inhibitor rapamycin and the dual ATP competitive PI3K and mTOR inhibitor Dactolisib which were consistent with the results of GSEA. In other words, hypoxia in AML was found to activate the mTOR signaling pathway, with high-risk groups characterised by high expression of hypoxia-related genes, exhibiting higher sensitivity to mTOR inhibitors. Furthermore, the high-risk group also displayed lower resistance to protease receptor inhibitor bortezomib and protease inhibitor MG132, suggesting that these drugs are effective in patients with high expression of PSMD11 and PSMD14, components of the 26S proteasome complex. (p < 0.05) (Figure 6B). Consequently, these anti-tumor drugs may yield better therapeutic outcomes in the high-risk group. Collectively, the risk prognosis model holds biological significance.

To further verify the influence of hub genes PSMD11 and PSMD14 in the risk prognostic model on patient prognosis, their expressions in bone marrow blood samples of AML patients and normal donors were examined. Real-time PCR and Western blot results revealed that the mRNA and protein expression levels of PSMD11 and PSMD14 were highest in the relapse group and lowest in the normal donors group, with significantly higher levels observed in the relapse group compared to the newly diagnosed group (p < 0.05) (Figures 7A, B). Notably, within a subset of 11 patients providing both newly diagnosed and relapse samples, the mRNA expression of PSMD11 and PSMD14 were consistently higher in the relapse group than that in the newly diagnosed group of the same patient (p < 0.05) (Figure 7C). Then, we divided the newly diagnosed clinical samples into 3 risk categories as favorable, intermediate, adverse groups based on the updated ELN AML Guidelines for 2022 (Döhner et al., 2022), and calculated risk scores for these clinical samples. The analysis revealed that the risk scores of clinical samples in the adverse group was significantly higher than those of the favorable group (p < 0.05) (Figure 7D). Collectively, these results further confirm the reliability of the risk prognosis model, and highlighted the close association of PSMD11 and PSMD14 with the poor AML prognosis.

Finally, to confirm the association between the risk prognostic model and the resistance towards conventional chemotherapeutic drug, we assessed the viability of AML-sensitive cell lines U937 and THP-1, as well as drug-resistant cell lines U937R and THP-1R under varying concentrations of Ara-C was examined using CCK8. The IC₅₀ values of U937 and THP-1 were 2.152 and 2.544 μM, whereas those of U937R and THP-1R were 126.5 and 131.7 μM, respectively, indicating that U937R, THP-1R cells were 58.8 and 51.8 times more resistant than U937and THP-1 cells (Figure 8A). Real-time PCR and Western blot results showed that the mRNA and protein expression levels of PSMD11 and PSMD14 were higher in AML drug-resistant cell lines (p < 0.05) (Figures 8B, C). These findings provide evidence that elevated levels of PSMD11 and PSMD14 are closely related to AML drug resistance, indicating a potential role of hypoxia in contributing to poor AML prognosis through drug resistance.