The expression and regulation of HOX genes and membrane proteins among different cytogenetic groups of acute myeloid leukemia.

BACKGROUND: The cytogenetic aberrations were considered as markers for diagnosis and prognosis in acute myeloid leukemia (AML), while the expression and regulation under different cytogenetic groups remain to be fully elucidated.
METHODS: In this paper, for favorable, poor, and cytogenetically normal groups of AML patients, we performed comprehensive bioinformatics analyses including identifying differentially expressed genes (DEGs) and microRNAs (miRNAs) among them, functional enrichment and regulatory networks.
RESULTS: We found that DEGs were enriched in membrane-related processes. Eleven genes and two miRNAs were significantly differentially expressed among these three AML groups. In survival analysis, membrane-related genes and several miRNAs were significant on prognostic outcome. Notably, six HOXA and three HOXB genes were significantly in low expression and high methylation in AML with favorable cytogenetics. Meanwhile, the miRNA-HOX gene co-regulatory networks revealed that HOXA5 was a hub node and regulated an AML oncogene SPARC.
CONCLUSION: Our work may provide novel insights to the molecular characteristics and classification between AML with different cytogenetics.

INTRODUCTION

Acute myeloid leukemia (AML) is a heterogeneous disease characterized by expansion of undifferentiated myeloid precursors, resulting in impaired hematopoiesis (Döhner, Weisdorf, & Bloomfield, 2015). The prognostic relevance of cytogenetic abnormalities has led to the widespread adoption of risk stratification for AML (Slovak et al., 2000). Overall, the 5‐year survival rate varies drastically based on the cytogenetic risk classification: 55%, 24%, and 5% overall survival (OS) for favorable (PML‐RARA, RUNX1‐RUNX1T1, or MYH11‐CBFB fusions), intermediate, and poor risk group patients (monosomy karyotype or complex alterations), respectively (Islam, Mohamed, & Assenov, 2017). However, nearly half of AML patients are cytogenetically normal (CN‐AML) with an intermediate prognosis (Slovak et al., 2000). As an important prognostic factor, the karyotype may involve in regulating gene expression to affect prognosis. Researchers examined the gene expression profiles driven by chromosomal lesions (e.g., t(8; 21), t(15; 17), and inv (16)) in AML and concluded that those cases with favorable cytogenetics are predictable with high accuracy by available genome‐wide gene expression technology (Verhaak et al., 2009). Recently, Mughal et al. (2016) identified differential expression of MYC protein related to cytogenetic risk groups in AML patients. Nonetheless, the key expressed molecules and their regulatory mechanisms between CN‐AML and other AML cytogenetics is unclear, which is an important issue to be revealed.

In recent years, microRNAs (miRNAs) and transcription factors (TFs) have drawn extensive attention in cancer research and their dysregulations have gradually been confirmed to be important in hematopoiesis (Undi, Kandi, & Gutti, 2013). Many miRNAs such as miR‐15/16, miR‐17‐92, and miR‐155 were reported as important regulators with essential functions in differentiation of hematopoietic stem cells, as well as the occurrence of leukemia (Musilova & Mraz, 2015). The aberrations of TFs (e.g., Homeobox genes) have been found to contribute to the disease progression of leukemia (Abramovich & Humphries, 2005). Furthermore, TFs and miRNAs are able to co‐regulate the expression of targets in forms of feed‐forward loops (FFLs) and feedback loops (FBLs), which are vital and common regulatory motifs in leukemia. Zhang et al. identified that miR‐146b‐5p as a key regulator accelerated the transformation by targeting NUMB and other genes in CML (Zhang et al., 2016). Pulikkan et al. (2010) confirmed CEBPA can regulate miR‐223, which in turn targets E2F1 forming a FFL in AML. Through network analysis, the complex regulatory relationships in leukemia will be illuminated on systematic level and the key regulators may be identified.

Homeobox (HOX) genes are TFs playing crucial roles in embryonic development and differentiation of hematopoietic cells (Spencer et al., 2015) and are organized by four clusters: HOXA (7p15), HOXB (17q21), HOXC (12q13), and HOXD (2q31) (Rice & Licht, 2007). To date, the dysregulation of HOX genes in AML has been attributed to specific chromosomal abnormalities involving mixed‐lineage leukemia (MLL) (Rice & Licht, 2007). Specifically, HOXA9 and HOXD13 are dysregulated through the t(7;11) and t(2;11) translocations, respectively (Nakamura et al., 1996). However, the different expression and function of HOX genes in various AML cytogenetics has not yet been completely elucidated.

In 2013, a systematic study of genomic landscape on AML through The Cancer Genome Atlas (TCGA) project provides a comprehensive framework of important data for AML clinical and research (Cancer Genome Atlas Research, 2013). Emerging evidence suggests miRNAs and gene expression signatures can predict outcomes of leukemia by analyzing TCGA data (Chuang et al., 2015; De Leeuw et al., 2016; Zhu et al., 2017). In addition, Li, Liang, and Zhang (2014) developed an effective model to integrate TCGA AML data and related TF binding data to reveal AML‐specific regulation. Here, we analyzed the difference of expression and regulation among AML with different cytogenetic groups and identified HOX genes and membrane‐related genes as key variance. Our work provides insights to the molecular difference of different AML cytogenetic groups and will be useful to AML classification and prognosis prediction.

METHODS

Patient's clinical information and data description

To analyze the AML data, we downloaded the normalized miRNA and mRNA expression data of AML patients from the TCGA data portal (https://tcga‐data.nci.nih.gov/tcga/) and microarray data set GSE6891 for validation. Clinical information (e.g., sex, age at diagnosis, cytogenetics, FAB, and OS) of AML patients was listed in Table S1. Based on the cytogenetic risk stratification (Cancer Genome Atlas Research, 2013), the TCGA cohort of 161 AML patients was further divided into three subgroups: (a) favorable cytogenetics (t(15; 17), t(8; 21) or Inv(16), n = 34); (b) CN‐AML (cytogenetically normal, n = 79); (c) poor cytogenetics (monosomy, complex karyotype, n = 48). We excluded 17 AML patients marked with intermediate cytogenetics but not CN‐AML (Table 1).

Table 1

Main clinical information of AML patients in different cytogenetic groups

Classification	Cytogenetic stratification (n = 161)
	Favorable (n = 34)	Intermediate (n=79)	Poor (n = 48)
M0	0	4	9
M1	4	23	10
M2	6	20	9
M3	16	0	0
M4	8	19	6
M5	0	11	9
M6	0	0	1
M7	0	1	2
Unknown	0	1	1
Median age	50	57	73
Median OS (month)	46.2	17	9.2
WBC ≥50 × 109/L	4	30	6
WBC <50 × 109/L	30	49	41

Differential expression analysis

To identify the differentially expressed genes (DEGs) and miRNAs (DEMs) between different AML subgroups by cytogenetic risk stratification in the TCGA cohort, we used the R‐package “NOISeq” (Tarazona, García, Ferrer, Dopazo, & Conesa, 2012), and thresholds for DEGs and DEMs were FDR ≤ 0.05 and |FC| >1.5. We used the R‐package “limma” (Phipson, Lee, Majewski, Alexander, & Smyth, 2016) with FDR ≤ 0.01 and |FC| >2 for the microarray data set GSE6891. Gene/miRNA expression patterns were indicated by heatmaps, which were performed using R packages such as “gglot2” and “ComplexHeatmap.” The functional annotation and KEGG pathway analysis of DEGs screened in different AML subgroups were analyzed using the Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool. Pathways and GO terms with p ≤ .01 based on Fisher's Exact Test were considered as statistical significance. DNA methylation data were downloaded from the TCGA portal and the analysis was performed by in‐house R scripts used in GSCALite (http://bioinfo.life.hust.edu.cn/web/GSCALite/) (Liu et al., 2018).

Construction of FFL and FBL in AML subgroup

The TF list was extract from Animal TFDB2 (Hu et al., 2019). We applied the method described in our FFLtool webserver (http://bioinfo.life.hust.edu.cn/FFLtool#!/) (Xie et al., 2020) to merge the miRNA and TF targets by both predicted and experimentally verified targets. Thus, the miRNA‐gene/TF and TF‐gene/miRNA regulatory relations were clearly clarified for DEGs and DEMs for different AML subgroups. Then, we constructed the miRNA‐TF‐gene FFLs and miRNA‐TF feed back loops using our in‐house scripts based on their regulations. These networks were visualized using Cytoscape (version 3.5).

Statistical analysis and survival analysis

We subdivided AML patients into two subgroups based on median expression values of the specific gene or miRNA. The Kaplan–Meier method with log‐rank test was applied to assess OS among patients with high expression or low expression levels of each gene or miRNA by GraphPad Prism (www.graphpad.com.). A one‐way analysis of variance (ANOVA) was performed to determine subgroup differences with Tukey tests used for comparisons between groups. p value ≤ .05 was considered as significance. Principle component analysis (PCA) analysis was used to distinguish samples based gene abundance. Partial least squares discriminant analysis (PLS‐DA) was employed to calculate the Variable Importance in the Projection (VIP) score and screen out the important genes which contribute the most to distinguish the groups (http://www.metaboanalyst.ca/) (version 3.0) (Xia, Sinelnikov, Han, & Wishart, 2015).

RESULTS

Differences of membrane‐related genes and HOX genes among AML patients of different cytogenetic groups

We performed differential gene expression analysis among AML samples of favorable, normal, and poor cytogenetics groups, and obtained a list of DEGs and DEMs for each comparison (Figure 1a). Notably, there were much fewer DEGs and DEMs in poor versus CN‐AML comparison than the other two comparisons (Figure 1a). Besides, we found there were 11 DEGs (SPARC, FAM105A, IFITM1, LASS4, CD99, SLC2A5, C10orf54, NINJ1, CAT, HOXA5, and HOXA7) (Figure S1a) and two DEMs (miR‐99b‐5p and miR‐500a‐3p) overlapped in these three pairwise comparisons (Figure S1b), suggesting their importance to AML prognosis. Among them, survival analysis indicated that SLC2A5, C10orf54, and HOXA5 showed significance for predicting prognostic outcome in the cohorts of TCGA (n = 159) and GSE6891 (n = 519) data set (Figure S1c,d). We also calculated the ratio of top 30 mutated genes in three cytogenetic groups and found the gene mutation ratio in favorable cytogenetic group was generally low except for the FLT3 mutation (24%). Besides, compared to favorable cytogenetics and poor cytogenetics, the mutation ratio was higher in CN‐AML and the top three genes were NPM1 (51%), FLT3 (38%), and DNMT3A (37%) (Figure S2). In CN‐AML, expression levels of most HOX genes in FLT3 + (n = 30) or NPM1 + (n = 40) samples were higher than samples with wide type of these two genes (p < .05) (Figure 2a,b), while there was no significant difference in other karyotypes.

Figure 1

Summary of DEGs and DEMs in comparisons of different cytogenetic groups of AML. (a) The number and ratio for up/down‐regulated DEGs and DEMs in comparisons of F versus C, P versus C, and F versus P in AML. F: favorable cytogenetics, P: poor cytogenetics, C: CN‐AML; (b) and (c) Venn diagrams of the numbers of DEGs (b) and DEMs (c) in pairwise comparisons; (d): Functional enrichment of merged DEGs in three pairwise comparisons of AML. EnrichRatio represents the proportion of the enriched genes accounting for the GO terms or KEGG pathways. AML, acute myeloid leukemia; CN, cytogenetically normal; DEGs, differentially expressed genes; DEMs, differentially expressed miRNAs

Figure 2

The expression of HOX genes in CN‐AML cohort between (a) FLT3 + mutation and FLT3 − mutation and (b) NPM1 + mutation and NPM1 − mutation (p < .05). AML, acute myeloid leukemia

To explore the functions of these DEGs, functional enrichment analysis of the union set of DEGs in three pairwise comparisons were carried out. Results indicated that most DEGs were involved in membrane‐related terms including integral component of membrane, cell surface, focal adhesion etc. (Figure 1d). Next, we investigated the expressions of top ten upregulated and downregulated DEGs, and top five upregulated and downregulated DEMs in each comparison (Figure 3a,b). Among these DEGs, 11 of them were membrane protein genes. Obviously, nine HOX genes were markedly decreased in favorable cytogenetics (Figure 3a) compared to other cytogenetics, which was consistent with the previous report (Alharbi, Pettengell, Pandha, & Morgan, 2013). As for DEMs, compared with other cytogenetics, let‐7b‐3p, miR‐196b‐5p and miR‐10a‐5p were downregulated in favorable cytogenetics (Figure 3b), which were reported to be positively correlated with the expression of HOXA and HOXB genes (Debernardi et al., 2007).

Figure 3

The expression heatmaps of top DEGs (A) and DEMs (B) in AML based on cytogenetic classification. F: favorable cytogenetics, P: poor cytogenetics, C: CN‐AML. The blue to red in “Expression Z score” indicates the expression level from low to high. The gene names with asterisk represent the HOX genes. DEGs in orange are the membrane‐related genes. The check mark indicates the DEGs and DEMs significantly in which comparisons. Genes have carried out hierarchical clustering by calculating average distances. AML, acute myeloid leukemia; CN, cytogenetically normal; DEGs, differentially expressed genes; DEMs, differentially expressed miRNAs

HOX genes stood out in differential expression, survival, and regulatory network analyses in cytogenetic group comparisons

As the dysregulation of HOX genes was stood out in differential expression analysis and was important in leukemia (Drabkin et al., 2002), we focused on them for further investigation. Totally, we obtained nine, six, and four differentially expressed HOX genes in the comparison of favorable cytogenetics versus CN‐AML, poor cytogenetics versus CN‐AML and favorable versus poor cytogenetics, respectively. Meanwhile, the nine and four differentially expressed HOX genes in the comparison of favorable cytogenetics versus CN‐AML and favorable versus poor cytogenetics all showed significance for predicting prognostic outcome, respectively (Figure 4). Moreover, six HOXA and three HOXB genes showed significantly lower expression and higher methylation in AML with favorable cytogenetics (p < .05) (Figure 5). We also calculated the VIP score and found that nine differentially expressed HOX genes showed high contribution rates in the separation of AML with different cytogenetics (Figure S3).

Figure 4

The overall survival of differentially expressed HOX genes in the comparison of (a) favorable cytogenetics versus CN‐AML and (b) favorable versus poor genetics group (p < .05). AML, acute myeloid leukemia; CN, cytogenetically normal

Figure 5

Expression and methylation rate of HOX genes in AML based on cytogenetic classification. The boxplots and scatter plots represent the expression and methylation rate of HOX family, respectively. AML, acute myeloid leukemia

To explore the HOX family and their regulations with miRNAs, we carried out differentially expressed HOX genes, DEMs and constructed the miRNA‐HOX gene co‐regulatory network in each comparison as described in Methods (Figure S4). Of note, the network has more nodes and edges in the favorable cytogenetics versus CN‐AML, implying the regulation was more complex. Moreover, all networks shared two HOX genes (HOXA5 and HOXA7), indicating their abnormal expression has a crucial role in genetic rearrangement. HOXA5 as a hub in these networks with the most nodes and edges, which suggests it plays key roles in these groups. Meanwhile, HOXA5 was predicted to regulate SPARC and the HOXA5‐miR‐10a‐5p‐SPARC formed a FFL in the comparisons of favorable cytogenetics versus CN‐AML and poor cytogenetics versus CN‐AML. In addition, HOXA5 and miR‐92a‐3p could co‐regulate SPARC when compares favorable cytogenetics versus poor cytogenetics, indicating the importance in leukemia.

Membrane‐ and prognosis‐related genes enriched and regulated by HOX genes in comparison of favorable cytogenetics and CN‐AML

As the comparison of favorable cytogenetics and CN‐AML was notable (Figure S4), we focused on this comparison for further investigation. By performing gene/miRNA expression analysis, we identified 443 DEGs (221 upregulated and 222 downregulated) and 38 DEMs (22 upregulated and 16 downregulated) (Figure 1a). Functional enrichment analysis showed that 110 upregulated DEGs were enriched in membrane‐related processes and 21 downregulated DEGs were enriched in immune processes (p < .05) (Figure 6a). We also used another independent data set (GSE6891, n = 344) to validate the result and found that the DEGs in the comparison of favorable cytogenetics and CN‐AML group were also highly enriched in membrane‐related biological processes (Figure S5). Next, we performed OS analysis to assess the association between the expression of DEGs/DEMs and prognostic significance in the TCGA cohort of AML with favorable cytogenetics and CN‐AML. As a result, we obtained 20 DEMs and 174 DEGs that showed significance for predicting prognostic outcome (p < .05) (defined as OS‐related DEGs/DEMs). The top ten genes and top five miRNAs ranked by p‐value were shown in Figure 7a,b. Meanwhile, seven genes (MPO, PLXNC1, LSP1, RHOBTB3, PDIA6, CCND2, and CCND3) were validated as the OS‐related DEGs in another AML data set (GSE6891) (Figure S6). Genes such as CCND2 and CCND3 were reported to play oncogenic roles in multiple malignancies (Deshpande, Sicinski, & Hinds, 2005; Liu et al., 2018). Among the top five DEMs, the overexpression of miR‐196b‐5p and let‐7a‐3p was associated with poor outcome in AML patients (Li et al., 2013; Li, Huang, Chen, et al., 2012; Li, Huang, Li, et al., 2012).

Figure 6

The functional enrichment and regulatory network for DEGs in AML comparison of favorable cytogenetics versus CN‐AML. (a) GO enrichment of all DEGs in comparison of favorable cytogenetics versus CN‐AML. Immune response genes downregulated and membrane‐related genes upregulated. (b) The HOX‐miRNA‐gene co‐regulatory network of immune and membrane‐related terms. Round Rectangles: miRNAs. Circles: target genes. Genes in green were downregulated and red ones were upregulated. AML, acute myeloid leukemia; CN, cytogenetically normal; DEGs, differentially expressed genes; miRNAs, microRNAs

Figure 7

Top five DEMs (a) and top ten DEGs (b) with OS significance in validated data set survival analysis. DEGs, differentially expressed genes; DEMs, differentially expressed miRNAs; OS, overall survival

To further explore the potential roles of HOX family in biological process, we constructed the HOX‐miRNA‐gene co‐regulatory network which contained six differentially expressed HOX genes (HOXA3, HOXA5, HOXA7, HOXA9, HOXA10, and HOXB3) and OS‐related DEGs/DEMs involved in GO/KEGG pathway of membrane and immune response (Figure 6b). The network contained 128 edges and 62 nodes including 30 membrane‐related genes and six immune response‐related genes. Compared with other five HOX genes, HOXA5 regulates most of membrane‐related gene targets and some of them may function as tumor suppressors, such as DGKG (Kai et al., 2017) and SPARC (DiMartino et al., 2006), which showed implication in proliferation and invasion of cancer cell. Meanwhile, the subnetwork centering on HOXA9 revealed its key regulatory role because it showed interactions with major miRNAs and several gene targets (TPM4, CNN2, and PDIA6). Moreover, the top four miRNAs regulating most gene targets were miR‐335‐5p, miR‐10a‐5p, miR‐181a‐5p and miR‐181b‐5p, which may act as the important regulatory factors and serve as biomarkers for the prognostic prediction of AML patients (Havelange et al., 2014; Super, 2015). The immune‐related gene, PRKCD, which is involved in the PRKCD‐P38‐C/EBPα pathway of AML, was demonstrated to be directly targeted by miR‐181a‐5p (Su et al., 2015). Our results were not only in accordance with these reports, but also revealed some synergistic regulations.

DISCUSSION

A range of studies reported the cytogenetic aberrations were mostly considered as markers for diagnosis, prognosis, and classification in AML. However, the regulatory interactions between dysregulated molecules in AML based on different clinical classification remain fully elucidated. In this study, using TCGA data, we performed comprehensive bioinformatics analysis, and validated in microarray data set GSE6891 among AML patients based on cytogenetic risk classification and identified several dysregulated DEGs/DEMs. Functional enrichment demonstrated that most DEGs were involved in membrane‐related GO terms and many of them were significant in survival analysis. We focused on HOX family genes in the favorable cytogenetics versus CN‐AML group to constructed miRNA and TF/genes regulatory networks, and identified the key regulators.

Specific HOX genes can be disrupted through chromosomal translocation. In this study, nine HOX genes were downregulated in the AML with favorable cytogenetics (Figure 2a) compared to other cytogenetics and showed high contribution rates in the separation of AML with different cytogenetics. Our results were consistent with the reports that the highly expression of HOXA genes has a relatively shorter OS in AML with intermediate or poor cytogenetic risk (Li, Huang, Chen, et al., 2012; Li, Huang, Li, et al., 2012; Xia et al., 2015). Moreover, previously study suggests that the regulation of HOX genes can provide a conceptual framework for the actions of some oncogenic TFs (Look, 1997). As a key regulator of myeloid differentiation, loss function of HOXA5 may limit leukemia associated with specific chromosomal translocations (Boucherat, Guillou, Aubin, & Jeannotte, 2009). Our study showed that the expression tendency of HOXA5 was closely correlated with HOXA9, which calls for further research to reveal the underlying mechanism of the interaction.

The functional enrichment showed most DEGs were involved in membrane‐related processes including integral component of membrane, cell surface, focal adhesion, etc. (Figures 1d and 6a). Although leukemic cells experience less anchorage than solid tumor cells, they do grow and reside within the bone marrow surrounded by a large network of microenvironmental cells (Bajaj et al., 2016). Membrane‐related proteins represent essential components of the leukemia microenvironment that regulate cell survival by interacting with the extracellular matrix (Gabarra‐Niecko, Schaller, & Dunty, 2003), which is important for leukemia stem cells’ biological behavior (Ayala, Dewar, Kieran, & Kalluri, 2009). Collectively, the membrane‐related terms might be involved in the progression and worse prognosis of AML. Recent studies have shown that regulations by miRNAs and TFs are tightly co‐regulated, which may work as the “core” of the whole gene regulatory network.

CONCLUSIONS

In conclusion, we explored the difference of miRNA/mRNA expression and regulatory networks associated with prognosis based on the cytogenetic classification, and identified key genes including HOX genes and membrane‐related genes in the favorable cytogenetics versus CN‐AML group. The miRNA‐HOX gene co‐regulatory networks revealed that HOXA5 was a hub node in three comparisons and it regulated a shared target SPARC. Our work may provide novel insights to the molecular characteristics between AML with different cytogenetics and will be useful for a better understanding of AML classification and prognosis.

CONFLICT OF INTEREST

The authors claim no conflicts of interest.

AUTHOR’S CONTRIBUTIONS

H. Wang and H. Hu designed the research; H. Wang, S.Y. Lin, H. Hu, and F.F. Hu analyzed and interpreted data; S.Y. Lin, H. Wang, and A.Y. Guo wrote and/or revised the manuscript.

Table S1‐Fig S1‐S6

Click here for additional data file.