Characterization of somatic mutation-associated microenvironment signatures in acute myeloid leukemia patients based on TCGA analysis

Recurrent genetic mutations occur in acute myeloid leukemia (AML) and have been incorporated into risk stratification to predict the prognoses of AML patients. The bone marrow microenvironment plays a critical role in the development and progression of AML. However, the characteristics of the genetic mutation-associated microenvironment have not been comprehensively identified to date. In this study, we obtained the gene expression profiles of 173 AML patients from The Cancer Genome Atlas (TCGA) database and calculated their immune and stromal scores by applying the ESTIMATE algorithm. Immune scores were significantly associated with OS and cytogenetic risk. Next, we categorized the intermediate and poor cytogenetic risk patients into individual-mutation and wild-type groups according to RUNX1, ASXL1, TP53, FLT3-ITD, NPM1 and biallelic CEBPA mutation status. The relationships between the immune microenvironment and each genetic mutation were investigated by identifying differentially expressed genes (DEGs) and conducting functional enrichment analyses of them. Significant immune- and stromal-relevant DEGs associated with each mutation were identified, and most of the DEGs (from the FLT3-ITD, NPM1 and biallelic CEBPA mutation groups) were validated in the GSE14468 cohort downloaded from the Gene Expression Omnibus (GEO) database. In summary, we identified key immune- and stromal-relevant gene signatures associated with genetic mutations in AML, which may provide new biomarkers for risk stratification and personalized immunotherapy.


Relationship between immune/stromal scores and cytogenetic risk.
Overall, patients in the favorable cytogenetic risk group had significantly lower immune scores and a higher 2-year OS rate. However, there were no significant differences in the 2-year OS rate, immune scores or stromal scores between the intermediate and poor cytogenetic risk groups. Therefore, the intermediate and poor cytogenetic risk groups (n = 137) were grouped together for the subsequent analysis. After grouping AML patients with intermediate and poor cytogenetic risk by the median score, those with high immune scores still had significantly lower 2-year OS rates than those with the low immune scores (23.5% [95% CI 13.2-35.5%] vs 48.8% [95% CI 34.0-62.1%], P = 0.011, log-rank test, Fig. 1f). Furthermore, patients in the high stromal score group tended to have lower 2-year OS rates compared with patients in the low stromal score group (29.1% [95% CI 18.0-41.1%] vs 43.8% [95% CI 28.9-57.8%], P = 0.14, log-rank test, Fig. 1g).
Somatic mutation-associated immune/stromal scores in the intermediate and poor cytogenetic risk groups. To further explore the association between immune/stromal scores and mutations, patients in the intermediate and poor cytogenetic risk groups were divided into two subgroups based on whether the individual somatic mutation existed, and their immune and stromal scores are presented in Fig. 2g,h.
The above analysis reflected that there were distinct relationships between somatic mutations and immune/ stromal scores. RUNX1 and TP53 mutations were related to both higher immune scores and higher stromal scores, FLT3-ITD was related to both lower immune scores and lower stromal scores, ASXL1 mutation was only related to higher stromal scores, NPM1 was only related to lower stromal scores, and biCEBPA was not related to immune scores and stromal scores. The results of GO term and KEGG pathway enrichment analysis are shown in Fig. 3. In general, the enrichment analysis results were consistent with the immune/stromal scores. In other words, both immune-and www.nature.com/scientificreports/ stromal-related GO terms and KEGG pathways were enriched for RUNX1, TP53 and FLT3-ITD mutations, which corresponded to their association with immune and stromal scores. Only stromal-related GO terms were enriched for ASXL1 mutation, which corresponded to its only association with stromal scores. However, discrepancies were observed for NPM1 and biCEBPA mutations: both immune-and stromal-related GO terms and KEGG pathways were enriched for NPM1 mutation, although it was only related to lower stromal scores; stromal-related GO terms were enriched for biCEBPA mutation, despite its lack of a relationship to immune/ stromal scores.
Characteristics and prognostic significance of mutation-associated immune/stromal cell-relevant DEGs. Mutation-associated immune/stromal score-relevant DEGs were selected according to ESTI-MATE algorithm gene lists and are shown in Supplementary Table S1. Overall, the number of immune/stromal www.nature.com/scientificreports/ cell-relevant DEGs was consistent with immune/stromal scores (Table 1): for RUNX1 and TP53 mutations, the majority of both immune and stromal cell-relevant DEGs were upregulated, which corresponded to their associated higher immune and stromal scores; for the ASXL1 mutation, only several stromal cell-relevant DEGs were upregulated, which corresponded to its associated higher stromal scores; the majority of both immune and stromal cell-relevant DEGs for FLT3-ITD and the majority of stromal cell-relevant DEGs for NPM1 mutation were downregulated, which corresponded to their associated lower scores; for biCEBPA mutation, almost no immune and stromal cell-relevant DEGs were observed, which corresponded to its lack of a relationship with immune and stromal scores. The only exception was that the majority of immune cell-relevant DEGs for NPM1 mutations were downregulated, despite their lack of association with immune scores. There were overlaps among the individual mutation-associated immune/stromal cell-relevant DEGs (Supplementary Fig. S2). The common upregulated stromal cell-relevant DEGs among RUNX1, TP53 and ASXL1 mutations included DDR2 and FRZB. The common downregulated immune cell-relevant DEGs between FLT3-ITD and NPM1 mutations included CD3D, CD48, GBP1, and IL18RAP. The common downregulated stromal cell-relevant DEGs between FLT3-ITD and NPM1 mutations included BGN, CDH5, COL1A2, COL6A3, CXCL12, DCN, FRZB, ISLR, ITIH5, MXRA5, TRAT1, and VCAM1 (n = 12). Furthermore, FRZB was the only DEG associated with all 5 mutations (RUNX1, TP53, ASXL1, FLT3-ITD and NPM1), ITIH5 was associated with RUNX1, ASXL1, FLT3-ITD and NPM1 mutations, and ISLR was associated with TP53, ASXL1, FLT3-ITD and NPM1 mutations.
The intermediate and poor cytogenetic risk patients were grouped by the median transcription levels of the individual mutation-associated immune/stromal cell-relevant DEGs to evaluate their effects on OS. The following genes were found to have prognostic significance, which was consistent with that of the corresponding mutation: high expression of SCUBE2 (RUNX1 mutation-associated upregulation) was shown to be related to lower OS, high expression of SPON1 (TP53 mutation-associated upregulation) was shown to be related to lower OS, high expression of GREM1 (ASXL1 mutation-associated upregulation) was shown to be related to lower OS, low expression of COL3A1, CXCL12, EMCN, FRZB, ITIH5, KDR, MXRA5, TART1, and VCAM1 (FLT3-ITD mutation-associated downregulation) were shown to be related to lower OS, high expression of ADAMTS5 (NPM1 mutation-associated upregulation) was shown to be related to higher OS, and low expression of SPON1 (NPM1 mutation-associated downregulation) was shown to be related to higher OS ( Supplementary Fig. S3).

Validation in the GEO database and identification of hub genes.
To verify whether these immune and stromal cell-relevant DEGs identified from TCGA AML patients are also associated with mutations in an independent AML cohort, we analyzed the gene expression levels of 524 AML cases from GSE14468, for which FLT3-ITD, NPM1 and biCEBPA mutation data were available. For immune score-involved genes, 7/7 of FLT3-ITD-associated, 13/13 of NPM1 mutation-associated and 3/3 biCEBPA mutation-associated DEGs were confirmed to be significantly related to the individual gene mutation by GSE14468. Similarly, for stromal score-involved genes, 7/22 FLT3-ITD-associated, 20/32 NPM1 mutation-associated and 2/2 biCEBPA mutationassociated DEGs were also confirmed. The confirmed mutation-associated immune and stromal cell-relevant DEGs are shown in Table 2.

Discussion
A number of genetic mutations have presented immune microenvironment modulatory properties in solid tumors: EGFR mutations correlate with an immunosuppressive TME and may impact the antitumor immune response in NSCLC 22,23 ; TP53 and KRAS mutations in lung adenocarcinoma can regulate the immune microenvironment to affect PD-1 blockade immunotherapy 24,25 ; JAK1 or JAK2 mutations may lead to acquired resistance to PD-1 blockade immunotherapy in patients with melanoma 26 . Recurrent genetic mutations found in AML have been heavily studied to classify and predict the risk of relapse after treatment. According to the ELN and NCCN guidelines, RUNX1, ASXL1, TP53, FLT3-ITD, NPM1 and biCEBPA mutations have been involved in AML prognostic stratification 2,3 . Interactions between leukemic stem cells and other cells in the BM microenvironment are known to be vital for the maintenance and progression of chemotherapy-resistant AML 14 . LSCs can remodel the BM niche into a favorable environment for expansion or even induce leukemic transformation. Nonetheless, the relationships between recurrent genetic mutations and the immune microenvironment in AML have not been comprehensively described 27 .
In our study, we calculated the immune and stromal scores for AML patients from the TCGA database based on the ESTIMATE algorithm. Our results showed that immune scores were significantly associated with OS and cytogenetic risk; a high immune score was a significantly poor prognostic factor for both the entire cohort and patients in the intermediate and poor cytogenetic risk groups, and a high stromal score tended to be correlated with poor OS in the intermediate and poor cytogenetic risk groups. The prognostic significance of immune and stromal scores was different in solid tumors: high immune and stromal scores correlated with poor survival in glioblastoma 18 , clear cell renal cell carcinoma (ccRCC) 28 and gastric cancer 16 , whereas high immune and stromal scores correlated with better survival in cervical squamous cell carcinoma (CSCC) 29 30 . These studies demonstrated the varied effects of immune and stromal scores on prognosis, and these effects are related to tumor type. Due to the similar OS rate, immune scores and stromal scores between the intermediate and poor cytogenetic risk groups, we considered these patients as a single group to explore the characteristics of the somatic mutation-associated immune microenvironment in AML. We compared the immune/stromal scores, identified DEGs between the mutation and WT groups, conducted functional enrichment analysis of DEGs and selected somatic mutation-associated immune/stromal cell-relevant DEGs. We found that similar to the impact of immune and stromal scores on prognosis, distinct relationships existed between somatic mutations and immune/ stromal scores. In other words, RUNX1, ASXL1 and TP53 mutations were related to higher immune or stromal scores, whereas FLT-ITD mutation was related to lower immune and stromal scores, although they were all poor  Table 2. FLT3-ITD, NPM1 and biCEBPA mutations-associated immune and stromal cell-relevant DEGs. Genes that displayed in italics were confirmed to be associated with mutations by GSE14468 data. *These mutation-associated genes have prognostic significance in the intermediate and poor cytogenetic risk patients.  www.nature.com/scientificreports/ prognostic mutations. Furthermore, patients with NPM1 mutation had lower stromal scores, while patients with biCEBPA mutation showed similar immune and stromal scores. despite their favorable prognostic risks. The functional enrichment analysis and immune/stromal cell-relevant DEGs were generally consistent with the immune/stromal scores for individual genetic mutations. There are unique and common genes among mutations associated with immune/stromal cell-relevant DEGs. The results obtained in this study demonstrated that RUNX1, TP53 and ASXL1 mutation-associated characteristics of the microenvironment are similar. Reports have revealed the pro-inflammatory impact of RUNX1, TP53 and ASXL1 mutations on the immune microenvironment. RUNX1 mutation has been shown to activate NF-κB signaling and has been proposed to promote inflammatory signaling pathways in the bone marrow microenvironment 31 . Previous reports have shown that TP53 mutations induce pro-inflammatory effects on epithelial cells through NF-κB-mediated production of inflammatory cytokines. Moreover, TP53 mutations in CAFs are associated with pro-tumor and pro-inflammatory effects through enhanced production of cytokines and chemokines, including CXCL12, SDF-1 and IL-6, which notably affect the immune microenvironment 32,33 . ASXL1 mutation is one of the most frequently observed mutations leading to clonal hematopoiesis (CH), which have been known to show elevated inflammation, impaired tumor suppressor function, and risk of eventual hematological malignancy (HM) [34][35][36] . Patients with NPM1 mutation and FLT3-ITD mutation not only had similar lower scores but also had multiple common immune/stromal cellrelevant DEGs, which were not consistent with their opposite prognostic significance. Our results indicated that there might be a common mechanism on the impact of NPM1 and FLT3-ITD mutations on the bone marrow microenvironment, which remains to be explored.

Status of DEGs Immune cell-relevant DEGs Stromal cell-relevant DEGs
Notably, several mutations associated with immune/stromal cell-relevant DEGs were observed to have prognostic significance in intermediate and poor cytogenetic risk patients. The results implied that these specific genes may play an important role in the formation of a mutation-associated microenvironment and may affect the survival of AML patients. Moreover, the majority of immune/stromal cell-relevant DEGs were confirmed to be significantly correlated with FLT3-ITD, NPM1, and biCEBPA mutations in the GSE14468 database. PPI networks were subsequently built based on the verified FLT3-ITD/NPM1-associated immune and stromal cellrelevant DEGs, and the top 10 hub genes were subsequently identified by the degree of interaction.
LCK (lymphocyte-specific protein tyrosine kinase) was the most significant hub gene associated with the FLT3-ITD mutation. LCK plays an essential role in the selection and maturation of developing T cells in the thymus, the activation of mature T-cells and the initiation of T cell antigen receptor (TCR) signal transduction pathways 37 . Studies have indicated higher expression of LCK in leukemic cells from less differentiated cases of AML (AML-0 and AML-1) 38 . A recent report found that LCK is overexpressed and mutated in CTV-1 cells (AML-M5 cell lines) 39 . Nonetheless, the expression of LCK in FLT3-ITD-mutated cells has not been studied to date. In the present study, the downregulation of LCK was correlated with the FLT3-ITD mutation. In a study of a zebrafish model 40 , FLT3 was found to initiate definitive hematopoietic stem cells, and the knockdown of FLT3 reduced hematopoiesis. The expression of the FLT3-ITD mutation resulted in the expansion of myeloid cells and the reduction of T cells. These results suggest that the FLT3-ITD mutation decreases the expression of LCK and reduces the production of functional T cells.
VCAM1 (vascular cell adhesion molecule-1) was shown to be the most significant hub gene associated with NPM1 mutation. VCAM1 is a cell adhesion molecule primarily expressed on endothelial cells, and its expression is induced by pro-inflammatory cytokines, such as TNFα 41,42 . VCAM1 has been identified to regulate vascular adhesion and transendothelial migration by binding to VLA-4 (very late antigen-4, an α4β1 integrin) on leukocytes. VCAM1 binding to VLA-4 confers AML blast cell protection from chemotherapy-induced apoptosis 43,44 . In our study, the downregulation of VCAM1 was confirmed to be correlated with NPM1 mutation. Although no study to date has explored the function of VCAM1 in NPM1-mutated AML patients, we speculated that the downregulation of VCAM1 may reduce the stroma-mediated protection of leukemic cells, which might confer favorable outcomes to AML patients with NPM1 mutations.
In conclusion, we focused on the relationship between recurrent genetic mutations and the immune microenvironment in AML patients based on TCGA database by integrated bioinformatic approaches. Important immune and stromal cell-relevant DEGs that affected the immune landscape of patients with individual gene mutations were identified and validated. Considering the specific properties of the hematopoietic microenvironment of leukemia 15 , ESTIMATE may not accurately predict infiltrating stromal and immune cells for the AML microenvironment, and we need to develop a more suitable and accurate algorithm. Due to the limited patient numbers in mutation subgroups in the TCGA database, further investigation of these mutation-associated stromal and immune signatures in large clinical AML patient cohorts is warranted, which may provide new prognostic biomarkers to achieve precision tumor therapy. Our results may help to elucidate how AML genetic mutations modulate the immune microenvironment to better guide personalized immunotherapy in the era of precision medicine 45 .

Materials and methods
Database. The transcriptional profiles and clinical and overall survival (OS) data of 173 AML patients were downloaded from the TCGA database (https ://porta l.gdc.cance r.gov/). The gene expression profile was measured experimentally using the Illumina HiSeq 2000 RNA Sequencing platform. Log2 transformations were performed for all gene expression data. Immune and stromal scores were calculated by applying the ESTIMATE algorithm 15 to the mRNA expression data (https ://bioin forma tics.mdand erson .org/estim ate/). The definitions of cytogenetic risk and risk-related somatic mutations were based on NCCN guidelines 2 .