The expression matrices and corresponding clinical characteristics of LUAD samples were obtained from the Gene Expression Omnibus (GEO) and The Cancer Genome Atlas (TCGA) databases. We excluded patients without survival, information, or incomplete clinicopathological characteristics from further assessment. A total of 936 patients were enrolled, including the cohorts GSE68465 (N = 438) and TCGA-LUAD (N = 498). The four validation databases were downloaded from the GEO database, including GSE11969, GSE13213, GSE37745, and GSE50081. The expression matrix data of the TCGA-LUAD cohort (FPKM format) were downloaded from the Genomic Data Commons platform, and FPKM units were converted to transcripts per kilobase million (TPM) units. The “Normalized Between Arrays” function of the R package “Limma” was performed for data standardization. Genome mutation data of TCGA-LUAD (including somatic mutation and copy number variation (CNV)) were downloaded from the UCSC Xena platform. Based on previous studies, we collected 24 m⁶A regulators, including 10 writers (CBLL1, ELAVL1, METTL3, METTL14, METTL16, KIAA1429, RBM15, RBM15B, WTA, and ZC3H13), two erasers (ALKBH5 and FTO), and 12 readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, FMR1, HNRNPA2B1, HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3, LRPPRC). Clinicopathological information and clinical immunotherapy scores (IPS) were obtained from the TCIA database.

Based on mRNA expression levels, 19 m⁶A regulators were extracted from the TCGA-LUAD and GSE68465 cohorts, and the samples were divided into diverse subtypes based on transcriptome data with the R package “Consensus Cluster Plus.” A thousand repetitions were performed to guarantee the stability of the classification (44).

To investigate the difference between m⁶A modification patterns in the biological process, we explored the variation in signaling pathways between each of the two subtypes of m⁶A regulators by using “Gene set variation analysis (GSVA)” R packages (45). We downloaded the gene set file “c2.cp.kegg.v7.4.symbols” from the MSigDB database for the GSVA analysis. An adjusted p-value of less than 0.05 was considered statistically significant.

To estimate the relative abundance of 23 immune cell types in the tumor microenvironment of LUAD, single-sample gene set enrichment analysis (ssGSEA) was used to calculate the enrichment scores and represent the relative abundance of each tumor-infiltrating immune cell type in each sample. The set of genes used to label each tumor-infiltrating immune cell type was obtained from the study by Charoentong (46).

The patients were divided into three subtypes according to the expression level of 19 m⁶A regulators. We identified differentially expressed genes (DEGs) among different subtypes using the empirical Bayesian method in the R package “Limma,” and selected p-values less than 0.001 as DEG candidates for further analysis (47).

To quantify the level of m⁶A modification in a single tumor, we established an m⁶A score model by performing principal component analysis (PCA). The procedure for the established m⁶A score model was as follows: first, overlapping DEGs were identified from different m⁶A clusters, and significant prognosis-related genes were identified by univariate Cox regression analysis; second, PCA was performed on the gene expression profile, and the principal components 1 and 2 were extracted as feature scores. This method minimized the deprivation of information contained in the original index and reduced the indicators to be analyzed, thereby allowing a comprehensive analysis of the collected data; lastly, the m⁶A score was defined by performing a formula similar to that used in previous studies (48, 49). m⁶A score = ∑ (PC1i + PC2i), where i is the expression of m⁶A phenotype-related genes.

All statistical analyses were performed with R software (version 4.0.3). The R package “Limma” was used for differential gene expression analysis. Spearman’s correlation analysis was used to calculate the correlation coefficients between the levels of different tumor-infiltrating immune cell types and the expression of m⁶A regulators. Kruskal-Wallis test and one-way analysis of variance were used to compare differences between more than two groups. The Kaplan-Meier (KM) method was used to draw the survival curve (5-year survival rate), and univariate Cox regression analysis was used to calculate the hazard ratios (HRs) for m⁶A regulators and m⁶A phenotype-related genes. LUAD samples were divided into high and low m⁶A score subgroups using the “surv-cutpoint” function in the R package “survival.”

Twenty-four m⁶A RNA methylation regulators were selected for LUAD according to previous studies, including 10 writers, 12 readers, and two erasers ( Table S1 ). We summarized the dynamic reversible epigenetic modification behavior of these m⁶A regulators and their biological functions in RNA, including mRNA export, mRNA translation, mRNA decay, and mRNA degradation/stability ( Figure 1A ). Somatic mutations in m⁶A regulators were found in 151 (26.63%) of 567 samples. ZC3H13 had the highest mutation frequency, followed by KIAA1429 and IGF2BP1, while no CBLL1 mutations were found in the samples ( Figure 1B ). In addition, YTHDC1 was significantly positively correlated with ZC3H13, YTHDC2, FMR1, and HNRNPA2B1 ( Supplementary Figure S1B ). There was widespread CNV in the 24 regulatory factors; METTL16, RBM15B, METTL14, ELAVL1, and RBM15 had copy number losses, while YTHDF1, KIAA1429, FMR1, IGF2BP2, and METTL3 showed numerous copy number gains (gene amplification) ( Figure 1C ). The locations of the CNVs in the 24 m⁶A regulators were labeled on the chromosomes ( Figure 1D ). Tumor samples were distinguished from normal samples by three-dimensional PCA (3D-PCA) of the 24 m⁶A regulators. The results showed that the two groups were completely separated from each other ( Figure 1E ). We then compared mRNA expression in normal and tumor samples to explore whether the expression levels of the m⁶A regulators were affected by the above gene variation. Seventeen of the 24 m⁶A regulators showed significant overexpression or downregulation in LUAD samples. YTHDF1, KIAA1429, HNRNPC, and METTL3 were most significantly upregulated in tumor samples, and METTL16 and METTL14 were markedly downregulated ( Figure 1F ). These results showed that the CNV was an important factor in controlling the expression of m⁶A regulators. Most m⁶A regulators underwent remarkable expression changes in LUAD, suggesting that the abnormal status of m⁶A regulators is involved in the development of LUAD.

Meta-analysis was performed using TCGA-LUAD and GEO (GSE68465) datasets to further elucidate the modification characteristics of the m⁶A regulators in LUAD. This algorithm finally identified 19 m⁶A regulators through data consolidation. To examine whether the 19 m⁶A regulators could be used as prognostic markers for LUAD, we used univariate Cox regression analysis ( Supplementary Figures S1C and S2A–S ). Eight of the 19 genes were significantly correlated with prognosis, including HNRNPA2B1, HNRNPC, IGF2BP2, IGF2BP3, LRPPRC, RMB15, WTAP, and ZC3H13 (HR >1). The prognostic significance of 19 m⁶A regulators in LUAD, the connection to the regulator, and the interactions are shown in the m⁶A regulator network ( Figure 2A ). The results showed a remarkable correlation among writers, readers, and erasers, and this crosstalk was essential for the generation of different m⁶A modification modes. The landscape also suggested that the occurrence and progression of LUAD are related to m⁶A regulators. Effective immune infiltration in tumors is considered a key factor in carcinogenesis and prognosis (50, 51). The correlations between individual regulators and each tumor-infiltrating immune cell type were analyzed using Spearman’s correlation analysis ( Figure 2B ). All m⁶A regulators were significantly correlated with some types of tumor-infiltrating immune cells, suggesting that m⁶A regulators are critical for immune infiltration in tumors. Most m⁶A regulators were significantly positively correlated with immune infiltration, while WTAP, an m⁶A methyltransferase, was negatively correlated.

Based on the expression of the 19 m⁶A regulators, the R package “Consensus Cluster Plus” classified patients into various m⁶A modification patterns. Consequently, the unsupervised consensus clustering algorithm revealed that patients were well defined when k = 3 ( Supplementary Figure S3A ). Thus, three distinct patient clusters were identified based on m⁶A modification patterns, including 253 cases of cluster A, 324 cases of cluster B, and 359 cases of cluster C. The Kaplan-Meier analysis of the three clusters revealed that cluster A had the worst prognosis, while cluster B had a significant survival advantage ( Figure 2C ). Next, the expression levels of the m⁶A regulators in the three clusters were analyzed; the expression profiles of most m⁶A regulators significantly differed among the three clusters. Cluster A was upregulated compared with the other two clusters, including CBLL1, ELAVL1, FMR1, HNRNPA2B1, HNRNPC, IGF2BP2, IGF2BP3, LRPPRC, METTL3, RBM15, RBM15B, WTAP, YTHDF1, and ZC3H13 ( Supplementary Figure S3B ). The expression levels of IGF2BP2 and IGF2BP3 were significantly downregulated in cluster B but upregulated in cluster A ( Figure 2D ). These results indicate that IGF2BP2 and IGF2BP3 are the dominant risk factors for malignant progression. Based on the KEGG gene set, GSVA enrichment analysis was used to explore the biological processes of the two m⁶A modification patterns ( Figure 2E ; Supplementary Figures S3C, D ). We found that cluster A was clearly enriched in the process of DNA repair, such as nucleotide excision repair, mismatch repair, DNA replication, and base excision repair; cluster B was prominently enriched in metabolism-associated pathways, such as fatty acid metabolism, amino acid metabolism, arachidonic metabolism, drug metabolism cytochrome P450, and primary bile acid biosynthesis. Cluster C was enriched in pathways associated with immune activity, including intestinal immune network production, antigen processing and presentation, and cytokine receptor interaction.

Subsequently, the correlation between 23 tumor-infiltrating immune cell types and three m⁶A cluster subsets was examined by ssGSEA analysis. The results showed that cluster C was associated with more adaptive and innate immune cell infiltration, including CD8 cells, NK cells, Tregs, MDSCs, macrophages, and B cells ( Figure 3A ). Consistently, cluster C exhibited a comprehensively elevated expression of MHC molecules ( Figure 3B ). Immune cell infiltration in the three m⁶A cluster subsets was assessed using the ESTIMATE algorithm. The results showed that cluster C exhibited higher immune and stromal scores, indicating that cluster C had significantly increased immune cell infiltration ( Figures 3C, D ), which was mainly due to increased immune infiltration and high MHC expression.

To further elucidate the underlying biological processes and functional annotation of a single m⁶A modification pattern, seventy-three m⁶A phenotype-related DEGs were identified using the R package “Limma” and represented on a Venn diagram ( Figure 4A ). GO enrichment analysis for DEGs was performed using the R package “Cluster Profiler.” Surprisingly, these genes were remarkably related to nuclear division and organelle fission ( Supplementary Figure S3E ). To further validate this regulatory mechanism, we obtained 68 prognosis-related DEGs using univariate Cox regression analysis. Unsupervised cluster analysis was performed on the 68 prognosis-related DEGs, which were divided into three subgroups: gene clusters a, b, and c. We performed PCA on the expression of the 68 DEGs, showing that the three groups were completely separated from each other ( Figure 4B ). The results showed 226, 193, and 311 patients in gene clusters a, b, and c, respectively. Patients in gene cluster a were associated with a worse prognosis, patients in gene cluster b had a better prognosis, and patients in gene cluster c had a moderate prognosis ( Figure 4C ). m⁶A regulator expression was significantly different between the three gene clusters, which was the same as that of the m⁶A modification patterns ( Figure 4D ).

To predict the immune status and prognosis in a single patient, we sought to develop an m⁶A score based on the 68 DEGs identified. Therefore, patients were divided into high- or low-m⁶A score groups based on the cutoff value. The high m⁶A score group had a better clinical survival profile ( Figure 5A ). In addition, we externally verified m⁶A score model to predict the prognosis of patients with lung adenocarcinoma from GSE11969, GSE13213, GSE37745, and GSE50081 datasets ( Supplementary Figures S4A–D ). An alluvial diagram was used to illustrate the workflow of the m⁶A score construction and to visualize the attribute changes in individual patients ( Figure 5B ). The results indicated that gene cluster b was associated with a high m⁶A score, whereas gene cluster c was associated with a lower m⁶A score. Notably, most patients who were still alive were included in the high m⁶A score group. Most patients with m⁶A cluster A were defined as low m⁶A, while patients with m⁶A cluster B had a high m⁶A score. The m⁶A cluster C was similarly distributed ( Figure 5C ). The group with a statistically low m⁶A score had more advanced patients ( Figure 5D ). We also examined the relationship between the above subtypes and the m⁶A score. Kruskal−Wallis analysis revealed that m⁶A cluster B and gene cluster b exhibited the highest m⁶A score, while m⁶A cluster A and gene cluster a showed the lowest score ( Figures 5E, F ). Multimodel crossvalidation suggested that the m⁶A score could serve as a prediction model for LUAD. Based on Spearman’s analysis, a heat map demonstrated a positive correlation between the m⁶A score and tumor-infiltrating immune cells ( Figure 5G ). Furthermore, the high m⁶A score group also exhibited a comprehensively elevated expression of MHC molecules ( Figure 5H ). These results indicated that the m⁶A score may be used as a model to predict the immune status in LUAD.

Several studies have confirmed correlations between tumor somatic mutations, genomic alterations, and immunotherapeutic effects. Therefore, we evaluated the distribution of the tumor mutation burden (TBM) in different m⁶A score groups. A box and scatter diagram showed that the high m⁶A score group had a lower TMB and that the m⁶A score was negatively correlated with TMB ( Figures 6A, B ). K-M survival analysis showed that TMB alone was not sufficient to accurately predict the prognosis of LUAD ( Figure 6C ). To further understand the relationship between TMB, m⁶A score, and survival outcomes, a K-M survival analysis based on the combination of m⁶A score and TMB was performed. The results revealed that the low TMB and low m⁶A score subgroups were associated with poor prognosis. The combination of the m⁶A score and TMB could accurately predict the quality of life of patients with LUAD ( Figure 6D ). These data emphasize the impact of the m⁶A score and TMB on cancer development. A list of the somatic mutations in the high- and low-m⁶A score groups is shown ( Figures 6E, F ). The low m⁶A score group presented more extensive TMB than the high m⁶A score group, except for KRAS (22% vs. 31%).

We examined the relationship between the m⁶A score and clinicopathological characteristics, and observed that the m⁶A score significantly differed between patients by stage, node (N), but not tumor (T) and metastasis (M) ( Figures 7A–D ). In this study, the prognostic value of the m⁶A score was determined for different clinicopathological characteristics, and it was found that patients with high m⁶A score had significantly longer overall survival than patients with low m⁶A score in stages I/II, T1/2, N0, and M0 ( Figures 7E–L ). Therefore, we considered the m⁶A score to be a possible prognostic factor for LUAD.

The prognostic value of the m⁶A score was investigated for patients with LUAD who accepted radiotherapy. Surprisingly, only in the m⁶A cluster C group did patients with radiotherapy had a better quality of life than those without radiotherapy ( Figures 8A–C ). Interestingly, patients with radiotherapy had a better quality of life than those without radiotherapy in the high m⁶A score group. The survival advantage from radiotherapy in the low m⁶A score group was virtually zero ( Figures 8D, E ).

Furthermore, the expression levels of PD-1 and CTLA4 were examined, and a remarkable elevation was observed in the low m⁶A score group ( Figures 8F, G ). The levels of other immune checkpoints were then compared in the high and low m⁶A score groups. The high m⁶A score group exhibited higher expression of CD27, CD28, LGALS9, TNFSF14, and TNFSF18, while the low m⁶A score group had higher expression of CD70, IDO1, LAG3, PDCD1LG2, and PVR ( Figure 8H ). In addition to the well-known TMB and checkpoint, IPS is one of the newly discovered predictive factors that is widely used and is strongly suggested to assess patients’ reaction to immunotherapy (52, 53). Compared with patients with low m⁶A score, patients with high m⁶A score exhibited significant clinical benefits from anti-PD-1/CTLA4 immunotherapy ( Figures 8I–L ). These findings confirmed that the levels of tumor m⁶A modification modes play a significant role in the regulation of the expression of immune molecules.