Three pairs of LUAD samples and tumor-adjacent normal tissues were obtained from the South University of Science and Technology Hospital. The LUAD samples and tumor-adjacent normal tissues were immediately collected and separated into centrifuge tubes after surgery. These samples were examined by experienced pathologists who confirmed the diagnosis of disease samples. All tissues were stored at −80°C until RNA isolation. This study was approved by the Ethics Committee of the South University of Science and Technology Hospital.

Total RNA was isolated and purified using TRIzol™ Reagent (Invitrogen™, cat. no15596018) following the instructions of the manufacturer. Total RNA was chemically broken into 200nt fragments. The product was added to anti-N6-Methyladenosine (m6A) antibody (Sigma-Aldrich, cat. no ABE572), protein A-magnetic beads (Invitrogen™, cat. no 10002D), protein G-magnetic beads (Invitrogen™, cat. no 10004D), which was then mixed and incubated overnight. RNA was extracted by phenol-chloroform lysate to afford the purified product. A portion of the initial fragmented RNA was used as the input RNA-seq library for MeRIP-seq. Ribosomal RNA was removed from the products. The first-strand cDNA was synthesized by the SMART principle, and enriched library fragments were amplified by PCR. Ultrafine RNA methylation m6A detection library fragments were obtained by DNA purification of magnetic bead library fragments. Bioptic Qsep100 Analyzer was used to perform a quality inspection on the library and detect whether the size distribution of the library is consistent with the theoretical size. NovaSeq high-throughput sequencing platform and PE150 sequencing mode were used.

Lung adenocarcinoma cell line A549 cells were cultivated in 1640 medium containing 10% fetal bovine serum (FBS) (ExCell Bio Inc., Shanghai, China) and 100 U/mL penicillin-streptomycin mixture (Gibco, Grand Island, NY, USA). The cells were maintained at 37°C in a humidified atmosphere of 95% air and 5% CO2. YTHDF1 siRNA was designed and constructed by Genechem Bio (Shanghai, China). In transfection, siRNAs were transfected by LipoRNAiMAX Transfection Reagent (Invitrogen) according to product instructions. Briefly, siRNAs were diluted in the 100 μL Opti-MEM (Gibco) at 10 nmol/L and mixed with 2.5 μL Lipofectamine RNAi Max. After incubation at room temperature for 18 min, the mixture was added to the cell culture medium, and the cells were collected after 48 h. The sequences of the targets were as follows: YTHDF1-Homo-370 sense 5’- GGAUACAGUUCAUGACAAUTT-3’, antisense 5’-AUUGUCAUGAACUGUAUCCTT-3’; YTHDF1-Homo-1076 sense 5’-CCUCCACCCAUAAAGCAUATT-3’, antisense UAUGCUUUAUGGGUGGAGGTT-3’; YTHDF1-Homo-1503 sense 5’-GCUCCAUUAAGUACUCCAUTT-3’, antisense 5’-AUGGAGUACUUAAUGGAGCTT-3’.

Total RNA was extracted from A549 cells using Trizol (Takara, Japan) according to the manufacturer’s instructions. For RT-qPCR, RNA was reverse transcribed to cDNA by using a reverse transcription kit (Takara, Japan). The levels of RNA transcripts were analyzed by the Jena qTOWER3 RT-qPCR system (Jena, Germany). The sequences of the primers were as follows: GAPDH forward primer, 5’-CTCCTCCTGTTCGACAGTCAGC-3’, reverse primer 5’-CCCAATACGACCAAATCCGTT-3’; YTHDF1 forward primer, 5’-GCACACAACCTCCATCTTCG-3’, reverse primer 5’-AACTGGTTCGCCCTCATTGT-3’; GPRIN1 forward primer 5’-AAAGCAGGCCGATTCCACTTC-3’, reverse primer 5’-TCCTTCCTCGGTGACACTGTA-3’.

The collected cell precipitates were coprecipitated by EZ-Magna RIP™ RNA-Binding Protein Immunoprecipitation Kit (Millipore 17-701) following the manufacturer’s procedure. Cells were lysed in a moderate volume RIP lysis buffer on ice for 5 min, and then the lysates were moved to -80°C overnight. Lysates were centrifuged at 4°C, and 12,000 rpm for 30 min to obtain the supernatant. 5μg YTHDF1 antibody or rabbit IgG was conjugated to protein A/G magnetic beads on a shaker for 30 min at room temperature. Magnetic beads containing antibodies were added to the supernatant, followed by rotating overnight at 4°C. TRIzol LS reagent (Invitrogen) was used to extract the RNA according to the manufacturer’s instructions. RT-qPCR follows the above method.

The Database for Annotation, Visualization and Integrated Discovery (DAVID) database (https://david.ncifc rf.gov/) was used in Gene Oncology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. 594 LUAD patients’ RNA-seq data were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). The expression levels of genes were analyzed using the R package DESeq2 [version 1.26.0]. Then, 526 LUAD patients including primary tumor, lymph node metastasis, distant metastasis, and clinical stage were selected for the further study. A univariate Cox regression model was adopted to calculate the hazard ratios (HRs) of the genes’ overall survival. Kaplan-Meier analysis was used to construct survival curves, and the Cox regression model was utilized to estimate the significance of the differences. ROC curves were constructed to estimate the accuracy of genes in predicting outcomes. YTHDF1-RIP-seq data was downloaded from the NCBI Gene Expression Omnibus (GEO) database under the accession number GSE136433.

MeRIP-seq analysis of tumor tissue and tumor-adjacent normal tissue from three LUAD patients was performed. The sequence logo showed the top three m6A motifs enriched from altered m6A peaks and the “GGAC” sequence was ranked in the first place ( Figure 1A ). The distribution of m6A peaks in normal and tumor tissues was further investigated. Figures 1B–D shows that the distribution of m6A signals around mRNA and lncRNA were comparable in the two classes of tissue samples. Also, m6A peaks in mRNA both in tumor and normal tissues were mainly enriched in the coding sequence near the stop codon. It was found that their m6A mRNA methylation was reduced globally in the tumor compartment compared with adjacent, normal control tissues ( Figure 1E ). In Figure 1F , by comparing with normal tissues, tumor tissues had 1192 significantly upregulated m6A peaks, corresponding to transcripts of 1092 genes, as well as 2849 significantly downregulated m6A peaks, representing transcripts of 2415 genes. The top 20 altered m6A peaks are listed in Table 1 . By analyzing the distribution of m6A peaks per mRNA, one m6A peak in the majority of mRNAs was found ( Figure 1G ). All altered m6A peaks were mapped to human chromosomes. Dysregulated m6A peaks were also observed in all chromosomes, particularly in chr1, chr2, chr16, chr17, and chr19 ( Figure 1H ).

To investigate the biological significance of m6A modification in LUAD, GO terms and KEGG pathway analyses of differential methylated genes were performed. For GO analyses, biological processes (BP), cellular components (CC), and molecular functions (MF) were taken into consideration. Figure 2A demonstrates the top 10 significantly enriched BP, CC, and MF of genes with downregulated m6A peaks, while GO analysis of genes goes along with upregulated m6A peaks as shown in Figure 2B . KEGG pathway analysis showed that genes with downregulated m6A peaks in LUAD were significantly associated with the pathways in cancer, Rap1 signaling pathway, and insulin resistance ( Figure 2C ). Genes with upregulated m6A peaks were significantly associated with the ascorbate and alternate metabolism, pentose and glucuronate interconversions, and sphingolipid signaling pathway ( Figure 2D ). Additionally, in terms of the fold change, differential methylated genes between tumor and normal were associated with the negative regulation of protein phosphorylation, glycoprotein biosynthetic process, and regulation of GTPase activity ( Figure 3 ).

The transcriptome profiles of tumor tissues versus normal tissues from three LUAD patients were detected using RNA-seq (MeRIP-seq input library). Compared to normal tissues, tumor tissues had 1726 up-regulated genes and 2011 down-regulated genes ( Figures 4A, B ). The top 20 differently expressed genes are listed in Table 2 . Then, all differential methylated m6A peaks with differential RNA levels (612) were divided into four groups by conjoint analysis for the MeRIP-seq and RNA-seq data. It was found that 177 hypermethylated m6A peaks in RNAs were significantly upregulated (128; hyper-up) or downregulated (49; hyper-down), while 435 hypomethylated m6A peaks in RNAs that were significantly upregulated (88; hypo-up) or downregulated (347; hypo-down) ( Figure 4C ). GO and KEGG pathway analyses for the biological significance of those genes (612) with differential methylated m6A peaks and synchronously differential expression were performed. GO analysis indicated that these genes were mainly enriched in the vasculogenesis (BP), an integral component of the plasma membrane (CC), and GTPase activator activity (MF) ( Figure 4D ). KEGG pathway analysis revealed that these genes were predominately enriched in the cGMP-PKG signaling pathway, focal adhesion, and Rap1 signaling pathway ( Figure 4E ). By contrast, GSEA enrichment of genes with significant changes in both the m6A modification and RNA levels showed that these genes were enriched in cell cycle phase transition, cell division, cellular response to DNA damage stimulus, and chromosome organization ( Figures 4F–I ).

To evaluate the clinical significance of m6A modification regulating gene expression, the TCGA database was explored. 428 differential expressed genes commonly expressed in this study and the TCGA database were identified ( Figure 5A ). The top 20 differently expressed genes are listed in Table 3 . The univariate Cox regression analysis of the expression of the top 20 differential expressed genes was performed to identify prognostic genes. The results demonstrated that the expression of 7 candidate genes was significantly related to the prognosis (p < 0.05) of LUAD patients ( Figure 5B ). Furthermore, we identified 2 genes (FCRL5 and GPRIN1) correlated with the prognosis of patients with LUAD. The overall survival curves of these 2 genes are shown in Figures 5C, D . This result showed that high expression of FCRL5 was associated with a better prognosis for LUAD patients, while high expression of GPRIN1 was associated with a worse prognosis for LUAD patients. The relationship between these 2 genes and clinical features was also analyzed in this study. GPRIN1 was found to be significantly high in patients with lymph node metastasis ( Figure 5E ), while FCRL5 was low expressed in patients with the clinical stage ( Figure 5F ), lymph node metastasis ( Figure 5G ), and distant metastasis ( Figure 5H ). In addition, the relationship between these 2 genes and clinical diagnosis showed that GPRIN1 had high accuracy in predicting the outcome, while FCRL5 had certain accuracy ( Figures 5I, J ). Previously, Li et al. (20) also reported GPRIN1 as a factor for poor prognosis of non-small cell lung cancer by analyzing RNA-seq data in TCGA and GEO databases. This indicates that GPRIN1 is a vital oncogene in LUAD, and m6A modification plays an important role in the expression level of GPRIN1.

To explore the potential roles of m6A modification in GPRIN1 expression, we conducted a bioinformatics analysis using the GEPIA2 and CPTAC databases. It was found that GPRIN1 was strongly correlated with m6A reader YTHDF1 (r = 0.41) ( Figure 6A ). Also, both GPRIN1 and YTHDF1 proteins were highly expressed in LUAD ( Figures 6B, C ). Further, the m6A RIP-seq data showed that the m6A peaks were enriched in the 3’UTRs of GPRIN1 transcript in LUAD tissues. Compared to the normal tissues, m6A methylation was up-regulated in the vicinity of the 3’UTR of GPRIN1. The region with a gain of m6A in LUAD is the same that gains enrichment of YTHDF1 in A549 cells ( Figure 6D ). It was reported that YTHDF1 regulated gene expression via control of mRNA translation efficiency in an m6A-dependent manner (21). To evaluate whether YTHDF1 involved the expression of GPRIN1, siRNAs targeting YTHDF1 were transfected efficiently into A549 cells. As shown in Figure 6E , knockdown YTHDF1 significantly down-regulated the expression level of GPRIN1. Furthermore, we verified that YTHDF1 could bind GPRIN1 mRNA by DF1-RIP-qPCR assay, in which the fold enrichment of the YTHDF1-IP group was significantly higher than that of the IgG-IP group ( Figure 6F ). Therefore, we hypothesized that YTHDF1 regulates the expression of GPRIN1 by recognizing the m6A modification on GPRIN1 mRNA.