The RNA-seq data in transcripts per million (TPM) format and relevant clinical information from the TCGA AML cohort were obtained from the TCGA data portal (https://portal.gdc.cancer.gov/projects) (Colaprico et al., 2016). A total of 173 patients with AML and 70 normal people (The GTEx Consortium atlas of, 2020) were included in the TGCA LAML program, and Log2 transformation of the TPM RNA-seq data was performed before further intrasample comparison. 153 AML patients were eligible and extracted for survival analysis after excluding patients without clinical survival outcome.

Bone barrow mononuclear cells (BMNC) were Isolated from six AML patients and three healthy donors by Histopaque®-10771 (Sigma-Aldrich) standard procedure. All participants provided written informed consent. The study was approved by ethics review boards of our institution. The clinical information of enrolled patients was listed in Supplementary Table S1. Total RNA of BMNC was isolated with TRIzol reagent (Life Technologies) and then reverse transcribed to complementary DNA (cDNA) using the Evo M-MLV RT Kit (AG) with gDNA Clean. The cDNAs were then mixed with SYBR reagent, Applied Biosystems real-time PCR system were used for gene expression analysis, and the expression level were normalized to GAPDH. The primer sequences for UNC93B1 and GAPDH are shown in Supplementary Table S2.

Patients were divided into two groups based on the median cutoff value of UNC93B1 expression, and the power of UNC93B1 expression to distinguish AML from healthy individuals was estimated using the receiver operating characteristic (ROC) curve drawn by the pROC software (Schlicker et al., 2007). Various alternative clinical events were analyzed using univariate Cox analysis to evaluate their potential impact on overall survival (OS) in AML patients. Furthermore, a multivariate Cox regression model with independent prognostic factors was applied to construct the final prognostic predictors. Independent factors were then recruited to build the nomogram prognostic model (Motwani et al., 2018). Additionally, calibration (Karapanagiotis et al., 2018) was used to estimate the predictive power of the nomogram model. The nomogram and calibration curve were constructed using the “RMS” R package. Finally, the time-dependent receiver operating characteristic (ROC) analysis (Krischer et al., 2019) was performed to assess the predicted value of UNC93B1 in the prognostic model.

Based on the median expression of UNC93B1 mRNA, AML patients were divided into two groups (UNC93B1-lower and UNC93B1-higher groups). The differentially expressed genes (DEGs) between the high-UNC93B1 and low-UNC93B1 groups were identified by DESeq2R (Love et al., 2014) package. DEGs were considered here as an adjusted p value < 0.05 and an absolute logarithmic 2-fold change (|log2 FC|) >2.0.

Functional Gene Ontology annotations were performed using R-package clusterProfiler (Yu et al., 2012), including biological process (BP), cellular component (CC), and molecular function (MF) (Gene Ontology Consortium, 2021), as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis (Kanehisa et al., 2022). To further evaluate the gene function, we also conducted the Gene Set Enrichment Analysis (GSEA) analysis (Szklarczyk et al., 2019) using hallmark gene sets (h.all.v7.2. symbols.gmt) from the MsigDB and it was defined as a statistically significant item with the p value was less than 0.05.

Protein–protein interaction (PPI) networks can discover and annotate the biological interactions between cellular proteins. In this study, potential protein–protein interactions were predicted by the website of STRING (https://string-db.org/) (Szklarczyk et al., 2019). The PPI for up-regulated DEGs were constructed and visualized by Cytoscape (version 3.7.1) (Shannon et al., 2003). Then hub genes were found using three different algorithms. Finally, the ggplot2 package was used to map the correlation of UNC93B1 and the expression of the hub genes in AML.

To explore the association between UNC93B1 expression levels and abundance of infiltrating immune cells, we used GSVA (version 1.34.0) to perform the ssGSEA algorithm and identified the qualified correlation of infiltrated immune cells and UNC93B1 expression by Pearson correlation coefficient (Bindea et al., 2013; Hänzelmann et al., 2013). A total of 24 subtypes of infiltrating immune cells were involved, including dendritic cells (DCs), activated dendritic cells (aDCs), immature dendritic cells (iDCs), plasmacytoid dendritic cells (pDCs), B cells, T cells, CD8 positive (CD8⁺) T cells, T helper cells, T central memory (Tcm) cells, T effector memory (Tem) cells, T follicular helper (Tfh) cells, T gamma delta (Tgd) cells, Th1 cells, Th17 cells, Th2 cells, regulatory T (Treg) cells, cytotoxic cells, eosinophils, macrophages, mast cells, neutrophils, NK cells, NK CD56 (bright) cells, and NK CD56 (dim) cells. Next, we draw a scatterplot in R by ggplot2 to show a quantitative correlation of UNC93B1-expressed level and well-recognized surface markers of specific immune cells (CD14, CD68), which were widely discussed in monocytes/macrophages.

Statistical analysis was carried out based on the R package (https://www.r-project.org/, version 3.6.3). Statistically significant was defined as a two-tailed p-value of <0.05. Differences in clinical features based on UNC93B1 expression were defined by the Kruskal–Wallis tests and the Wilcoxon rank sum tests. The Pearson χ2 test was used to detect differences in clinical features between low and high UNC93B1-expressed AML. The survival curve was plotted using the Kaplan-Meier (KM) method, and the logarithmic rank p-value < 0.05 was considered statistically significant. Univariate and multivariate Cox regression analysis were selected to determine prognostic factors. Risk ratios (HR) and 95% confidence intervals (CI) for survival-related genes or events were calculated.

To explore the expression landscape and potential indication of UNC93B1, the mRNA levels of this molecular in tumors and normal tissues were extracted from TCGA database. Pan-cancer analysis revealed that UNC93B1 was widely and highly expressed in most malignancies (Figure 1A), especially in AML (p < 0.001, Figure 1B). We validated the high expression level of UNC93B1 in another GEO datasets (GSE13159, p < 0.001, Figure 1C).

Intriguingly, UNC93B1 expression distinguished AML from healthy individuals with a predictive power of 0.998 (95% confidence interval, CI = 0.994–1.000), as shown by AUC values analyzed by the ROC curve (Figure 1D). Therefore, we doubted whether UNC93B1 expression in AML patients is clinically relevant. We surprisingly observed UNC93B1 expression is linked to final OS events (p < 0.01, Figure 1E) and other clinical characteristics (Table 1). The UNC93B1 expression level was higher in AML patients with older age (age > 60, p < 0.01, Figure 1F), non-M3 FAB classification (p < 0.05, Supplementary Figure S1A), especially M4/M5 FAB subtype (p < 0.001, Figure 1G), in which monocytes were prominent. Unexpectedly, unfavorable cytogenetic risk including intermediate or poor cytogenetic risk, which indicated inferior survival in 2017 European Leukemia Net (ELN) genetic risk stratification had higher level of UNC93B1 (p < 0.05, Figure1H). Moreover, among the favorable cytogenetic risk, upregulated UNC93B1 was associated with inv (16) (p < 0.01, Supplementary Figure S1B), which results in the fusion oncoproteins CBFβ-SMMHC. As we all known, AML subtype M4-with eosinophilia (M4-eo) is classified by AML blasts with inv (16) that have a myelomonocytes with immature basophilic granules (Sangle and Perkins, 2011), which is consistent with our results that UNC93B1 is higher in M4/M5 FAB subtypes leukemia. Furthermore, we performed quantitative real-time PCR (qRT-PCR) of our clinical samples to confirm that UNC93B1 were high-expressed in AML-M4/M5 FAB subtypes (p < 0.005, Figure 1I). Whereas, there is no association with gender, WBC counts, FLT3, IDH1, RAS and NPM1 mutation status (Table1). Taken together, high UNC93B1 expression was closely linked to unfavorable clinical characteristics and deserved further exploration in AML.

To explore the prognostic value of UNC93B1 in AML, we draw the overall survival (OS) curve using the Kaplan-Meier method. Patients with higher levels of UNC93B1 than the median expression value had a worse overall survival (HR = 2.08, p = 0.001, Figure 2A). Furthermore, we rationally proposed UNC93B1 level can also distinctly reflect survival outcomes of different subgroups, and the analysis showed higher UNC93B1 level had a more pronounced shortens of survival in AML patients with older age (age > 60, HR: 10.98, p = 0.021, Figure 2B), non-M3 FAB subtype (HR: 1.90, p = 0.004, Figure 2C), WBC counts in peripheral blood ≤ 20 x10^9/L (HR: 2.76, p = 0.001, Figure 2D), and unfavorable cytogenetic risk (including Intermediate cytogenetic risk and Poor cytogenetic risk, HR = 2.49, p = 0.001, Figure 2E). Additionally, higher expression of UNC93B1 also a poor survival factor in AML patients with normal karyotype (HR: 2.29, p = 0.008, Figure 2F), which is predominant in intermediate-risk group (Lin et al., 2017). Interestingly, we found high expression of UNC93B1 just linked to worse prognosis in AML patients without FLT3 (p < 0.05, Supplementary Figure S2A), IDH1 (p < 0.05, Supplementary Figures S2B–D) or NRAS (p < 0.05, Supplementary Figure S2H) mutation, once acquired mutations of these genes, UNC93B1 expression level had no effects on patient’s survival outcome (p > 0.05, Supplementary Figures S2E–G). Meanwhile, there is a significantly negative correlation between UNC93B1 level and OS of both subgroups divided by BM blasts and NPM1 (p < 0.05, Supplementary Figures S2I–L).

To further identify the indicated prognosis of UNC93B1 in AML and other survival-related factors, we performed univariate and multivariate Cox regression analyses (Supplementary Table S3). Univariate analysis indicated that older age (age > 60, HR: 3.33, 95% CI: 2.16–5.13, p < 0.001), unfavorable cytogenetic risk (including Intermediate cytogenetic risk and Poor cytogenetic risk, HR: 3.21, 95% CI: 1.65–6.24, p < 0.001), non-M3 FAB subtype (HR: 3.93, 95% CI: 1.43–10.81, p = 0.008) and higher UNC93B1 levels (HR: 2.08, 95% CI: 1.35–3.20, p < 0.001) predicted poor survival (Figure 2G). Multivariate Cox regression analyses showed that except for age > 60 (HR: 2.78, 95% CI: 1.77–4.38, p < 0.001), expression level of UNC93B1 (HR: 2.18, 95% CI: 1.38–3.42, p < 0.001) was the only factor which have clinically meaningful effects to overall survival (Figure 2H). Thus, UNC93B1 is an independent prognostic factor found in both univariate and multivariate analysis.

Additionally, a nomogram model was built using age, FAB subtype, cytogenetic risk and UNC93B1 gene expression, which were independent prognostic factors in univariable analysis, to provide the model’s quantitative predictive power for survival outcomes in AML patients (Figure 2I). The predictive OS concordance index (C-index) was 0.729 (95% CI: 0.701–0.756, p < 0.001, Figure 2I), indicating that these four univariable factors had a certain predictive accuracy for OS. Calibration analysis was performed to further assess the predictive value of the nomogram model, in which the calibration curves showed the agreement between the predicted OS and the actual proportion of OS at different time points, and the 5-year median survival rate is sharply decreased compared to 2-year median survival rate (Figure 2J), probably in part due to the higher relapse rate of AML within 2 years clinically. Thus, this nomogram model is well calibrated and had an accuracy to predict the patients’ clinical outcomes. Finally, time-dependent analysis of sensitivity, specificity, and ROC curve was constructed to evaluate UNC93B1 for predicting survival outcome at 2 years or 5 years, and the cumulative 2-year survival rate is 43.7% (2-year AUC = 0.70, 95% CI: 0.60–0.80), while the cumulative 2-year survival rate was dramatically decreased to 21.9% (5-year AUC = 0.76, 95% CI:0.59–0.94, Figure 2K), which is consistent with clinical data. In summary, the survival-predictive model constructed here is significant for AML and UNC93B1 served as an independent prognostic factor for poorer outcome.

To explore the potential mechanisms of UNC93B1 in leukemogenesis, we first identified the DEGs between two AML populations divided by high and low UNC93B1 expression. In total, 614 DEGs (|log2(FC)| >2 and adjusted p value < 0.05) were obtained and shown in volcano plots (Figure 3A; Supplementary Table S4), including 342 upregulated genes and 272 downregulated genes. To elucidate the intracellular signaling pathways influenced by the DEGs, we performed GO and KEGG analysis. Firstly, we investigated the biological function of the up-regulated DEGs, the top five of BP enrichment items, CC enrichment items, MF enrichment items as well as top 5 KEGG pathways were obtained from GO annotation (Figure 3B; Supplementary Table S5). We can find two main signaling sets were enriched, one set is immune-related signaling, especially the innate immune associated pathways, such as antigen processing and presentation, macrophage activation, pattern recognition receptor activity, Phagosome and Toll-like receptor signaling pathway were significantly enriched (Figure 3B; Figure 3D; Supplementary Table S5). Another notable signaling cascade was metabolism-related set, including NADPH oxidase complex, carbohydrate binding, lipopeptide binding, superoxide-generating NADPH oxidase activity (Figure 3B; Figure 3E; Supplementary Table S5). Moreover, to better understand the mutual connection among genes and the aforementioned two main signaling sets, we performed an interactive analysis of two main enriched signaling sets and its corresponding genes (Figure 3D-E). We discovered that CD14, TLR8, TLR7, ITGAM, NCF2, TREM2, CYBB, CLEC7A, CD209, NCF1, SLC11A1 and MARCO shared by at least two pathways in innate immune-related sets (Figure 3D); while NCF2, CYBB, NCF1B, NCF1C and NCF1 were commonly enriched in metabolism-related sets (Figure 3E).

To validate the activation of innate immune signaling in UNC93B1-high expressed AML, we next use GSEA analysis to investigate the enrichment pathway, and finally, we found four significantly activated pathways, which were innate immune system (NES = 2.67, p = 0.01, Figure 3C), immunoregulatory interactions between a lymphoid and a non-lymphoid cell (NES = 2.60, p = 0.01, Supplementary Figures S3A), neutrophil degranulation (NES = 2.23, p = 0.01, Supplementary Figures S3B) and adaptive immune system (NES = 2.94, p = 0.01, Supplementary Figure S3C; Supplementary Table S6). So, we assumed that immune system, especially innate immune signaling might play an important role in AML with high level of UNC93B1.

To determine possible proteins directly interact with UNC93B1, protein-protein interaction (PPI) analyses were used. We calculated all of the 342 up-regulated genes and the top 20 hub genes were obtained by the maximal clique centrality (MCC, Supplementary Figures S4A), maximum neighborhood component (MNC, Supplementary Figures S4B), and density of maximum neighborhood component (DMNC, Supplementary Figures S4C) algorithms respectively. Clearly, four common hub genes (S100A9, CCR1, MRC1 and CD1C) shared from the aforementioned three gene lists (Figure 4A), and all of these four genes are related to innate immune response. To gain further deep insight in these molecules, we next turn to interrogate the expression and survival outcome of these four genes. All hub genes were upregulated in AML patients compared to healthy individuals (p < 0.05, Figure 4B). Furthermore, increased expression of S100A9, CCR1 and MRC1 were also associated with poor survival in AML (p < 0.05, Figure 4E). Additionally, UNC93B1 expression positively correlates with S100A9, CCR1, MRC1 and CD1C (p < 0.001, Figure 4C), and the detailed co-expression pattern were shown in Figure 4D.

To investigate the association between UNC93B1 expression and AML-infiltrating immune cells, we studied a total of 24 immune cell subtypes. The results displayed that AML with high expression of UNC93B1 had a significant trend to be infiltrated by dendritic cells (DC), immature dendritic cells (iDC), Macrophages, Neutrophils, Eosinophils, NK CD56(dim) cells, Effective Memory T Cell (Tem), Th17 cells and Treg cells (Figure 5A), most of which were responsible for innate immune system instead of adaptive immune. The refined details of a quantified analysis of UNC93B1 and infiltrated innate immune cells (DC, iDC, macrophages, neutrophils, Eosinophils and NK CD56dim cells) by Spearman’s correlation are shown in Figure 5B. More strikingly, we found the level of UNC93B1 expression was positively corelated with macrophages (r = 0.572, p < 0.001), and regulatory T cells (Tregs, r = 0.41, p < 0.001) in AML (Figure 5E), which were responsible for suppressing functional T cell activation and regulating immune systems. We also discovered UNC93B1 expression levels positively correlated with CD14 (r = 0.672, p < 0.001, Figure 5C), a surface antigen that is preferentially expressed on monocytes and a molecule associated to shorter overall survival and lower complete remission rate in AML (Solary et al., 1992; Bradstock et al., 1994). CD68, generally represented macrophages and associated with a lower complete remission in AML (Liu et al., 2022), were also found to be higher in UNC93B1-high expressed AML patients (r = 0.618, p < 0.001, Figure 5C), which further validated the role of innate immune cells in AML pathogenesis and therapy response. TLRs, the most important family of receptors in the early to middle stages of innate immunity response, play a crucial role in recognizing microorganisms and augmenting the inflammatory mediators (Shey et al., 2010; Takeuchi and Akira, 2010; Trejo-de la et al., 2014; Schnetzke et al., 2015), were all extremely upregulated in AML (p < 0.05, Figure 5D), consistent with dysfunctional innate immune response in AML. Unsurprisingly, the level of UNC93B1 was positively and significantly correlated with CD14, CD68, and almost all of TLRs except for TLR3 (p < 0.001, Figure 5C; Figure 5F).

To refine the potential mechanism of UNC93B1 on leukemogenesis, we quantified the correlation of expression between UNC93B1 and the genes related to innate immune response-activating signal transduction (Figure 6B) and cytokine metabolic process (Figure 6C), both were significantly enriched in GO/KEGG (Figure 3B). Gene lists related to each pathway were shown in Figure 6A, and CLEC7A, TLR7, TLR8 played important roles in both signaling. C-Type Lectin Domain Family 7 Member A (CLEC7A), encoding a glycoprotein that is a small type II membrane receptor with extracellular C-type lectin-like domain folding and a cytoplasmic domain with an activated motif based on immune receptor tyrosine. Gene Ontology (GO) annotations associated with this gene include carbohydrate binding and pattern recognition receptor activity, which is consistent with the mutually shared genes in aforementioned immune and metabolic pathway. What’s more, TLR7 and TLR8 predominantly expressed in dendritic cells, TLR8 is also expressed in monocytes, macrophages, neutrophils, and Tregs (Gorden et al., 2005; Peng et al., 2005). TLR7 activation triggers the production of interferon-α (Gorden et al., 2005), and TLR8 activation leads to the production of immunomodulatory cytokines, such as interleukin-12, interleukin-18, and pro-inflammatory cytokines, such as tumor necrosis factor alpha and interleukin-1β (Cristescu et al., 2018), which is associated with leukemic microenvironment and immune tolerance (Pyzer et al., 2017). Therefore, the high expression of TLRs on leukemia cells and the consequences of TLR activation might be a potentially therapeutic target.

There have been multiple studies that indicate innate immune cells undergo substantial metabolic reprogramming following activation (Angiari et al., 2020). Meanwhile, various metabolic enzymes involved in glycolysis and mitochondrial metabolic pathways have been found to play a vital role in influencing innate immune cell function (Li et al., 2019). So, we supposed the coupling of metabolic process and innate immune responses play an essential role in the AML development, especially in the patients with high-expressed UNC93B1. And we draw the correlation heat maps of the key regulators in metabolic process and innate immune pathways (Figure 6D), which demonstrated obvious co-expression patterns. Therefore, we conclude that immunometabolism has taken on a different meaning in AML by unknown immune-metabolism interactions, for example, contributions of key metabolic factors contribute to immune cell development or fate; in turn, responses of immune cells resulted in metabolic reprogramming. Altogether, metabolic system and innate immune affected AML patients with aberrantly high-expressed UNC93B1 through combined strategies.

As we show in this work, UNC93B1 is high-expressed in AML, especially in M4/M5 (Figure 1G). Recent studies have shown that AML hierarchy composition as a determinant of response to targeted therapy (Zeng et al., 2022), and Monocytic AML (M5) is intrinsically resistant to venetoclax + azacytidine (Pei et al., 2020). We then further analyzed the expression of UNC93B1 in AML and found that the level of UNC93B1 is significantly higher in FAB-M5 patients than any other FAB-subtype patient (p < 0.05, Figure 7A-B), which is consistent with the positive correlation expression of UNC93B1 and CD14, CD68 shown in Figure 5C. As we demonstrated in Figure 3B, another signaling enriched by UNC93B1 were metabolic pathways, and AML blasts might switch from BCL2-dependence to MCL1-dependence to drive energy metabolism under specific circumstances (Pei et al., 2020), we next focused on expression of BCL2 and MCL-1 in AML. Contrary to the increasing expression of UNC93B1 in FAB-M5 AML, BCL2 expression decreased remarkably in FAB-M5 AML (p < 0.05, Figures 7C–D). However, MCL1 presented similar expression pattern to UNC93B1 in AML, tended to be high-expressed in FAB-M5 (p < 0.05, Figures 7E–F). To clear the quantitative expression correlation of UNC93B1, BCL2 and MCL1, we performed Spearman Rank Correlation Analysis, and as expected, we revealed significantly negative correlation of UNC93B1 expression and BCL2 expression (r = −0.345, p < 0.001, Figure 7G), meanwhile, consistent with the previous results, MCL-1 expression is positively related to UNC93B1 level in AML (r = 0.347, p < 0.001, Figure 7H). Given that venetoclax (VEN) is an inhibitor specific to BCL2 and some studies have discovered that the level of BCL2 expression is closely related to VEN response in vitro (Souers et al., 2013; Pan et al., 2014), we hypothesized high-UNC93B1 AML might be resistant to venetoclax, and selective MCL1 inhibitor, such as VU661013 (Ramsey et al., 2018), alone or in combination with AZA might be effective in venetoclax-resistant AML and high-UNC93B1 AML. Together, these data indicated that high-UNC93B1 expressed AML cells shows greater dependence on MCL1 than BCL2 for survival, and UNC93B1 might be a biomarker to predicting drug response in acute myeloid leukemia.