Potential Impact of ALKBH5 and YTHDF1 on Tumor Immunity in Colon Adenocarcinoma

Background ALKBH5 and YTHDF1 are regarded as the eraser and reader, respectively, in N6-methyladenosine (m6A) modification. Recently, immune contexture has been drawing increasing attention in terms of the progression and treatment of cancers. This study aimed to determine the relationship between ALKBH5/YTHDF1 and immunological characteristics of colon adenocarcinoma (COAD). Methods Expression of ALKBH5 and YTHDF1 was investigated across TCGA and GEO validated in our study. Patients with COAD were divided into two clusters using consensus clustering based on the expression of ALKBH5 and YTHDF1. We then compared their clinical characteristics and performed gene set enrichment analysis (GSEA) to identify the functional differences. Immune infiltration analyses were conducted using ESTIMATE, CIBERSORT, and ssGSEA. In addition, we evaluated the expression of the targets of immune checkpoint inhibitors (ICIs) and calculated the tumor mutation burden (TMB) of the tumor samples. Weighted gene co-expression network analysis (WGCNA) was used to identify the genes related to both ALKBH5/YTHDF1 expression and immunity. GSE39582 was utilized for external validation of immunological features between the two clusters. Results Cluster 2 had high expression of ALKBH5 and lesser so of YTHDF1, whereas Cluster 1 had just the reverse. Cluster 1 had a higher N stage and pathological stage than Cluster 2. The latter had stronger immune infiltration, higher expression of targets of ICIs, more TMB, and a larger proportion of deficiency in mismatch repair-microsatellite instability-high (dMMR-MSI-H) status than Cluster 1. Moreover, WGCNA revealed 14 genes, including PD1 and LAG3, related to both the expression of ALKBH5/YTHDF1 and immune scores. Conclusions ALKBH5 and YTHDF1 influence immune contexture and can potentially transform cold tumors into hot tumors in patients with COAD.


INTRODUCTION
Colorectal cancer (CRC) ranks third in incidence, and its mortality ranks second in both sexes worldwide (1). Approximately 70-80% of patients with early stage CRC are eligible for surgery, and the five-year survival rate of these patients is approximately 90%. However, the five-year survival rate of patients with distant CRC (stage IV) is as low as 10-15% (2,3).
Colon adenocarcinoma (COAD) is the main type of CRC that originates from adenomatous lesions and evolves into cancer due to the accumulation of genetic mutations (4). Recently, a study has shown that mutations can generate new antigens, which can be recognized by the immune system (5). Moreover, immunotherapy has proven to be effective in treating advanced carcinomas (6). Nevertheless, immunotherapy has limitations in some patients with microsatellite instability (MSI) and most patients with microsatellite stability (MSS) (7). Therefore, identification of novel immunotherapy markers and uncovering the underlying mechanisms of immune checkpoints would be important.
N6-methyladenosine (m6A) is the most abundant RNA modification that occurs in both coding and non-coding RNAs, and is a crucial post-transcriptional regulator in various cancers (8)(9)(10)(11). Proteins involved in m6A modification may be divided into three categories: writer (catalyzes the occurrence of m6A modification), eraser (catalyzes the removal of m6A modification), and reader (recognizes and binds m6A modification) (12). Although tumor-intrinsic carcinogenic processes are vital, the impact of m6A modification on tumor and immunity is also worth attention. In recent years, some studies have suggested that targeting of dysregulated m6A regulators with small molecule inhibitors has potential in treating cancer. Given the functional importance of m6A modification in various cancers, targeted treatment against m6A regulators may be applicable in the clinic, in combination with chemotherapy or immunotherapy, to improve cancer therapy (13).
ALKBH5 plays the role of eraser in m6A modification and has been proven to regulate suppressive immune cell accumulation in melanoma (14,15). In pancreatic ductal adenocarcinoma, ALKBH5 inhibits tumorigenesis by decreasing m6A modification of WIF-1 RNA and mediating the Wnt signaling pathway (16). On the other hand, YTHDF1 is a reader in m6A modification that can improve the efficiency of mRNA translation (17). It regulates the expression of lysosomal proteases in an m6A-dependent manner to control anti-tumor immunity and improve the efficacy of immunotherapy. Deficiency of YTHDF1 can enhance the therapeutic effect of PD-L1 checkpoint blockade (18).
To comprehensively understand the role of ALKBH5 and YTHDF1 in tumor immunity, we analyzed the transcriptome profiling data of colon adenocarcinoma from The Cancer Genome Atlas (TCGA) and Gene Expression Omnibus (GEO). We divided patients into two clusters based on the expression of ALKBH5 and YTHDF1 and compared their differences in clinical characteristics, biological pathways, immune infiltration, immune checkpoint expressions, and mutational landscapes using bioinformatics methods. Our results indicated that Cluster 2, with high expression of ALKBH5 and low expression of YTHDF1, had more immune infiltration, immune checkpoint inhibitor expression, and tumor mutation burden than Cluster 1, hence suggesting that Cluster 2 might respond better to immunotherapy. ALKBH5 and YTHDF1 may remarkably influence the immune contexture of colon adenocarcinoma.

Data Sources and Preprocessing
Transcriptome profiling data (HTSeq-Counts and HTSeq-FPKM) with clinical information were downloaded from TCGA-COAD project by R (version 4.0.2) with R package TCGAbiolinks (19). Cases that contained intact clinical information (age, sex, T stage, N stage, M stage, and prognostic information) were included. Level 3 HTSeq-FPKM of 435 primary solid tumor samples were treated by log2(FPKM+1) transformation for further analyses, and HTSeq-Counts were used for differential analysis.
Simple nucleotide variation data (MuTect2) of 376 patients with COAD were collected using R package maftools (20). Due to the lack of mutation information for some patients with COAD, we only included 376 patients in the analysis involving mutational landscape. Waterfall plots were used to show the genetic mutation of patients using the R package ComplexHeatmap (21). Tumor mutation burden (TMB) was calculated based on simple nucleotide variation, defined as the number of mutations per megabase.
Expression profiling by GSE39582 array was downloaded from the Gene Expression Omnibus (GEO) database (https:// www.ncbi.nlm.nih.gov/gds/) (22). The dataset with 566 colon cancer tissues was used to verify the immune characteristics of patients with colon cancer.

Immune Infiltration Analysis
ESTIMATE is a method that determines the fractions of stromal and immune cells based on gene expression signatures in tumor samples. It was applied to evaluate the tumor microenvironment (TME) of each patient with COAD, along with stromal score (stromal content), immune score (extent of immune cell infiltration), ESTIMATE score (synthetic mark of stroma and immune), and tumor purity by R package estimate (23).
CIBERSORT is a means of computing cell composition based on the expression profiles. This deconvolution algorithm was used to calculate the proportion of 22 immune cells in each patient with COAD (24). The sum of 22 immune cell type fraction in each sample was 1.
By applying the single-sample gene set enrichment analysis (ssGSEA) method from R package GSVA (25), we calculated the extent of infiltration of 28 immune cell types according to the expression levels of genes in 28 published gene sets for immune cells (26).

Consensus Clustering Based on ALKBH5 and YTHDF1
Expression of ALKBH5 and YTHDF1 was extracted and clustered coherently using the R package ConsensusClusterPlus (27). The samples were divided into two clusters. We used the R package CMScaller to identify the consensus molecular subtypes (CMS) of each sample (28). CMS is a robust classification system for CRC. Every CMS has distinct features: CMS1 (immune), CMS2 (canonical), CMS3 (metabolic), and CMS4 (mesenchymal) (29). A Sankey diagram was used to indicate the relationship between the two clusters and CMS.

Gene Set Enrichment Analysis
GSEA was performed using the R package clusterProfiler to discover the significant functional difference between the two clusters (30). Significant pathway enrichment was identified by the normalized enrichment score (|NES| >1), P value <0.05, and FDR q value <0.05.

Differential Expressed Genes
Expression profiling data (HTSeq-Counts) were compared to identify the DEGs of two clusters using the R package DESeq2 (31). The threshold values were |log2FoldChange | >1 and adjusted P value <0.05.

Weighted Gene Co-Expression Network Analysis
We performed WGCNA on DEGs using the R package WGCNA (32). To ensure that the constructed co-expression network approached scale-free distribution, we chose 5 as the soft power. We obtained nine modules and calculated their correlation with cluster, stromal score, immune score, ESTIMATE score, and tumor purity. Subsequently, we acquired 14 genes according to the calculation of module membership (MM) and gene significance (GS).

Functional Enrichment Analysis
Gene Ontology (GO) analysis was applied to understand the functions of 14 selected DEGs using the R package clusterProfiler (30). We then constructed a protein-protein interaction (PPI) network using the STRING database (33). Next, we analyzed the Spearman's correlation of gene-gene, gene-ESTIMATE, and gene-ssGSEA using the R package corrplot.

Specimen Collection and Real-Time Quantitative PCR
Twelve pairs of CRC tissues and their adjacent tissues were collected from the First Hospital of China Medical University with informed consent and approval from the Institutional Ethics Board of the First Hospital of China Medical University.
Total RNA was extracted using TRIzol reagent (Invitrogen). cDNA was synthesized using the PrimeScript RT Reagent Kit (TaKaRa). The SYBR Prime Script RT-PCR kit (TaKaRa) was used to perform RT-qPCR according to the manufacturer's protocol.

Statistical Analysis
All statistical analyses were conducted by R (4.0.2) and SPSS (25.0) software. Figure panels were pieced together by Adobe Illustrator (CC 2019). Box plot analyses were performed using the Wilcoxon rank-sum test. Correlation analysis was performed using the Spearman's coefficient. Chi-square test was used to compare the clinical characteristics between the two clusters (Fisher's exact test was used when required). Multivariate logistic regression analysis was used to evaluate the clinical characteristics affecting the clusters. Survival curves were constructed using the Kaplan-Meier method (log-rank test). All hypothetical tests were two-sided, and a P value <0.05 was considered significant.

Identification of Immune-Related m6A Regulators and Consensus Clustering of Patients
First, we calculated four indices of ESTIMATE in each sample to assess the fractions of stromal and immune cells. In order to explore the role of m6A modification in tumor immunity of patients with COAD, 21 m6A regulators were identified, and correlation between the expression of m6A regulators and results of ESTIMATE was evaluated ( Figure 1A). Considering the highest absolute value of correlation with immune score, ALKBH5 and YTHDF1 were included in subsequent analyses. Next, Wilcoxon rank-sum test was conducted between tumor and normal tissues using RNA-seq data of TCGA-COAD; the tumor tissues were found to have lower ALKBH5 expression and higher YTHDF1 expression than normal tissues ( Figures 1B, C). We thereafter performed a consensus clustering on 435 TCGA-COAD samples according to the expression matrix of ALKBH5 and YTHDF1, and divided the samples into two clusters ( Figure 1D and Supplementary Figure 1). The heatmap shows Cluster 1 (n = 217) to have low expression of ALKBH5 and high expression of YTHDF1 while Cluster 2 (n = 218) had low expression of YTHDF1 and high expression of ALKBH5. Since the two genes showed opposite trends in the two clusters, Spearman's correlation between ALKBH5 and YTHDF1 was investigated, and a weak, negative correlation (R = −0.30, P = 1.34e-10) was found ( Figure 1E). To understand the features of the two clusters better, we evaluated the CMS of each sample and drew a Sankey diagram to indicate their relationship (Supplementary Figure 2).

Evaluation of Clinical Characteristics
To identify the differences in clinical characteristics between the two clusters, we drew a survival curve first, and found no significant difference in prognosis between the two clusters (Supplementary Figure 3). Associated analysis showed Cluster 1 to have higher N stage, higher pathological stage, and lesser age than Cluster 2 ( Table 1). Moreover, multivariate logistic regression analysis showed age, T stage, and N stage to be independent factors affecting clustering ( Table 2).

Identification of Immune-Related Pathways by GSEA
GSEA was performed to understand the functional differences between the two clusters. All differentially expressed genes (Cluster 2 vs. Cluster 1) were included in the GSEA. We identified many significant pathways related to immunity in the enrichment of MSigDB Collection (c5.cp.v7.0.symbols.gmt), including adaptive immune response, cell killing, humoral immune response, positive regulation of cytokine production, T cell activation, and T cell proliferation (Figure 2A).

Comparison of Immune Infiltration
ESTIMATE, CIBERSORT, and ssGSEA were performed to understand the differences in immunological function better. Cluster 2 had higher stromal, immune, and ESTIMATE scores, and lower tumor purity than Cluster 1 in the ESTIMATE analysis ( Figures 2B-E). Furthermore, CIBERSORT analysis demonstrated that Cluster 2 to have a higher proportion of CD8 T cells ( Figure 2F). ssGSEA showed 25 immune cell subtypes (such as activated B cells, activated CD4 T cells, activated CD8 T cells, activated dendritic cells, natural killer  cells, and natural killer T cells) to be highly expressed in Cluster 2 ( Figure 2G). A heatmap was prepared to show the overall conditions of the 28 immune cell subtypes in the two clusters (Supplementary Figure 4). Results indicated that Cluster 2 tended to have a stronger immune infiltration than Cluster 1, especially regarding CD8 T cells.

Evaluation of Sensitivity to Immunotherapy
To evaluate the sensitivity of patients with COAD to immunotherapy, we identified some targets of immunomodulatory drugs in clinical trials for metastatic colorectal cancer.   We then compared the expression of these immunomodulatory targets between the two clusters and found most of the immunomodulatory targets (PD1, PDL1, PDL2, CTLA4, CD80,  CD86, LAG3, TIM3, TIGHT, OX40, GITR, 4-1BB, ICOS, CD27, and CD70) to be expressed significantly higher in Cluster 2 ( Figures  3A-D). The results suggested that Cluster 2 may show better response to immunotherapy than Cluster 1.

Comparison of Genetic Mutation
Different genetic mutations can influence the efficacy of immunotherapy differently; therefore, we evaluated the mutational conditions in COAD. Landscapes of mutation profiles between the two clusters are shown in Figures 4A, B. Cluster 2 had higher TMB and more numbers of MLH1, MSH2, MSH6, PMS2, POLE, and POLD1 mutations than Cluster 1 ( Figures 4C, D), which once again indicated that the effect of immunotherapy may be better in Cluster 2.

Functional Enrichment of Hub Genes and Their Correlation With Immune Infiltration
To determine the biochemical functions of the 14 hub genes, we performed a GO enrichment analysis ( Figure 6A). The most significant GO term was negative regulation of immune system. We also conducted PPI analysis and checked correlations across the genes (Figures 6B, C). Spearman's correlation analysis between genes and immune infiltration (ESTIMATE and ssGSEA) showed most of the genes to be significantly correlated with immunity ( Figures 6D, E).

GEO Validation of Immune Characteristics Between Two Clusters
First, we divided 566 colon cancer samples from GSE39582 into two clusters in the same way as performed in TCGA ( Figure 7A) and found the distribution of ALKBH5 and YTHDF1 expression in the two clusters to be quite similar to that in TCGA. We also calculated the Spearman's correlation coefficient for ALKBH5 and YTHDF1 (R = −0.36, P = 1.24e-18) ( Figure 7B).  The expression of immunomodulatory targets and extent of immune infiltration (ESTIMATE, CIBERSORT, and ssGSEA) were evaluated in the same manner ( Figures 7C-L). Cluster 2 was clearly more active regarding the immune system than Cluster 1.

Verification of Expression Levels of ALKBH5 and YTHDF1 in CRC and Adjacent Tissues
We tested the expression levels of ALKBH5 and YTHDF1 in 12 CRC tissues and paired adjacent tissues using RT-qPCR. Results indicated CRC tissues to have lower expression of ALKBH5 and higher expression of YTHDF1 than the paired adjacent tissues (Figure 8).

DISCUSSION
ALKBH5 is expressed at low levels in colon cancer; its overexpression inhibits cell metastasis in vivo and cell invasion in vitro, thus suggesting it as a tumor suppressor (34). A recent study has reported that the expression and function of ALKBH5 in different types of cancer are variable (35). Although ALKBH5 has been proven to correlate with the response to anti-PD1 therapy in melanoma, the association between ALKBH5 expression and response to immunotherapy in patients with COAD still remains unclear (15). YTHDF1 is highly expressed and enhances stem cell-like activity in CRC (36). The high expression of YTHDF1 could lead to low immune cell abundance, since high stemness indices represent a low immune cell fraction and low PD-L1 expression (37). The expression of ALKBH5 and YTHDF1 in patients with COAD in our study was consistent with those in previous studies. In addition, we found a negative correlation between ALKBH5 and YTHDF1, suggesting that their functions may have a cross-talk or interaction upon m6A modification. A previous study had reported that ALKBH5 suppresses tumor progression in nonsmall cell lung cancer in a YTHDF1-dependent manner (38). Moreover, their relevance to the immune score in the ESTIMATE analysis was just the opposite. Therefore, both ALKBH5 and YTHDF1 may participate in the regulation of m6A modification, which can in turn influence immune infiltration and response to immunotherapy in patients with COAD. We applied consensus clustering to divide patients with COAD from TCGA into two clusters: Cluster 1 (ALKBH5: low expression; YTHDF1: high expression) and Cluster 2 (ALKBH5: high expression; YTHDF1: low expression). Moreover, we investigated the relationship between the two clusters and the CMS. We found CMS2 to be mostly classified into Cluster 1, whereas CMS1 and CMS3 were mostly classified into Cluster 2. There was no difference in CMS4. CMS1 represented microsatellite instability immune type with hypermutation, MSI, and strong immune activation; CMS2 represented canonical type with remarkable WNT and MYC activation; CMS3 represented metabolic type with epithelial and evident metabolic dysregulation (29). From the perspective of CMS, we inferred Cluster 2 to possibly possess more immune infiltration than Cluster 1. We then evaluated their clinical characteristics and discovered Cluster 1 to have a higher N stage than Cluster 2, which may have resulted from the inhibition of metastasis by ALKBH5.
To further investigate the functional differences between the two clusters, we used TCGA-COAD data for GSEA and found some immune-related pathways, such as adaptive immune response, cell killing, cytokine production, and T cell activation to be enriched in Cluster 2. This suggested that Cluster 2 may act more actively in immune response than Cluster 1.
Next, we compared the immune characteristics of the two clusters using the ESTIMATE, CIBERSORT, and ssGSEA methods. In the ESTIMATE analysis, Cluster 2 was proven to possess higher stromal, immune, and ESTIMATE scores than Cluster 1, thereby implying that Cluster 2 had a vibrant tumor immune microenvironment. In the CIBERSORT analysis, we found the proportion of CD8 T cells and M1 subtype macrophages to be significantly elevated in Cluster 2. Previous studies had demonstrated CD8 T cells to have the strongest effect on patient prognosis in most tumor-infiltrating immune cell subtypes (39). In the ssGSEA analysis, 25 immune cell subtypes showed significantly higher expression in Cluster 2, including

CD8 T cells, T helper cells (CD4), dendritic cells (DCs), natural killer (NK) cells, natural killer T (NKT) cells, and macrophages.
Tumor-infiltrating T cells have a major impact on the clinical attributes of CRC. High infiltration of CD8 T cells can predict the response to drugs and improve survival in patients with CRC and hepatic metastases (40). A previous study had illustrated that patients with high expression of Th1 have a prolonged prognosis, whereas those with high expression of Th17 have poor prognosis in CRC. In addition, the effect of Th1 seemed to surpass the effect of Th17 on survival (41). DCs have been reported as key antigenpresenting cells that promote anti-tumor immunity by activating T cells (42). Moreover, conventional type 1 dendritic cells are known to be recruited into the tumor microenvironment following stimulation by NK cells (43). The latter have cytotoxic capacity in anti-tumor immunity, and their extensive infiltration leads to a favorable outcome in CRC (44). NKT cells could reinvigorate the exhausted CD8 T cells in an anti-PD-1resistant tumor model, hence playing a pivotal role in anti-tumor immunity (45). A previous study had shown that high NKT cell infiltration to be an independent favorable prognostic factor in CRC (46). Macrophages are conventionally divided into M1 (proinflammatory; anti-tumor) and M2 (anti-inflammatory; tumor-promoting) subtypes. According to the results of CIBERSORT analysis, Cluster 2 had a higher proportion of M1 subtype macrophages than Cluster 1, which suggested that Cluster 2 could easily achieve anti-tumor Th1-type responses while Cluster 1 tended to establish a tolerogenic microenvironment (47). Based on our study of immune contexture, Cluster 2 had more extensive immune cell infiltration than Cluster 1. Therefore, Cluster 2 may have more immunological competence and be more likely to benefit from immunotherapy. Previous studies reported that programmed cell death 1 (PD1), programmed cell death 1 ligand 1 (PDL1), and cytotoxic T lymphocyte antigen 4 (CTLA4) are approved as targets of immune checkpoint inhibitors (ICIs) by the FDA (48). In addition, lymphocyte activation gene-3 (LAG3), T cell immunoglobulin-3 (TIM3), and T cell immunoglobulin and The P values are labeled using asterisks (ns, no significance, *P < 0.05, **P < 0.01, ***P < 0.001). ITIM domain (TIGIT) are regarded as co-inhibitory receptor targets (49). In this study, we compared two clusters of immunomodulatory drugs, which have been included in clinical trials for metastatic CRC (48). Most of these targets were found to be significantly high in expression in Cluster 2. Next, we analyzed the mutational landscapes of the two clusters, and found a remarkable difference between them. We found Cluster 2 to have more TMB than Cluster 1. TMB may affect the generation of immunogenic peptides and thereby influence the response to immunotherapy (50). Furthermore, CRC can be categorized into two groups based on microsatellite instability (MSI) and mismatch repair (MMR), namely dMMR-MSI-H and pMMR-MSI-L. The signature of dMMR-MSI-H in patients with CRC is a specific biomarker for evaluating the response to immunotherapy. In addition to the hypermutation caused by dMMR-MSI-H, POLE proofreading domain mutation also leads to a remarkable hypermutation, which may result in excellent prognosis (48). Our study suggested that Cluster 2 possesses a higher mutation rate of MMR-related genes (MLH1, MSH2, MSH6, and PMS2) and POLE/POLD1 than Cluster 1, hence implying that Cluster 2 would show better effect of immunotherapy.
We next performed WGCNA to identify the blue module that differed between the two clusters in relation to both ALKBH5/ YTHDF1 and immune score. Based on MM and GS, we obtained 14 genes from this module, including PD1 (PDCD1) and LAG3. An explicit synergistic interaction between PD1 and LAG3 has already been reported, and they have been shown to mediate T cell exhaustion together (51); this probably occurs in Cluster 2. Therefore, we speculated that Cluster 2 could acquire a better response to immunotherapy than Cluster 1 by evaluating the extent of infiltration extent of immune cells, expression of ICIs, and mutational burden; this difference might result from m6A modification mediated by ALKBH5 and YTHDF1.
Based on the immunological features, Cluster 2 was considered to have a hot tumor, whereas Cluster 1 tended to be have a cold tumor (52). Since our RT-qPCR results verified the expression of ALKBH5 and YTHDF1 in 12 pairs of cancer and normal tissues, we inferred Cluster 1 to represent the overall immunological characteristics of COAD, thus showing poor immune response. ALKBH5 and YTHDF1 could possibly play a potential role in the transformation of cold to hot tumor in COAD.
Although the comprehensive analysis improved our understanding of the relationship between ALKBH5/YTHDF1 and immunity, and we used 566 patients with GSE39582 as the external validation set, there are still some limitations in the current study. First, it was a retrospective study. Therefore, a prospective study should be conducted in future in order to avoid analysis bias associated with retrospective studies. Moreover, the study was performed based on TCGA and GEO; we could not illustrate the expression of ALKBH5 and YTHDF1 from the protein level or demonstrate the direct mechanisms of ALKBH5/ YTHDF1 in anti-tumor immunity. So further studies to unravel the direct mechanisms should be performed.

CONCLUSION
By clustering patients of TCGA-COAD and GSE39582 based on the expression of ALKBH5 and YTHDF1, we demonstrated Cluster 2 (ALKBH5: highly expressed; YTHDF1: lowly expressed) to have more infiltration of immune cells, expression of ICI targets, TMB, and proportion of dMMR-MSI-H than Cluster 1 (ALKBH5: lowly expressed; YTHDF1: highly expressed), thereby suggesting that Cluster 2 acquired better response to immunotherapy. Our findings illustrated that ALKBH5 and YTHDF1 have potential in tumor immunity and provided novel insights into the relationship between m6A modification and immunity.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding authors.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the institutional ethics board of the First Hospital of China Medical University. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
GY designed the study and wrote the manuscript. YA performed the literature search and collected data for the manuscript. BX analyzed the data. NW edited the figures and tables. MS and XS FIGURE 8 | Verification of ALKBH5 and YTHDF1 expression in CRC tissues using RT-qPCR. The P values are labeled using asterisks (ns, no significance, *P < 0.05, ***P < 0.001).
revised the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by the National Natural Science Foundation of China (81670506).