Samples of paracancerous and LUAD tissues used for immunohistochemistry were obtained from 10 patients admitted to the Cancer Center of Guangzhou Twelfth People's Hospital (Guangzhou, China) from February 2019 to May 2021. The cohort comprised 6 men and 4 women, with a median age of 57 years (range, 31–76 years). Clinical data for the patients are shown in Table SI. All patients provided their written informed consent for study participation.

Gene expression profiles and clinical information for a cohort of 585 patients with LUAD from The Cancer Genome Atlas (TCGA-LUAD) were downloaded to serve as the training set. Notably, all 585 samples belong to this single TCGA-LUAD cohort, with data obtained from two distinct sources: TCGA database (http://portal.gdc.cancer.gov/), which provided 284 tumor and 37 adjacent normal samples, and the UCSC Xena database (https://xena.ucsc.edu/), which provided 242 tumor and 22 adjacent normal samples. The dataset used in the analysis is provided in Table SII. In addition, gene expression datasets [accession nos. GSE135222 (19) and GSE126044 (20)] were downloaded from the Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) series to analyze the expression levels of m6A-related genes.

Based on the method previously described (21), the ‘DESeq2’ R package (v 3.5.1; http://bioconductor.org/packages/release/bioc/html/DESeq2.html) was used to identify genes that were differentially expressed between the 526 LUAD and 59 adjacent non-tumor samples, with a threshold of log2 FC >1 and P<0.05 was considered to indicate a statistically significant difference.

The Tumor IMmune Estimation Resource (TIMER) database (http://timer.cistrome.org/) was used to evaluate the correlation between m6A regulators and the levels of immune cell infiltration [including cancer-associated fibroblasts, bone marrow dendritic cells (BMDCs), CD4⁺ T cells, neutrophils, regulatory T cells (Tregs), CD8⁺ T cells and macrophages]. The Tumor and Immune System Interaction Database (TISIDB) was used to assess the correlation between the expression of m6A regulators and immunoregulators [including immunosuppressant, immunostimulator and major histocompatibility complex (MHC) molecules].

cBio Cancer Genomics Portal (cBioPortal; http://www.cbioportal.org/), an open-access website that explores, visualizes and analyzes multidimensional cancer genomics data, was used to analyze the genetic alterations of m6A regulators in LUAD.

Mutation annotation format files were analyzed using the R package ‘Maftools’ (v 2.25.10; Bioconductor; http://bioconductor.org/packages/devel/bioc/html/maftools.html) to summarize, visualize and interpret somatic mutations in m6A regulator genes across the LUAD cohort.

The correlation between m6A regulator aberrations and the survival time of patients with cancer was determined using the cBioPortal database. The correlation between the overall survival time and the expression of m6A regulators was evaluated using Kaplan-Meier curves.

First, the correlation between m6A gene expression and biomarkers of T cell exhaustion [TNF, IL-2, IFN-γ and cytotoxic T lymphocyte (CTL)] was analyzed using the R software package ‘psych’ (version 2.4.3; http://www.rdocumentation.org/packages/psych/versions/2.4.3), based on research datasets in a previous study (22).

Next, the expression levels of m6A genes in patients with LUAD treated with anti-programmed cell death protein-1(PD-1)/PD-L1 were analyzed the GSE135222 and GSE126044 datasets. In these datasets, patients with LUAD were grouped according to whether they had responded to anti-PD-1/PD-L1 treatment [19 non-responders and 8 responders in GSE135222 (23); 11 non-responders and 5 responders in GSE126044 (20)]. Furthermore, analyzing the correlation between the expression levels of m6A-related genes and immune cell marker genes was analyzed using Spearman's rank correlation to clarify the role of m6A-related genes in predicting immune infiltration responses.

Samples of LUAD and paracancerous tissue were collected at Cancer Center of Guangzhou Twelfth People's Hospital. Immunohistochemical procedures were performed on the tissues to assess gene expression according to the methods reported in a previous study (24). In brief, the tissues were fixed in 10% formalin, embedded in paraffin, and cut into 4-µm sections. Subsequently, the sections were dewaxed using gradient ethanol and subjected to antigen retrieval in citrate buffer (cat. no. G1202; Wuhan Servicebio Technology Co., Ltd.) at 95°C for 20 min. The slides were then treated with 3% H₂O₂ at 25°C for 30 min to block endogenous peroxidase activity, followed by blocking with 3% bovine serum albumin (cat. no. GC305006; Wuhan Servicebio Technology Co., Ltd.) at room temperature for 30 min. Next, the sections were incubated with anti-heterogeneous nuclear ribonucleoprotein C (HNRNPC; cat. no. ab75822; dilution 1:800), anti-insulin-like growth factor 2 mRNA binding protein (IGF2BP)1 (cat. no. ab184305; dilution 1:4,000), anti-IGF2BP3 (cat. no. ab179807; dilution 1:1,000), anti-PD-L1 (cat. no. ab205921; dilution 1:1,000), anti-CD4 (cat. no. ab133616; dilution 1:500), anti-CD8 (cat. no. ab217344; dilution 1:2,000) and anti-Forkhead box P3 (FOXP3; cat. no. ab20034; dilution 1:500), respectively, followed by incubation with goat anti-rabbit IgG (HRP; cat. no. ab7090; dilution 1:1,000) (all from Abcam) and DAB chromogen. The tissue sections were then counterstained with hematoxylin, observed and captured under a light microscope.

Statistical analyses were performed using SPSS (version 26.0; IBM Corp.) and GraphPad Prism (version 10.1.2; Dotmatics). Differential expression analysis of m6A methylation regulators between tumor and paired adjacent normal tissues was assessed using the paired Student's t-test for parametric data or the Wilcoxon signed-rank test for non-parametric data, with 59 paired tumor and normal tissue samples. Univariate Cox proportional hazards regression analysis was used to evaluate the association between individual m6A regulator expression levels and overall survival time. Kaplan-Meier survival curves were generated and compared using the log-rank test. The correlation between T-cell exhaustion markers and other variables of interest was assessed using Spearman's rank correlation coefficient (ρ), with statistical significance determined by the corresponding P-value. Differences in gene expression across ordinal clinical variables [tumor, lymph node and metastasis (TNM) stages) were assessed using one-way ANOVA with Tukey's post hoc test for multiple comparisons. Mutation analysis of m6A genes was performed. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed on genes identified through Pearson's correlation analysis (r>0.5, P<0.001) of m6A-related genes from TCGA expression matrix that were co-expressed with HNRNPC, IGF2BP1 and IGF2BP3. P<0.05 was considered to indicate a statistically significant difference for all analyses, unless otherwise indicated.

To evaluate the mutation load in LUAD, the tumor mutation load in 33 TCGA tumor types was compared. The present study data demonstrated that LUAD had a high mutation load, which was lower compared with the mutation loads of bladder urothelial carcinoma, lung squamous cell carcinoma and skin cutaneous melanoma (Fig. 1).

The differential expression levels of m6A genes in the tumor and control groups and the correlation patterns were analyzed. Excluding YTHDF3, YTHDC protein (YTHDC)1, YTHDC2 and IGF2BP2, the expression levels of all m6A genes significantly differed between the tumor and control groups. (P<0.01; Fig. 2A).

A correlation analysis indicated a moderate positive correlation between m6A genes, among which, YTHDF3 had the strongest correlation with vir like m6A methyltransferase associated (KIAA1429) and YTHDC1 a slightly weaker correlation with METTL14 and Zinc finger CCCH domain-containing protein 13 (ZC3H13). Furthermore, FTO was negatively correlated with METTL5 (r=0.19; Fig. 2B).

To identify somatic mutations in patients with LUAD, mutation data were analyzed using the R software package ‘Maftools’. Results indicated that 110 of 561 patients with LUAD had mutated m6A-related genes. As shown in Fig. 3A, the mutated m6A-related genes included KIAA1429 (3%), ZC3H13 (3%), IGF2BP1 (3%), YTHDC1 (2%), YTHDC2 (2%), RNA binding motif protein 15 (RBM15; 2%), IGF2BP2 (1%), IGF2BP3 (1%), HNRNPC (1%), YTHDF1 (1%), METTL14 (1%), YTHDF3 (1%), METTL14 (1%), FTO (1%), Wilms' tumor 1-associating protein (1%), bacterial alkane hydroxylase homolog 5, RNA demethylase (1%) and METTL116 (1%).

Furthermore, the somatic interactions between m6A genes and the status of single nucleotide polymorphisms (SNPs) and highly mutated genomic regions were further explored. C>A and C>T were the two main mutation types. The proportion of each sample is shown in a stacked histogram (Fig. 3B). The rainfall plot demonstrated high mutation genomic regions based on different SNP mutation types (Fig. 3C). As shown in Fig. 3D, IGF2BP1 had significant coexpression frequencies with RBM15 and IGF2BP2, and KIAA1429 had significant coexpression frequencies with IGF2BP2 and YTHDC1 (P<0.05).

To further evaluate the prognostic value of m6A RNA methylation regulators in LUAD, the relationship between their expression levels and overall patient survival in TCGA database was determined by a univariate Cox regression and Kaplan-Meier analysis. The Cox regression results revealed that five genes were associated with LUAD prognosis. High expression levels of HNRNPC, IGF2BP1, IGF2BP2, IGF2BP3 and METTL5 were with significantly associated with a poor prognosis in LUAD (P<0.001; Fig. 4A-E).

The associations between HNRNPC, IGF2BP1 and IGF2BP3 expression and clinical parameters, including tumor and TNM stages (25), were evaluated. The present study results indicated that there was a significant difference in HNRNPC expression between the normal group and the four stages (P<0.001; Fig. 5A). In addition, a statistical difference was observed between stage 1 and 4 (P=0.031). There were significant differences between the normal group and the four stages in the T stage (all P<0.001), and there was a significant difference between the normal group and the two stages in the M stage (both P<0.001). Significant differences were also observed between the normal group and the three stages in the N stage (all P<0.001).

For IGF2BP1 expression, results indicated that there was a statistically significant difference between the normal group and the four stages in stage (P<0.001; Fig. 5B). Significant differences were also found between the normal group and the four stages in the T stage (P<0.01). Furthermore, a significant difference was observed between T stages 1 and 2 (P=0.044). There were significant differences between the normal group and M0 and M1 stages in the M stage (P<0.001 and P<0.01, respectively). Significant differences were also found between the normal group and the three stages in the N stage (P<0.001).

For IGF2BP3 expression, results indicated that there was a statistically significant difference between the normal group and the four stages (P<0.001; Fig. 5C). Significant differences were observed between the normal group and the four stages in the T stage (all P<0.001). There was also a significant difference between the normal group and M0 and M1 stages in the M stage (both P<0.001). Significant differences were also found between the normal group and the three stages in the N stage (P<0.001 for all). Next, the expression levels of HNRNPC, IGF2BP1 and IGF2BP3 were detected in LC tissues by IHC. The results demonstrated that these three genes were all upregulated in LUAD when compared with paracancerous tissues (Fig. 6).

The TIMER database was used to evaluate the correlation between m6A genes and immune cell infiltration level in LUAD. Results revealed that BMDCs (ρ=0.494; P=8.96×10⁻³²) and Tregs (ρ=−0.181; P=5.02×10⁻⁵) were significantly positively and negatively correlated with HNRNPC, respectively (Fig. 7A). Cancer-associated fibroblasts (ρ=0.155; P=5.02×10⁻⁴), BMDCs (ρ=0.374; P=6.84×10⁻¹⁸) and macrophages (ρ=0.155; P=5.54×10⁻⁴) were significantly positively correlated with IGF2BP1 (Fig. 7B). IGF2BP3 indicated significant positive correlations with cancer-associated fibroblasts (ρ=0.208; P=3.09×10⁻⁶), BMDCs (ρ=0.389; P=2.65×10⁻¹⁹), neutrophils (ρ=0.16; P=3.56×10⁻⁴), CD8⁺ T cells (ρ=0.224; P=4.66×10⁻⁷), and macrophages (ρ=0.239; P=7.39×10⁻⁸); however, there was a negative correlation with CD4⁺ T cells (ρ=−0.094; P=3.74×10⁻²) (Fig. 7C). Furthermore, IHC staining was performed to assess the expression levels of PD-L1, CD4, CD8 and FOXP3 in clinical samples. Results indicated that the expression levels of PD-L1, CD4, CD8 and FOXP3 were elevated in LUAD cancer tissues (Fig. 7D).

First, the correlation between m6A gene expression and biomarkers of T cell exhaustion (TNF, IL-2, IFN-γ and CTL). The results indicated that HNRNPC was negatively correlated with TNF expression (P<0.05), IL-2 expression (P<0.001) and CTL expression (P<0.05), IGF2BP3 was positively correlated with IFN-γ expression (P<0.01), while IGF2BP1 (P>0.05) was not correlated with any of the biomarkers (Fig. 8A). Next, the expression levels of m6A genes in patients treated with anti-PD-1/PD-L1 were analyzed based on GEO datasets (dataset nos. GSE135222 and GSE126044). However, the expression levels of HNRNPC (P=0.163; P=0.583), IGF2BP1 (P=0.69; P=0.743) and IGF2BP3 (P=0.621; P=0.583) were not significantly different between the non-response group and response group (Fig. 8B and C).

The TISIDB database was used to analyze the correlation between m6A genes and the expression of immunomodulators. Furthermore, to explore the effect of m6A regulatory factors on tumor immune response, the correlation between m6A regulatory factors and the expression of immunoregulatory factors was assessed. The present results demonstrated that HNRNPC and IGF2BP1 were negatively correlated with immunosuppressants (Fig. 9A), immunostimulators (Fig. 9B) and MHC molecules (Fig. 9C), whereas IGF2BP3 was negatively correlated with immunostimulants and MHC molecules and positively correlated with immunosuppressants.

To identify the genes co-expressed with HNRNPC, IGF2BP1 and IGF2BP3, a Pearson's correlation analysis of m6A genes from TCGA expression matrix was performed. r>0.5 and P<0.001 were used for screening in the subsequent analysis. A total of 276 genes exhibiting co-expression with HNRNPC were identified, along with 22 for IGF2BP1 and 80 for IGF2BP3. To investigate the downstream pathways of hub m6A regulators in LUAD, GO and KEGG (26–28) analyses were performed using coexpression genes of the three m6A regulators. The present study demonstrated that HNRNPC was mainly enriched in the ‘ribonucleoprotein complex biogenesis’ pathway (Fig. 10A), IGF2BP1 in the ‘mitotic cell cycle checkpoint’ pathway (Fig. 10B) and IGF2BP3 in the ‘nuclear division’ pathway (Fig. 10C).