Identification of aberrantly methylated and underexpressed genes
The clinical and genetic features of patients of GSE15061 and GSE58477 were described in the original article(9). In comparison with the control group, 9954 DEGs and 4536 underexpressed DEGs were identified in the AML group of GSE15061. A total of 1139 exonic hypermethylated DMRs were revealed in AML group of GSE58477. 198 shared genes between DEGs and DMRs were found by Venn analysis (the shared genes are list in Supplementary table 1).
Overrepresentation analysis of aberrantly methylated genes
The shared genes were revealed to be significantly enriched in the following biological processes (BP): antigen processing and presentation, positive regulation of NK cell cytokine product, etc. Cell component (CC) analysis demonstrated that shared genes were predominantly located in the MHC class I protein complex, recycling endosome membrane, etc. Molecular function (MF) analysis indicated that shared genes were enriched in peptide antigen binding, CD8 receptor binding, etc. Based on KEGG pathway analysis, the shared genes were significantly enriched in the following signaling pathways: endocytosis, cell adhesion molecules (CAMs), etc. In Reactome analysis, the shared genes were enriched in neutrophil degranulation, signaling by NOTCH3, etc. The detail results of overrepresentation analysis (ORA) were listed in Supplementary Table 2, meanwhile the top enriched pathways are shown in Figure 2.
Survival analysis
The assumption for the proportional hazards model was tested for each variable. The age variable did not meet the P-H model and a time-dependent Cox hazards model was used to adjust for age. In univariable analysis, the methylation level of 6 genes (out of 198 shared genes) significantly predicted unfavorable prognosis (Table 2): CORO1A (cg06038367, CpG shore), EHD1 (cg01151584, open sea), MPO (cg14619064, CpG island), SLC11A1 (cg07015784, CpG shore), SYNJ2 (cg10977795, CpG island), and GAS2L1 (cg04929173, CpG island). While expression level of MPO and SYNJ2 were also significantly associated with OS (Table 2). Notably, we summarized the expression profiles of the 6 genes in previous studies by Oncomine online tools (https://www.oncomine.org/resource/)(Figure 3), indicating that all 6 genes were predominantly underexpressed in leukemia patients compared with their expression in normal controls.
In multivariable analysis, age, race, mutation count, induction type, FAB subtype (M3 or non-M3), risk stratification by cytogenetics, risk stratification by molecular mutations, transplant or not, WBC count, methylation level (beta value) of 6 genes, expression level (M value) of CORO1A/MPO/SLC11A1/SYNJ2 were included, according to the results of univariable analysis. The stepwise forward linear ranking (LR) selection was used to screen independent prognostic factors to establish the prediction model. Finally, age, race, risk stratification of molecular mutations, induction type, transplant or not, WBC count, and methylation level of SYNJ2 (cg10977795) were significantly associated with OS (Table 2). The hazard ratio (HR) of methylation of SYNJ2 was 2.130 (95%CI 1.010-4.492) and the p value was 0.047, indicating hypermethylation of SYNJ2 is an independent unfavorable prognostic factor.
In the following subgroup analysis, the samples were classified as “HMA group” and “chemotherapy group” according to the determined induction type. The clinical features of the HMA and chemotherapy groups were compared and shown in Table 1. The patients receiving HMA had a significantly higher median age, lower ratio of AML-M3, more mutations, inferior risk stratification, lower transplant ratio and lower WBC count than the chemotherapy group.
For the HMA group, univariable/multiple analysis was shown in Table 3. Mutation count, risk stratification of molecular mutations and cytogenetics, expression and methylation level of EHD1 and methylation level of CORO1A were included in multivariable analysis, based on the results of univariable analysis. The Cox hazards model identified methylation of CORO1A as the only significant favorable factor (p value = 0.012, HR = 0.004), associating with longer OS. For the chemotherapy group, univariable analysis revealed that age, race, FAB subtype (M3 or non-M3), risk stratification of molecular mutations and cytogenetics, WBC, expression and methylation level of CORO1A/MPO/SYNJ2, methylation level of EHD1/SLC11A1/GAS2L1 were significant prognostic factors and included in multivariable analysis. Notably, the expression of CORO1A was reported to be downregulated in CEBPA-mutated AML by Federzoni et al(26). Therefore, we downloaded corresponding CEBPA mutation data from the TCGA database and used CEBPA mutation status as a stratifying factor in the following multivariable Cox hazards analysis. The age (p value = 0.008, HR = 1.004), risk stratification of molecular mutations (p < 0.05), transplant (p value = 0.012, HR = 0.512), WBC count (p value = 0.011, HR = 1.006), and SYNJ2 methylation (p value = 0.041, HR = 2.544) were assessed. Notably, the methylation of CORO1A was not a significant predictor for chemotherapy patients.
Moreover, we classified the samples as “CORO1A hypermethylated group” and “CORO1A hypomethylated group” inside the HMA and chemotherapy groups, in comparison with median beta value of each group. Kaplan-Meier plot displayed the opposite impact by methylation of CORO1A in the HMA and chemotherapy groups (Figure 4A&B). Therefore, CORO1A methylation level may serve as a biomarker for HMA treatment. The AML patients harboring hypermethylated CORO1A may get superior response and survival benefit in HMA treatment instead of chemotherapy.
Furthermore, to explore the association of CORO1A and HMA, GSE40871 dataset was employed, including expression and methylation data of primary AML samples with decitabine or mock treatment. The transcription of CORO1A in decitabine treated samples was upregulated compared with that in the control group (p value = 0.0484, log2foldchange = 0.668), while the expression in cytarabine treated samples didn’t vary compared with that in control group, indicating CORO1A may be a substantial target demethylated by HMA. While the methylation analysis found no methylation changes among target genes by comparing decitabine and mock treated group. Lacking direct methylome evidence, the relation of HMA and CORO1A still need validation.
Association between genome-wide methylation profile and SYNJ2 methylation
74 differentially methylated regions (DMRs) were identified to be associated with SYNJ2 methylation significantly (FDR-adjusted p value < 0.05, the list of DMRs are shown in Supplementary Table 3). The DMRs included the following vital genes: 1) HOXA9, an oncogenic transcription factor(27); 2) CALR, the mutation of which is often reported in myeloid neoplasms (MPN), but the methylation of which was rarely studied for AML; 3) ABI1, identified as a new partner gene for MLL in an AML patient(28) ;4) STAT3, involving in leukemogenesis and immune evasion in AML(29). All abovementioned genes were hypomethylated in SYNJ2 hypermethylated group. The DMRs are associated with SYNJ2 were enriched in Chemokine signaling pathway, Hematopoietic cell lineage, Graft-versus-host disease, etc.
Association between genome-wide gene/microRNA expression profile and SYNJ2 methylation
322 upregulated DEGs and 814 downregulated DEGs were identified, comparing SYNJ2 hyper- and hypo-methylated groups (adjusted p value < 0.05, |log2foldchange| >1, Figure 5). The heatmap of top SYNJ2 methylation associated DEGs is shown in Figure 6. Gene set enrichment analysis (GSEA) using MSigDB database identified 3 suppressed and 30 activated pathways, significantly correlated with SYNJ2 methylation level (q value < 0.05, Figure 7 & Table 4). The activated pathways included PI3K-Akt signaling, PD-L1 expression and PD-1 checkpoint pathway, JAK-STAT signaling, NF-kappaB signaling, etc. while the suppressed pathways included ascorbate and aldarate metabolism, Pentose and glucuronate interconversions, etc.
26 upregulated differentially expressed microRNAs (DEmiRs) and 75 downregulated DEmiRs were identified as significantly associated with SYNJ2 methylations (adjusted p value < 0.05, |log2foldchange| > 1). We selected the AML-related DEmiRs for further exploration, according to previously validated evidence by experiment in vitro or vivo. The low level of miR-217 expression was reported to be associated with unfavorable prognosis in AML(30). The miR-485-3p was reported to regulated expression of DNA topoisomerase and induce chemoresistance via targeting nuclear factor YB, in leukemic cells(31). The miR-889-3p may inhibit myeloid cell nuclear differentiation antigen expression, acting as potentially targets in AML(32). The potential target mRNAs were predicted by miRWalk 2.0 online tools. An integrated network of DEmiRs, target genes and DEGs was displayed in Figure 8.