key: cord-0069655-5ggt53ip
authors: Zhong, Jinman; Wu, Hang; Bu, Xiaoyin; Li, Weiru; Cai, Shengchun; Du, Meixue; Gao, Ya; Ping, Baohong
title: Establishment of Prognosis Model in Acute Myeloid Leukemia Based on Hypoxia Microenvironment, and Exploration of Hypoxia-Related Mechanisms
date: 2021-10-26
journal: Front Genet
DOI: 10.3389/fgene.2021.727392
sha: 98ec793051df030ca918522dc8e962ffec475b5b
doc_id: 69655
cord_uid: 5ggt53ip

Acute myeloid leukemia (AML) is a highly heterogeneous hematologic neoplasm with poor survival outcomes. However, the routine clinical features are not sufficient to accurately predict the prognosis of AML. The expression of hypoxia-related genes was associated with survival outcomes of a variety of hematologic and lymphoid neoplasms. We established an 18-gene signature-based hypoxia-related prognosis model (HPM) and a complex model that consisted of the HPM and clinical risk factors using machine learning methods. Both two models were able to effectively predict the survival of AML patients, which might contribute to improving risk classification. Differentially expressed genes analysis, Gene Ontology (GO) categories, and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed to reveal the underlying functions and pathways implicated in AML development. To explore hypoxia-related changes in the bone marrow immune microenvironment, we used CIBERSORT to calculate and compare the proportion of 22 immune cells between the two groups with high and low hypoxia-risk scores. Enrichment analysis and immune cell composition analysis indicated that the biological processes and molecular functions of drug metabolism, angiogenesis, and immune cell infiltration of bone marrow play a role in the occurrence and development of AML, which might help us to evaluate several hypoxia-related metabolic and immune targets for AML therapy.

Acute myeloid leukemia (AML) is a clonal malignant aggressive hematological tumor, resulting in the accumulation of acquired chromosomal, genetic, and epigenetic abnormalities in highly heterogeneous myeloid precursors. It is the most common acute leukemia and accounts for approximately 80% of cases in adults. In the United States, the age-adjusted incidence of AML is 4.3 per 100,000 population annually (Shallis et al., 2019) , which has a high mortality rate and variable prognosis. In recent years, the incidence rate of AML is getting increasingly serious and poses an enormous threat to human health. In the research of AML, there remain several challenges, advances in treatment for AML have remained quite limited, and the current prognostic evaluation system cannot completely distinguish the prognosis of AML patients. One of the hotspots and critical points for medical research is to identify specific prognostic factors that may help predict the outcomes.

Hypoxia is a common condition in the solid-tumor microenvironment (Kim et al., 2007) , playing an important role in various biological processes, such as metabolic alteration, angiogenesis, and metastasis (Cosse and Michiels, 2008; Godet et al., 2019) . However, the role of the pathophysiological implications of hypoxia in AML remains controversial, and the mechanism is still not clear. In bone marrow (BM), the low oxygen partial pressure (pO 2 ) is physiological (Harrison et al., 2002) . And the hypoxic microenvironmental niches within leukemic BM compared with those of the normal BM were expanded, accompanied by leukemia stem cell (LSC) proliferation (Keith and Simon, 2007; Benito et al., 2011; Benito et al., 2013; Zhou et al., 2016) . The hypoxic BM microenvironment has also been shown to contribute to acute leukemic progression, resistance to chemotherapy, and minimal residual disease (MRD) (Frolova et al., 2012; Matsunaga et al., 2012; Tabe and Konopleva, 2015) . In response to hypoxia, cells change their hypoxia-related gene expression, which was proved to be correlated with prognosis for various solid tumors.

However, the European LeukemiaNet (ELN) 2017 risk classification (ELN 2017) (Dohner et al., 2017) , an important AML risk stratification standard that has been widely used to estimate prognosis of AML, is based on cytogenetic and molecular features. Mutations in the FMS-like tyrosine kinase 3 gene (FLT3-ITD) are quite common in AML and have been associated with poorer overall survival (OS) (Kayser et al., 2009) . Nucleophosmin (NPM1) gene mutations have been associated with improved outcomes in patients with AML (Becker et al., 2010) . Mutations of the CCAAT/enhancer binding protein alpha (CEBPA) gene have been associated with a favorable outcome in patients with AML, but mainly in those patients with cytogenetically normal AML (Renneville et al., 2009; Rockova et al., 2011) . Rare studies to date have developed a hypoxia-related prognosis model (HPM) of AML based on gene expression profiles. AML-suitable hypoxia gene signatures still need to be developed.

To evaluate the potential utility of hypoxia-related gene expression profiles in AML prognosis, The Cancer Genome Atlas (TCGA) and MsigDB databases were analyzed, and the clinical features of patients were considered to construct an 18gene-based hypoxia risk classifier. The model could be useful for the prognostic evaluations and development of novel therapeutic modalities aimed at interfering with hypoxia-sensing pathways and modifying the hematopoietic microenvironment.

Acquiring and Pre-Processing of Sample Data and Primary Screening of Acute Myeloid Leukemia Hypoxia-Related Genes

In this study, three public accessible transcriptome datasets of BEATAML1.0 (Cohort 1, https://portal.gdc.cancer.gov/projects/ BEATAML1.0-COHORT; Supplementary Data Sheet S1), TARGET-AML (Cohort 2, https://portal.gdc.cancer.gov/ projects/TARGET-AML; Supplementary Data Sheet S2), and TCGA-LAML (Cohort 3, https://portal.gdc.cancer.gov/projects/ TCGA-LAML; Supplementary Data Sheet S3) were used throughout the training and validation stages. The latest clinical follow-up information was also obtained from Vizome (http://www.vizome.org/; Supplementary Data Sheet S4) (Tyner et al., 2018) or TCGA database (https://cancergenome.nih.gov/; Supplementary Data Sheet S5, S6). A total of 315 hypoxiarelated genes defined in the Molecular Signatures Database (http://www.gsea-msigdb.org/gsea/msigdb/; Supplementary Table S1) were used as the initial candidates for Cox and least absolute shrinkage and selection operator (LASSO) survival analysis. We applied strict quality control (QC) for these datasets on sample and gene levels, respectively. On the sample level, we removed patients diagnosed with myelodysplastic syndromes (MDSs), myeloproliferative neoplasms (MPNs), or other non-AML diseases; we also filtered out individuals with no survival information. On the gene level QC, we kept genes having expression information in all three transcriptome datasets; genes with low expression quantity in all samples [reads per kilobase of transcript per million mapped reads (RPKM) < 1] were removed from downstream analysis. The principal component analysis (PCA) was also performed to identify the outlier; we removed individuals who deviated from the study samples. The ComBat method was used to correct the potential batch effects of RNA sequencing (Supplementary Figure S1 ). At this point, a total of 419 samples in Cohort 1, 156 samples in Cohort 2, and 151 samples in Cohort 3 (315 hypoxia-related genes) were kept for survival analysis (Supplementary Data Sheet S7-S9).

Univariate Cox proportional hazards regression was first used to preliminarily screen the AML prognostic genes (p < 0.05). Next, Cohort 1 was randomly divided into a training set of 293 cases and a test set of 126 cases (7:3 ratio). To narrow down the prognostic genes for prediction, a Cox proportional hazards regression model combined with the LASSO (Goeman, 2010) using the "glmnet" package was applied to select the most important hypoxia-related signatures, and the optimal values of the penalty parameter λ were determined by 10-fold cross-validations at which the minimal mean squared error (MSE) is achieved in the training set (Simon et al., 2011) . Afterward, the multivariate Cox regression analysis was performed to estimate independent prognostic factors associated with patient survival. Finally, the stepwise method was employed to select the best subset of predictors in a risk score model. To this, a hypoxia-related prognostic risk (HRS) score model was built, with the regression coefficients (β) weighted by the multivariate Cox proportional hazards regression model in the training set. The HRS model formula was as follows: We calculate the HRS of the study samples and use the median HRS of the training group as the cutoff point to label the low-risk or high-risk individual of the three cohorts. The trained model was then tested using the test and validation sets (other independent cohorts: Cohort 2 and Cohort 3). The "survminer" package was used for the Kaplan-Meier (KM) survival analysis of patients in the high-risk and low-risk groups, while the "timeROC" package was used to construct time-dependent receiver operating characteristic (ROC) curves and calculate the area under the ROC curve (AUC) at 1-, 3-, and 5-years OS.

We explored the relationship between HPM and other clinical parameters for AML patient outcomes. Univariable Cox analysis and multivariate Cox analysis were performed with all patients' clinical covariates in the BEATAML1.0 cohort by the "rms" package. Samples with incomplete data about potential prognostic factors were excluded from the multivariable Cox analyses. Pearson's correlation between HPM and different clinical characteristics was calculated by "stats" package and plotted by "corplot" package. A nomogram was formulated using the "rms" r package based on the results of the multivariate analyses; and calibration plots and time-dependent ROC plots were performed to assess the prognostic accuracy of the nomogram. The predicted outcomes and observed outcomes of the nomogram were plotted in the calibration curve to evaluate the degree of fitting of the nomogram, and the 45°line represented the best prediction.

We procured 3,695 human metabolic genes concerning 145 metabolic subsystems from the Recon 3D (http://vmh.life) (Brunk et al., 2018) . Among them, 3,224 metabolic genes were matched to our Cohort 1 data. The empirical Bayes algorithm of the R package "limma" (Ritchie et al., 2015) was used to identify differentially expressed genes (DEGs) between the top 25% samples of the high-and low-risk groups (totally 210 samples). All gene expression values were log2 transformed to identify the metabolic genes with significant differential expression during hypoxia stress [logarithmically transformed fold change (log2(FC)) ≥ 1 or (log2(FC)) ≤ −1 and p-value < 0.05]. Besides, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) enrichment analyses of DEGs with clusterProfiler package (Yu et al., 2012) . KEGG pathway enrichment analysis utilized the KEGG database (http://www.genome.jp/kegg) while GO enrichment was utilized (http://www.geneontology.org).

Regarding the association between the hypoxia risk score and immune cells in the BM microenvironment, CIBERSORT algorithm (Newman et al., 2015) was used to estimate the relative immune cell fractions in the BM samples of Cohort 1, based on the standard LM22 leukocyte signature matrix that distinguishes 22 immune cell subtypes and 1,000 permutations (CIBERSORT R script v1.03 is available on http://cibersort.stanford.edu/). We performed the following analyses in CIBERSORT: B cells, CD4 + T cells, CD8 + T cells, Tregs, NK T cells, γδ-T cells, lymphocytes, and macrophages. Further, xCell (http://xcell.ucsf.edu/) was used to validate the result.

A flowchart overviewing the procedures of this study is presented in Figure 1 .

We analyzed clinical characteristics of patients with AML from Cohort 1, including age, gender, ELN 2017, the mutations of NPM1 and FLT3, and CEBPA Biallelic status ( Table 1) . Differences in general clinical information of two sets are not statistically significant. Among the 419 AML patients with complete clinical information in BEATAML1.0, the mean diagnosis age was 56.22 ± 18.23 years, while the proportion of males was 55.61% (233/419).

A univariate Cox regression was performed to identify prognostic hypoxia-related genes associated with OS in BEATAML1.0 dataset. A total of 33 genes (ALDH1A1, ALDOC, BACE2, BATF3, CA9, CALD1, COL5A1, DR1, EGLN3, ELOB, HBP1, HK1, ID2, KRT14, LRP8, NOS1, NOS2, PDK3, PLOD2, PSMA2, PSMA7, PSMB6, PSMC1, PSMC4, PTGS1, RPS27A, SIAH2, SLC16A1, SORL1, TGM2, THBS1, TPD52, and UBA52) were identified to have a significant prognostic value in patients with AML (p < 0.05). Then, LASSO-penalized Cox analysis with 10-fold cross-validation (Figures 2A,B ) was performed for further screening, and 23 genes were left. Finally, a total of 18 hub genes (ALDOC, BATF3, COL5A1, DR1, ELOB, HBP1, HK1, KRT14, NOS2, PSMA2, PSMA7, PSMB6, PSMC1, PTGS1, SIAH2, SORL1, THBS1, and UBA52) were identified from the stepwise multivariate Cox regression ( Figure 2C ), and the formula to calculate the hypoxia risk score was as follows:

October 2021 | Volume 12 | Article 727392

Prognostic Value of the Hypoxia Risk Score

The KM survival curves revealed that patients in the high-risk group exhibited a significantly lower OS rate than the low-risk group in all cohorts (training set, test set, Cohort 2, and Cohort 3;

p < 0.0001, p 0.04, p < 0.001, and p < 0.01, respectively; Figure 3 ). In the training dataset of BEATAML1.0 Cohort, the median OS of low-risk patients was 1.970 years (95% CI: 1.734-3.775), whereas the median OS of high-risk patients was Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 727392 5 median OS of patients with high-risk scores was 2.227 years (95% CI: 1. 689-5.195) , and the HR is 2.283 (95% CI: 1.455-3.582). In TCGA-LAML cohort (validation Cohort 3), the median OS of patients with low-risk scores was 2.170 years (95% CI: 1.581-3.838), the median OS of patients with high-risk scores was 0.833 years (95% CI: 0.586-1.003), and the HR is 1.773 (95% CI: 1.190-2.642). Furthermore, the prognostic accuracy of the hypoxia risk score was assessed with time-dependent ROC analysis for OS at 1, 3, and 5 years in all datasets, and the results were as follows. In the training set, the AUC was 0.813 at 1 year, 0.788 at 3 years, and 0.899 at 5 years ( Figure 3E ); 0.675 at 1 year, 0.767 at 3 years, and 0.753 at 5 years ( Figure 3F ) in the test set; 0.616 at 1 year, 0.684 at 3 years, and 0.690 at 5 years ( Figure 3G ) in Cohort 2; and 0.712 at 1 year, 0.657 at 3 years, and 0.640 at 5 years ( Figure 3H ) in Cohort 3, which revealed that the hypoxia risk score was a valuable predictor. We next studied whether the HRS could improve prognostic assessment in AML patients on the 2017 ELN genetic risk stratification (Dohner et al., 2017) basis. In each ELN2017 stratification, the HPM high-risk group had a poorer prognosis than the HPM low-risk group ( Figure 3I) . Furthermore, we found that 37 out of 108 ELN2017 favorable patients (51.7%) were high-risk for the HPM and had significantly worse survival. In this patient subgroup (ELN favorable-HPM high), representing 9.3% (37/397, only contains patients who had clear ELN2017 stratification) of the BEATAML1.0 cohort, their survival was similar to that of ELN intermediate/ adverse-risk patients with HPM low-risk. To further evaluate the performance of the HPM, we compared our HPM with other gene expression based AML prognostic models, which were published within the last 5 years, including PMID29138577 (Huang et al., 2017) , PMID32268820 (Zhang and Xiao, 2020) , PMID29956722 (Zhao et al., 2018) , and PMID34282207 (Jiang et al., 2021) . Risk score was calculated based on formulas from the corresponding literatures. We used the p-value of the KM survival analysis to reflect the discrimination and the AUC value to evaluate the accuracy. The survival analysis showed that the HPM was significantly correlated with the survival of the patients in all three datasets. In contrast, other models could only perform well in at most two of the three datasets (Supplementary Figure S2) . The AUC values of the HPM were more stable, which means a wider range of suitable population (Supplementary Figure S3 ; Supplementary Table S2 ).

We ranked the risk scores of patients in the training and test sets, and the distributions associated with gene expression, survival time, and status are shown in Figures 4A,B . Hypoxia risk score increased with higher age, worse 2017 ELN genetic risk stratification, NPM1 wild type, and FLT3-ITD wild type and have no CEBPA double mutation ( Figure 4C ).

To study the relationship between HRS and metabolic flux, we performed differential gene expression analysis using "limma" package, 115 DEGs were shortlisted from the raw dataset, and 93 genes were upregulated in the high-risk group while 22 genes were downregulated. Results were visualized into volcano plot ( Figure 5A ) and heatmap ( Figure 5B ). To evaluate the molecular mechanisms of DEGs, KEGG metabolic pathways enrichment analyses, and GO functional annotation were conducted (Table 2; Figure 5 ). Unsurprisingly, pathways associated with oxidoreduction activity were enriched. Unexpectedly, we found 12 immune-related pathways such as "leukocyte mediated immunity" and "cell activation involved in immune response," which were significantly overrepresented. Besides, two pathways of angiogenesis were significantly enriched in the upregulated genes for the high-risk group. And the enrichment of "heme binding" and "tetrapyrrole binding" might be associated with oxygen-carrying capacity. Moreover, drug catabolic process pathways were markedly enriched, such as "xenobiotic metabolic process," "drug catabolic process," "Drug metabolism-other enzymes," "Drug metabolism-cytochrome P450," and "Metabolism of xenobiotics by cytochrome P450." Moreover, the enrichment of the pathway "Chemical carcinogenesis" suggests different disease susceptibility between the high-and low-risk groups.

Based on the foregoing GO enrichment analysis results, we conjecture that hypoxia is related to immunity. Using CIBERSORT, we found the proportion of resting mast cells and plasma cells in the high-risk group significantly to be lower than that of the low-risk group, whereas the proportion of neutrophils, monocytes, M0 macrophages, γδ-T cells, and regulatory T cells (Tregs) increased in high-risk patients ( Figure 6A ; Supplementary Data Sheet S10). Cell types were also predicted and visualized using xCell, giving similar results. Further, xCell also shows that the Immune Score and Microenvironment Score of the high-risk group were significantly higher than of the low-risk group ( Figure 6B ). An overview of the predictive result of mechanism exploration is presented in Figure 7 .

To explore the independent prognostic factors for AML, univariate and multivariate Cox analyses were sequentially performed in the BEATAML1.0 dataset ( Figure 8A ), including the hypoxia risk score and other available clinical characteristics, such as age, gender, ELN 2017, the mutations of NPM1 and FLT3, and CEBPA Biallelic status. The HRS, ELN 2017, and age remained statistically significant (p < 0.05) in both the univariate and multivariate Cox analyses, indicating that HRS, ELN 2017, and age were independent prognostic factors. NPM1 status was an independent factor of prognosis after adjustment for other clinical factors. The correlation between clinical factors was analyzed and visualized in Figure 8B . HRS was positively correlated with survival status, age, and ELN2017 while negatively correlated with survival time and NPM1 status. Furthermore, Figures 8C -E indicate that the HRS had a higher AUC than other clinical factors in 1-, 3-, and 5-years survival prediction. Overall, these results demonstrated that the HPM can predict the AML prognosis independently and effectively. To reveal the prognostic value, maximize practicability, and facilitate clinicians' usage of our model, we constructed a nomogram that was composed of both the hypoxia risk score and available clinical risk factors based on BEATAML1.0 cohort ( Figure 8F) . A combination of HRS and clinical risk factors was found to improve its prognostic value with a markedly better AUC than the 2017 ELN genetic risk stratification ( Figure 8G ). It is shown that the complex model can predict more accurately in the long term with the increasing tendency of AUC. To validate the predicted and actual probabilities at 1, 3, and 5 years, calibration plots were constructed ( Figures 8H-J) , and the nomogram performs well. These findings demonstrated that the nomogram is an optimal model for predicting the survival probability of AML patients than individual prognostic factors.

AML is a highly aggressive and heterogeneous hematologic malignancy (Döhner et al., 2015) . Hypoxia is a significant outcome factor of leukemia patients, which could be reflected by the changed expression of related genes . Here, in our analysis, we identified an 18-hypoxia Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 727392 11 related gene-based prognosis model that could independently predict the survival probability of AML patients. To facilitate clinical application, we combined our HPM, 2017 ELN genetic risk stratification, and clinical risk factors to construct a complex model and nomogram, which outperformed 2017 ELN genetic risk stratification in the prediction of survival rate. Thus, the comprehensive nomogram could be utilized by clinicians in the near future. Further, we surmise that hypoxia may affect the prognosis of AML patients by affecting the angiogenesis, drug metabolism, and BM immune microenvironment.

The risk model created in this study consisted of 18 hypoxiarelated genes, many of which were reported in cancer. ALDOC encodes a member of the class I fructose-biphosphate aldolase gene family and is involved in HIF-1 signaling pathway. ALDOC was identified as activators of Wnt signaling, a signaling pathway involved in cancer genesis and progression when it was over-activated (Caspi et al., 2014) . Meanwhile, ALDOC was overexpressed in gallbladder carcinoma (Fan et al., 2020) , melanoma (Izraely et al., 2020) , and lung cancer (Yuan et al., 2021) , associated with their growth or pathogenesis. BATF3 is an AP-1 family transcription factor that controls the differentiation of CD8(+) thymic conventional dendritic cells in the immune system. According to the KEGG database, it was involved in PD-L1 expression and PD-1 checkpoint pathway in cancer. Immunotherapy with immunomodulatory monoclonal antibodies targeting PD-1 or CD137 requires Batf3dependent dendritic cells (Sanchez-Paulete et al., 2016) . COL5A1 is a member of the fibrous subfamily of collagen. Overexpression of COL5A1 may promote metastasis of lung adenocarcinoma (Liu et al., 2018) and the progression of muscle-invasive bladder cancer (Ewald et al., 2013) and may increase the risk of hematogenous and lymphatic metastasis in serous ovarian cancer (Yue et al., 2019) . COL5A1 is also overexpressed in gastric cancer, which may regulate the proliferation of gastric cancer cells by affecting the tumor microenvironment and is associated with poor prognosis (Wei et al., 2020) . HBP1 plays a role in the regulation of the cell cycle and is a tumor suppressor (Bollaert et al., 2019) . Downregulating HBP1 promotes the migration and invasion of oral squamous cell carcinoma and breast cancer (Li et al., 2011) . HK1 encodes hexokinase 1, which is the first rate-limiting enzyme in glycolysis, is related to the progression of ovarian cancer (Li Y. et al., 2020) and colorectal cancer (Li S. et al., 2020) . PSMA7 interacts with proteins such as HIF-1α, EMAP II, c-Abl, and Arg tyrosine kinases, which participated in tumorigenesis. Studies reported that PSMA7 expression was elevated in testicular, liver, breast, prostate, cervical, gastric, and colorectal cancers (Xia et al., 2019) , while UBA52 is overexpressed in the colon (Barnard et al., 1995) , and renal cancers (Kanayama et al., 1991) . THBS1 was also implicated in the development of several cancers, including breast, gastric, melanoma, and cervical cancers and glioblastoma (Qi et al., 2021) .

To probe into the mechanism of hypoxia in leukemia, differential gene expression analysis and enrichment analysis were implemented. The enriched biological process like "leukocyte mediated immunity" and "cell activation involved in immune response" suggested that hypoxia might affect cellmediated immunity against cancer cells. Angiogenesis, which was particularly important for tumor survival in the hypoxic condition (Hanahan and Weinberg, 2011) , was significantly enriched in the upregulated genes for the high-risk group.

Although there are few studies about angiogenesis in hematological malignancy, BM angiogenesis in AML patients has been observed and may play a role in the pathogenesis (Hussong et al., 2000; Padro et al., 2000; Testa et al., 2020) . BM microvessel density of AML patients is higher than that of healthy individuals at the time of diagnosis and decreases after remission (Padro et al., 2000; Song et al., 2015) . The vascular niches could support the survival of leukemic cells and protect AML by regulating AML cell cycle through paracrine secretion and adhesive contact with endothelial cells, helping to resist chemotherapy (Cogle et al., 2016) . Moreover, higher microvessel density at the time of diagnosis was associated with poor prognosis (Kuzu et al., 2004; Rabitsch et al., 2004) . In addition, drug metabolism is significantly affected by changes in pharmacokinetics, expression, and function of drug metabolic enzymes and transporters under hypoxia (Donovan et al., 2010) . Hypoxia affects the transcription and function of cytochrome P450 (CYP450) through HIF-1α (Min et al., 2019) . CYP enzymes and other drug-metabolizing enzymes are expressed in BM stroma (Alonso et al., 2015) . CYP450 is involved in the stroma-mediated resistance of AML cells to chemotherapy (Alonso et al., 2015) . Drug catabolic process pathways, such as "xenobiotic metabolic process," "drug catabolic process," "Drug metabolism-other enzymes," "Drug metabolism-cytochrome P450," and "Metabolism of xenobiotics by cytochrome P450," were markedly enriched, which might correlate with chemotherapeutic drugs concentration and sensitivity change in a hypoxic microenvironment of AML. Hypoxia might promote angiogenesis, disturb cell-mediated immunity balance, and affect pharmacokinetics to result in a bad prognosis.

To explore hypoxia-related changes in the BM immune microenvironment, immune cell composition analysis was performed. We found that the proportion of resting mast cells in the high-risk group was significantly lower than in the low-risk group, while activated mast cells in the high-risk group were higher than in the low-risk group, even though the difference did not reach statistical significance (mean proportion 0.12 vs. 0.02, p 0.22). The present study shows that mast cells can promote cancer growth by stimulation of neoangiogenesis and modulation of the immune response (Dyduch et al., 2012) . This might indicate that hypoxia affects angiogenesis and immune response by affecting the proportion of mast cells. The proportion of monocytes and M0 macrophages in the high-risk group were significantly higher than the low-risk group, which contribute to the poor prognosis of high-risk group patients. Tumor-infiltrating monocytes and macrophages have well-recognized facilitative roles in the initiation, migration, and invasion of solid tumors (Condeelis and Pollard, 2006; Qian and Pollard, 2010; Richards et al., 2013) in vitro (Lee et al., 2012) . Al-Matary et al. showed an increase of monocytes/macrophages in the BM of AML patients and AML mouse models, which might support the proliferation of AML cells in vitro, and it was also observed that the grade of macrophage infiltration was correlated with the survival of AML mice (Al-Matary et al., 2016) . The proportion of monocytes and M0 macrophages in the high-risk group was significantly higher than in the low-risk group, which contributes to the poor prognosis of highrisk group patients. Tregs, which were regarded as suppressor T cells preventing autoimmunity, were found to be aberrantly accumulated in some types of tumor, playing a crucial role in dampening antitumor immunity and establishing an immunosuppressive microenvironment (Wang et al., 2017) . Higher Treg percentages could indicate poor prognosis in a variety of cancer types (Overacre-Delgoffe and Vignali, 2018), and AML is no exception. Williams et al. found that an increased amount of Tregs in peripheral blood of AML patients is also associated with an increased risk of relapse (Williams et al., 2019) . These studies are consistent with our findings that Treg proportion is higher in high-risk patients than lowrisk patients. The detailed mechanism needs further investigation in the future. And the scRNA-seq transcriptome data can provide precision and details of the interaction between the tumor cells and the microenvironment. Zhang et al. developed a novel scRNA-seq data-based approach to reconstruct a multilayer signaling network that contains pathways from intercellular ligand-receptor interactions, intracellular transcriptional factors, and their target genes . Meanwhile, the singlecell RNA-sequencing data based on multilayer network method (scMLnet) (Cheng et al., 2021) also help to resolve tumor-microenvironment interactions and dissect the microenvironment-mediated intercellular and intracellular signaling pathways of tumor cells, which might help to investigate the influence of microenvironment on the tumor growth, drug resistance, and patient prognosis.

Concerning treatment, Tregs have been targeted in the clinic, although the efficacy is limited (Overacre-Delgoffe and Vignali, 2018). Targeting macrophages with bisphosphonate could reduce angiogenesis and tumor growth in melanoma-bearing mice (Gazzaniga et al., 2007) . In a mouse model of AML, Tregs accumulate at the site of disease and suppress the function of adoptively transferred cytotoxic T cells (CTL), and depletion of Tregs restored CTL function and reduced leukemia progression (Zhou et al., 2009) . Anti-CD47 monoclonal antibody can inhibit the immune escape of AML leukemic stem cells to macrophages and play an antileukemic role by phagocytosis of leukemic stem cells through macrophages (Majeti et al., 2009) . Clinical trials of CD47 mAb magrolimab (Hu5F9-G4) as a single agent or in combination for the treatment of relapsed refractory AML have shown promising efficacy (magrolimab, Phase I, ClinicalTrials.gov ID: NCT02678338; magrolimab + atezolizumab, Phase I, ClinicalTrials.gov ID: NCT03922477; magrolimab + azacitidine, Phase Ib, ClinicalTrials.gov ID: NCT03248479). Antiangiogenic therapy may be an effective method in AML patients. Reduced BM angiogenesis may help to restore drug sensitivity of drug-resistant AML (Lin et al., 2019) . Antiangiogenic therapy may be an effective method in AML patients. Lin et al. demonstrated that wogonoside, one of the metabolites of traditional Chinese medicine Huangqin, could inhibit the BM angiogenesis and tumor progression of AML in vivo and in vitro (Lin et al., 2019) . Based on the results of clinical trials, some antiangiogenic drugs that inhibit vascular endothelial growth factor (VEGF), such as bevacizumab (Karp et al., 2004) , cediranib (Fiedler et al., 2010) , AG-013736 (Giles et al., 2006) , and SU5416 (Fiedler et al., 2003) , could be an effective treatment for AML, either alone or in addition to chemotherapy that works independently on different targets. Hypoxia is becoming an emerging target in AML. Hypoxia-activated prodrug, a new class of anti-cancer agents, selectively deliver cytostatic or cytotoxic agents to hypoxic subregions, uncloak at low oxygen pressure, and release the active drug. Small-scale clinical trials about hypoxia-activated prodrugs PR-104 and TH-302 treating patients with relapsed and/or refractory AML were conducted, showing a definite antileukemia activity Badar et al., 2016) (TH-302, Phase 1, clinicaltrials.gov ID: NCT01149915; PR-104, Phase 1, ClinicalTrials.gov ID: NCT01037556). BCL-2, a proapoptotic protein, could be overexpressed in hypoxia conditions; its inhibitors can reduce oxidative phosphorylation and eradicate quiescent chemoresistant AML stem cells (Ashton et al., 2018) . Echinomycin, a hypoxia-inducible factor HIF-1α inhibitor, can selectively kill the leukemia-initiating cell without affecting host HSCs in relapsed AML mice (Wang et al., 2014) . The preliminary pathophysiological observation of this study may provide a perspective for further investigation and a potential therapeutic target in the future.

Although large cohorts were utilized to establish our model, there are still certain limitations in our study. The study was based on retrospective cohorts, lacking prospective cohorts, and experimental evidence. Data from different centers and various platforms are necessary to validate the performance of our model. Further studies including animal experiments and in vitro cellular experiments will be needed to confirm our findings and delineate the pathophysiological mechanism.

In summary, our HPM, complex model, and nomogram had excellent predictive power for clinical applications, helping clinicians in making clinical decisions. Our hypoxia-related immune and metabolic alterations might help to find a potential therapeutic target.

The original contributions presented in the study are included in the article material. The R Code used for this analysis is available as an additional file (Supplementary Data File 11). Further inquiries can be directed to the corresponding authors.

Conception and design were contributed by BP and YG. Collection and assembly of data were contributed by XB, WL, Frontiers in Genetics | www.frontiersin.org October 2021 | Volume 12 | Article 727392 SC, MD, JZ, and HW. Data analysis and interpretation were contributed by JZ and HW. Article writing was contributed by JZ and HW. All authors read and approved the final article.

Acute Myeloid Leukemia Cells Polarize Macrophages towards a Leukemia Supporting State in a Growth Factor independence 1

Human Bone Marrow Niche Chemoprotection Mediated by Cytochrome P450 Enzymes

Oxidative Phosphorylation as an Emerging Target in Cancer Therapy

Phase I Study of Evofosfamide, an Investigational Hypoxia-Activated Prodrug, in Patients with Advanced Leukemia

Ubiquitin Fusion Proteins Are Overexpressed in colon Cancer but Not in Gastric Cancer

Favorable Prognostic Impact ofNPM1Mutations in Older Patients with Cytogenetically Normal De Novo Acute Myeloid Leukemia and Associated Gene-and MicroRNA-Expression Signatures: A Cancer and Leukemia Group B Study

Pronounced Hypoxia in Models of Murine and Human Leukemia: High Efficacy of Hypoxia-Activated Prodrug PR-104

Targeting Hypoxia in the Leukemia Microenvironment

Hypoxia-Activated Prodrug TH-302 Targets Hypoxic Bone Marrow Niches in Preclinical Leukemia Models

The HMG Box Transcription Factor HBP1: a Cell Cycle Inhibitor at the Crossroads of Cancer Signaling Pathways

Recon3D Enables a Three-Dimensional View of Gene Variation in Human Metabolism

Aldolase Positively Regulates of the Canonical Wnt Signaling Pathway

Inferring Microenvironmental Regulation of Gene Expression from Single-Cell RNA Sequencing Data Using scMLnet with an Application to COVID-19

Acute Myeloid Leukemia in the Vascular Niche

Macrophages: Obligate Partners for Tumor Cell Migration, Invasion, and Metastasis

Tumour Hypoxia Affects the Responsiveness of Cancer Cells to Chemotherapy and Promotes Cancer Progression

Acute Myeloid Leukemia

Diagnosis and Management of AML in Adults: 2017 ELN Recommendations from an International Expert Panel

Mast Cells and Cancer: Enemies or Allies?

Expression Microarray Meta-Analysis Identifies Genes Associated with Ras/MAPK and Related Pathways in Progression of Muscle-Invasive Bladder Transition Cell Carcinoma

Terminal Binding with ALDOC Disrupts the Ability of ALDOC to Sense Glucose and Promotes Gallbladder Carcinoma Growth

A Phase 2 Clinical Study of SU5416 in Patients with Refractory Acute Myeloid Leukemia

An Open-Label, Phase I Study of Cediranib (RECENTIN) in Patients with Acute Myeloid Leukemia

Regulation of HIF-1α Signaling and Chemoresistance in Acute Lymphocytic Leukemia under Hypoxic Conditions of the Bone Marrow Microenvironment

Targeting Tumor-Associated Macrophages and Inhibition of MCP-1 Reduce Angiogenesis and Tumor Growth in a Human Melanoma Xenograft

The Anti-angiogenesis Agent, AG-013736, Has Minimal Activity in Elderly Patients with Poor Prognosis Acute Myeloid Leukemia (AML) or Myelodysplastic Syndrome (MDS)

Fatemapping post-hypoxic Tumor Cells Reveals a ROS-Resistant Phenotype that Promotes Metastasis

L1 Penalized Estimation in the Cox Proportional Hazards Model

Hallmarks of Cancer: the Next Generation

Oxygen Saturation in the Bone Marrow of Healthy Volunteers

Identification and Validation of Potential Prognostic Gene Biomarkers for Predicting Survival in Patients with

Evidence of Increased Angiogenesis in Patients with Acute Myeloid Leukemia

The Melanoma Brain Metastatic Microenvironment: Aldolase C Partakes in Shaping the Malignant Phenotype of Melanoma Cells -a Case of Inter-tumor Heterogeneity

A Hypoxia Risk Signature for the Tumor Immune Microenvironment Evaluation and Prognosis Prediction in

Changes in Expressions of Proteasome and Ubiquitin Genes in Human Renal Cancer Cells

Targeting Vascular Endothelial Growth Factor for Relapsed and Refractory Adult Acute Myelogenous Leukemias

Insertion of FLT3 Internal Tandem Duplication in the Tyrosine Kinase Domain-1 Is Associated with Resistance to Chemotherapy and Inferior Outcome

Hypoxia-inducible Factors

Effects of Hypoxia on Tumor Metabolism

Phase I/II Study of the Hypoxia-Activated Prodrug PR104 in Refractory/relapsed Acute Myeloid Leukemia and Acute Lymphoblastic Leukemia

Bone Marrow Microvessel Density (MVD) in Adult Acute Myeloid Leukemia (AML): Therapy Induced Changes and Effects on Survival

Protumoral Role of Monocytes in Human B-Cell Precursor Acute Lymphoblastic Leukemia: Involvement of the Chemokine CXCL10

miR-17-5p Promotes Human Breast Cancer Cell Migration and Invasion through Suppression of HBP1

miR-146b Functions as an Oncogene in Oral Squamous Cell Carcinoma by Targeting HBP1

lncARSR Sponges miR-34a-5p to Promote Colorectal Cancer Invasion and Metastasis via Hexokinase-1-mediated Glycolysis

Prognostic Significance and Related Mechanisms of Hexokinase 1 in Ovarian Cancer

Wogonoside Impedes the Progression of Acute Myeloid Leukemia through Inhibiting Bone Marrow Angiogenesis

COL5A1 May Contribute the Metastasis of Lung Adenocarcinoma

CD47 Is an Adverse Prognostic Factor and Therapeutic Antibody Target on Human Acute Myeloid Leukemia Stem Cells

Elevated HIF-1α Expression of Acute Myelogenous Leukemia Stem Cells in the Endosteal Hypoxic Zone May Be a Cause of Minimal Residual Disease in Bone Marrow after Chemotherapy

Modulation of Drug-Metabolizing Enzymes and Transporters under Hypoxia Environment

Robust Enumeration of Cell Subsets from Tissue Expression Profiles

Treg Fragility: A Prerequisite for Effective Antitumor Immunity?

Increased Angiogenesis in the Bone Marrow of Patients with Acute Myeloid Leukemia

Thrombospondin-1 is a Prognostic Biomarker and is Correlated with Tumor Immune Microenvironment in Glioblastoma

Macrophage Diversity Enhances Tumor Progression and Metastasis

Bone Marrow Microvessel Density and It's Prognostic Significance in AML

The Favorable Impact of CEBPA Mutations in Patients with Acute Myeloid Leukemia Is Only Observed in the Absence of Associated Cytogenetic Abnormalities and FLT3 Internal Duplication

Monocytes and Macrophages in Cancer: Development and Functions

Limma powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies

Risk Stratification of Intermediate-Risk Acute Myeloid Leukemia: Integrative Analysis of a Multitude of Gene Mutation and Gene Expression Markers

Cancer Immunotherapy with Immunomodulatory Anti-cd137 and Anti-PD-1 Monoclonal Antibodies Requires BATF3-dependent Dendritic Cells

Epidemiology of Acute Myeloid Leukemia: Recent Progress and Enduring Challenges

Regularization Paths for Cox's Proportional Hazards Model via Coordinate Descent

Levels of Bone Marrow Microvessel Density Are Crucial for Evaluating the Status of Acute Myeloid Leukemia

Role of Microenvironment in Resistance to Therapy in AML

Functional Genomic Landscape of Acute Myeloid Leukaemia

Echinomycin Protects Mice against Relapsed Acute Myeloid Leukemia without Adverse Effect on Hematopoietic Stem Cells

Metabolic Regulation of Tregs in Cancer: Opportunities for Immunotherapy

LncRNA HOTAIR Promotes the Growth and Metastasis of Gastric Cancer by Sponging miR-1277-5p and Upregulating COL5A1

The Distribution of T-cell Subsets and the Expression of Immune Checkpoint Receptors and Ligands in Patients with Newly Diagnosed and Relapsed Acute Myeloid Leukemia

Overexpression of PSMA7 Predicts Poor Prognosis in Patients with Gastric Cancer

clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters

EPB41 Suppresses the Wnt/β-Catenin Signaling in Non-small Cell Lung Cancer by Sponging ALDOC

Gene Signature Characteristic of Elevated Stromal Infiltration and Activation Is Associated with Increased Risk of Hematogenous and Lymphatic Metastasis in Serous Ovarian Cancer

Identification and Validation of a Prognostic 8-gene Signature for Acute Myeloid Leukemia

Single-cell Transcriptome-Based Multilayer Network Biomarker for Predicting Prognosis and Therapeutic Response of Gliomas

A Novel Scoring System for Acute Myeloid Leukemia Risk Assessment Based on the Expression Levels of Six Genes

Depletion of Endogenous Tumor-Associated Regulatory T Cells Improves the Efficacy of Adoptive Cytotoxic T-Cell Immunotherapy in Murine Acute Myeloid Leukemia

Bone Marrow Niche-Mediated Survival of Leukemia Stem Cells in Acute Myeloid Leukemia: Yin and Yang

The authors thank teachers in the Department of Hematology, Huiqiao Medical Center, Nanfang Hospital and Laboratory Department, Zhujiang Hospital, for the excellent technical assistance.

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2021.727392/ full#supplementary-material Conflict of Interest: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.Publisher's Note: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.Copyright © 2021 Zhong, Wu, Bu, Li, Cai, Du, Gao and Ping. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.