The process of the present study is presented in Supplementary Figure S1 . Public genetic expression information and complete clinical annotations were retrieved from the GEO, TCGA, and GTEx databases. Patients with no survival data were excluded. Data from the eligible GSE71014 GEO dataset were used for deeper assays. For the microarray performed using this dataset, we acquired the normalized matrix data directly through download. For the TCGA dataset, RNA sequencing data of gene expression were acquired from the GDC via the R Package-Age TCGA biolinks (29), which was developed particularly for comprehensive assays of GDC data. Somatic mutation information was obtained from the TCGA database. Datasets were obtained for copy number variant (CNV) assays. The results were assayed using R 3.6.1 and R Bioconductor software.

Twenty-three modulators from the integrated GEO dataset were extracted to identify the different m6A modifications mediated by m6A regulators. These 23 m6A modulators included eight writers (METTL3, METTL14, METTL16, WTAP, VIRMA, ZC3H13, RBM15, RBM15B), two erasers (ALKBH5 and FTO), and 11 readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, HNRNPC, FMR1, LRPPRC, HNRNPA2B1, IGFBP1, IGFBP2, IGFBP3, RBMX). Based on the expression of 23 m6A regulators, an unsupervised cluster assay was employed to determine different m6A modification patterns and to categorize patients for more in-depth assays. ConsensusClusterPlus software was employed to implement the aforementioned procedures.

ssGSEA algorithms were used to quantify the comparative abundance of each infiltrative cell activity within the TME. The gene markers for each TME infiltrative immune cytotype were acquired from the research of Charoentong, which identified a variety of human immunocyte subgroups such as stimulated CD8 T cells, stimulated dendritic cells, macrophages, and NKT and modulatory T cells (30). The enrichment values computed via ssGSEA assay were used to calculate the comparative abundance of each TME infiltrative cell in each specimen.

To identify the m6A-associated genes, the patients were divided into three different m6A modification pattern groups based on their expression of the 23 m6A modulators. The experiential Bayesian method of the limma R software was employed to identify differently expressed genes (DEGs) between different modification patterns (31). Statistical significance was set at p < 0.001 for identifying DEGs. Functional enrichment analysis of DEGs was performed using clusterProfiler (32) R software to determine GO classes, such as BP, MF, and CC. A pathway enrichment assay using the KEGG database was also completed using this software. Statistical significance was set at p < 0.05.

To quantify the m6A modification pattern of a single tumor, we established a scoring system to assess the m6A modification features of each patient with AML, i.e., the m6A genetic signature. The DEGs detected in distinct m6A clusters were normalized, and overlapping genes were extracted. Patients were divided into groups for more in-depth analysis using unsupervised clustering methods to analyze the overlapping DEGs. Consensus clustering algorithms were used to define the quantity of genetic clusters and the relevant steadiness. We then completed a prognosis assay for each gene in the signature using a univariable Cox regression pattern. Genes with remarkable prognostic outcomes were extracted for deeper assays. Our team subsequently performed PCA to establish m6A-associated genetic signatures.

Plasma samples from 10 patients with AML and 10 healthy controls were collected from the People’s Hospital of Guizhou Province. Ethical approval was granted by the Human Research Ethical Board of GPPH. The mRNA expression levels of the four immune genes in the specimens were determined via qRT-PCR. Overall RNA from AML and normal control plasma specimens was prepared using Trizol reagent (Servicebio, Wuhan, China) as per the specification. RNA was then converted to cDNA by reverse transcription using the RevertAid First-Strand cDNA Preparation Kit (Thermo Fisher Scientific, Waltham, MA, USA). The genetic expression was normalized to that of GAPDH. FastStart Universal SYBR Green Master (Roche) was used to quantify the real-time PCR assay, using StepOne (Applied Biosystems, Waltham, MA, USA). The primer sequences are presented in Supplementary Table S1 . Each RNA specimen was prepared three times. To compare the expression levels of the different specimens, the comparative expression of genes associated with inflammatory responses was calculated using the 2^−ΔΔCt method.

The sections were dewaxed two times with xylene and gradients of different concentrations of anhydrous ethanol (100%, 95%, 90%, 80%, and 70% alcohol), placed into citrate buffer solution (0.01 mol/L, pH 6.0), and then boiled for 3 min in an autoclave. The sections were then sealed with 5% BSA and placed in an incubator at 37°C for 0.5 h. The slices were cultivated with the first antisubstance at 4°C nightlong. The sections were cultivated with a second antibody (horseradish peroxidase-labeled) for 60 min at room temperature and stained with DAB. Eventually, the sections were restained with hematoxylin, desiccated, and then fixed. The expression of PD-L1, PD-L2, MRP1, and MRP2 proteins in bone marrow tissue was detected by immunofluorescence. In the paraffin-embedded bone marrow pathology sections, the aforementioned identification process was completed according to the specifications provided by the supplier. The nuclei were stained with DAPI, and the relevant images were examined using an inverted fluorescence microscope.

One-way ANOVA and Kruskal-Wallis tests were used to compare the diversities among ≥3 groups (33). The cutoff points for each dataset subgroup were identified via the Survminer R package based on the association between m6A score and patient survival. Applying the “surv-cutpoint” function, all potential cutoff points were iteratively tested to identify the maximal rank statistic to dichotomize the m6A score. Patients were then separated into high and low m6A score groups based on the maximum chosen log rank statistic to reduce the batch effect of the computation. Survival curves for progression assays were produced using the KM approach, and log-rank tests were used to analyze the significance of differences. The waterfall function of the maftools package was employed to present the mutations in patients in the high and low m6A score groups from the TCGA-LAML cohort. The R package from RCircos was used to map CNVs in the 23 m6A modulators. All statistical p results were bilateral, and p < 0.05 was considered statistically significant. The entire processing was completed using R3.6.1 software.

We first summarized the incidence of CNVs and somatic mutations in 23 m6A regulators in AML, including eight writers, two erasers, and 13 readers ( Supplementary Table S2 ). Mutations in the m6A regulators occurred in two of the 134 samples, with a frequency of 1.49%. The assay of AML specimens found that only WTAP and RBM15 showed any mutation frequency, while the other genes did not ( Figure 1A ). Wild-type WTAP was more significantly expressed ( Supplementary Figure S2 ).

Investigation of CNVs showed prevalent CNV alterations in the 23 regulators, most of which were deletions in copy number, with RBM15B, YTHDC2, and METTL16 frequently demonstrating CNV deletion ( Figure 1B ). The site of CNV variation in the m6A regulators on the chromosome is shown in Figure 1C . To determine if these gene variants affect the expression of m6A modulators in patients with leukemia, we investigated the expression levels of regulator mRNA in healthy and tumor specimens and discovered that CNVs may be a major factor contributing to the disruption of m6A regulator expression. The expression of CNV-amplified m6A regulators YTHDF1, YTHDF2, YTHDF3, and FTO was remarkably higher in AML specimens than in healthy specimens ( Figure 1D ). Survival analysis of high and low m6A regulator expression is shown in Supplementary Figure S3 . The survival analysis demonstrates a remarkable level of inhomogeneity in the gene and expression of m6A modulators in healthy and AML samples, suggesting that imbalances in m6A regulator expression play a key role in the development of AML.

An overview of m6A modulator mutual effects, modulator association, and prognosis of AML patients is described via the m6A modulator net ( Figure 2A ). The R software of Consensus Cluster Plus was employed to categorize patients with distinct m6A modification patterns via quantification of the expression of the 23 m6A modulators, and three different modification patterns were ultimately determined via unsupervised clustering ( Figure 2B ), including 117 patients with pattern A, 89 patients with pattern B and 25 patients with pattern C. We named these patterns m6A clusters A–C, respectively ( Supplementary Table S3 ). The prognosis assay for the three primary m6A modification subgroups revealed specific survival benefits for the m6A cluster C modification pattern ( Figure 2C ). To explore the biological activities of the different m6A modification patterns, a KEGG enrichment assay was performed, and it was found that the significant enrichment of m6A cluster A in stromal and carcinogenesis stimulation pathways such as T-cell receptor, B-cell receptor, and JAK-STAT signaling pathways, cell apoptosis, and NK-mediated cell toxicity ( Figure 2D ). m6A cluster B displayed enrichment of pathways related to pyruvate, fructose, mannose, glycerolipid, and galactose metabolic activities, as well as glycosaminoglycan degradation ( Figure 2E ). m6A cluster C was significantly associated with B-cell receptor, chemokine, cell attachment, MTOR signaling pathways, and immune suppression biological processes ( Figure 2F ). We found that the three m6A modification patterns exhibited distinctly different TME infiltrative cell profiles. Cluster A was categorized as an immunorejection phenotype characterized by congenital immunocyte infiltrative activity and stromal stimulation. Cluster B was categorized as immune metabolism, regulating immune activation through the metabolic programs of different cell subpopulations, and cluster C was characterized by immunosuppression.

The following assay of TME infiltrative cell results revealed that m6A cluster A had significant enrichment of congenital immunocyte infiltrative activity, including NK, macrophages, eosinophils, mastocytes, MDSCs, and PDCs ( Figure 3A ). However, significant differences occurred in the m6A transcriptional profile among the three distinct m6A modification patterns ( Figure 3B ). To explore the underlying biological behavior of each m6A modification pattern, the limma package was used to identify 31 DEGs associated with m6A phenotypes ( Figure 3C and Supplementary Table S4 ). Cluster Profiler software was employed to complete the GO enrichment assay of the DEGs ( Supplementary Table S5 ). The biological processes with significant enrichment were denoted as KEGG pathways ( Supplementary Table S6 ). Interestingly, the genes that were enriched were significantly associated with m6A modification and immune activity, which confirms that m6A modification has a nonnegligible effect on the immunoregulation of the TME ( Figures 3D, E ).

To explore the effects of these m6A modification phenotypes on diverse clinical characteristics and biological behaviors, unsupervised clustering of the AML dataset was performed, and different patterns of m6A modification was found in the cohorts ( Supplementary Figure S4 ). An unsupervised cluster analysis was then completed based on the 31 m6A phenotype-associated genes, which was then used to classify patients into different subgroups. Consistent with the cluster group categorization of the m6A modification patterns, the unsupervised cluster algorithm also determined three different m6A modification phenotypes, which we termed m6A gene clusters A–C ( Figure 4A ). This suggests that three different patterns of m6A methylation exist in AML cancers. Of the 233 AML patients, 72 were clustered into the B phenotype, which was associated with a poorer prognosis. In contrast, the patients in cluster C (34 patients) displayed better prognostic outcomes ( Figure 4B ). In these three clusters, obvious differences in m6A regulator expression were discovered, which is consistent with the predicted outcomes of the m6A methylated modification patterns ( Figure 4C ). The obvious differences in m6A gene expression across these three m6A genetic clusters have a nonnegligible regulatory effect on the formation of distinct TME conditions. Nevertheless, these assays are only based on AML patient populations and cannot precisely forecast the m6A methylation modification patterns of individual patients. Given the inhomogeneity and intricacy of m6A modifications, based on these phenotypically associated genes, we established a scoring system to quantify the pattern of m6A modifications in each AML patient, which we termed the m6A score. Alluvial plots were used to visualize the changes in various parameters of a single patient ( Figure 4D ). Genetic cluster A displayed the lowest median and the highest survival rate, while cluster C had the highest median score and increased mortality. The patients with greater m6A scores exhibited significantly improved survival ( Figure 4E ). To better characterize m6A, known characteristics and the m6A score were correlated ( Figure 4F ). The Kruskal-Wallis test showed significant differences in m6A scores among the m6A gene clusters. Genetic cluster A had the lowest median, whereas genetic cluster C displayed the greatest median, suggesting that a low m6A score may be strongly associated with immune stimulation-associated traits, while a high m6A score may be associated with stromal stimulation-associated traits ( Figure 4G ). In addition, m6A score was significantly higher in m6A cluster C and lowest in m6A cluster A ( Figure 4H ). These results suggest that low m6A scores are negatively correlated with immune stimulation, and high m6A scores are also negatively correlated with stromal stimulation. The m6A score can be used to evaluate the m6A modification pattern of individual tumors and further assess the TME infiltrative cell characteristics of tumors.

Subsequently, we determined the m6A score and tumor mutation load and how this was related to patient prognosis. Patients were separated into high or low m6A score groups ( Figure 5A and Supplementary Table S7 ). m6Ascore and TMB showed a positive correlation ( Figure 5B ). Significantly improved survival existed in the high mutational load group, and there was an even more pronounced survival advantage in patients with high mutational load and high m6A score ( Figures 5C, D ). These results indicate that the m6A score could also be utilized to assess the clinical prognosis of patients. Our team subsequently analyzed the differences in the distribution of certain mutations between the high and low m6A score groups from the TCGA database using the maftools package. As shown in Figures 5E, F and Supplementary Table S8 , the low m6A score group displayed a broader tumor mutational load compared with the high m6A score group, with a ratio of 80.65% and 66.7% mutated genes, respectively. Accumulating evidence suggests that patients with high TMB status exhibit persistent clinical responses to anti-PD-1/PD-L1 immune therapy. Thus, the differences in tumor m6A modification patterns may be key factors mediating clinical reactions to anti-PD-1/PD-L1 immunotherapies. These results will provide new perspectives to explore the mechanisms of m6A methylation patterns in tumor somatic mutations, the formation of TME status, and their role in ICB therapies.

Immunotherapy, comprising PD-L1 and PD-1 blockade, has become an important breakthrough in tumor treatment. We tested whether m6A modifications could forecast patient reactions to ICB treatment and discovered that patients with high m6A scores demonstrated significant clinical benefits and significantly prolonged survival ( Figures 6A, B ). To further elucidate the correlation between m6A-modified gene expression patterns and host antitumor immune responses, we assessed the expression of immunoregulatory genes in m6A score patients and found that PD-L1, PD-L2, MRP1, and MRP2 expression was significantly higher in patients with low m6A scores, in contrast to patients with high m6A scores, suggesting a possible response to antitumor immune responses ( Figures 6C–F ). Furthermore, we evaluated the proportion of different genes in the peripheral blood of healthy controls (NC), and patients with AML, identified by qRT-PCR. As shown in Figure 7A , compared with the control group, the levels of PD-L1, PD-L2, MRP1, and MRP2 in the peripheral blood of the patients with AML was significantly increased (p < 0.01). Moreover, IHC staining demonstrated the same pattern ( Figure 7B ), suggesting that m6A modifications were crucial for the antitumor immune response in AML. In conclusion, the present research suggests that m6A methylation patterns are closely associated with tumor immunophenotypes and antitumor immune responses, which will help to provide new therapeutic guidance for antitumor immunotherapies.