Comprehensive Analysis of Differentially Expressed Profiles of mRNA N6-Methyladenosine in Colorectal Cancer

Aim: To comprehensively profile the landscape of the mRNA N6-methyladenosine (m6A) modification in human colorectal cancer (CRC). Methods: Methylated RNA immunoprecipitation sequencing (MeRIP-seq) was explored to compare the difference in mRNA N6-methyladenosine (m6A) methylation between CRC tissues and adjacent normal control (NC) tissue. RNA-sequencing (RNA-seq) was performed to transcribe differentially expressed mRNAs. Conjoint analysis of MeRIP-seq and RNA-seq data was conducted to predict RNA-binding proteins (RBPs). Results: MeRIP-seq identified 1110 differentially m6A methylated sites (DMMSs) and 980 differentially m6A methylated genes (DMMGs) in CRC, with 50.13% of all modified genes showing unique m6A-modified peaks in CRC. RNA-seq showed 915 upregulated genes and 1463 downregulated genes in CRC. QRT-PCR verified the RNA-seq results by detecting the expression of some mRNAs. Conjoint analysis of MeRIP-seq and RNA-seq identified 400 differentially m6A methylated and expressed genes (DEGs), and pathway analysis detected that DMMGs and DEGs were closely related to cancer. After analyzing these DMMGs and DEGs through the GEPIA database, we found that the expression of B3GNT6, DKC1, SRPK1, and RIMKLB were associated with prognosis, and the expression of B3GNT6 and RIMKLB were associated with clinical stage. 17 RBPs were identified based on the DMMGs and DEGs, among which FXR1, FXR2, FMR1, IGF2BP2, IGF2BP3, and SRSF1 were obviously highly expressed in CRC, and FMR1, IGF2BP2, and IGF2BP3 were closely related to methylation, and might be involved in the development of CRC. Conclusion: This study comprehensively profiled m6A modification of mRNAs in CRC, which revealed possible mechanisms of m6A-mediated gene expression regulation.

RNA m 6 A modification is not only involved in normal metabolism in the body, but also associated with the development of different types of tumors, such as colorectal tumors (CRC), acute myeloid leukemia, non-small cell lung cancer, hepatocellular carcinoma, cervical squamous cell carcinoma, and breast cancer (He et al., 2019). A number of studies display that RNA m 6 A methylation of CRC is involved in the activity regulations of cancer Stem Cell, as well as the growth, proliferation, anti-chemotherapy and immunotherapy of tumour cells (Bai et al., 2019;Wang et al., 2020;Yang et al., 2021). Therefore, to further reveal the pathogenesis of CRC, it is worthwhile to study the modification patterns of CRC methylation.
Methylated RNA immunoprecipitation sequencing (MeRIPseq) technology have become the main detection method for m 6 A modification (Dominissini et al., 2012). In this study, we conducted MeRIP-seq and RNA-sequencing (RNA-seq) in CRC tissues and the corresponding normal control (NC) tissues to investigate the m 6 A modification patterns in colorectal cancer, the involvement of m 6 A modification in the development of CRC, the relationship between mRNA m 6 A modification and mRNA expression, together with the predicted RBPs closely related to mRNA m 6 A methylation for a comprehensive understanding of m 6 A methylation in CRC.

Tissue Samples
All five specimens of patients diagnosed with CRC were obtained at the Chengdu Medical College from March 2017 to June 2018. None of these patients received radiation or chemotherapy before the specimens were collected. Clinicopathological informations of 5 CRC patients were displayed in Supplementary Table S1. Three to four specimens with diameters of 0.5-1.0 cm were collected from the center of the CRC tissues. At the same time, three to four pieces of the NC tissues of a similar size were collected 5 cm away from the edge of the cancer tissue. Informed consent was obtained from patients regarding the collection and use of these specimens. The collected specimens were frozen sectioned and stained with hemAtoxylin-eosin to confirm that the collected specimens were colon cancer (COAD) or rectal cancer (READ) tissues and corresponding NC tissues. This study had been approved by the Ethics Committee of the First Affiliated Hospital of Chengdu Medical College.

RNA Isolation
RNA was extracted using TRIzol Reagent (Invitrogen, Carlsbad, CA, United States) according to the manufacturer's instructions. The NEBNext rRNA Depletion Kit (New England Biolabs, Inc., Ipswich, MA, United States) was used to remove ribosomal RNA from the total RNA. The concentration and purity of RNA were detected using a NanoDrop ND-1000 spectrophotometer (Thermo, Waltham, MA, United States), and the integrity of RNA was evaluated using denaturing gel electrophoresis.

MeRIP-Seq and RNA-Seq
MeRIP-seq and RNA-seq were provided by Cloudseq Biotech Inc. (Shanghai, China). RNA was fragmentated using RNA fragmentation reagents (Thermo, Waltham, MA, United States). A total of 5 μg of fragmented mRNA was saved as input control for RNA-Seq, while 500 μg of fragmented mRNA was used to perform m 6 A RNA immunoprecipitation with GenSeqTM m 6 A-MeRIP Kit (GenSeq, Beijing, China). Both the input samples without immunoprecipitation and the m 6 A IP samples were used to generate libraries for RNA sequencing using the NEBNext ® Ultra II Directional RNA Library Prep Kit (New England Biolabs, Inc., Ipswich, MA, United States). Library quality was evaluated using a BioAnalyzer 2100 system (Agilent Technologies, Inc., Palo Alto, CA, United States). Libraries were sequenced on an Illumina HiSeq instrument 4000 with 150 bp paired-end reads. The main reagents and instruments of the research were shown in Supplementary Tables S2, S3 respectively.

Real-Time Quantitative PCR
The cDNA primer sequences used in the study designed and synthesized by QingKe Biotechnology Co., Ltd. (Chengdu, China). The primer sequences of the target genes, β-actin and GAPDH used for RT-qPCR were displayed in Supplementary  Table S4. Total RNA Extraction KIT (Solarbio, Beijing, China) was used to extract total RNA from tissues according to the manufacturer's instructions. The concentration and purity of the isolated RNA were determined by NanoDrop ND-1000 spectrophotometer and the iScript cDNA Synthesis Kit (Bio-Rad, Hercules, CA, United States) was used for reversal. RT-qPCR was performed using SsoAdvanced Universal SYBR Green Supermix (Bio-Rad, Hercules, CA, United States). All qPCR analyses were conducted in triplicate and the average value was calculated. The expression values of genes were normalized β-actin and GAPDH, and then log2-ΔΔCt for the gene expression.

Data Analysis
Reads were harvested from Illumina HiSeq sequencer, and were quality controlled by Q30. After 3' adaptor-trimming and low quality reads removing by cutadapt software (v1.9.3). First, clean reads of input libraries were aligned to reference genome (UCSC HG19) by Hisat2 software (v2.0.4). After that, The MACS (v1.4.2) software compared the clean reads with the human genome to identify the methylation sites (peaks) on RNAs. Differentially m 6 A methylated sites (DMMSs, fold-change ≥ 2 and p-value < 0.05) were identified by diffReps (v1.55.3). Sequence motifs were identified using Homer (v3.0). Under the guidance of the Gene Transfer Format gene annotation file, gene expression was calculated by Cufflinks (v2.2.1) using FPKM values (Fragments per kilobase of exon per million fragments mapped). Fold change ≥2, p-value <0.05 and FPKM ≥0.1 were acted as the threshold for differentially m 6 A methylated genes (DMMGs) screening. These peaks identified by both software with overlapping mRNA exons were identified by home-made scripts. Genes of interest were visualized in the Integrative Genomics Viewer (IGV) software (v2.8.3). Screening of differentially expressed genes (DEGs) in CRC group and NC group using Cuffdiff (v2.2.1) software.
Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed by DAVID (v6.7) database. p < 0.05 was used as the threshold of significant enrichment.
The combined analysis of MeRIP-seq and RNA-seq screened out the genes that were differentially methylated and expressed. The differential methylation regions of CRC vs. NC were divided into 9 portions from −5 to 4 according to the methylation fold change. The binding region of RBP intersected with the differentially methylated region of the selected gene to find the RBP that can bind to the methylated region via RMbase (v2.0) database. [columns (Log2FC): the ratio of differentially methylated areas of CRC/NC; rows: RBPs]. The color in the figure represented the RBPs binding rate. The RBPs binding rate was calculated through the number of RBPs binding divided by the number of differential methylation regions within the interval multiplied by the number of binding information of the RBP in the database. The darker the color, the greater the binding rate, which meant that the differentially methylated regions in the range were more likely to be bound by the RBPs. The main software information used in the research was shown in Supplementary Table S5.

Statistical Analysis
Data from three or more independent experiments were presented as mean ± standard deviation. Paired student's t-tests were performed between cancer tissue samples and NC tissue samples. One-way analysis of variance was used to access the differences among five or more groups. Differences with p < 0.05 were considered statistically significant (*p < 0.05, **p < 0.01, ***p < 0.001).

General Features of mRNA m 6 A Modification Patterns in CRC and NC Tissues
The MeRIP-seq technology was used to perform methylation sequencing on samples from CRC and NC. Reads data and quality testing were shown in Supplementary Table S6. A total of 23181 m 6 A peaks were identified in the CRC tissues, representing 13890 m 6 A-modified transcripts. In NC tissues, 20183 m 6 A peaks were identified, representing 12903 m 6 A-modified transcripts. Compared to NC tissues, CRC tissues had 9456 individual peaks and 3873 m 6 A-modified transcripts indicating a significant difference in global m 6 A modification patterns between CRC tissues and NC normal tissues ( Figures 1A,B). The identified m 6 A peak was consistent with the most common m 6 A methylation motif GG/A (m 6 A) CH ( Figure 1C) . It also confirmed the reliability of the obtained MeRIP-seq results. In addition, it was found that 50.13% of all modified genes had the unique m 6 A-modified peak in CRC tissues, and 50.33% in NC tissues. A relatively small number of m 6 A-modified genes contained two or more peaks ( Figure 1D).

Distribution of Differentially Methylated m 6 A Peaks
With respect to the m 6 A distribution patterns in the entire transcriptome, the results displayed that the m 6 A peaks were mainly enriched in the vicinity of the CDs and at the beginning of the 3'untranslated region (3′UTRs) in both groups ( Figures  1E,F). IGV was used to visually analyze the differential m 6 A peak of AQP5 (HGNC:638), KISS1 (HGNC:6341) and MUC16 (HGNC:15582) in CRC and NC. The differential hypermethylated peak of AQP5 on outside of the 5′UTRs was detected. The differential hypermethylated peaks of KISS1 and MUC16 were both enriched in the CDS regions ( Figure 1G, Supplementary Figure S1).
A total of 1110 DMMSs and 980 differentially m 6 A methylated genes (DMMGs) were identified in CRC tissues. Out of the 980 DMMGs, 43.78% (429/980) were significantly hypermethylated, and 56.22% (551/980) were significantly hypomethylated (Supplementary Figure S2A). The top 10 genes of m 6 A hypermethylated or hypomethylated with the highest foldchange values were shown in Table 1. There was no significant difference in the distribution pattern of hypermethylated sites. The distribution of hypomethylation sites in the 5′UTR was significantly less than the 3′ UTR and CDS regions (Supplementary Figure S2B, Supplementary Table  S7). All DMMSs were mapped to chromosomes to analyze their distribution profiles. DMMSs were mostly on chromosomes 1, 2, 11, and 19 (Supplementary Figure S2C).   Frontiers in Cell and Developmental Biology | www.frontiersin.org January 2022 | Volume 9 | Article 760912 5

Differentially Methylated RNAs Were Involved in Important Biological Pathways
To explore the biological significance of m 6 A modification in CRC tissues, GO analysis and KEGG pathway analysis were performed for DMMGs. Go analysis revealed that the hypermethylated genes were significantly enriched in the cell cycle process, intracellular organelle lumen, and nucleic acid binding (Figure 2A). The hypomethylated genes were significantly involved in the regulation of ion homeostasis, intracellular, and ion channel binding ( Figure 2B). Cell cycle process, intracellular, and nucleic acid binding played an important role in tumor development (Dominissini et al., 2012).
KEGG pathway analysis revealed that hypermethylated genes were enriched in signaling pathways regulating the pluripotency of stem cells and the tumor protein p53 (TP53, HGNC:11998) signaling pathway ( Figure 2C). Hypomethylated genes were enriched in the calcium signaling pathway and the nuclear factor kappa B (NF-κB, HGNC: 9954) signaling pathway ( Figure 2D). Stem cell pluripotency regulation, the TP53 signaling pathway, the NF-κB signaling pathway, and calcium ion regulation were closely related to the formation, metastasis, and drug resistance of multiple tumors (Ha et al., 2019;Liu et al., 2020b). Functional Analysis of Differentially Expressed Genes in CRC RNA-seq was used to detect DEGs (fold-change ≥ 2 and p-value < 0.05) in CRC and NC tissues, which showed 915 upregulated genes and 1463 downregulated genes ( Figure 3A). Reads data and quality testing were shown in Supplementary Table S8. The top 10 upregulated and downregulated genes were listed in Table 2.
To explore the physiological and pathological significance of DEGs in CRC, GO analysis and KEGG pathway analysis were performed. GO analysis revealed that the upregulated genes were significantly enriched in chromosome organization, extracellular organelles, and structural constituent of ribosomes, while the downregulated genes were significantly enriched in the developmental process, cell periphery, and calcium ion binding ( Figure 3C). KEGG pathway analysis revealed that the upregulated genes were enriched in DNA replication and the TP53 signaling pathway, while downregulated genes were involved in the calcium signaling pathway and the MAPK signaling pathway ( Figure 3D). The TP53 signaling pathway, calcium signaling pathway, and MAPK signaling pathway were important for the malignant progression of tumors (Nucera et al., 2011).

Differentially Methylated and Expressed mRNAs Were Related to Cancer
Using conjoint analysis of MeRIP-seq and RNA-seq data, 400 genes with differential m 6 A methylation and differential expression were identified. Among the 118 m 6 A hypermethylated genes, 117 genes had upregulated mRNA expression, and one gene had downregulated mRNA expression ( Figure 4A). The top 10 hypermethylated and upregulated genes with the highest methylation fold changes were listed in Table 3. Additionally, there were 281 downregulated genes and one upregulated gene among the 282 m 6 A hypomethylated genes ( Figure 4A). The top 10 hypomethylated and downregulated genes with the highest methylation fold changes were listed in Table 4. Almost all genes showed the same trend of methylation and expression. Hyper-methylation was accompanied by upregulation of mRNA expression, while hypo-methylation was accompanied by downregulation of mRNA expression. The results showed that the methylation of mRNA in CRC tissues may maintains the stability of mRNA and promotes post-transcriptional translation of mRNA.
To explore the role of differentially methylated and expressed mRNAs in the biological process of CRC, GO analysis and KEGG pathway were performed on hypermethylated-upregulated (hyper-up) genes and hypomethylated-downregulated (hypodown) genes.
GO analysis indicated that hyper-up genes were mainly involved in the processes of microtubule motor activity, spindle and nuclear division ( Figure 4B). The hypo-down genes were mainly involved in actin binding, cation channel complex, cytosolic calcium ion transport ( Figure 4C). KEGG pathway analysis showed that hyper-up genes were enriched in the cell cycle and protein processing endoplasmic reticulum pathways ( Figure 4D). Hypo-down genes were mainly enriched in calcium signaling pathway, vascular smooth muscle constraction and purine metabolism pathways ( Figure 4E).

Correlation Between m 6 A-Regulated Gene Expression and Clinical Parameters
To assess the clinical significance of m 6 A-regulated gene, UALCAN and GEPIA databases were explored. These databases were used to investigate the relationship between the   expression of 400 differentially methylated-expressed genes with tumor staging, lymphatic metastasis, and survival period. B3GNT6 (HGNC:24141), RIMKLB (HGNC:29228) in COAD, and DKC1 (HGNC:2890), SRPK1 (HGNC:11305) in READ were analyzed. Survival analysis showed that patients with high B3GNT6, DKC1, SRPK1, and low RIMKLB expression had a longer overall survival period ( Figure 5). Compared with NC, the expressions of B3GNT6, RIMKLB, DKC1, and SRPK1 in COAD or READ were significantly different, which consistented with our sequencing datas. There was a significant difference in the expression of B3GNT6 in COAD stage 2 and stage 3. The expression of RIMKLB was significantly lower in COAD stage 1 than stage 2, 3, and 4. Meawhile, the expression of RIMKLB in COAD N1 stage was different from that of N0 and N2, respectively ( Figure 6).

Conjoint Analysis of RBPs Based on Differentially Methylated and Differentially Expressed mRNA
Many proteins encoded by m 6 A methylated genes played an important role in CRC tissues. The biological functions of RNAs modified by m 6 A methylation were closely related to RBPs. In order to further studied the RBPs related to m 6 A modification in the occurrence and development of colon cancer, 400 different m 6 A methylation modified and expressed genes were screened out to predict the RBPs through the RMBase database. According to the ratio of potential RBPs bound to m 6 A sites in all differentially methylated and differentially expressed mRNAs. The overview of the RBPs was presented as a heatmap ( Figure 7A). In order to further verify the expression level of RBPs in colorectal cancer, we searched the protein expression in 100 NC tissues and 97 COAD tissues through the CPTAC database. Consistent with our conclusion, FXR1 (HGNC:4023), FXR2 (HGNC:4024), FMR1 (HGNC:3775), IGF2BP2 (HGNC: 28867), IGF2BP3 (HGNC:28868), and SRSF1 (HGNC:10780) were obviously highly expressed in COAD (Figures 7B-G).
Among the 17 RBPs, FXR1/2, FMR1, PUM2 (HGNC:14958), LIN28A (HGNC:15986), IGF2BP2/3, HNRNPH1 (HGNC:5041), and SRSF1 were closely related to mRNA m 6 A methylation in CRC. FXR1/2, FMR1, PUM2, LIN28A, and IGF2BP2/3 were involved in the process of hyper-methylation. HNRNPH1 and SRSF1 were involved in the process of hypo-methylation. Specifically, FXR1/2, FMR1, IGF2BP2/3, and SRSF1 were closely related to the occurrence and malignant progress of CRC. The corresponding mRNAs of FMR1, IGF2BP2, and SRSF1 were presented in Table 5, and a diagram showing the pattern of m 6 A methylation regulation in CRC was shown in Figure 8. FXR1/2, FMR1, IGF2BP2/3, and SRSF1 were RBPs that shuttled between the nucleus and cytoplasm, which could increase the stability of targeted mRNA translation. In this study, we detected many hypermethylated and highlyexpressed mRNAs related to these RBPs, which were closely related to the occurrence and malignant progression of CRC. For example, HOXD10 (HGNC:5133) was a major factor that negatively regulates tumor metastasis, and the expression of the HODX10 gene in the CRC tissue with lymphatic metastasis was lower than that in tissues without lymphatic metastasis (Wang et al., 2016). ORC6 (HGNC:17151) was related to the survival outcome of CRC (Hu et al., 2019). CDC20 (HGNC:1723), a cell cycle regulator, can bind to and activate the APC (HGNC:583) complex and participate in the development of CRC (Yu, 2007). Overexpression of the BLM (HGNC:1058) gene increased the risk of CRC (Laitman et al., 2016). Patients with a high expression of ASPM (HGNC:19048) were at more advanced clinical stages and were more prone to lymph node metastasis .

DISCUSSION
N 6 -methyladenosine methylation is the most abundant modification in eukaryotic mRNAs. Compared to the normal physiological process, there were cases of both hyper-and hypomethylations in m 6 A methylation in tumors, and they both affected the development of tumors by regulating oncogene expression, cancer cell growth, survival, and invasion. METTL14 inhibited hematopoietic stem cell/progenitor cell differentiation and promotes leukemia by enhancing the m 6 A modification of MYB (HGNC:7545)/MYC (HGNC:7553) . YTHDF1 augmented the translation of EIF3C (HGNC:3279) in an m 6 A-dependent manner by binding to m 6 A-modified EIF3C mRNA and concomitantly promoted the overall translational output, thereby facilitating tumorigenesis and metastasis of ovarian cancer (Liu et al., 2020a). M 6 A demethylase ALKBH5 affected the proliferation and invasion of lung adenocarcinoma cells by downregulating m 6 A modification on FOXM1 (HGNC:3818) mRNA and promoting FOXM1 expression (Chao et al., 2020). Therefore, the study of methylation contributed to further understand the pathogenesis of tumors and provide directions for treatment.
In this study, MeRIP-seq and RNA-seq analyses were conducted in CRC tissues and corresponding NC tissues to profile m 6 A modification patterns in CRC. MeRIP-seq analysis revealed a total of 23181 m 6 A peaks in the CRC group, representing 13890 m 6 A-modified transcripts. Compared to the NC group, the CRC group had 9456 individual peaks and 3873 m 6 A-modified transcripts. A total of 1110 DMMSs and 980 DMMGs were identified in the CRC group, and 50% of the modified genes had the unique m 6 A-modified peak in CRC. DMMSs were more abundant on chromosomes 1, 2, 11, and 19, and predominantly distributed in the CDs and 3′ UTR.
To understand the relationship between m 6 A methylation and mRNA expression, conjoint analysis of MeRIP-seq and RNA-seq datas was conducted, and 400 genes with differential m 6 A methylation and differential expression were identified. There was a positive correlation between methylation and expression that hyper-methylation was found to be accompanied by Frontiers in Cell and Developmental Biology | www.frontiersin.org January 2022 | Volume 9 | Article 760912 10 upregulation of mRNA expression, and hypo-methylation was accompanied by downregulation of mRNA expression. Notably, m 6 A methylation can either reduce or enhance the stability, cleavage, and transport of target mRNAs (Taketo et al., 2017;Lin et al., 2019).
GO and KEGG analysis showed that the differentially methylated and expressed genes mainly enriched in Spindle and Nuclear division, calcium signaling pathway, cell cycle and protein processing endoplasmic reticulum pathways etc. Parts of pathways were closely related to cancer. Previous studies showed that ASPM was highly expressed in hepatocellular carcinoma (HCC), and the high expression of ASPM was correlated with poor HCC prognosis. Mechanism studies showed that METTL3 promoted HCC cells growth and metastasis via m 6 A modification of ASPM (Wang et al., 2021). Mettl14-mediated m 6 A played a critical role in liver regeneration by regulating the expression of polypeptide-processing proteins and maintaining endoplasmic reticulum homeostasis (Cao et al., 2021). The prediction of these pathways provided clues for further research on the methylation mechanism of colorectal cancer.
The regulation of m 6 A methylation was mediated by key enzymes and RBPs. RBPs prediction showed that 17 RBPs may be closely related to the differentially methylated and expressed genes in CRC. Among the 17 RBPs, FXR1, FXR2,  Frontiers in Cell and Developmental Biology | www.frontiersin.org January 2022 | Volume 9 | Article 760912 FMR1, PUM2, IGF2BP2, IGF2BP3, and SRSF1 were closely related to mRNA m 6 A methylation in CRC. FXR1 was a new driver in the 3q26-29 amplicon and predicted poor prognosis in human cancers (Qian et al., 2015). It shuttled between the nucleus and cytoplasm, highly expressed in lung and ovarian cancer (Le Tonquese et al., 2016), and may act as a tumor promoter to increase the proliferation, migration, and invasion of colorectal cancer (Jin et al., 2016). FMR1 facilitated the nuclear export of N 6 -methyladenosine-containing mRNAs (Hsu et al., 2019). The expression of PUM2 was significantly increased in glioblastoma, and the downregulation of PUM2 can significantly inhibit the proliferation and migration of glioblastoma cells . METTL3 stabilized the expressions of HK2 (HGNC:4923) and SLC2A1 (HGNC:11005) in CRC through a m 6 A-IGF2BP2/3dependent mechanism to regulate glycolytic pathways and promote colorectal cancer progression (Shen et al., 2020).
Research demonstrated that IGF2BP3 bound to the mRNA of CCND1 (HGNC:1582), and reduced its mRNA stability via reading m6A modification in the CDs region. Overexpression of CCND1 in IGF2BP3 down-regulated cells completely rescued the inhibited percentage of S phase in cell cycle as well as cell proliferation. . SRSF1 was identified as a proto-oncogene that shuttled between the cytoplasm and nucleus and bound to the GGAG base sequence on mRNAs to regulate the processes of mRNA transcription, stabilization, nucleation, and translation (Long and Caceres, 2009;Das and Krainer, 2014). Long non-coding RNA AGAP2-AS1 (HGNC:48633) accelerated cell proliferation, migration, invasion and the EMT process in CRC via regulating the miR-4668-3p/SRSF1 axis . The involvements of FMR1, IGF2BP2, IGF2BP3 in methylation have been preliminarily explored. However, the researchs on the mechanism of more RBPs involvement in CRC methylation deserve to be investigated further.

CONCLUSION
This study established the mRNA m 6 A methylation modification map of CRC. Genes of differential m 6 A methylation and differential expression were identified and the RBPs of these genes were predicted. The possible molecular mechanism of m 6 A methylation in the development of CRC was revealed. This study provided a foundation for in-depth mechanistic analysis of m 6 A methylation in CRC.

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

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the First Affiliated Hospital of Chengdu Medical College. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.

AUTHOR CONTRIBUTIONS
YZ and X-AL were involved in the conception and design of the study. QG and NL acquired the data and drafted the manuscript. QZ conducted data analysis. B-JC collected the samples. YZ critically revised the manuscript for intellectual content, approved the final version. All authors read and approved the final manuscript.

FUNDING
The present study was supported by the Seedling Engineering Cultivation Foundation of Science and Technology (grant No. 2018075).