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 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 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⁶A RNA immunoprecipitation with GenSeqTM m⁶A-MeRIP Kit (GenSeq, Beijing, China). Both the input samples without immunoprecipitation and the m⁶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.

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.

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⁶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⁶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.

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).

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⁶A peaks were identified in the CRC tissues, representing 13890 m⁶A-modified transcripts. In NC tissues, 20183 m⁶A peaks were identified, representing 12903 m⁶A-modified transcripts. Compared to NC tissues, CRC tissues had 9456 individual peaks and 3873 m⁶A-modified transcripts indicating a significant difference in global m⁶A modification patterns between CRC tissues and NC normal tissues (Figures 1A,B). The identified m⁶A peak was consistent with the most common m⁶A methylation motif GG/A (m⁶A) CH (Figure 1C) (Wu et al., 2019). 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⁶A-modified peak in CRC tissues, and 50.33% in NC tissues. A relatively small number of m⁶A-modified genes contained two or more peaks (Figure 1D).

With respect to the m⁶A distribution patterns in the entire transcriptome, the results displayed that the m⁶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⁶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⁶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⁶A hypermethylated or hypomethylated with the highest fold-change 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).

To explore the biological significance of m⁶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).

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.

Among these differentially expressed mRNAs, qRT-PCR was used to verify the expression of 4 mRNAs in 8 pairs of CRC and NC tissues. The results of qRT-PCR showed same expression trend with the RNA-Seq results. ESM1 (HGNC:3466), SALL4 (HGNC:15924), TOP2A (HGNC:11989) were upregulated, and MYH11 (HGNC:7569) was downregulated in CRC (Figure 3B).

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).

Using conjoint analysis of MeRIP-seq and RNA-seq data, 400 genes with differential m⁶A methylation and differential expression were identified. Among the 118 m⁶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⁶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 (hypo-down) 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).

To assess the clinical significance of m⁶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).

Many proteins encoded by m⁶A methylated genes played an important role in CRC tissues. The biological functions of RNAs modified by m⁶A methylation were closely related to RBPs. In order to further studied the RBPs related to m⁶A modification in the occurrence and development of colon cancer, 400 different m⁶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⁶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⁶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⁶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 highly-expressed 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 (An et al., 2017).