Comprehensive Analysis of the Transcriptome-wide m6A Methylome in Lung Adenocarcinoma by MeRIP Sequencing

N6-methyladenosine (m6A) is the most abundant internal modification on eukaryotic mRNAs. There is increasing evidence that m6A plays a key role in tumor progression, so it is important to analyze m6A modifications within the transcriptome-wide in lung adenocarcinoma (LUAD). Three pairs of LUAD samples and tumor-adjacent normal tissues were obtained from the South University of Science and Technology Hospital. And then methylated RNA immunoprecipitation sequencing (MeRIP-seq) and RNA sequencing (RNA-seq) were used to identify differential m6A modifications between tumor and tumor-adjacent normal tissues. We identified 4041 aberrant m6A peaks, of which 1192 m6A peaks were upregulated and 2849 m6A peaks downregulated. It was found that genes with the dysregulated m6A peaks were enriched in the pathways in cancer, Rap1 signaling pathway, and insulin resistance. Additionally, 612 genes with abnormal regulation of m6A peaks and RNA expression were identified by combining MeRIP-seq and RNA-seq data. Through KEGG analysis, the 612 genes were enriched in cancer-related signaling pathways, such as the cGMP-PKG signaling pathway, and the Rap1 signaling pathway. What’s more, GSEA enrichment analysis showed these genes were enriched in cell cycle phase transition, cell division, cellular response to DNA damage stimulus, and chromosome organization. To further explore the relationship between differential m6A modified genes and clinical parameters of LUAD patients, we searched The Cancer Genome Atlas (TCGA) and identified 2 genes (FCRL5 and GPRIN1) that were associated with the prognosis and diagnosis of LUAD patients. Furthermore, we found a positive correlation between GPRIN1 and m6A reader YTHDF1 in the GEPIA2 database. It was verified that YTHDF1 binds to GPRIN1 mRNA and regulates its expression. Our study results suggest that m6A modification plays important role in the progression and prognosis of LUAD and maybe a potential new therapeutic target for LUAD patients in the future.

Lung cancer is the most common malignant tumor, with the highest mortality rate in the world, and only 21% of the 5-year survival rate (15). LUAD is the most common histological type of lung cancer, accounting for about 40% (16). In recent years, some studies have reported that m6A modification plays an important role in the progression of lung adenocarcinoma (17). Li et al. (16) revealed that YTHDF2 was involved in the development of LUAD via promoting the attenuation of AXIN1 mRNA to activate the Wnt/b-catenin signaling pathway, which may be a potential therapeutic target for LUAD. In addition, miR-600 induces degradation of METTL3 mRNA by targeting the 3'UTR of METTL3 mRNA and leads to inactivation of the PI3K/AKT pathway related to cell growth and survival (18). What's more, METTL3 promotes the translation of YAP mRNA and increases YAP activity by regulating the MALAT1-miR-1914-3p-YAP axis, inducing drug resistance and metastasis of A549 cells (19). Many studies have described the important role of m6A modification in LUAD. However, the transcriptome-wide m6A methylome of LUAD has not been determined.
In this study, we identified differential methylated peaks and expressed genes through MeRIP-seq and RNA-seq. Furthermore, 2 genes with significant changes in both the m6A modification and expression associated with prognosis and diagnosis of LUAD were discovered by additional TCGA clinical data analysis.

Patients and Samples
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.

MeRIP Sequencing and RNA Sequencing
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 highthroughput sequencing platform and PE150 sequencing mode were used.

RIP-qPCR
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. 5mg 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.

Public Databases and Analysis
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.

Overview of the m6A Methylation Map in LUAD
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).

m6A-Containing Genes are Involved in Important Biological Processes and Pathways
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).

Overview of Transcriptome Profiles and Conjoint Analyses of MeRIP-Seq and RNA-Seq Data
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).

Identifying the 2 genes Containing m6A Modification Correlates with Clinical Parameters of LUAD Patients
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 nonsmall 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.

DISCUSSION
m6A is the most abundant internal modification on eukaryotic mRNAs. Recently, the significance of m6A modification is revealed that it involves almost all aspects of mRNA metabolism, including pre-mRNA processing, splicing, nuclear export, translation, and degradation (22). m6A modification also plays critical roles in physiology and pathology such as embryonic development (23), neurodevelopment (24), immunoregulation (25), and tumorigenesis (26). It is also reported that m6A modification plays a crucial role in the development and progression of LUAD (16,27). However, the transcriptome-wide m6A methylome of LUAD has not been characterized. In this study, MeRIP-seq and RNA-seq were performed to reveal the m6A mapping in LUAD. We identified 4041 anomalously regulated m6A peaks in tumor tissues, of which 1192 m6A peaks were upregulated and 2849 m6A peaks were downregulated. Moreover, GO and KEGG pathway analyses revealed the potential functions of differential methylated transcripts. Additionally, by combining MeRIP-seq and RNA-seq data, genes with differential methylated m6A peaks and synchronously differential expression were enriched in cancer-related pathways. Finally, 2 genes were further proven to be associated with the prognosis and diagnosis of LUAD patients in the TCGA database. In this study, it was detected that differential methylated RNAs related to important biological pathways via KEGG pathway analysis. It was found that upregulated m6A modification genes were associated with the ascorbate and alternate metabolism, pentose, and glucuronate interconversions, and sphingolipid signaling pathway. Moreover, downregulated m6A modification genes were significantly associated with the pathways in cancer, Rap1 signaling pathway, and insulin resistance. Combining MeRIP-seq and RNA-seq data, 612 genes with differential methylated m6A peaks and synchronously differential expressed in RNA were identified. We performed GO and KEGG pathway analyses to investigate the biological significance of those genes. KEGG pathway analysis revealed that these genes were primarily enriched in the cGMP-PKG signaling pathway, focal adhesion, and Rap1 signaling pathway. It was revealed that these biological pathways are closely implicated in the onset and development of tumors, including tumor growth, and metastasis (28,29). Rap1 is activated in response to upstream signalings, such as growth factors, cytokines, and chemokines that act on receptor tyrosine kinases and G-protein coupled receptors (30). Rap1 plays an important role in cell-matrix adhesion and supports the key role of the RAP1/RAC1 signaling axis in head and neck squamous cell carcinoma (HNSCC) cell migration via inducing a5b1 integrin via the extracellular matrix molecule fibronectin (31). Rap1 activity is increased in progressively metastatic prostate cancer cell lines and promotes metastasis in in vivo models through the integrins (32). Xiang et al. found that H. pylori infection enhances PRTG expression by promoting the stabilization of transcription factor ZEB1 and recruitment of PRTG promoters, and subsequently activates the cGMP/PKG signaling pathway to further promote proliferation, metastasis, and chemotherapy resistance of gastric cancer cells (33). The effects of EMX2OS and FUS on proliferation, invasion, and migration of human cell lines DU145 and PC3 mediated by the cGMP-PKG pathway were investigated by functional gain and functional loss experiments (34). The attenuation of cellular b-catenin accumulation through the activated cGMP-PKG pathway and the increased phosphorylation of b-catenin were observed in exisulindtreated colonic tumor cells (35). In further investigation, 2 genes (FCRL5 and GPRIN1) were identified as most correlated with the development of LUAD using TCGA databases. Overexpression of GPRIN1 served as a factor for poor prognosis, while FCRL5 favors prognosis. Besides, high expression of GPRIN1 is highly involved in the lymph node metastasis, and FCRL5 with not only lymph node metastasis but also distant metastasis and the clinical stage of LUAD patients. Further analysis revealed that the region with a gain of m6A in LUAD is the same that gains enrichment of YTHDF1 in A549 cells. What's more, YTHDF1 binds GPRIN1 mRNA by YTHDF1-RIP-qPCR assay, in which the fold enrichment of the YTHDF1-IP group was significantly higher than that of the IgG-IP group.
However, our research still has some shortcomings. Firstly, LUAD is a kind of highly heterogeneous cancer, and the tissue samples of clinical patients vary greatly. Although the three pairs   of clinical samples met the statistical requirements, they were too small a sample size to reflect the differential m6A modification and gene expression between LUAD and adjacent tissue in the LUAD population, and the number of cases still must be further increased. Secondly, although we have performed some experiments to verify that YTHDF1 regulates GPRIN1 expression, more experiments are needed to strengthen the credibility, which is also the direction of our future indepth research. In summary, the different m6A methylome in LUAD relative to the corresponding normal controls was analyzed, which demonstrated a strong association between m6A modifcation and LUAD progression. 2 genes (FCRL5 and GPRIN1) with differential methylated m6A peaks and synchronously differential expression were associated with the prognosis and diagnosis of LUAD patients. Further one step research, GPRIN1 was found to act as an downstream target of YTHDF1. Targeting the YTHDF1-m6A-GPRIN1 axes may offer a promising therapeutic approach in curing LUAD.

DATA AVAILABILITY STATEMENT
The data presented in the study are deposited in the GEO repository, accession number GSE198288.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee of the South University of Science and Technology Hospital. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
WM was a major contributor to all of the experimental work, data analysis, and manuscript writing. KW, YZ, and WL were involved in the experimental work. QY, QM, NW, and GZ were involved in the data analysis. YW acquired funding and assisted with the manuscript development. The final manuscript was reviewed and approved by all of the authors.

FUNDING
This work was supported by the Stability Support Plan for Higher Education Institutions in Shenzhen (20200809124551001).