key: cord-1019992-7gsiyj0z
authors: Cai, Changjing; Long, Jie; Huang, Qiaoqiao; Han, Ying; Peng, Yinghui; Guo, Cao; Liu, Shanshan; Chen, Yihong; Shen, Edward; Long, Kexin; Wang, Xinwen; Yu, Jian; Shen, Hong; Zeng, Shan
title: M6A “Writer” Gene METTL14: A Favorable Prognostic Biomarker and Correlated With Immune Infiltrates in Rectal Cancer
date: 2021-06-17
journal: Front Oncol
DOI: 10.3389/fonc.2021.615296
sha: a41ad42a79b8f211da1668f29c4e29f4c68dd65e
doc_id: 1019992
cord_uid: 7gsiyj0z

Rectal cancer (RC) is the leading cause of tumor-related death among both men and women. The efficacy of immunotherapy for rectal cancer is closely related to the immune infiltration level. The N6-methyladenosine (m6A) modification may play a pivotal role in tumor-immune interactions. However, the roles of m6A-related genes in tumor-immune interactions of rectal cancer remain largely unknown. After an evaluation on the expression levels of m6A-related genes and their correlations with the prognosis of rectal cancer patients, we found that METTL14 was the only gene to be significantly correlated with prognosis in rectal cancer patients. Therefore, we further observed the impact of METTL14 expression and m6A modification on the immune infiltration in rectal cancer. Our study indicates that low expression of the m6A “writer” gene METTL14 in rectal cancer may lead to the downregulation of m6A RNA modification, thus reducing the level of immune cell infiltration and resulting in poor prognosis. METTL14 expression level is an independent prognostic factor in rectal cancer and is positively correlated with the immune infiltration level. Our study identified METTL14 as a potential target for enhancing immunotherapy efficacy in rectal cancer.

Rectal cancer (RC) is the leading cause of tumor-related death among both men and women (1) . Thus, clinicians are calling for more effective treatments for rectal cancer. Immune checkpoint inhibitors have achieved successful responses in few rectal cancer patients (2, 3) . Mismatch repair (MMR) status functions as a major predictor of the efficacy of immune checkpoint inhibitor therapy, which has many limitations: MMR deficiency is found in 10% to 15% of colorectal cancer (CRC) patients, among which only 30% to 50% of them can benefit from the immunotherapy (4) (5) (6) (7) . Therefore, there is an urgent need to identify new biomarkers that can accurately predict the immunotherapy response of rectal cancer patients, to reveal resistance mechanisms, and to seek potential targets for enhancing immunotherapy efficacy.

The efficacy of immunotherapy for rectal cancer is closely related to the immune infiltration level in RC (2, 8) . N6methyladenosine (m6A) modification, the most common internal modification of messenger RNAs (mRNAs) in eukaryotes (9) , is a reversible event modulated by "writers" (WTAP, KIAA1429, RBM15, RBM15B, METTL3, METTL16, and METTL14), "erasers" (FTO and ALKBH5), and "readers" (HNRNPA2B1, HNRNPC, YTHDF1, YTHDF2, YTHDF3, and YTHDC1) (10). Li et al. (11) reported that m6A modification controls T cell homeostasis by targeting the IL-7/STAT5/SOCS pathways. Han et al. (12) reported that YTHDF1-dependent m6A mRNA methylation controls antitumor immunity mediated by dendritic cells (DCs). These findings suggest that m6A modification plays a crucial role in tumor-immune interactions. Meanwhile, several studies indicated altered expression of m6A-related genes in gastrointestinal cancers (13) . However, the roles of m6A-related genes in tumorimmune interactions of rectal cancer remain largely unknown.

In this study, we assessed the expression levels of m6A-related genes and their correlations with the prognosis of rectal cancer patients. We observed that METTL14 expression is positively correlated with overall survival (OS) and tumor immune cell infiltration. Our work suggests an important role of METTL14 in regulating tumor and immune microenvironment interactions, and identifies METTL14 as a prognostic biomarker as well as a potential target for enhancing the immunotherapy effect in rectal cancer.

The expression levels of m6A-related genes in rectal cancer were obtained from four GEO data sets (GSE123390, GSE87211, GSE60331, GSE68204) and other 3 TCGA data sets from TIMER (14) , GEPIA (15) and UCSC Xena (16) for further study (data set A, data set B, and data set C).

The relationship between m6A-related gene expression levels and the prognosis of rectal cancer patients was first explored in the TIMER (14) database (https://cistrome.shinyapps.io/timer/) and GEPIA (15) database (http://gepia.cancer-pku.cn/). TIMER is a web resource that includes 10,897 samples from 32 cancer types from TCGA, and users can explore the clinical effect of genes in the "Survival" module. GEPIA is a web resource that contains 9,736 tumors and 8,587 normal samples from TCGA and the GTEx database, and users can conduct overall survival analysis of an input gene in specific cancers. The threshold for splitting the high-expression and low-expression cohorts can be adjusted. In our study, 50% of patients had a higher expression level than the threshold in both TIMER and GEPIA.

The "GENE" module of TIMER allows users to easily explore the association between the expression of a certain gene and the level of infiltration of multiple immune cell types, including macrophages, neutrophils, dendritic cells, B cells, CD4+ T cells, and CD8+ T cells, in a given cancer type. Because tumor purity is an important confounding factor, the first panel of this analysis displayed METTL14 expression levels against tumor purity. GEPIA provides an interface to conduct gene correlation analysis by using Pearson, Spearman, or Kendall methods; thus, we further explored the correlation between METTL14 expression and markers of diverse immune cells in the GEPIA database. We chose the Spearman method and used TCGA tumor and TCGA normal data sets. TISIDB (17) (http://cis. hku.hk/TISIDB) is a user-friendly platform to investigate the role of a certain gene in tumor-immune interactions. Tumor-infiltrating lymphocyte (TIL), immunoinhibitor, immunostimulator, chemokine, chemokine receptor, major histocompatibility complex (MHC), immune subtype, and molecular subtype analyses were performed.

The rectal cancer cell line cannot be acquired in China amid the pandemic of COVID-19, and human colorectal cancer cell line (HCT116) is generally considered to be a representative colorectal cancer cell line, and we found several studies focusing on rectal cancer, which employed HCT116 to conduct research in vitro (Published in Annals Of Surgery (18), British Journal Of Cancer (19) . We have tried to find normal rectal cell lines when conducting the research as well, and the only normal rectal cell line in ATCC named Hs 680.Rec was not available for international distribution. We found several studies focusing on rectal cancer which employed human normal colorectal epithelial cell line (FHC) to conduct research in vitro [Published in Biomedicine & Pharmacotherapy (20) , Journal of cellular biochemistry (21) ]. As far as we know, rectal cancer and Abbreviations: RC, rectal cancer; m6A, N6-methyladenosine; CNV, copy number variation; MMR, mismatch repair; CRC, colorectal cancer; DC, dendritic cells; OS, overall survival; READ, rectum adenocarcinoma; TIL, tumor infiltrating lymphocyte; Th, T helper; Tfh, T follicular helper cell; Treg, regulatory T cell; MHC, major histocompatibility complexes; iDC, immature dendritic cell; TGF-b1, transforming growth factor-b1; TGF-bR1: transforming growth factor-b receptor 1; TNFSF15: tumor necrosis factor superfamily-15; TNFRSF25, tumor necrosis factor receptor superfamily 25; CCL, C-C motif chemokine ligand; CXCL, C-X-C motif chemokine ligand; CCR, C-C chemokine receptor; CXCR, C-X-C motif chemokine receptor; HLA, human leucocyte antigen; TAP, antigen peptide transporter; TAPBP, tapasin binding protein; B2M, b-2 microglobulin; CIN, chromosomal instability; GS, genome stable; HM-SNV, hypermutated single nucleotide variants; HM-indel, hypermutated insertion-deletion. colon cancer cells have similar biological behaviors (22), so we used HCT116 and FHC in our study. HCT116 and FHC were obtained from the Cell Bank of the Chinese Academy of Sciences and were cultured in RPMI 1640 (GIBCO, Grand Island, NY) supplemented with 10% fetal bovine serum, 100 U/mL penicillin, and 100 µg/mL streptomycin (GIBCO). Cells were maintained in a humidified incubator with a 5% CO 2 atmosphere at 37°C.

Arraystar Human m6A-mRNA Epitranscriptomic Microarray Analysis Total RNA was extracted from HCT116 cells and FHC cells using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer's instructions and was immunoprecipitated with anti-N6-methyladenosine antibody (Synaptic Systems, 202003). The modified RNAs were eluted from the immunoprecipitated magnetic beads as the "IP." The unmodified RNAs were recovered from the supernatant as "Sup." The IP and Sup RNAs were labeled with Cy5 and Cy3, respectively, as cRNAs in separate reactions using the Arraystar Super RNA Labeling Kit (Arraystar, AL-SE-005). The cRNAs were combined and hybridized onto an Arraystar Human mRNA Epitranscriptomic Microarray (8x60K, Arraystar). After washing the slides, the arrays were scanned in two-color channels by an Agilent Scanner G2505C. Agilent Feature Extraction software (version 11.0.1.1) was used to analyze acquired array images. Differential m6A-methylated mRNAs between two samples were identified through fold change filtering and were analyzed for Gene Ontology (GO) and pathway enrichment by using GO (http://www.geneontology. org) and KOBAS (http://kobas.cbi.pku.edu.cn/index.php).

Tissue samples of 150 patients diagnosed with rectal cancer from July 2015 to January 2018 in Xiangya Hospital of Central South University were collected to establish a tissue bank. The demographic characteristics, cancer stages, and pathological reports were obtained from the electronic medical records (EMR) system. A retrospective cohort was performed, as of July 31, 2019, a total of 89 patients were included. Survival analysis, multivariate Cox regression, and the clinical features analysis were conducted after follow-up. This study was reviewed and approved by the Xiangya Hospital Medical Ethics Committee of Central South University.

Formalin-fixed, paraffin-embedded tissue array slides were used to detect METTL14 protein expression. Briefly, after deparaffinization and rehydration, citrate buffer (ZLI-9064; ZSGB-BIO, Beijing, China) was used for heat-induced epitope retrieval. Endogenous peroxidase activity was inhibited for 10 min with 3% hydrogen peroxide (reagent 1; PV-9000; ZSGB-BIO). Nonspecific binding was blocked with normal goat serum (ZSGB-BIO) for 15 to 20 min at room temperature. The slides were then incubated overnight at 4°C with METTL14 rabbit polyclonal antibody at a dilution of 1:100 (catalog no. 26158-1-AP; Proteintech, Chicago, USA). Next, the slides were incubated in polymer helper (reagent 2; PV-9000; ZSGB-BIO) for 20 min at 37°C and then incubated with polyperoxidase-anti-mouse/rabbit IgG (reagent 3, PV-9000, ZSGB-BIO) for 20 min at 37°C. 3-3'-diaminobenzidine was used for coloration, and hematoxylin was used for counterstaining.

METTL14 staining was defined as positive when rectal cells showed nuclear staining. A METTL14 staining score was defined by adding the staining intensity score and the positive staining percentage score. The staining intensity was categorized into three grades: score 1, yellow; score 2, light brown; score 3, brown. Positive staining percentage patterns were categorized into four groups: score 1, 0% to 25% staining of rectal cells; score 2, 25% to 50% staining of rectal cells; score 3, 50% to 75% staining of rectal cells; score 4, 75% to 100% staining of rectal cells. The percentage and intensity scores were added, and the results were classified into a high expression group (total scores >4) and a low expression group (total scores ≤4) by using a mean score of 4 as cutoff.

Analysis of differentially expressed genes (DEGs) of METTL14 expression level (high vs. low) was based on an empirical Bayesian approach using the limma R package (23) . An adjusted P value less than 0.05 and log2 fold change (log2FC) greater than 1 indicated statistical significance for further gene ontology (GO) enrichment.

The pathway activation levels (PAL) calculated from gene expression data between METTL14-High vs METTL14-Low groups were performed. Both the single sample gene set enrichment analysis (ssGSEA) algorithm (24) and Oncobox library (oncoboxlib) (25) were used for PAL calculation. ssGSEA was performed based on "c2.cp.hallmark.v7.1.symbols" gene sets downloaded from MSigDB database using the "ssGSEA" R package. Oncoboxlib calculates PAL according to Sorokin et al. (25) , and it takes a file that contains gene symbols in HGNC format, their expression levels for one or more samples (cases and/or controls) and calculates PAL values for each pathway in each sample (Table S1 ).

Statistical data were analyzed by using SPSS 20.0 and presented by using GraphPad Prism 7.0. Kaplan-Meier curves and log-rank tests were performed to evaluate the prognostic value of METTL14 in rectal cancer. Overall survival (OS) was measured from the date of diagnosis to either the date of death or the date of the final follow-up, with final evaluation on July 31, 2019. Multivariate Cox regression analysis was performed to identify whether METTL14 was an independent prognostic factor of survival in rectal cancer. Multivariate logistics regression analysis was performed to identify the impact factors of METTL14 expression level (high vs. low) in rectal cancer. The results with a P value less than 0.05 were considered statistically significant if not specified.

We used four GEO data sets (GSE123390, GSE87211, GSE60331, GSE68204) and other two TCGA data sets (data sets A and B) to determine the expression level of m6A-related genes in rectal cancer and normal tissues. Besides the data sets without the statistically significant difference, the results showed that compared to normal tissues, the expression levels of KIAA1429, METTL3, METTL16, and HNRNPA2B1 increased in rectal cancer, and the expression of METTL14 and ALKBH5 in rectal cancer was significantly decreased. The expression levels of other nine genes showed different alterations in the tumors using different data sets (Figures 1 and S1 ).

The mRNA expression levels of m6A-related genes suggested that they might play an important role in rectal cancer. To determine whether the expression level of m6A-related genes is correlated with the prognosis in rectal cancer patients, we selected the latest TCGA data set (data set C), to conduct retrospective studies. In the univariate cox analysis, RBM15, YTHDF2, and METTL14 were significantly correlated with prognosis in rectal cancer patients ( Figure 2A) . Nonetheless, in the log-rank test of RBM15, YTHDF2 and METTL14, only METTL14 were correlated with OS in these patients ( Figures 2B-D) . The hazard ratio (HR) for patients with high METTL14 expression levels was 0.464 (p = 0.0266, Figure 2A ), which showed that patients with high METTL14 expression had longer OS time and better prognosis than those with low METTL14 expression. These results indicated that only METTL14 expression was significantly correlated with prognosis in rectal cancer patients, therefore we further observed the impact of METTL14 expression and m6A modification on the immune cell infiltration in rectal cancer.

The tumor immune infiltration level plays a crucial part in affecting the prognosis of rectal cancer patients (26) . Therefore, we evaluated whether METTL14 expression was related to immune infiltration levels in rectal cancer using the TIMER and GEPIA databases. The results from the TIMER database showed that METTL14 expression level had a significant positive correlation with the levels of B cell (cor. = 0.244, p=0.00377), CD8+ T cell (cor. = 0.56, p<0.0001), macrophage (cor. = 0.244, p=0.00379), neutrophil (cor. = 0.301, p=0.000331), and dendritic cell infiltration (cor. = 0.213, p=0.0119) in rectum adenocarcinoma (READ), while METTL14 expression was not significantly correlated with CD4+ T cell infiltration (cor. = 0.042, p = 0.619) ( Figure 3A) . To verify the results of the TIMER database analysis, we then investigated the correlation between METTL14 expression level and immune marker sets of diverse immune cells in the GEPIA database. As expected, the results were basically consistent with those of the TIMER database analysis. Immune markers of CD8+ T cells, B cells, M1 macrophages, M2 macrophages, neutrophils and dendritic cells were positively correlated with METTL14 expression in READ. In addition, in the GEPIA database, the immune markers of CD4+ T cell, T cell (general), monocyte, tumor-associated macrophage (TAM), T helper 1 (Th1) cell, Th2 cell, Th17 cell, T follicular helper cell (Tfh), regulatory T cell (Treg), natural killer cell, and T cell exhaustion also showed positive correlations with METTL14 expression levels ( Table 1) .

After confirming that METTL14 expression was positively correlated with immune activation signature, we further used the TISIDB database to analyze the correlations between FIGURE 1 | The relative expression of m6A-related genes in rectal cancer. Heatmap showing the alterations in the mRNA expression of m6A-related genes in the TCGA and GEO data sets. The red color indicates upregulated expression; the blue color indicates downregulated expression; the blank indicates no significant changes, and the black color indicates that the related gene is absent in the data sets. The data were statistically analyzed by Student's t test (unpaired, two-tailed). METTL1 4 expression and TILs, immunoinhibitors, immunostimulators, chemokines, chemokines-receptors, and major histocompatibility complexes (MHCs). P-values less than 0.05 are considered to be significant. The results of TIL analysis were basically consistent with those of the TIMER and GEPIA database analyses. The levels of CD4+ T cell, Th2 cell, and immature dendritic cell (iDC) infiltration were positively correlated with METTL14 expression, and the CD4+ T cell infiltration level had the strongest correlation (cor. = 0.4, p < 0.0001). In contrast to the GEPIA results, the level of CD56 dim nature killer cell infiltration (cor.=-0.377, p<0.0001) and monocyte infiltration (cor.= -0.185, p=0.017) was negatively related with METTL14 expression ( Figure 3B ). The results of immunoinhibitor analysis revealed that METTL14 expression negatively correlated with four of nine immunoinhibitors' expression. Among them, transforming growth factor-b1 (TGF-b1) had the strongest correlation coefficient (cor. = −0.314, p < 0.0001). Additionally, five out of nine immunoinhibitors were positively correlated, and transforming growth factor-b receptor 1 (TGF-bR1) had the strongest correlation (cor. = 0.304, p < 0.0001) ( Figure 3C ). Immunostimulator analysis showed that METTL14 expression was negatively correlated with most immunostimulators' expression (9/15), and tumor necrosis factor receptor superfamily 25 (TNFRSF25) had the most significant negative correlation (cor. = −0.468, p < 0.0001). However, a few immunostimulators (6/15) had positive correlations, and tumor necrosis factor superfamily-15 (TNFSF15) had the strongest correlation (cor. = 0.349, p<0.0001) ( Figure 3D) .

Chemokine analysis showed that all C-C motif chemokine ligand (CCL) family members and METTL14 expression were negatively correlated. CCL23 had the strongest correlation (cor. = −0.378, p < 0.0001). Three of four C-X-C motif chemokine ligand (CXCL) family members and METTL14 expression were positively correlated, and CXCL13 had the strongest correlation (cor. = 0.271, p=0.000424) ( Figure 3E) . As for the C-C chemokine receptor (CCR) family, chemokine receptor analysis revealed that CCR10 had a negative correlation with METTL14 expression (cor. = −0.474, p < 0.0001), and CCR5 had a positive correlation with METTL14 expression (cor. = 0.202, p = 0.00892). In the C-X-C motif chemokine receptor (CXCR) family, there is a negative correlation between CXCR3 and METTL14 (cor. = −0.343, p < 0.0001) ( Figure 3F) . MHC analysis showed that all human leucocyte antigen (HLA) family members and METTL14 expression were negatively correlated, and HLA-G had the strongest correlation (cor. = −0.386, p < 0.0001). Antigen peptide transporter (TAP) family members were also negatively correlated with METTL14 expression, and tapasin binding protein (TAPBP) had the strongest correlation (cor. = −0.364, p < 0.0001). However, b-2 microglobulin (B2M) was positively correlated with METTL14 expression (cor. = 0.211, p = 0.00618) ( Figure 3G) . We also assessed whether METTL14 expression was related to immune subtypes (27) and molecular subtypes (28) in READ. The results showed that there was no significant difference in METTL14 expression among the six immune subtypes (p = 0.393), while there were differences among molecular subtypes (p = 0.00145). The expression level of METTL14 in the chromosomal instability (CIN) subtype was the lowest, while METTL14 expression in genome stable (GS), hypermutated single nucleotide variants (HM-SNV) and hypermutated insertion-deletion (HM-indel) subtypes was higher than that in CIN ( Figures 3H, I) .

To verify the impact of METTL14 on immune infiltration, we cultured HCT116 cells, a human colorectal cancer cell line, and FHC cells, a human normal intestinal epithelial cell line, and conducted Arraystar Human m6A-mRNA Epitranscriptomic microarray analysis. Since METTL14 functions as an m6A writer gene, its low expression in rectal cancer suggests that the m6A level in rectal cancer cells may be lower than that in normal cells. Our results showed that the expression level of METTL14 in HCT116 is significantly lower than that in FHC. In HCT116, there are 1103 genes with reduced m6A methylation levels and 509 genes with increased m6A methylation levels. About two thirds of the genes take on reduced m6A methylation levels ( Figure 4) .

We analyzed the microarray results and selected all 1103 genes with reduced m6A methylation levels in HCT116 cells for KEGG and GO analysis. The results of GO analysis showed that these genes were significantly enriched in multiple immune functions, and the five functions with the highest enrichment were immune system process, immune response, immune effector process, activation of immune response, and regulation of immune response ( Figure 5A ). KEGG analysis showed that these genes were enriched in multiple immunerelated signaling pathways, and the pathways with the highest enrichment were inflammatory mediator regulation of TRP channels, human T-cell leukemia virus-1 infection, Toll-like receptor signaling pathway, and natural killer cell mediated cytotoxicity. Enriched genes are highlighted in red in the figure  (Figures 5B-F) .

To validate above findings in silico analysis, we conducted a retrospective study. There were 150 rectal cancer tissues and 150 adjacent nontumor tissues in our tissue bank, upon which immunohistochemistry staining for METTL14 was performed using tissue microarrays. We found that 69 patients had high METTL14 expression and 81 patients had low METTL14 expression. Compared with adjacent nontumor tissues, the expression levels of METTL14 protein in rectal cancer tissues significantly decreased (p < 0.0001), which was consistent with the results from the public database ( Figures 5G-I) . A retrospective cohort was performed, and the clinical information of 89 patients was collected. After follow-up, we conducted survival analysis of the cohort. The results showed that patients with high METTL14 expression had longer OS (p = 0.004), and Cox survival analysis identified the expression level of METTL14 as an independent prognostic factor of patients with rectal cancer. The HR value of the patients with high expression levels was 0.077, and the P value was 0.0164. Advanced clinical stages were associated with decreased survival time (HR: 8.392, p=0.0439), and age, gender, and differentiation level were not found to be statistically significant related to survival ( Figure 5J and Table 2 

To investigate the correlation between clinical features and METTL14 expression level, in TCGA data set C, we analyzed the correlation between MSI status, TMB, neoantigen, k-Ras mutations and METTL14 expression level, and there was no statistically significant difference ( Figures 6A-D) . According to the results of IHC, the expression level of METTL14 were classified into a high expression group and a low expression group. We further conducted a multivariate logistics regression analysis in our validation cohort to quest for the relationship between age, gender, clinical stage, differentiation level, and METTL14 expression level (high expression group vs. low expression group), and no statistically significant differences were found ( Table 3) .

To further explore the mechanisms of METTL14, in TCGA (data set C), we performed DEG analysis for METTL14-High vs METTL14-Low groups and identified a total of 11214 DEGs. There were 1080 DEGs with Log FC >1 and P<0.05. GO analysis showed that the top 20 functions of these DEGs were mostly tied up with tumor-related functions and cell development-related functions ( Figures 6E-F) .

To further analyze the impact of the expression level of METTL14 on signal pathways, we used ssGSEA algorithm (24) and Oncobox library (oncoboxlib) (25) to compare PAL calculated from gene expression data between METTL14-High and METTL14-Low groups. The ssGSEA results showed that the top 20 activated pathways were basically tumor-related pathways, including immunity pathways (TGF-b), cell junctions pathway (adherens junction), and metabolism pathway (mTOR signaling pathway) ( Figure 7A ). The results in oncoboxlib included more than 3,000 signal pathways, which was different from ssGSEA. In METTL14-High group, the top 20 pathways with high PAL were primarily immune-related, migration-related and apoptosis-related pathways ( Figure 7B) . In METTL14-Low group, the top 20 pathways with high PAL were mainly related to the negative regulation of apoptosis, cell migration, and immunity ( Figure 7C) . Compared with ssGSEA algorithm, the results from oncoboxlib reflect the biological function of METTL14 more properly, which is consistent with the results of our clinical survival cohort. In sum, METTL14 is a good prognostic factor in rectal cancer.

Notably, decreased METTL14 expression was found to be significantly correlated with poor prognosis in rectal cancer patients. Furthermore, our study indicates that the low expression of the m6A writer gene METTL14 in rectal cancer may lead to a decrease in m6A RNA modification, thus reducing the level of immune cell infiltration and resulting in poor prognosis. This mechanism for our findings needs further validation.

The analysis of the expression of m6A-related genes in rectal cancer showed that the expression levels of KIAA1429, METTL3, METTL16, and HNRNPA2B1 were higher and the expression levels of METTL14 and ALKBH5 in rectal cancer were significantly lower than those in normal tissues, showing that aberrant expression of m6A regulatory genes commonly occurs in rectal cancer.

Although many m6A-related genes were dysregulated in rectal cancer, only METTL14 was related to the prognosis of rectal cancer patients. The TCGA data sets and our retrospective cohort confirmed that METTL14 had a lower expression level in rectal cancer tissues than in normal adjacent tissues, and that patients with low METTL14 expression had shorter OS times than those with high METTL14 expression, which is consistent with previous study (29) . Multivariate survival analysis (COX) of our cohort revealed METTL14 as an independent prognostic factor in rectal cancer patients. These findings strongly suggest that METTL14 is a prognostic biomarker in rectal cancer.

Immunogenomic analysis of over 10,000 tumors in TCGA data sets revealed that tumor immune landscapes differ greatly between and within cancer types (30, 31) . CRC development is driven by Wnt/Myc hyperactivation, KRAS/BRAF mutation, genetic instability and accompanied by progressive immunosuppressive tumor microenvironment (TME) (32) . CRCs have poor response to immune checkpoint inhibitors, expect for those with MMR deficiency or high levels of T-cell infiltration (32). Han D. et al. (12) reported that m6A-related genes can affect tumor progression through the tumor immune microenvironment, but their roles in rectal cancer remain unknown. Our data support a potential role of METTL14 in rectal cancer with increased antitumor immunity. Using the TIMER, GEPIA, and TISIDB databases, our results revealed the positive correlation of the expression level of METTL14 with the immune infiltration level in rectal cancer, including gene signature of CD8+ T cell, CD4+ T cell, B cell, DCs, macrophages, neutrophils and Th2, which were consistent in at least two databases mentioned above. It is well established that increased immune infiltration is conducive to more favorable prognosis in many types of cancers such as melanoma, head and neck, breast, bladder, urothelial, ovarian, renal, prostatic, lung, chordoma, and colorectal cancer (26, (33) (34) (35) , and infiltrating immune cells mainly consist of the T lymphocyte family, such as CD8+ T cell and CD4+ T cell (36) (37) (38) (39) . In addition to T cells, B cells also play an important role in antitumor immunity. Some studies have shown that B cells, the key cells in the humoral immune response, can act on cells and molecules in the tumor We expanded our in silico analysis to examine the correlation of METLL14 expression with several classes of immune modulators. METTL14 expression showed a negative association with TGFb1, which is known to promote tumor immune tolerance by upregulating PD1 expression on T cells (46) . The expression level of TGF-bR1, the receptor of TGF-b1, was positively correlated with METTL14 expression, which might be explained by a compensatory increase caused by the decrease in TGF-b1 expression levels. Immunostimulators analysis revealed that METTL14 expression was positively correlated with the expression level of some immunostimulators, and TNFSF15 had the strongest correlation coefficient. Zhao et al. (47) found that TNFSF15 had great importance in the DC-involved Th9 differentiation process, and Th9 cells contributed to antitumor immunity. However, as an immunostimulator, TNFRSF25 had a negative correlation with METTL14 expression. TNFRSF25, also known as death receptor 3, is mainly expressed on the surface of T cells. Interestingly, TNFRSF25 plays a complex dual role in tumors: on the one hand, it binds to its ligand TL1A to activate T cells and inhibit tumor progression (48) ; on the other hand, through the PI3K/NF-kB pathway, it can reduce the apoptosis of colon cancer cells and promote tumor proliferation and metastasis (49) . Chemokine analysis showed that CCL family members, especially CCL23, showed strong negative correlation with METTL14 expression, whereas most CXCL family members, especially CXCL13, showed strong positive correlation. Hannah H. Yan et al. (50) found that CCL23 was mainly secreted by CD33+ myeloid cells and could act on the TGF-b signaling pathway together with CCL9 to promote the progression and metastasis of breast cancer. The role of CXCL13 in tumors is controversial. It can not only promote the invasion and metastasis of tumor cells, but also activate immune cells, increase the immune infiltration level in cancer, and inhibit tumor proliferation (51, 52) . Chemokine receptor analysis showed that CCR5 had a weak positive correlation with METTL14 expression, and CCR10 had a negative correlation with METTL14 expression. CCR5 is widely involved in tumor proliferation and metastasis, and anti-CCR5 therapy has made some progress in various tumors (53) , suggesting that CCR5 and METTL14 play the opposite role in cancer. However, due to the weak correlation (cor. = 0.202), their relationship needs further verification. Hao-yu Lin et al. (54) found that CCR10 promoted the metastasis and invasion of breast cancer through the ERK1/2/ MMP-7 signaling pathway. In the CXCR family, CXCR3 had a negative correlation with METTL14. The role of CXCR3 in tumors is also under dispute. Some studies pointed out that CXCR3 promoted tumor proliferation and invasion through autocrine mechanisms, while other studies showed that CXCR3 inhibited tumor growth by promoting the differentiation of immune cells and activating immune cells (55) .

The results of MHC analysis showed that all HLA family members and METTL14 were negatively correlated in terms of expression, and HLA-G had the highest correlation coefficient. HLA-G has been reported to be a key molecule in tumor immune tolerance and is associated with poor prognosis in cancer patients (56) . TAPBP, a member of the TAP family, has the strongest negative correlation with METTL14. Inconsistent with our results, previous studies have shown that TAPBP participates in the antigen presentation pathway, mediates the immune process, and inhibits tumor proliferation (57), suggesting complex interaction mechanism between METTL14 and TAPBP. B2M was the only MHC gene that was significantly positively correlated with METTL14. Several studies have reported that the loss or mutation of the B2M gene was one reason for tumor immune escape. Considering that B2M loss is a mechanism for resistance to anti-PD-1 therapy, METTL14 may be associated with anti-PD-1 efficacy (58) .

The expression level of METTL14 in rectal cancer varies according to molecular subtype. The chromosomal instability (CIN) subtype had a low expression level of METTL14, which was consistent with previous studies. CIN can cause tumor immune escape, drug resistance, and metastasis, leading to poor clinical outcomes (59) .

In conclusion, the impact of METTL14 on the immune infiltration level of rectal cancer may be related to TGFb1, TNFSF15, CCL23, CCR10, HLA-G, B2M, and CIN. These findings suggest an interesting possibility that METTL14 expression is associated with anti-PD-1 efficacy. Further studies are needed to determine the cause of reduced METTL14 expression, and how it contributes to suppressive immune TME through abovementioned mediators.

As an m6A writer gene (60) , METTL14 may affect the survival of rectal cancer patients by downregulating m6A methylation levels. Human cell m6A-mRNA epitranscriptomic microarray analysis showed that m6A methylation levels of over 1,000 genes were reduced in rectal cancer cells. Further analysis showed that GO terms and pathways were significantly enriched in immunerelated functions and pathways. In the human T-cell leukemia virus-1 infection signaling network, the TGF-b pathway, NF-kB pathway, TNF pathway, P53 pathway, and PI3K/AKT pathway were enriched, which was consistent with the results we mentioned earlier. Macrophages and mast cells are the main cells involved in the inflammatory mediator regulation of TRP channels, and T cells are the major activated cells in the Toll-like receptor signaling pathway, which is consistent with the infiltrating immune cells mentioned before. However, in the natural killer cell-mediated cytotoxicity signal network, we found that METTL14 was negatively correlated with the CD56 dim nature killer cell infiltration level; therefore, a negative regulatory relationship may exist between METTL14 and NK cells. Since the m6A microarray results are based on the single experiment, so the further validation is needed. Our research preliminarily explores the mechanism by which METTL14 regulates the immune infiltration of rectal cancer and provides inspiration for further research.

The mechanisms of METTL14's role in rectal cancer are not only associated with immune system. The agnostic analysis between METTL14-High and METTL14-Low groups showed that besides immunity, the expression of METTL14 was also widely associated with cell apoptosis, adhesion, migration, and cell development process, suggesting that METTL14 plays a multifaceted role in rectal cancer, which needs further validation. Some studies have reported that METTL14 inhibited tumorigenicity, CRC cells growth, invasion, migration, and metastasis, which is consistent with our findings (61-63).

METTL14 expression level is an independent prognostic factor in rectal cancer and is positively correlated with the immune infiltration level. Furthermore, we identified METTL14 as a potential target for enhancing immunotherapy efficacy in rectal cancer.

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

The studies involving human participants were reviewed and approved by Xiangya Hospital Medical Ethics Committee of Central South University. 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.

SZ, HS, and JY designed the study. SZ and HS supervised the study. CC and JL performed the experiments. CC, JL, QH, YP, SL, KL, CG, XW, YC, and YH analyzed and interpreted the data. ES and CC performed the statistical analysis. CC, JL, SZ, HS wrote the manuscript. All authors contributed to the article and approved the submitted version. 

Cancer Statistics, 2020

Neoadjuvant Immunotherapy Leads to Pathological Responses in MMR-Proficient and MMR-Deficient Early-Stage Colon Cancers

Immunotherapy in Colorectal Cancer: Rationale, Challenges and Potential

PD-1 Blockade in Tumors With Mismatch-Repair Deficiency

Nivolumab in Patients With Metastatic DNA Mismatch Repair-Deficient or Microsatellite Instability-High Colorectal Cancer (Checkmate 142): An Open-Label, Multicentre, Phase 2 Study

Durable Clinical Benefit With Nivolumab Plus Ipilimumab in DNA Mismatch Repair-Deficient/Microsatellite Instability-High Metastatic Colorectal Cancer

Tumor Mutational Burden Is Predictive of Response to Immune Checkpoint Inhibitors in MSI-High Metastatic Colorectal Cancer

Comprehensive Analyses of Tumor Immunity: Implications for Cancer Immunotherapy

N6-Methyladenosine-Dependent Regulation of Messenger RNA Stability

Functions of N6-Methyladenosine and Its Role in Cancer

Mrna Methylation Controls T Cell Homeostasis by Targeting the IL-7/STAT5/ SOCS Pathways

Anti-Tumour Immunity Controlled Through mRNA Ma Methylation and YTHDF1 in Dendritic Cells

The Potential Role of RNA N6-Methyladenosine in Cancer Progression

TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells

GEPIA: A Web Server for Cancer and Normal Gene Expression Profiling and Interactive Analyses

Visualizing and Interpreting Cancer Genomics Data Via the Xena Platform

TISIDB: An Integrated Repository Portal for Tumor-Immune System Interactions

Understanding and Resetting Radiation Sensitivity in Rectal Cancer

Comprehensive Molecular Exploration Identified Promoter DNA Methylation of the CRBP1 Gene as a Determinant of Radiation Sensitivity in Rectal Cancer

Long Noncoding RNA LINC00461 Mediates Cisplatin Resistance of Rectal Cancer Via Mir-593-5p/CCND1 Axis

LncRNA NKILA Correlates With the Malignant Status and Serves as a Tumor-Suppressive Role in Rectal Cancer

Cancers of the Colon and Rectum: Identical or Fraternal Twins?

Limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies

Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles

Algorithmic Annotation of Functional Roles for Components of

Prognostic and Predictive Values of the Immunoscore in Patients With Rectal Cancer

The Correlation Between Immune Subtypes and Consensus Molecular Subtypes in Colorectal Cancer Identifies Novel Tumour Microenvironment Profiles, With Prognostic and Therapeutic Implications

Comparative Molecular Analysis of Gastrointestinal Adenocarcinomas

Expression Patterns and Prognostic Value of M(6)a-Related Genes in Colorectal Cancer

Cancer-Cell-Intrinsic Mechanisms Shaping the Tumor Immune Landscape

The Immune Landscape of

Immunogenic Cell Death in Colon Cancer Prevention and Therapy

Pattern of Programmed Cell Death-Ligand 1 Expression and CD8-Positive T-Cell Infiltration Before and After Chemoradiotherapy in Rectal Cancer

Clinical Impact of the Immune Microenvironment in Spinal Chordoma: Immunoscore as an Independent Favorable Prognostic Factor

Intratumoral CD8 + Cytotoxic Lymphocyte Is a Favorable Prognostic Marker in Node-Negative Breast Cancer

The Interplay Between Neutrophils and CD8 T Cells Improves Survival in Human Colorectal Cancer

The Evolving Role of CD8+CD28-Immunosenescent T Cells in Cancer Immunology

Tissue Resident CD8 Memory T Cell Responses in Cancer and Autoimmunity

Telomerase and CD4 T Cell Immunity in Cancer

CXCL10, CXCL11/CXCR3 Axis for Immune Activation -A Target for Novel Cancer Therapy

The Role of Dendritic Cells in Cancer

Dendritic Cells and Cancer Immunity

The Prognostic Role of Macrophage Polarization in the Colorectal Cancer Microenvironment

CD16 Expression on Neutrophils Predicts Treatment Efficacy of Capecitabine in Colorectal Cancer Patients

The Role of Sex in the Innate and Adaptive Immune Environment of Metastatic Colorectal Cancer

Tgfbeta1-Mediated SMAD3 Enhances PD-1 Expression on Antigen-Specific T Cells in Cancer

Dectin-1-Activated Dendritic Cells Trigger Potent Antitumour Immunity Through the Induction of Th9 Cells

Triggering of TNFRSF25 Promotes CD8+ T-Cell Responses and Anti-Tumor Immunity

Survival Advantages Conferred to Colon Cancer Cells by E-Selectin-Induced Activation of the PI3K-Nfkb Survival Axis Downstream of Death Receptor-3

CCL9 Induced by Tgfbeta Signaling in Myeloid Cells Enhances Tumor Cell Survival in the Premetastatic Organ

CXCL13-Producing TFH Cells Link Immune Suppression and Adaptive Memory in Human Breast Cancer

CXCL13 and its Receptor CXCR5 in Cancer: Inflammation, Immune Response, and Beyond

Recent Advances Targeting CCR5 for Cancer and Its Role in Immuno-Oncology

CCR10 Activation Stimulates the Invasion and Migration of Breast Cancer Cells Through the ERK1/2/MMP-7 Signaling Pathway

Cell and B Cell-Related Pathways for Novel Cancer Treatments

Heterogeneity of HLA-G Expression in Cancers: Facing the Challenges

Downregulation of Antigen Presentation-Associated Pathway Proteins Is Linked to Poor Outcome in Triple-Negative Breast Cancer Patient Tumors

Genetic Mechanisms of Immune Evasion in Colorectal Cancer

The Multifaceted Role of Chromosomal Instability in Cancer and Its Microenvironment

Structural Basis of N (6)-Adenosine Methylation by the METTL3-METTL14 Complex

METTL14 Suppresses Proliferation and Metastasis of Colorectal Cancer by Down-Regulating Oncogenic Long Non-Coding RNA XIST

METTL14-Mediated N6-Methyladenosine Modification of SOX4 Mrna Inhibits Tumor Metastasis in Colorectal Cancer

METTL14 Suppresses CRC Progression Via Regulating N6-Methyladenosine-Dependent Primary Mir-375 Processing

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2021. 615296/full#supplementary-material Conflict of Interest: The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.Copyright © 2021 Cai, Long, Huang, Han, Peng, Guo, Liu, Chen, Shen, Long, Wang, Yu, Shen and Zeng. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.