The N⁶-methyladenosine Epitranscriptomic Landscape of Lung Adenocarcinoma

Lung cancer remains the most commonly diagnosed cancer globally and the leading cause of cancer-related death in both men and women (1). Non–small cell lung cancer (NSCLC), which accounts for approximately 80% to 85% of all cases, is predominantly represented by lung adenocarcinoma (LUAD). LUAD is particularly notable for its poor survival rates (1), underscoring the urgent need for improved understanding of the molecular mechanisms and therapeutic strategies. Recent advancements in molecular technologies, including genomics, transcriptomics, epigenomics, and proteomics, have unveiled some of the genetic, genomic, and epigenetic alterations underlying LUAD lethality (2, 3). In contrast, posttranscriptional mechanisms contributing to tumor biology remain underexplored. Epitranscriptomics, often referred to as the “epigenetics” of RNA, encompasses myriad RNA modifications that refine gene expression (4). N⁶-methyladenosine (m⁶A), the most prevalent modification in eukaryotic RNAs, governs mRNA export, splicing, translation, and degradation (5, 6). A large number of m⁶A modifications are dynamically catalyzed (i.e., added or removed) by writer complex METTL3/METTL14/WTAP and removed by erasers FTO and ALKBH5 (7). These modifications occur in more than 25% of human transcripts (8), often near the stop codons of mRNA transcripts with the consensus motif RRACH (R = A or G; H = A, U, or C; refs. 5, 8), influencing RNA fate and gene expression through recruitment of m⁶A readers like YTHDF1/2/3, YTHDC1/2, and IGF2BP1/2/3 (5, 6, 9).

There is growing evidence that underscores the significant dysregulation of m⁶A modifications in tumorigenesis and metastasis (7, 10). m⁶A influences cancer hallmarks such as proliferation, stemness, metastasis, and immune evasion. Most m⁶A regulators (writers, erasers, and cofactors) and readers have been implicated in modulating tumor biology via regulating the RNA stability or protein translation of key oncogenes and tumor suppressors during cancer onset and progression (7, 10). Notably, while many of these studies have explored the role of m⁶A by modulating the abundance of one or multiple m⁶A regulators in vitro or in vivo, the composite m⁶A methylome in primary tumors, a product of the orchestrated dysregulation of m⁶A regulators, remains uncharacterized.

m⁶A immunoprecipitation followed by high-throughput sequencing (m⁶A MeRIP-seq) offers a robust method for analyzing transcriptome-wide m⁶A modifications (8, 11). Employing a previously developed low-input m⁶A MeRIP-seq protocol (12), we profiled the m⁶A epitranscriptomes of 10 nonneoplastic lung (NL) tissue samples and 51 primary LUAD tumors. This comprehensive dataset allowed us to assess the overarching consequences of the disrupted m⁶A regulators on m⁶A epitranscriptome in LUAD tumors. Our analysis discerned three gene clusters, each of which was characterized by unique genomic features and underwent methylation via either canonical or noncanonical manners. By integrating this epitranscriptomic data with clinicopathologic and patient outcome metrics, we identified distinctive m⁶A patterns and site-specific modifications that were exclusively associated with patient outcome, bypassing any mRNA and protein correlation. We pinpointed hypermethylation of EML4 as a new driver of metastasis. Pharmacologic inhibition of METTL3 profoundly reduced EML4 m⁶A modification and protein abundance, leading to a significant reduction in tumor metastases. This study not only provides an extensive m⁶A profile resource for LUAD research but also lays the groundwork for further studies on m⁶A regulatory dynamics and novel therapeutic targets.

Applying an optimized m⁶A MeRIP-seq approach previously developed by our group (12), we profiled 10 NL tissues and 53 primary LUAD tumors, using 2 μg of total RNA for each sample (Fig. 1A; Supplementary Fig. S1A; Supplementary Table S1). Primary lung tumors were sourced from a subset of patients accrued to the University Health Network Patient-derived Lung Cancer Model Establishment (UHN-PLCME) project, with transcriptomic and proteomic data and rich clinical annotation available (13–15). NL tissues were histologically nonneoplastic lung tissues collected from normal-appearing lung of lobectomy specimens from patients undergoing lung cancer resection. Escherichia coli (E. coli) K-12 spike-in RNA was used for library size normalization. On average, more than 80% of the MeRIP IP reads were uniquely mapped, with 85.8% originating from protein-coding genes and 10.5% from long intergenic noncoding RNA (lincRNA) genes (Supplementary Table S1). Notably, only 0.06% of the uniquely mapped IP reads were from ribosomal RNA (rRNA)/small nuclear RNA (snRNA)/small nucleolar RNA (snoRNA). Principal component analysis conducted on read counts from annotated genes revealed distinct separation between NL tissues and tumors in both m⁶A immunoprecipitated RNA (IP) and the total RNA (INPUT) libraries (Fig. 1B). Two tumor samples were identified as outliers with high percentages of spike-in reads in the INPUT libraries (Fig. 1B; Supplementary Table S1), indicating either poor RNA quality or inaccurate RNA quantification. Consequently, these two tumor samples were excluded from downstream analysis. m⁶A peak numbers in NL and tumor samples varied between 14,170 and 15,806 and between 8,141 and 15,878, respectively. While m⁶A peaks were relatively evenly distributed across the gene bodies of lincRNAs (Supplementary Fig. S1B), m⁶A peaks in protein-coding genes predominantly clustered near the stop codon for both NL and tumor samples (Fig. 1C). The motif “RRACH” was consistently enriched in the majority of samples (Supplementary Table S1), in accordance with previous reports (8, 11). The quantity and distribution of peaks, along with the detection of the consensus motif, underscore the effectiveness of the low-input m⁶A MeRIP-seq method in profiling the m⁶A epitranscriptome from patient specimens.

A gene was defined as m⁶A methylated if it contains an m⁶A peak(s) in at least two samples, either NL or tumor. In total, 11,409 such genes were detected across NL and tumor samples, with an average of 1.8 peaks per gene (Fig. 1D). To reliably quantify m⁶A level, we further excluded genes with an abundance below Reads Per Kilobase Million of 1 in both INPUT and IP libraries for more than half of NL or tumor samples. This filtration yielded 6,140 unmethylated and 8,030 methylated genes, inclusive of 184 methylated lincRNAs (Fig. 1D). We then quantified the m⁶A level of each gene with three distinct methods: (i) the maximum ratio of each peak’s IP signal to the corresponding INPUT signal (max peak), (ii) the ratio of the sum of all peaks’ IP/INPUT (sum peak), and (iii) the ratio of whole-transcript IP/INPUT (sum gene). We cross-referenced each method using an external dataset with matched m⁶A MeRIP-seq and m⁶A-LAIC-seq in human embryonic stem cell line H1-hESC (16, 17). m⁶A-LAIC-seq, a technology that extracts both m⁶A methylated and unmethylated transcripts simultaneously (17), serves as a benchmark for m⁶A quantification at the gene level.

All three methods were strongly correlated in H1-hESC cells (Supplementary Fig. S1C) and exhibited similar trends in both NL and tumor samples (Supplementary Fig. S1D). The whole-transcript IP/INPUT ratio showed the strongest correlation with the percentage of methylated transcripts as determined by m⁶A-LAIC-seq (Supplementary Fig. S1C; r = 0.67; P < 2.2 × 10⁻¹⁶), making it the preferred method for subsequent analyses. Whole-transcript IP/INPUT ratios (m⁶A level hereafter) of the methylated genes were significantly higher compared with the unmethylated genes in both NL and tumor samples (Fig. 1E), further supporting its efficacy as the optimal method for m⁶A quantification from m⁶A MeRIP-seq data. In both NL and tumor samples, m⁶A levels were highly correlated with those observed in H1-hESC (Supplementary Fig. S1E), in agreement with previous findings that the m⁶A level for a given gene is highly conserved across varied cell types (8, 18). This consistency allowed us to project m⁶A levels in our dataset to the percentage of methylation as characterized by m⁶A-LAIC-seq. Notably, log₂(m⁶A level) ranges of −6 to 0, 0 to 2, and 2 to 6 corresponded to methylation percentages of 20% to 40%, 40% to 60%, and 60% to 80%, respectively (Fig. 1F).

The m⁶A level is dynamically regulated by m⁶A writers and erasers (7). Thus, we sought to systematically assess the correlation between the mRNA abundance of m⁶A regulators and m⁶A levels for each methylated gene in our cohort. We included three genes (METTL3, METTL14, and WTAP) in the catalytic core components and three essential regulatory adaptors (VIRMA, CBLL1, and ZC3H13) in the writer complex, along with two erasers (ALKBH5 and FTO) in our analysis. First, we examined the mRNA abundance of these regulators in our cohort and the TCGA LUAD cohort (Supplementary Fig. S2A). Both METTL3 and VIRMA showed significant upregulation, while WTAP and FTO were significantly downregulated in tumors compared with NL tissues in both cohorts (Supplementary Fig. S2A). This observation suggests a complex dysregulation of m⁶A in tumors.

Next, we conducted clustering analysis based on the correlations between mRNA abundance of the regulators and m⁶A levels of the 8,030 methylated genes in tumor samples. The writer complexes and erasers were clearly separated and clustered into two distinct groups (Fig. 2A). These 8,030 methylated genes formed three groups (G1, G2, and G3) based on their correlation patterns (Fig. 2A; Supplementary Table S2). The m⁶A levels of most G1 genes displayed a positive correlation with writer complexes and a negative correlation with erasers (Fig. 2A). G2 genes were featured by positive correlation with METTL3 and METTL14 but a negative correlation with CBL11, VIRMA, and the two erasers. The majority of G3 genes exhibited a negative correlation with writer complexes and positive correlation with erasers (Fig. 2A). This diverse association between m⁶A levels and regulator abundances aligns with previous reports in which m⁶A levels in certain genes increased, while others decreased after knocking down m⁶A writer complexes (19). In A549 cells, the proportions of the binding target genes of m⁶A regulators in the G3 group were significantly lower than in the G1 and G2 groups (Fig. 2B; Supplementary Table S2; Pearson’s χ² test). Similarly, G3 genes exhibited significantly decreased proportions of target genes of m⁶A regulators compared with G1 genes in HeLa and HEK-293 cells (Supplementary Fig. S2B and S2C; Pearson’s χ² test; see Supplementary Methods). About the binding targets of m⁶A readers, we collected the public data of the target genes of YTHDF, YTHDC2, and IGF2BP proteins identified in A549, HeLa, and HEK-293 cells (see Supplementary Methods). We found that the targets of YTHDC2 and IGF2BP proteins were significantly depleted in the G3 group, whereas the proportions of the binding target genes of YTHDF proteins in the G2 group were significantly higher than those in the G1 and G3 groups (Supplementary Fig. S2D and S2E; Pearson’s χ² test). Collectively, G3 genes exhibited distinct binding patterns of m⁶A regulators and readers, with evident lower binding by m⁶A regulators and readers. The m⁶A levels of G3 genes were likely influenced by indirect effects or other uncharacterized regulatory mechanisms.

Among the three groups, G1 genes have the lowest m⁶A levels (Wilcoxon rank-sum one-sided test: G1 vs. G2 and G1 vs. G3: P values < 2.2 × 10⁻¹⁶; G2 vs. G3: P = 0.49) and highest mRNA abundance (Wilcoxon rank-sum one-sided test: P values < 2.2 × 10⁻¹⁶; Fig. 2C). Notably, G1 genes exhibited a canonical peak distribution pattern, primarily enriched in 3′-UTR and near stop codon (Fig. 2D). In contrast, G2 and G3 genes had a relatively higher number of m⁶A peaks in their coding sequence (CDS) region compared with G1 (Fig. 2D; Kolmogorov–Smirnov test: G1 vs. G2: P = 1.06 × 10⁻⁴; G1 vs. G3: P = 4.38 × 10⁻³; G2 vs. G3: P = 6.38 × 10⁻²). We further assessed the genomic characteristics of these genes because they might influence the m⁶A signal distribution. Specifically, because fragmented RNA for m⁶A MeRIP is approximately 200 nucleotides (nt), peak distribution in meta-gene analysis over the genomic regions will shift toward the 5′ end of the gene in cases with extremely short CDS. There were negligible differences in the length of 5′-UTR, whole CDS region, and 3′-UTR between the three groups (Supplementary Fig. S2F), indicating that the noncanonical peak distribution of G2 and G3 genes is not a technical artifact. G1 genes contained the highest number of exons with the shortest length of the last CDS, whereas G3 genes harbored the fewest number of exons with the most extended length of the last CDS (Fig. 2E). This was exemplified by genes RBM22, PLA2G15, and GIMAP4 in the G1, G2, and G3 groups, respectively (Fig. 2F). The Kyoto Encyclopedia of Genes and Genomes pathway analysis revealed that G1 genes were linked to terms related to cancer cell survival and metastasis, such as “focal adhesion,” “adherens junction,” “hypoxia-inducible factor-1 (HIF-1),” and “NOTCH signaling pathways” (Fig. 2G; Supplementary Table S2). Interestingly, the top enriched signature in both G2 and G3 genes was herpes simplex virus 1 (HSV-1) infection, which mainly consisted of zinc finger (ZNF) genes (Fig. 2G; Supplementary Table S2). The ZNF genes are a large class of putative transcription factors featured by tandem ZNF structures (20), displaying a more profound shift of peak distribution toward the 5′ of CDS (Supplementary Fig. S2G). There were no statistically significant differences in CDS length between ZNF genes and the remaining G2 and G3 genes (Supplementary Fig. S2H), further reinforcing the notion that the noncanonical distribution was not a technical artifact. The unique association pattern with m⁶A regulators, combined with the noncanonical m⁶A peak distribution, suggests that this gene group might undergo a yet-to-be-identified regulation by m⁶A.

To explore the heterogeneity of m⁶A modification in lung patient tumors, we performed consensus clustering based on m⁶A levels from the top 20% of methylated genes, ranked by interquartile range (IQR). This led to the identification of five m⁶A subtypes (P1–5; Fig. 3A). We applied a similar consensus clustering approach to RNA and protein levels. Consistent with the positive correlations between RNA and protein abundance for a significant proportion of genes, there was significant overlap between the RNA and protein subtypes (Supplementary Fig. S3A). However, no significant overlap was found between the m⁶A subtypes and those of RNA or protein (Supplementary Fig. S3A; Supplementary Table S3). Furthermore, two RNA expression-based LUAD classification methods were compared with the m⁶A subtypes (21, 22). In the four-cluster (C1–C4) LUAD classification system, C1 and C2 were characterized by high expression of proliferation-associated genes and similarities with small-cell lung cancer. In contrast, the C4 subclass was akin to normal lung and associated with better survival. The C3 subclass displayed features of both C4 and normal lung (21). In the “bronchioid/magnoid/squamoid” classification system, the bronchioid subclass was associated with early-stage tumors and favorable outcome (22).

As shown in Supplementary Fig. S1A, the UHN-PLCME cohort comprises 71.7% early-stage tumors. In line with the characteristics of the aforementioned classifications, the majority of the tumors in the UHN-PLCME cohort was categorized into the C3/C4 subclasses (92.2%; 47 out of 51) and the bronchioid subclass (58.8%; 30 out of 51; Supplementary Fig. S3B and S3C; Supplementary Table S3). These classifications demonstrated significant associations with our RNA and protein subtypings (Supplementary Fig. S3B and S3C). In contrast, when comparing these classifications with the m⁶A subtypes, we observed no significant overlap (Supplementary Fig. S3B and S3C). This lack of correlation highlights the uniqueness of m⁶A subtyping, underscoring potential distinct molecular mechanisms that might be overlooked by conventional mRNA expression-based classifications.

Notably, patients in the P2 group demonstrated the lowest m⁶A levels compared with other groups. This observation held true whether considering the top 20% of methylated genes by IQR ranking (Fig. 3A and 3B) or all methylated genes (Supplementary Fig. S3D). The P2 subtype exhibited relatively low RNA and protein abundances of the core writer complex component WTAP and relatively high protein levels of eraser ALKBH5, along with multiple reader complex proteins (Supplementary Fig. S3E). Intriguingly, the P2 subtype was associated with a significantly better survival rate compared with other groups (Fig. 3C; Supplementary Fig. S3F and S3G). A comparative analysis of the RNA expression profiles between the P2 subtype and the rest of the groups revealed 246 upregulated and 199 downregulated genes (Fig. 3D). A functional enrichment analysis of these differentially expressed genes identified “cell cycle” as the top enriched Gene Ontology/Biological Process (GO/BP) and Kyoto Encyclopedia of Genes and Genomes pathways for the upregulated genes in the P2 group (Supplementary Table S3), whereas the downregulated genes were enriched in immune-related pathways (Fig. 3E; Supplementary Table S3). The majority of the downregulated genes contributing to the enrichment of immune-related pathways were immunoglobulin molecules (Supplementary Table S3). Taken together, we identified a unique m⁶A subtype associated with a favorable clinical outcome.

To understand the clinical implications of individually methylated genes, we examined the correlations between mRNA, protein, and m⁶A levels and various clinicopathologic features. These features, which include sex, smoking status, and tumor stage, are documented for all the 51 tumors studied. Each of the transcriptomic, proteomic, and epitranscriptomic datasets contained genes uniquely correlated with a specific feature (Supplementary Fig. S4A and S4B). The m⁶A levels of TBRG4, SEC24B, and TXLNA were most significantly associated with sex, smoking status, and tumor stage, respectively (Supplementary Fig. S4A and S4B). These associations were generally independent of the mRNA and protein levels of corresponding genes (Supplementary Fig. S4B). We then extended the analysis to clinical outcomes. The cohort included 45 patients with clinical follow-up data, 82% of whom had early-stage cancer (stages I and II; Supplementary Table S1). A total of 635 methylated genes with complementary RNA and protein data available in all the 45 patients were included for the analysis. The m⁶A levels of 40 genes were associated with overall survival (P < 0.05), with BLVRA being the most significant (Fig. 4A and 4B; log-rank test: P = 3.4 × 10⁻⁴). Notably, neither the RNA nor the protein abundances of BLVRA served as prognostic indicators in this cohort (Fig. 4B). Consistently, BLVRA mRNA level was not associated with patient survival in the TCGA LUAD cohort (Supplementary Fig. S4C). In line with these findings, BLVRA m⁶A level was not correlated with its mRNA or protein abundance (Fig. 4C). However, BLVRA protein and mRNA levels were highly correlated (Supplementary Fig. S4D), indicating that the steady-state abundance of BLVRA protein was mainly driven by its mRNA.

BLVRA m⁶A profile demonstrated a canonical signal distribution, with a single prominent peak near the stop codon (8), suggesting that this singular peak likely influences BLVRA m⁶A level as a prognostic factor. Indeed, the m⁶A peak intensity around the stop codon (peak region IP/INPUT ratio) was highly correlated with the m⁶A level of the entire transcript (Fig. 4D), and this association was also significant about survival (Supplementary Fig. S4E). Near the peak summit, we identified a consensus motif (Fig. 4E). To validate this site, we selected three tumor samples representing low, medium, and high m⁶A levels, to quantify the m⁶A level at the “A” site within the motif using the single-base elongation- and ligation-based quantitative polymerase chain reaction (qPCR) amplification method (termed “SELECT”; Fig. 4D–4F; ref. 23). Methylation levels at this “A” site aligned well with that at both peak and whole-transcript levels in these tumor samples (Fig. 4F).

Using the SELECT method, we further assessed m⁶A levels at the “A” site in a validation cohort of 80 patients with early-stage LUAD with transcriptomic profiling and clinical follow-up data available (24). Consistent with the UHN-PLCME cohort, only the high levels of m⁶A, but not the mRNA levels of BLVRA, predicted poor survival in the validation cohort (Fig. 4G). No significant correlation was observed between the m⁶A and mRNA levels of BLVRA (Fig. 4H). In a multivariate Cox regression analysis, after adjustment for age, sex, smoking status, and tumor stage, high BLVRA m⁶A level remained significantly associated with increased risk of cancer-related death [UHN-PLCME cohort: hazard ratio (HR) (95% confidence interval (CI)): 10.47 (1.02–107.7), P = 0.048; validation cohort: HR (95% CI): 2.48 (1.08–5.70), P = 0.033; Supplementary Table S4]. Therefore, BLVRA’s m⁶A level, being the first of its kind identified, emerges as a potential prognostic factor for patients with LUAD.

To understand why high BLVRA m⁶A levels correlate with shorter survival in patients with LUAD, we conducted a differential RNA expression analysis on tumor samples dichotomized by high and low BLVRA m⁶A levels in our UHN-PLCME cohort. Gene set enrichment analysis identified that epithelial-to-mesenchymal transition (EMT) was the most significantly enriched signature in genes upregulated in tumors with high BLVRA m⁶A levels (Supplementary Fig. S5A). This is particularly relevant as EMT is known to be associated with drug resistance, immune evasion, metastasis, and poor outcome during cancer progression (25–27). Next, we analyzed BLVRA m⁶A levels under two well-established EMT stimuli, including transforming growth factor beta (TGF-β) and hepatocyte growth factor (HGF; ref. 28). BLVRA m⁶A modification was significantly elevated upon HGF treatment with negligible changes in its mRNA or protein levels in NSCLC cells, whereas TGF-β had no impact on BLVRA m⁶A levels (Supplementary Fig. S5B and S5C). The dCas13b-ALKBH5 fusion protein system was utilized to specifically remove the m⁶A modification from BLVRA transcripts in HGF-treated A549 cells (29). Despite a significant reduction in m⁶A levels (Supplementary Fig. S5D), we did not observe substantial changes in the abundance of BLVRA mRNA or protein (Supplementary Fig. S5D and S5E), suggesting that while BLVRA m⁶A level may serve as a biomarker for EMT stimuli and clinical outcome in LUAD, it does not substantially influence RNA stability or translation efficiency.

To identify the unique characteristics of m⁶A modifications in LUAD between tumor and NL samples, we analyzed the m⁶A MeRIP-seq data from both sample types. The median number of m⁶A peaks detected in NL tissues was 15,442, significantly higher than the 11,618 in tumor samples [Fig. 5A (left); Wilcoxon rank-sum one-sided test: P = 1.16 × 10⁻⁶]. Every RNA category, including protein-coding RNAs and lincRNAs, displayed significantly higher numbers of m⁶A peaks in NL tissues (Supplementary Fig. S6A). A much higher variation in peak number was observed in tumors compared with NL tissues [Fig. 5A (left); Levene’s test: P = 3.5 × 10⁻³]. To further investigate the variation, a pair-wise intragroup comparison was performed in NL and tumor samples. High overlapping rate (>68.6%) of m⁶A peaks was observed in NL tissues, whereas that of tumor samples was significantly lower, ranging from 38.4% to 81.4% [Fig. 5A (right); Wilcoxon rank-sum one-sided test: P < 2.2 × 10⁻¹⁶], indicating high inter-tumor heterogeneity. The majority of m⁶A peaks (12,562) were detected in both NL and tumor samples, which was defined as “common peaks” (Fig. 5B; hypergeometric test: P < 2.2 × 10⁻¹⁶). There were 1,607 and 6,519 m⁶A peaks exclusively detected in NL and tumor samples, which were termed as “NL-specific peaks” and “tumor-specific peaks,” respectively (Fig. 5B). The preponderance of tumor-specific peaks likely stems from the larger sample size in tumor and the noted high intertumor heterogeneity. Indeed, only 10% of the tumor-specific peaks appeared in >20% of tumor samples, whereas half of the NL-specific and 80% of common peaks were observed in >20% of their respective samples (Supplementary Fig. S6B). A negative correlation between m⁶A and mRNA levels was observed consistently within each sample, with correlation coefficients in NL ranging from −0.45 to −0.39 (mean = −0.42; Supplementary Fig. S6C, left), similar to that observed in human embryonic stem cells (17). In contrast, correlations in tumor samples were less consistent, showing greater variability [Supplementary Fig. S6C, left; −0.33 (−0.45, −0.19); Wilcoxon rank-sum one-sided test: P = 2.44 × 10⁻⁵]. Inverse correlations were evident in NL samples when analyzing individual genes across the cohort, but this pattern was nearly absent in tumor samples (Supplementary Fig. S6C, right).

The m⁶A levels in tumors were significantly lower than those in NL tissues (one-sided Z-test: P = 3.04 × 10⁻⁸) and displayed greater variance (Kolmogorov–Smirnov test: P < 2.2 × 10⁻¹⁶; Fig. 5C; Supplementary Fig. S6D), consistent with observation at the peak level (Fig. 5A). Differential methylation analysis identified 652 differentially methylated genes (DMG), including 430 hypo-methylated (20 lincRNAs) and 222 hyper-methylated (4 lincRNAs) genes in tumors (Fig. 5D and 5E; Supplementary Fig. S6E and S6F; Supplementary Table S5). To ensure the robustness of these findings, we performed a validation analysis by randomly subsampling 10 tumors from the 51 LUAD samples 100 times and comparing these subsets with the 10 NL samples. This comparison demonstrated that, on average, more than 70% of the hyper- and hypo-methylated genes identified in the full dataset were consistently observed across the subsampled comparisons (Supplementary Fig. S6G). Additionally, targeted validation was performed on three representative hypo-methylated (AP1M1, KLRG2, and ZNF688) and three representative hyper-methylated (MAT2A, EML4, and R3HDM2) genes, using the m⁶A MeRIP qPCR assay in an independent cohort comprising 15 LUAD tumor and 15 NL samples. We observed significantly lower m⁶A levels of the hypo-methylated genes and higher levels of the hyper-methylated genes within tumor samples compared with NL samples (Supplementary Fig. S6H).

Notably, 91.2% of hypo-methylated genes contained common peaks, whereas the proportion of hypermethylated genes with common peaks was only 59% (Supplementary Fig. S6I). Across the 51 tumor samples, more than 2,600 tumor-specific methylated genes were identified, which were specifically methylated in a subgroup of tumors (Supplementary Table S5), in line with the high inter-tumor heterogeneity observed at peak levels (Fig. 5A). These genes were most strongly enriched for pathways related to immunology, including “immunoglobulin receptor binding” and “antigen binding” (Supplementary Fig. S6J; Supplementary Table S5). Although the lung tissue samples used in this study consisted mainly of lung epithelial cells with high purity (>60% tumor cellularity), the bulk tissue samples were infiltrated with other cell types, such as immune cells. To test whether the enriched immunology signatures were related with tumor-specific immune infiltration, we assessed the abundance of six major immune cell subtypes (B cells, CD4 T cells, CD8 T cells, macrophages, neutrophils, and dendritic cells) in both NL and tumor samples in the UHN-PLCME and the TCGA LUAD cohorts (Supplementary Fig. S6K). B cells were the only cell type that showed significantly higher infiltration in tumors compared with NL tissues in both cohorts (Supplementary Fig. S6K; Wilcoxon rank-sum test: UHN-PLCME: P = 3.2 × 10⁻⁵, TCGA LUAD: P = 3.5 × 10⁻²). Consistently, tumor samples with high B-cell infiltration had significantly higher proportion of B-cell-specific genes with m⁶A peaks (Supplementary Fig. S6L; Wilcoxon rank-sum one-sided test: P = 1.2 × 10⁻²). Considering that immunoglobulin molecules primarily originate from B lineage cells (30), the elevated B-cell infiltration in the tumor microenvironment likely influences the observed tumor-specific m⁶A peaks.

Hypomethylation tends to occur in transcripts with low RNA abundance but high m⁶A levels, whereas hyper-methylation predominantly emerges in transcripts with high RNA abundance but low m⁶A levels (Supplementary Fig. S6M; Wilcoxon rank-sum one-sided test: P < 2.2 × 10⁻¹⁶). Neither hypo- nor hypermethylated genes showed notable RNA abundance variations between NL and tumor samples (Supplementary Fig. S6M; Wilcoxon rank-sum two-sided test: hypo: P = 0.18; hyper: P = 0.17). Specifically, 68.1% of hypomethylated (293/430) and 69.4% of hyper-methylated (154/222) genes did not show a significant difference at RNA levels between NL and tumor samples (Fig. 5D). These data indicate that the vast majority of the methylation differences between NL and tumor samples do not result in significant alteration in RNA abundance in tumors.

Among the 447 DMGs with no significant changes at RNA levels, EML4 is of particular interest: (i) It ranked second in hypermethylated genes (Fig. 5F) and has a single m⁶A peak near the stop codon (Fig. 5G), ideal for further functional validation, and (ii) EML4 is well-known for its role in the EML4-ALK fusion gene, which produces a transforming fusion tyrosine kinase, an oncogenic driver mutation present in approximately 5% of NSCLC cases (31). However, the independent function of EML4 in cancers, including LUAD, has largely been overlooked. Consistent with the UHN-PLCME cohort, there was no obvious upregulation in EML4 mRNA abundance in tumors compared with adjacent normal samples in the TCGA LUAD cohort [Wilcoxon rank-sum one-sided test: log₂(FC (tumor/normal)) = 0.53, P < 1.55 × 10⁻⁸; Supplementary Fig. S7A]. Aligning with the m⁶A level of the entire transcript, EML4 m⁶A peak intensity (peak region IP/INPUT ratio) in tumor samples was markedly elevated compared with NL tissues (Fig. 5H; Wilcoxon rank-sum one-sided test: P = 2.03 × 10⁻⁴).

While knocking down the key writer of m⁶A, METTL3, had a minimal effect on EML4 mRNA abundance, it resulted in a notable decline in methylation levels around the peak region and at the “A” site within the motif predicted to be closest to the peak summit, along with reduced EML4 protein abundance and translation efficiency (Fig. 5I). Accordingly, the site-specific removal of m⁶A from this “A” site on EML4 transcript, facilitated by the dCas13b-ALKBH5 fusion protein system, led to a significant reduction in EML4 protein abundance and translation efficiency, with a negligible effect on its mRNA level compared with the controls (Fig. 5J and 5K; ref. 29). Proteomic and transcriptomic data of 99 paired LUAD tumors and adjacent normal tissues from the Clinical Proteomic Tumor Analysis Consortium (CPTAC) LUAD cohort revealed a more pronounced upregulation of EML4 protein abundance in tumor samples than of mRNA abundance [protein: FC (tumor/normal) = 2.78, Wilcoxon rank-sum one-sided test: P = 1.27 × 10⁻¹²; mRNA: FC (tumor/normal) = 1.40, Wilcoxon rank-sum one-sided test: P = 7.35 × 10⁻⁹; Fig. 5L]. Moreover, 77.8% of LUAD tumors (77 of 99) exhibited an upregulation of EML4 protein in comparison to their corresponding adjacent normal tissues (Supplementary Table S6; ref. 3). Immunohistochemistry (IHC) also underscored the enhanced expression of EML4 protein in primary LUAD (n = 55, paired Student t test: P < 0.0001; Fig. 5M; Supplementary Fig. S7B).

We further set out to assess which m⁶A reader(s) mediate EML4 protein translation. Given that translation occurs in the cytoplasm, we examined six major m⁶A readers with cytoplasmic localization, including YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, and IGF2BP3 (5). Knockdown of IGF2BP1 and IGF2BP3 led to a significant decrease in EML4 protein abundance and translation efficiency, with negligible changes on EML4 mRNA levels (Supplementary Fig. S7C–S7E). In contrast, the depletion of other m⁶A readers did not significantly impact the abundance or translation efficiency of EML4 (Supplementary Fig. S7C–S7E). RNA-binding protein immunoprecipitation (RIP) qPCR assay confirmed the binding of IGF2BP1 and IGF2BP3 to EML4 mRNA (Supplementary Fig. S7F). Similar regulation patterns in EML4 translation efficiency, as measured by ribosome profiling (Ribo-seq; see Supplementary Methods), were observed in HeLa cells. Specifically, EML4 translation efficiency was significantly reduced upon METTL3 knockdown, whereas knockdown of YTHDF family proteins had no obvious effect [log₂(METTL3 knockdown/control FC) = −1.07; Supplementary Fig. S7G].

The protein-to-mRNA ratio served as a proxy for the translation rate, under the assumption that translation efficiency is the dominant posttranscriptional mechanism regulating protein abundance (32). In line with the translation-promoting function of EML4 m⁶A in cell line models, both EML4 protein and protein-to-mRNA ratio levels were significantly higher in tumors with high EML4 m⁶A levels compared with those with low levels in the UHN-PLCME cohort (Supplementary Fig. S7H). This reinforces the role of EML4 m⁶A in enhancing protein translation in LUAD. Taken together, these data suggest that m⁶A hypermethylation of EML4, catalyzed by classic METTL3-METTL14 methyltransferase complex and recognized by m⁶A readers IGF2BP1 and IGF2BP3, leads to enhanced protein translation and pervasive overexpression in primary LUAD.

We sought to delineate the functional roles of EML4 in LUAD tumorigenesis. Knockdown of EML4 using siRNA substantially decreased NSCLC cell (A549 and H358) migration and invasion in vitro, as assessed by a Transwell assay, with no notable changes in cell proliferation (Fig. 6A and 6B; Supplementary Fig. S8A and S8B). Similarly, compromised migration abilities were evident in EML4 CRISPR knockout (KO) cells, confirming the specificity of the effects of siRNAs and shRNAs (Fig. 6A–6C; Supplementary Fig. S8A, S8C, and S8D). Statistically significant yet modest effects on cell proliferation were observed with EML4 overexpression, with increased growth in H358 cells and decreased growth in A549 cells (Supplementary Fig. S8A, right). Nevertheless, significantly enhanced tumor cell migration abilities were observed in both EML4-overexpressing A549 and H358 cells (Fig. 6A and 6D; Supplementary Fig. S8E). Wound healing assay and live cell imaging further confirmed the influence of EML4 on cell motility (Fig. 6D and 6E; Supplementary Fig. S8F; Supplementary Videos S1–S5). We thus further investigated the role of EML4 in metastasis in vivo using mouse lung metastases models, established by tail vein injection of A549 and H358 cells into immunodeficient mice. EML4 knockdown strikingly reduced lung metastasis and prolonged survival (Fig. 6E–6G), whereas its overexpression resulted in increased tumor metastasis in the lungs of these mice (Fig. 6H). Collectively, these in vitro and in vivo results indicate that EML4 primarily promotes tumor cell migration, invasion, and ultimately metastasis in LUAD.

In line with its role in regulating cell motility, EML4 knockdown, KO, and overexpression induced pronounced cellular morphology changes in NSCLC cells (Fig. 6I; Supplementary Fig. S8G and S8H), suggesting a role in cytoskeleton regulation. Specifically, EML4 knockdown and KO resulted in impaired tumor cell spreading [Fig. 6I (left); Supplementary Fig. S8G], whereas cells with EML4 overexpression exhibited more expansive shapes compared with control cells [Fig. 6I (right)]. Interestingly, axon-like cellular protrusions were induced upon EML4 overexpression in a large proportion of H358 cells [Fig. 6I (right); Supplementary Fig. S8H]. As a member of the echinoderm microtubule-associated protein-like family, EML4 has been reported to bind directly to microtubules, colocalizing with microtubules in interphase in HeLa and U2OS cells and throughout all mitosis phases in HeLa cells (33–35). Yet, immunofluorescence (IF) staining revealed diffuse EML4 localization in the cytoplasm of interphase NSCLC cells (Fig. 6J; Supplementary Fig. S8I). EML4 was only colocalized with α-tubulin in the mitotic spindle during mitosis and in the midbody during cytokinesis in A549 cells (Supplementary Fig. S8J). This suggests that EML4 promotes lung tumor cell motility without directly binding to microtubules. Moreover, in a portion of NSCLC cells, EML4 showed membrane localization, which increased upon its overexpression (Fig. 6J; Supplementary Fig. S8I and S8K).

To further elucidate the molecular mechanism underlying the pro-metastatic role of EML4, we performed RNA-sequencing (RNA-seq) in EML4-modified and control A549 and H358 cells (Fig. 6K; Supplementary Fig. S8L; Supplementary Table S7). We also compared the gene expression profiles of LUAD tumors with high and low EML4 protein expression levels from the CPTAC cohort (Fig. 6K; Supplementary Table S7). Pathway enrichment analysis revealed that “membrane” and “cell periphery” were common top ranked terms. They were depleted upon EML4 knockdown and enriched upon EML4 overexpression in both cell lines tested and also enriched in LUAD tumors with high EML4 protein expression (Fig. 6K; Supplementary Fig. S8L; Supplementary Table S7). Other commonly enriched terms identified in both cell line and clinical samples related to cellular protrusions and neuronal synapses (e.g., synapse, axon, and dendrite; Fig. 6K; Supplementary Fig. S8L; Supplementary Table S7). Those terms were closely associated with cell shape control and cell motility and migration, suggesting that the transcriptomic profiles mirror the phenotypic changes induced by EML4. However, an overlapping analysis of differentially expressed genes from the RNA-seq data revealed only a few common differentially expressed genes (Supplementary Table S7), such as synaptotagmin 12 (SYT12) and serpin family E member 1 (SERPINE1; Supplementary Fig. S8M)–both known to promote tumorigenesis and metastasis in multiple types of cancers including NSCLC (36–38). The consistency at pathway level, but not at individual gene level, suggests that EML4 may promote cell motility through mechanisms other than direct transcriptional regulation.

To further understand the mechanisms underlying EML4 function, we performed co-immunoprecipitation coupled with mass spectrometry (MS) (Co-IP/MS) in A549 cells, identifying 26 potential interacting proteins unique to EML4 and not found in IgG control precipitates (Supplementary Table S8). Of particular interest is APRC1A, a component of the actin-related protein 2/3 (Arp2/3) complex known for its crucial role in cell motility and metastasis and its cytoplasmic and membrane localization as EML4 (39). The Arp2/3 protein complex is a central modulator of actin cytoskeleton and controls the nucleation of actin filament branches and actin polymerization, thereby regulating a variety of cellular structures and functions that involve the actin cytoskeleton (39, 40). The Arp2/3 complex is crucial for the formation and dynamics of lamellipodia, a major type of membrane protrusions and the crucial force leading to mesenchymal cell migration (41). Overactivation of the Arp2/3 complex drives metastasis and is associated with aggressiveness and tumor progression in varieties of cancers (40, 41). Western blot analysis against ARPC1A using immunoprecipitated Flag-tagged EML4-associated protein complex confirmed the interaction between EML4 and ARPC1A in NSCLC cells (Fig. 6L). IF staining demonstrated colocalization of EML4 and ARPC1A in the cytoplasm as well as at the edge of lamellipodia (Fig. 6M). Functional studies further revealed that overexpression of EML4 significantly enhanced ARPC1A enrichment at the edge of lamellipodia (Fig. 6N). In line with the known function of the Arp2/3 complex, lamellipodia formation was significantly reduced in EML4 KO cells (Supplementary Fig. S9A), whereas EML4 overexpression significantly increased lamellipodia formation (Supplementary Fig. S9B). Reciprocally, silencing of ARPC1A phenocopied the effect of EML4 loss, inducing similar morphology changes and inhibiting both cell migration and lamellipodia formation (Supplementary Fig. S9C–S9F). Furthermore, the cell migration and lamellipodia formation enhancements induced by EML4 overexpression were significantly reduced upon ARPC1A depletion (Fig. 6O–6Q; Supplementary Fig. S9G). These data collectively suggest that EML4 enhances the formation of lamellipodia and the migration capabilities of tumor cells through interacting with ARPC1A in LUAD.

The m⁶A-mediated widespread overexpression of EML4 protein in primary LUAD tumors, combined with the pro-metastatic functions of EML4, prompted us to explore the clinical correlation between EML4 protein abundance and pathological characteristics related to metastasis. Elevated EML4 protein expression in primary LUAD tumors significantly correlated with lymph node status in both the CPTAC cohort (Pearson’s χ² test; P = 0.049) and our internal cohort (Pearson’s χ² test; P = 0.029; Supplementary Table S6). In contrast, no significant association between EML4 mRNA expression in primary LUAD tumors and lymph node status was observed (Pearson’s χ² test; P = 0.331; Supplementary Table S6). Furthermore, EML4 protein abundance in metastases (lymph node metastases: n = 36; distant metastases: n = 3) was significantly higher than that in matched primary LUAD tumors (n = 39, paired Student t test; P < 0.0001; Fig. 7A). An increase in EML4 protein level was evident in 79.5% of metastases (31 of 39) compared with their matched primary tumors (Supplementary Table S6). EML4 protein abundance in both primary tumors and lymph node metastases was significantly correlated with the number of positive lymph nodes and the number of positive lymph node stations (Supplementary Fig. S10A and S10B). In addition, EML4 protein expression in widespread lymph node metastases (N2) was significantly higher than that in limited lymph node metastases (N1; Supplementary Fig. S10C). However, METTL3 protein abundance was not significantly elevated in metastases compared with primary tumors as EML4 (Supplementary Fig. S10D and S10E; Supplementary Table S6). There was also no significant correlation between METTL3 and EML4 protein abundance either in primary LUAD tumors (the CPTAC cohort: Supplementary Fig. S10F, left; the internal cohort: Supplementary Fig. S10F, middle) or in metastases (Supplementary Fig. S10F, right), suggesting that EML4 m⁶A hypermethylation is not caused by METTL3 protein upregulation. Taken together, these analyses in clinical specimens underscore EML4 as a novel driver of LUAD metastatic progression.

Our findings suggest EML4 as a promising therapeutic target for preventing LUAD metastases. Because silencing METTL3 efficiently reduced EML4’s m⁶A modification and protein abundance in vitro (Fig. 5I), we hypothesized that suppressing METTL3 would reduce EML4 m⁶A modification and protein abundance, thereby inhibiting metastasis in LUAD. Consistent with our expectations, METTL3 knockdown significantly decreased tumor cell migration in vitro (Supplementary Fig. S10G) and metastasis in vivo (Fig. 7B and 7C; Supplementary Fig. S10H). Similarly, the METTL3 inhibitor STM2457 markedly reduced both EML4 m⁶A level and its protein expression (Fig. 7D and 7E; Supplementary Fig. S10I; ref. 42). This reduction coincided with a decline in tumor cell migration (Fig. 7E; Supplementary Fig. S10J). Importantly, overexpression of EML4 significantly restored the cell migration capability in METTL3-knockdown and STM2457-treated cells (Fig. 7E and 7F; Supplementary Fig. S10J). These data suggest that the cell migration-inhibitory effects of METTL3 knockdown and STM2457 were primarily mediated through the modulation of EML4 m⁶A and protein abundance.

Next, we assessed the efficacy of STM2457 in preventing metastases in vivo. We developed a mouse lung metastasis model by intravenously injecting 1 × 10⁶ A549 cells into each nonobese diabetic/severe combined immunodeficiency (NOD/SCID) mouse. Starting from the day of tumor cell injection, mice were randomly assigned to receive either a control (vehicle) or STM2457 treatment. This treatment was administered daily for an initial five consecutive days, followed by every 3 days for a duration of 36 days (Fig. 7G). We evaluated mouse lung metastases 2 months post the tail vein injection. Our data showed that STM2457 treatment significantly alleviated lung metastases compared with the vehicle group (P < 0.01; Fig. 7H). STM2457 treatment demonstrated no obvious toxicity as evaluated by mice body weight and hematology parameters (Fig. 7H; Supplementary Fig. S11A and S11B). Thus, targeting EML4 through METTL3 inhibition may offer a feasible and promising therapeutic strategy for suppressing LUAD metastases.

Utilizing the low-input m⁶A MeRIP-seq protocol, we generated the epitranscriptome of 10 NL tissues and 51 LUAD tumors from a well-characterized clinical cohort. This dataset afforded an unparalleled avenue to investigate genome-wide m⁶A regulation and in relation to clinicopathologic characteristics. Notably, tumors displayed a global decrease and higher variance of m⁶A levels compared with NL tissues. While the m⁶A methyltransferase METTL3 is significantly upregulated in LUAD, the methyltransferase core subunit WTAP is significantly downregulated. Previous studies indicated that disruption of WTAP impacts the m⁶A levels to a greater extent compared with that of METTL3 and METTL14 (19). Thus, the downregulation of WTAP might contribute to the global hypomethylation in tumor samples. Another study reported that histone modification H3K36me3 recruits the m⁶A methyltransferase complex to deposit m⁶A on nascent RNAs co-transcriptionally in the nucleus. SET domain containing 2 (SETD2), the histone H3K36me3 methyltransferase, is mutated in 22% of LUAD in the TCGA cohort (2, 43), which might further contribute to m⁶A hypomethylation in tumors.

m⁶A methylation is a finely tuned process governed by a robust mechanism. The cornerstone of this mechanism is the m⁶A methyltransferase complex, comprising the core component METTL3/METTL14/WTAP and adaptors like VIRMA, CBLL1, and ZC3H13). This complex is responsible for the deposition of m⁶A, while erasers such as ALKBH5 and FTO remove these modifications. A hallmark of m⁶A is its typical enrichment near the stop codon. While the G1 gene group aligns well with this model, group G3 deviates, suggesting secondary effects or even uncharacterized methylation/demethylation mechanisms. Deepening the complexity, various transcription factors, such as ZFP217, SMAD2/3, and CEBPZ, and exon junction complex can either recruit or repel the m⁶A methyltransferase complex to specific chromatin sites under specific conditions, thereby shaping m⁶A topologies (44–46). In addition, writers and erasers may have crosstalk. A previous study discovered a mutual regulatory relationship between METTL14 and ALKBH5 (47). Furthermore, because FTO preferentially demethylates nuclear m⁶A, the spatial regulation may result in different accessibility to its target m⁶A (48). Thus, these identified and potentially uncharacterized regulatory pathways further complicate the regulatory network, leading to disparities in m⁶A levels in specific genes when compared with the abundance of their regulators.

By integrating multiomics and survival data, we discovered BLVRA m⁶A level is associated with patient survival independent of the gene’s mRNA and protein levels. We showed that BLVRA m⁶A was induced by HGF treatment an EMT-associated stimuli in LUAD cells. Other environmental stimuli may also exist to induce BLVRA m⁶A. However, BLVRA m⁶A does not substantially influence RNA stability or translation efficiency, two major functions of m⁶A modifications. Further study is warranted to investigate the functional impact of BLVRA m⁶A.

We observed negative correlations between m⁶A and RNA levels at sample levels in both NL and tumor samples and at gene levels in NL samples. Whether this phenomenon is caused by the RNA-destabilizing function of m⁶A warrants further investigation (49). Our differential methylation analysis spotlighted EML4 as a notably hyper-methylated gene, underlining its significant pro-metastatic capacity in lung cancer cells. While EML4 is well-known for its EML4-ALK translocation (31), emphasis has largely been on the role of ALK tyrosine kinase in driving the proliferation and growth of the tumor cells, and the function of EML4 remains largely unexplored. One reason is the minimal variation in EML4 mRNA levels observed between tumors and normal samples across various cancer types, including LUAD. Our study revealed an overexpression of EML4 protein in primary LUAD tumors, a process driven by robust m⁶A modification, promoting its translation. Data from both the CPTAC LUAD cohort and our in-house IHC revealed that EML4 protein levels in primary tumors were approximately three- to fourfold higher than in normal samples. Despite its varied expression in primary tumors, mirroring the inherent intratumoral heterogeneity in LUAD (50), metastases predominantly showed a consistent higher EML4 protein expression. Furthermore, EML4 protein abundance correlates with the degree of lymph node metastases, hinting at its role in tumor aggression. This pattern suggests that cells with elevated EML4 protein might be the primary source in invasion and metastatic spread. EML4 exhibited a distinct membrane localization in LUAD, which was enhanced upon EML4 overexpression along with boosted cell migration abilities. Transcriptomic profiling of EML4-modulated cells pinpointed pathways related to cellular membranes and membrane-formed structures, such as cell protrusions. We further characterized that the migration-promoting effect by EML4 was dependent on an EML4-interacting protein, ARPC1A, which is one of seven subunits of the Arp2/3 protein complex. The Arp2/3 protein complex governs the organization of actin cytoskeleton and drives mesenchymal cell migration by modulating the formation and dynamics of lamellipodia (39–41). The molecular mechanism of how EML4 regulates the activity of the Arp2/3 protein complex such as actin polymerization and nucleation of actin filaments warrants further investigation.

Metastasis is the leading cause of cancer-related deaths. Therefore, preventing tumor progression and metastasis is a pivotal strategy to bolster NSCLC outcome. We have pinpointed EML4 as a novel metastatic driver in LUAD, advocating for therapeutic interventions targeting EML4. Our in vitro and in vivo studies confirmed the viability of such an approach, in which METTL3 small-molecule inhibitor mitigated EML4 m⁶A marking, effectively suppressing its protein translation and potently reducing tumor metastasis in A549 cells with relatively high EML4 m⁶A modification. Whereas existing research underscores the therapeutic promise of METTL3 inhibitors for tumor growth reduction (51, 52), our findings suggest that these inhibitors could be particularly effective against LUAD metastasis, especially in cases with prominent EML4 m⁶A modification. Whereas silencing of EML4 has very minimal, if any, impact on cell proliferation, inhibition of METTL3 results in noticeable reduction of lung cancer cell growth. The profound reduction of tumor metastases upon METTL3 knockdown or pharmacologic inhibition is likely because of combinatorial effects on both tumor growth and metastasis. Our study provides strong rationale for developing strategies, such as small-molecule inhibitors, to directly target wild-type EML4 in patients with high EML4 m⁶A and protein levels to prevent or reduce tumor metastasis. In addition, METTL3 small-molecule inhibitor has been reported to regulate the differentiation and stemness of cancer cells in an m⁶A-dependent manner (42, 53). Both differentiation and stemness play critical roles in metastasis (54). Therefore, while we demonstrate EML4 as a primary target of METTL3 in driving tumor migration and metastasis in our experimental research models, the additional effects of METTL3 inhibition on tumor biology contributing to metastasis in LUAD warrant further investigation.

As an antibody-based immunoprecipitation method, m⁶A MeRIP-seq is sensitive to the quality of the antibody used. Cautions must be taken to mitigate potential batch effects. The utility of current nonantibody-based genome-wide m⁶A quantification methods, such as glyoxal- and nitrite-mediated deamination of unmethylated adenosines (GLORI), m⁶A-selective allyl chemical labeling and sequencing (m⁶A-SAC-seq), and evolved TadA-assisted N6-methyladenosine sequencing (eTAM-seq; ref. 55), has yet to be tested in patient tumors with limited RNA material. Therefore, there is significant interest in the continuous development of robust and precise methods for m⁶A measurement.

The tumor heterogeneity, including intratumoral heterogeneity, tumor microenvironment heterogeneity, and inter-tumor heterogeneity, of the epitranscriptomic marks is underexplored (56). In our study, the tumor-specific m⁶A methylated genes were enriched in immune-related pathways. Some of these genes, such as PAX5 and CD79B, were specifically expressed in immune cells, prompting the question of whether m⁶A modifications in those genes are due to increased immune cell infiltration in LUAD. On the other hand, some of these genes, such as HLA genes, were strongly expressed in both LUAD and immune cells. Whether those tumor-specific m⁶A modifications originated from LUAD or immune cells could not be adequately addressed with bulk-level m⁶A profiling. Furthermore, it is possible that LUAD subclones exist that possess unique pro-tumorigenic or pro-metastatic epitranscriptomic modifications and are more prone to malignant transformation, metastasis, drug resistance, and immune evasion. These epitranscriptomic modifications were probably underestimated in bulk-level analysis when the quantity of these malignant subclones is small. With the aid of single-cell technologies for the analyses of the epitranscriptomic marks (57, 58), these questions could be addressed in the near future.

In summary, our work pioneered an m⁶A epitranscriptomic blueprint in LUAD. By aligning it with corresponding transcriptomic, proteomic, and clinical data, we offer an enriched dataset that unearths the intricate dynamics of m⁶A modifications in a clinical milieu, paving the way for novel therapeutic insights.

The NSCLC A549 (Cat. #CCL-185) and H358 (Cat. #CRL-5807) cell lines were obtained from the ATCC and maintained in RPMI 1640 medium (Cat. #350-000-CL, Wisent), supplemented with 10% FBS (Cat. #080150, Wisent) and 1% penicillin-streptomycin (Cat. #450-201-EL, Wisent). All cell lines were at early passages (n < 10) and tested for mycoplasma contamination once a month by the MycoAlert Mycoplasma Detection Kit (Cat. #LT07-118, Lonza). All cell lines were authenticated by short tandem repeat every 3 months by Wuhan Pricella Biotechnology Co., Ltd. (China).

For EML4 knockdown, METTL3 knockdown, and EML4 overexpression assays, A549 or H358 cells (overexpression: 5 × 10⁵ cells per mouse; knockdown: 1 × 10⁶ cells per mouse) were injected intravenously through the tail vein into each 4- to 6-week-old male NOD/SCID (Strain #: 001303, The Jackson Laboratory) or NOD/SCID gamma (NSG) mouse (Cat. NO. NM-NSG-001, Shanghai Model Organisms Center, Inc., China). The humane intervention points were closely monitored. The mice were sacrificed ∼8 weeks after tail vein injection. For survival assay, EML4-knockdown and control H358 cells were injected intravenously through the tail vein into each 5-week-old male NSG mouse (1 × 10⁶ cells per mouse). The health conditions of the mice were closely monitored. Mice were euthanized when the humane intervention points were met (e.g., suppressed activity, being unresponsive to touch, marked hunched posture, dehydration, rough hair coat, dyspnea, and weight loss of 30% or more from maximum body weight). For STM2457 treatment assay, the mouse lung metastasis model was first established by tail vein injection of 1 × 10⁶ A549 cells into each ∼5-week-old male NOD/SCID mouse. On the same day of tumor cell injection, the mice were randomly assigned for treatment with vehicle control or 50 mg/kg STM2457 (Cat. #S9870, Selleck) daily for five consecutive days and then every 3 days for 36 days. STM2457 was freshly prepared at each time of treatment. To prepare STM2457 working solution, 10 mg STM2457 was dissolved by adding each solvent one by one: 100 μL DMSO, 800 μL PEG300, 100 μL Tween-80, and 1 mL H₂O. The mice were intraperitoneally injected with 50 mg/kg of STM2457 or the same volume of the solvent mix (vehicle group). The mice were sacrificed 8 weeks after tail vein injection. All the five lobes of mouse lungs were harvested and embedded in paraffin. Lung metastasis was assessed either by the number of pulmonary tumor nodules or the proportion of tumor metastases in the mouse lung when there were too many tumor nodules so that single nodule was indistinguishable based on hematoxylin and eosin staining. All mice experimental procedures were approved by the Animal Ethics Committee of the University Health Network (UHN, Canada) or Sun Yat-sen University (China).

LUAD snap-frozen specimens were collected and banked from patients who underwent surgical resection at the UHN using study protocols approved by the UHN Human Research Ethics Board (protocols 08-0230 and 09-0510). Informed written consent was obtained from each patient. Histology slides were reviewed, and cases were selected for >60% tumor cellularity. The 53 tumors were collected from the UHN-PLCME project during 2005 to 2010 (annotated as PLC samples; Supplementary Table S1; refs. 13–15). The validation cohort included a subset of 80 cases from a prognostic gene classifier validation project collected during 1996 to 2005 (annotated as MB samples; Supplementary Table S4; ref.24); these included only cases with residual mRNA still available after previous study. All patients with LUAD received no prior treatment for their disease (targeted therapy, immunotherapy, chemotherapy, or radiotherapy) before surgical resection. Quality assurance protocol included tumor histology review prior to RNA isolation. Diagnoses were based on the WHO classification for lung cancer (59). Patient and tumor demographic features and survival were retrieved from electronic patient records, as per approved protocols mentioned above. NL tissues were collected from patients with early-stage LUAD (IA) with >1 cm resection margin distance from the tumor edge. Fifteen NL tissue samples and 15 primary LUAD tissue samples used for DMG validation were obtained from the First Affiliated Hospital of Sun Yat-sen University. LUAD IHC specimens were obtained from patients with LUAD who underwent surgery at the First Affiliated Hospital of Sun Yat-sen University (Supplementary Table S6). The study protocol was approved by the Clinical Research Ethics Committee of the First Affiliated Hospital of Sun Yat-sen University. Written informed consent was obtained from each patient. All studies were conducted in accordance with the Declaration of Helsinki.

To generate stable Cas9-expressing cells, A549 cells were transduced with Cas9 lentivirus, which was generated by transfecting 293T cells with lentiCas9-Blast (Cat. #52962, Addgene) and packaging plasmids. After blasticidin selection, Cas9-expressing A549 cells were transduced with doxycycline-inducible sgRNA lentivirus and enriched by puromycin selection. EML4 knockdown was induced by adding 1 μg/mL doxycycline to Cas9-inducible sgRNA-expressing cells. After 72 hours, cells were harvested for Western blot analysis and cell migration assays.

Two EML4 sgRNAs (EML4 sgRNA-1: 5′-CUGUUAAGAUAUGGUUAUCG-3′, Cat. #LM2301206-A; EML4 sgRNA-2: 5′-CUAUUUUCCCGGUCGGAAGA-3′, Cat. #LM2301206-B; Supplementary Table S9) were purchased from Beijing Tsingke Biotech Co., Ltd. (China). Cas9 protein (GenCrispr NLS-Cas9-NLS Nuclease, Cat. #Z03469) was purchased from GenScript Biotech (China). Cas9/sgRNA ribonucleoprotein complexes were prepared by mixing Cas9 protein with two EML4 sgRNAs at 1:2.5:2.5 ratio and transfected into A549 and H358 cells by electroporation using Neon NxT Electroporation System (Thermo Fisher Scientific). Two EML4 KO clones (EML4-KO-1 and EML4-KO-2) were then isolated from the transfected cells by single-cell sorting and validated by Western blot analysis.

Total RNAs were purified with RNeasy Mini Kit (Cat. #74106; QIAGEN) or RNAprep Pure Cell/Bacteria Kit (Cat. #4992235/DP430, TIANGEN, China). Approximately 5 μg of total RNA was used to deplete rRNA using Ribo-Zero Gold rRNA Removal Kit (Cat. #MRZG12324, Illumina). rRNA-depleted RNAs were fragmented using NEBNext Magnesium RNA Fragmentation Module (Cat. #E6150S, NEB) and then reverse-transcribed by SuperScript II Reverse Transcriptase (Cat. #18064071, Thermo Fisher Scientific). Next, U-labeled second-stranded DNAs were synthetized with DNA Polymerase I (E. coli; Cat. #M0209L, NEB), RNase H (Cat. #M0297L, NEB) and dUTP Solution (Cat. #R0133, Thermo Fisher Scientific). After adding adapters, the fragments were purified with AMPure XP beads, were treated with Uracil-DNA Glycosylase (Cat. #M0280, NEB), and were then amplified by PCR. RNA-seq libraries were subjected for 2 × 150 bp paired-end sequencing (PE150) in duplicates at ∼40 million reads per library on an Illumina NovaSeq 6000 (LC-Bio Technology CO., Ltd., Hangzhou, China) following the vendor’s recommended protocol.

Cells were seeded on glass coverslips coated with fibronectin (Cat. #F0895, MilliporeSigma) in six-well plates. Approximately 16 to 24 hours after seeding, cells were fixed with 4% paraformaldehyde (Cat. #P6148, MilliporeSigma)/PBS for 20 minutes and then permeabilized by 0.2%Triton X-100 (Cat. #T9284, MilliporeSigma)/PBS for 5 minutes. Nonspecific binding sites were blocked with 3% FBS/PBS for 1 hour. Primary antibodies (EML4: Cat. #12156S, Cell Signaling Technology; α-tubulin: Cat. #3873S, Cell Signaling Technology; ARPC1A: Cat. #PA5-82149, Thermo Fisher Scientific; Flag: Cat. #F1804, MilliporeSigma) were applied to the cells at room temperature for 2 hours. After PBS wash, secondary antibodies with desired fluorescence probes were applied. Then cells were washed and mounted with ProLong Gold Antifade Mountant with DAPI (4′, 6-diamidino-2-phenylindole; Thermo Fisher Scientific). To stain F-actin, cells were fixed, permeabilized, blocked, and stained with Alexa Fluor Plus 555 phalloidin solution (Cat. #A30106; Thermo Fisher Scientific; dissolved in DMSO) for 30 minutes.

STM2457 was purchased from MedChemExpress (Cat. #HY-134836) and dissolved in DMSO as the stock solution. Twenty-four hours after seeding, A549 cells were treated with STM2457 (60 μmol/L) or DMSO for 3 days and harvested for Western blot analysis, cell growth analysis, and cell migration analysis.

To assess the correlation between genes’ m⁶A levels and mRNA abundances of m⁶A regulators, we calculated the Pearson correlation coefficient between log-transformed m⁶A levels and mRNA abundances of m⁶A regulators for individual genes (Fig. 2A). Hierarchical clustering was performed to group the methylated genes according to the correlation coefficient between genes’ m⁶A levels and mRNA abundances of m⁶A regulators. Consensus clustering was employed to identify subtypes of patients based on the m⁶A, RNA, and protein levels of the top 20% of the genes by IQR ranking after filtering out genes with missing measurements (Fig. 3A). The optimal number of clusters was chosen based on elbow and silhouette methods. The classification of previously reported RNA-based LUAD subtypes, including C1–C4 subtypes (21) and bronchioid/magnoid/squamoid subtypes (22), was determined by comparing the similarities between RNA levels of the reported marker genes from each UHN-PLCME sample and the centroid of the maker genes for each subtype (21, 22). The similarity was measured using the Pearson correlation coefficient.

To assess the correlations with clinicopathologic features, we focused on the genes with quantitative data for its mRNA, protein, and m⁶A levels. To minimize the sample size effects, especially for proteomic data, we used a more stringent cutoff that requires the protein to be detected in at least 40 samples, resulting in a dataset of 1,011 genes. Then, using the set of clinical characteristics including sex, smoking status, and tumor stage, an ANOVA was conducted for each gene and each of the three types of omics data. For survival analysis, we further filtered our dataset to 635 genes with all three types of omics data in 45 patient samples with documented survival information. We classified the patients into two even groups by median dichotomizing the abundance of mRNA, protein, and m⁶A levels and conducted a log-rank test for each gene.

All statistical tests were performed using R. P <0.05 was considered significant. Data are presented as mean ± SD. The Shapiro–Wilk test and Levene’s test were used to assess the normality and equality of variances, respectively. The Student t test (for small sample size), Z-test (for large sample size), and Wilcoxon rank-sum test (also known as the Mann–Whitney U test) were used to assess the difference between two samples. The ANOVA test and the Kruskal–Wallis test were used for comparison of more than two groups. The Kolmogorov–Smirnov (K-S) test was applied to compare two distributions. Pearson’s χ² test or Fisher’s exact test were used for enrichment analysis or independence analysis. Log-rank test was performed for survival analysis. The Benjamini–Hochberg correction was used to correct for multiple hypothesis testing.

m⁶A MeRIP-seq data is available upon request through EGA with accession number EGAS00001005524 (https://ega-archive.org/studies/EGAS00001005524). RNA-seq and RIP-seq data of NSCLC cells have been deposited in NCBI’s Gene Expression Omnibus (GEO) database under accession number GSE263727 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?&acc=GSE263727). The code and processed files to reproduce all analyses and figures are available at https://github.com/HansenHeLab/Epitranscriptomic_m6A_Profiling_in_LUAD_Data_and_Codes.

No disclosures were reported.

H.H. He, M.S. Tsao, J. Yu, S. Wang, and Y. Zeng designed studies; S. Wang performed the majority of the experiments; Y. Zeng and S. Wang performed data analysis; L. Zhu, M. Zhang, L. Zhou, W. Luo, L. Wang, Y. Liu, H. Zhu, X. Xu, P. Su, X. Zhang, M. Ahmed, W. Chen, M. Chen, S. Chen, M. Slobodyanyuk, Z. Xie, J. Guan, P.C. Boutros, Z. Ke, and Z. Cai assisted experiments or analysis; W. Zhang, A.A. Khan, S. Sakashita, N.Liu, N.-A. Pham, M.F. Moran, W. Yang, C. Cheng, and M.S. Tsao generated patient cohort coupled with transcriptome and proteomic dataset; S. Wang, Y. Zeng, H. Zhu, and H.H. He wrote the first draft of the manuscript; all authors revised and approved manuscript.

This work was supported by the Princess Margaret Cancer Foundation (886012001223 to H.H. He), Canada Foundation for Innovation and Ontario Research Fund (CFI32372 to H.H. He), NSERC discovery grants (498706 to H.H. He; 06305 to M.F. Moran), Canadian Cancer Society innovation grants (703800 to H.H. He), IMPACT grant (701595 to M.S. Tsao), CIHR operating grants (142246, 152863, 152864, and 159567 to H.H. He; 364778 to M.F. Moran), Foundation grant (FDN-148395 to M.S. Tsao), and Terry Fox New Frontiers Program Project Grant (1090 and 1124 to H.H. He). H.H. He was supported by TFRI New Investigator Awards and CIHR New Investigator Awards. H.H. He holds an OMIR Early Researcher Award. Y. Zeng was supported by a postdoctoral fellowship by Princess Margaret Cancer Foundation. We thank Jonathan R. Krieger (SPARC BioCentre, Hospital for Sick Children, Toronto Canada) for help with MS data analysis.