This study was conducted using The Cancer Genome Atlas (TCGA, https://cancergenome.nih.gov/) and Gene Expression Omnibus (GEO, National Center for Biotechnology Information, USA)) databases. We acquired transcript sequencing array data as the training cohort, measured in fragments per kilobase million (FPKM), for individuals identified as LUAD from TCGA. Additionally, we downloaded the transcriptome profile expression levels of two cohorts (GEO: GSE30219 and GSE50081) for external validation. The cohorts selected for analysis were based on the following criteria: 1. large-scale human samples of mRNA gene-expression patterns from untreated primary LUAD tissues with more than 30 samples; 2. assessed on the same technological platform that stores raw expression data and clinical information (such as survival times, censored information, and TNM stage); and 3. Peer reviews or publications in scientific journals proved or evaluated the data quality. A log2 scale was applied to all raw data after quantile standardization. Subsequently, we exercised prudence by excluding genes from the dataset with expression levels of 0 FPKM in at least 50% of the samples.

The human non–small-cell lung cancer (NSCLC) cell line (A549 and H460) was acquired from the American Type Culture Collection (ATCC, USA). The DMEM medium (BI, Israel) supplemented with 10% fetal bovine serum (BI, Israel) was used to cultivate the above cells in a temperature of 37°C under 5% CO2. To knock down METTL3, METTL14, and WTAP, all cells were transfected with small interfering RNA (siRNA) targeting these genes or control siRNA using Lipofectamine RNAiMAX (Invitrogen, USA) as per the manufacturer’s instructions. Untreated cells served as negative controls, and siRNA-targeting scrambled sequences were used as a transfection control. Transfection efficiency was validated using RT-qPCR, confirming significant knockdown of METTL3, METTL14, and WTAP, as shown in Supplementary Figure 4 . Supplementary Table 1 listed these siRNA sequences.

For METTL3 knockdown, lentiviral vectors harboring shRNA for knockdown and overexpression of METTL3 and negative control underwent syncretization and then cloned into pLKO.1 vector. The plasmids were transfected using lipofectamine LTX and Plus™ Reagent (Invitrogen, USA) into A549 cells according to the manufacturer’s protocol. The sequences are presented in Supplementary Table 1 . Briefly, stably transfected cells were selected with 10 μg/ml puromycin (MCE, USA) for 3 weeks.

The proteins of cells were extracted using the RIPA buffer, which was cooled before use (Beyotime, China). Identical protein samples were measured, loaded, separated on a 10% SDS-PAGE, and then transferred to 0.45 μm PVDF membranes (Beyotime, China). Following a 1.5-hour blocking step with 5% non-fat milk in TBST, the membranes were subjected to overnight incubation at 4°C with the primary antibodies ( Supplementary Table 2 ). Later on, the secondary antibodies were introduced and the concoction was subjected to incubation at room temperature for one hour. The proteins in the immunoblots were identified using the GelDoc XR System (BioRad, SA).

Human T cells were isolated from three NSCLC patients. After obtaining informed consent, 20 ml of peripheral venous blood was collected from the donor and T cells were then isolated via density gradient centrifugation using Ficoll-Paque solution (Biolegend, USA). The isolated T cells were subsequently cultured in Roswell Park Memorial Institute (RPMI)-1640 medium (BioInd, Israel), supplemented with 10% fetal bovine serum (Gibco, USA) and IL-2 (200 U/mL; SinoBiological, China). PBMCs were added to the A549 cell cultures at a ratio of 1:5 (A549 cells: PBMCs) in fresh complete medium, with or without 10 μg/ml atezolizumab (Genentech, USA), and incubated for 48 hours at 37°C in a humidified atmosphere containing 5% CO2. The supernatant was isolated from each group and the LDH release assay (Beyotime, China) was performed based on the instructions from the manufacturer. Absorbance was detected at 490 nm by Biotek microplate reader. In addition, the supernatant was also used for the enzyme-linked immune-sorbent assay (ELISA) for quantifying IFN-γ and IL-2 production (R&D System, USA) as per instructions given by the manufacturer. Processed data from plate readings taken at 450 nm.

Commercially available tissue microarray slides (HLugA180Su07) containing 93 histologically confirmed LUAD tissues were purchased from Biochip (Shanghai Biochip Co., Ltd., China) for immunohistochemistry (IHC) analysis. Immunohistochemical staining was performed on tissue microarrays (TMAs) incubating with the specific antibodies ( Supplementary Table 2 ) overnight at a temperature of 4°C. Subsequently, they were incubated with polyclonal peroxidase-conjugated anti-rabbit IgG (Boster Biological Technology co.ltd, USA) at room temperature for 20 min, as per the instructions provided by the manufacturer. Three expert pathologists blind to the clinical information independently graded each tissue sample. The intensity of staining was categorized as follows: 0 (negative), 1(weak), 2 (moderate), or 3 (strong). The extent of staining varied according to the proportion of positive cells (out of 200 cells examination): 0 (less than 5%), 1 (5%–25%), 2 (26%–50%), 3 (51%–75%), or 4 (>75%). The IHC expression scores were determined by multiplying the staining intensity by the staining extent, and then the scores were normalized by the z-score in order to calculate the risk scores.

RNA was isolated from the cells using the TRIzol Reagent (Invitrogen, USA) according to the manufacturer’s instructions. Reverse transcription of the isolated RNA was performed using HiScript II Q RT SuperMix for qPCR (+gDNA wiper) (Vazyme, China). RT-qPCR was performed using the AceQ qPCR SYBR Green Master Mix (Vazyme, China). The mRNA levels were standardized using glyceraldehyde-3-phosphate dehydrogenase (GAPDH) as the reference gene. The oligonucleotide sequences are listed in Supplementary Table 1 .

The Magna RIP Kit (Millipore, MA, USA) was used to conduct the RIP assay. Briefly, 5 μg anti-METTL3 (Abcam, USA) or anti-m⁶A (Millipore, Germany) and anti-rabbit IgG (Millipore, Germany) were incubated with 50 μL of magnetic beads before cell lysates were added (approximately 2 × 10⁷ cells per sample). Following six rounds of washing, the RNA–protein immunoprecipitation (IP) complexes were incubated in proteinase K digestion solution to extract the proteins. Finally, the RNA was purified for RT-qPCR analysis after being extracted using phenol–chloroform. Normalizing relative enrichment to the input was done as % input =1/10 × 2^{Ct [IP] – Ct [input]}.

To discern the relationship between m⁶A and immune status, the t-distributed Stochastic Neighbor Embedding (t-SNE) algorithm was performed (16). Moreover, signature gene sets exhibit consistent expression and serve as summaries of certain clearly defined biological states or processes. According to previous publications, we retrieved expression matrixes of m⁶A regulators (17). Additionally, 29 immune-related genes were selected, reflecting diverse array of immune cell types, roles, and pathways (IMMPORT) (18). We employed the Non-negative Matrix Factorization (NMF) technique to cluster the expression patterns of immune- or m⁶A-associated genes. We used k-means clustering because it allows for the identification of distinct subgroups based on m⁶A modification patterns, assuming that the data are relatively homogeneous within each cluster. Using Cox regression analysis, we assessed the correlation between each potential gene and overall survival (OS). Three clusters were determined to be the optimal number when the correlation coefficient decreased. Using the aforementioned immune/m⁶A gene mRNA expression data and the t-SNE algorithm, LUAD subtypes were identified. The differentially expressed genes (DEGs) were determined with a criterion of |logFC| > 0.585 and P<0.05, after controlling for false discovery rate (FDR). In addition, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) (https://www.kegg.jp/) and Gene Ontology (GO) pathway analyses (https://david.ncifcrf.gov//) using the R Cluster-Profiler. Pathways were considered significant when P and q values were below the 0.05 threshold.

Genes that overlapped between the m⁶A- and immune-associated DEGs were selected for further analysis. The risk scoring system was established using Cox regression analysis and the Least Absolute Shrinkage and Selection operator (LASSO) (19). Using five-fold cross-validation (CV), we were able to find the optimal parameters while reducing the bias that overfitting the training samples may induce. The prognostic roles of candidate genes were subsequently investigated using multivariate Cox regression analysis. Risk scores were generated by multiplying the multivariate Cox regression’s coefficient of gene expression. All patients were categorized as either high- or low-risk according to the median risk score. The data were analyzed using Kaplan-Meier survival analysis, and statistical significance was determined by log-rank tests. We tested the model’s prediction power using a receiver operating characteristic (ROC) curve.

To assess gene set enrichment in transcriptomes, we conducted Gene Set Variation Analysis (GSVA). This technique scores gene sets to convert gene expression levels to pathway levels, thus determining the biological function of samples (20). We evaluated the immunological state and m⁶A levels in risk groups using data from the Molecular Signature Database, with significant results at FDR q < 0.25 and P < 0.05. Immunocyte infiltration in the groups was analyzed via the CIBERSORT algorithm, focusing on 22 distinct immunocyte subunits. The TIMER (version 2.0) database provided data on tumor-infiltrating immune cell abundance (21). Spearman’s correlation analysis was employed to study the relationship between gene expression and immune cell concentrations. Additionally, using data from the Genomics of Drug Sensitivity in Cancer database (GDSC), we estimated the IC50 of chemotherapeutic medications through regression modeling.

Data analysis was carried out using the R (v.4.2.0) software. The Chi-square test was used for the analysis of qualitative variables. A statistically significant result was defined as a significance level of P < 0.05, unless specified otherwise, in specific conditions as outlined independently.

The processed original mRNA expression data for LUAD were acquired from the TCGA database. Subsequently, 1811 immunity genes and 1670 immune regulatory factors were identified from the IMMPORT database. We further determined three immune clusters and calculated the Euclidean distance using t-SNE for each patient ( Figure 1A and Supplementary Figures 1A, B ). Statistically significant disparities were observed across the three clusters using survival analysis. Notably, C3 exhibited longer median survival when comparing both C1 and C2 clusters ( Figure 1A ). Similarly, three m⁶A clusters were detected based on the m⁶A regulators expression matrices. The results indicated cluster C2 had a comparatively longer median survival time ( Figure 1B and Supplementary Figures 1C, D ). Furthermore, we identified 955 immune-associated DEGs and 282 m⁶A -associated DEGs. Among these, 145 genes were found to be co-expressed in both the immunological and m⁶A subtypes, making them potential candidate DEGs for further investigation.

Using pathway analysis, the possible roles of the co-expressed DEGs were determined (Fisher’s exact test, P <0.05). GO enrichment analysis yielded biological process (BP), molecular function (MF), and cellular component (CC) terms. For BP, DEGs were predominantly enriched in nuclear division, regulation of cell-cycle phase transition, regulation of immune effector processes, and immune responses. In the MF, the majority of DEGs were involved in antigen binding and cyclin-dependent protein serine/threonine kinase regulatory activity. The spindle, immunoglobulin complex, and midbody were most abundant in CC enrichment ( Figure 1C ). KEGG pathway enrichment analysis showed that DEGs had a significant enrichment mostly in the Cell cycle, p53 signaling, and T cell receptor signal transduction pathways ( Figure 1D ). Collectively, the pathway enrichment results suggested that the DEGs we chose have a tight connection to cell proliferation and immune response, which may not only serve as prognostic indicators for patients with LUAD but also contribute significantly to the interplay between immune infiltration and well-established signaling pathways associated with tumor development and invasion.

There were 18 immunologically and m⁶A-linked DEGs that contributed to OS after conducting the univariate Cox regression analysis ( Supplementary Table 3 ). We then used LASSO to identify seven signatures for a risk-based prognostic assessment model ( Figure 2A and Supplementary Figure 2A ). We used the following formula to calculate risk scores: Risk Score = SFTPC*(-0.061883497) + CYP24A1*0.086651739 + KRT6A*0.127175139 + PTTG1*0.146161987 + S100P*0.14825993 + FAM83A*0.164912979 + ANLN*0.169266477. All patients were classified as the high-risk or low-risk cohort based on the median risk score. In both the training and internal test set, the high-risk group exhibited a lower OS compared to the OS seen in the low-risk group (P < 0.001 and P=0.026, Supplementary Figure 2B and Figure 2B ) and the AUC values of the risk score for 3-year survival were 0.70 and 0.71 ( Supplementary Figure 2C and Figure 2B ). Subsequently, two additional independent cohorts were included for external validation ( Supplementary Tables 4 , 5 ). As expected, high-risk patients in both cohorts had shorter OS compared to low-risk individuals (P = 0.002 and 0.045, respectively, Figures 2C, D ). The AUCs with regard to 3-year OS were 0.65 and 0.68, respectively ( Figures 2C, D ). This further implied that the risk score system showed excellent repeatability and stability during validation. Furthermore, we performed Cox regression analyses and determined that the risk score served as a reliable and independent prognostic indicator ( Figures 2E, F ). Moreover, there was a substantial correlation between the risk scores and the tumor, node, metastasis (TNM), T, and N stages ( Figure 2G ). Taken together, the risk score model we constructed is not only capable in predicting the survival but also has a noteworthy connection with clinical characteristics.

To evaluate the risk cohorts and their correlation with immune status, we primarily programmed GSVA and suggested the risk cohorts were enriched in the IL2-STAT5, PI3K/AKT/MTOR, and Reactive oxygen species pathways ( Figure 3A and Supplementary Figures 3A, B ). Furthermore, assessment of immune cell infiltration found that the risk score was positively correlated with T cell CD4 + memory activation, Macrophages M0, Macrophages M1, activated NK cells, follicular helper T cells, CD8 + T cells, and activated Mast cells, but negatively associated with Macrophages M2, Monocytes, resting CD4 + memory T cells, dendritic cells, and resting mast cells in patients with different risk levels ( Figure 3B and Supplementary Figure 3C ). Interestingly, the expression of PD-L1, CTLA-4, and IDO1 increased as the risk score elevated ( Figure 3C and Supplementary Figure 3D ). Genetic mutation analysis revealed that individuals considered to be at high risk had higher incidences of TP53, TTN, and MUC16 mutations ( Figure 3D ). Quantitative analysis further confirmed that load of tumor mutational burden (TMB) and neoantigen levels were increased in the high-risk group ( Figure 3E ). In addition, we predicted the chemosensitivity of each tumor sample and the results indicated that low-risk individuals showed greater chemosensitivity to Bleomycin, Cytarabine, Paclitaxel and Docetaxel ( Figure 3F ). Taken together, the findings indicate that the low-risk group could exhibit increased sensitivity to conventional chemotherapy. However, patients with high-risk scores exhibited higher immune cell infiltration and were positively correlated with the expression of immune checkpoints, indicating that they may potentially receive therapeutic benefits from immunotherapy.

Due to the proteins associated with LUAD and the TIME perform important biological functions, we conducted further experimental verification using LUAD tissue microarrays and IHC ( Figure 4A ). Following the exclusion of invalid samples, the risk score for every LUAD patient was determined using the aforementioned formula and further classified as low- or high-risk according to the median. Survival analysis revealed that patients with high risk exhibited worse OS (P<0.001; Figure 4B ). The Cox regression analyses also showed that the risk score was the only independent prognostic indicator ( Figure 4C ). Consistent with prior findings, there was a substantial correlation between the risk scores and both TNM and T stages, which indicated a robust relationship between the risk score and tumor invasion ( Figure 4D ). Notably, PD-L1 protein expression levels were significantly higher in the high-risk group of patients (P = 0.002; Figure 4E ). Collectively, these results experimentally verified the stability and reliability in the scoring system at protein level.

As a crucial component for RNA m⁶A modification, METTL3 is essential for the regulation of TME and antitumor immunity in NSCLC (8, 22). To investigate regulatory role of METTL3 in the risk score model, we initially found a positive correlation between METTL3 expression and risk scores using Spearman correlation (R = 0.137, P<0.05; Figure 5A ). Additionally, m⁶A modification data for candidate signatures were retrieved from the m⁶A target database (23). We subsequently investigated whether METTL3 exerts regulatory effects on candidate signatures, given its role as a key m⁶A writer in lung cancer. In A549 cells, METTL3 knockdown resulted in significant inhibition of CYP24A1, KRT6A, S100P, FAM83A, PTTG1, and ANLN expression, while SFTPC expression was significantly enhanced ( Figure 5B ). Similar expression changes were observed in H460 cells, except for ANLN, which did not show significant alteration ( Figure 5C ). This discrepancy likely reflects the distinct genetic backgrounds and m⁶A regulatory landscapes of the two cell lines, highlighting the complexity of m⁶A modifications. To confirm the role of METTL3 as the potential “writer”, we evaluated the direct binding interaction between METTL3 and the mRNAs of the seven candidate genes using RIP-qPCR assays in both cell lines. There was significant METTL3 enrichment in the mRNAs of all seven candidate genes ( Figure 5D ), and this enrichment decreased upon METTL3 silencing ( Figure 5E ). These findings confirm that METTL3 directly binds to and potentially methylates the mRNAs of the candidate genes in an m⁶A-dependent manner.

Subsequently, we explored the functional consequences of METTL3 modulation. We constructed A549 cell lines with stable overexpression (oe-METTL3) and knockdown (sh-METTL3) of METTL3 ( Supplementary Figure 4A ). The results found that oe-METTL3 increased PD-L1 expression, while sh-METTL3 decreased PD-L1 expression at both the mRNA and protein levels ( Figures 5F, G ). Coculture experiments with A549 cells and PBMCs demonstrated that METTL3 knockdown significantly enhanced T cell-mediated cytotoxicity, whereas METTL3 overexpression inhibited these T cell functions; notably, the addition of atezolizumab partially restored T cell antitumor effects with METTL3 overexpression ( Figure 5H ). Furthermore, T cell proliferation, indicated by increased IL-2 and IFN-γ production, was enhanced following METTL3 knockdown, while METTL3 overexpression inhibited these functions and atezolizumab also mitigated the inhibitory effects of METTL3 overexpression on T cell activity to some extent ( Figures 5I, J ). Collectively, these results indicate that METTL3-mediated m⁶A modifications regulate candidate gene expression in a cell line-specific manner and modulate the antitumor immune response, underscoring the potential of targeting METTL3 in NSCLC therapy.