Exploration of Novel Potential Biomarkers for Cancer ImmunotherapyView this Special Issue
Characterization of m6A Methylation Modification Patterns in Colorectal Cancer Determines Prognosis and Tumor Microenvironment Infiltration
Cumulative studies have suggested that dysregulation of m6A regulators and immunity is highly linked to the prognosis of patients with cancer. However, the potential contribution of m6A modification patterns to the tumor microenvironment (TME) and the therapeutic efficacy of immunotherapy for colorectal cancer (CRC) remain elusive. A comprehensive analysis of the m6A modification profiles of 458 patients with CRC was performed by clustering 21 genes encoding m6A methylation regulators and linking the m6A modification pattern with TME characteristics. Using principal component analysis (PCA), a risk model was constructed to quantify individual m6A modification patterns in patients with CRC. The results indicated that the expression profiles and genetic mutations of 21 genes encoding m6A methylation regulators in CRC were characterized by a high degree of heterogeneity. Three m6A clusters had significant differences in prognosis, m6A modification patterns, and TME characteristics. Furthermore, a risk model, termed m6Ascore, was developed by PCA to quality m6A methylation patterns at an individual level. The m6Ascore could stratify patients into high- and low-m6Ascore groups. Further analyses demonstrated that the m6Ascore had a good predictive performance for overall survival and clinical efficacy of immunotherapy in patients with CRC. Finally, the predictive value of the model was validated by external cohorts. In conclusion, the comprehensive characterization of m6A methylation modification patterns might contribute to our understanding of the TME in CRC and the development of personalized antitumor immunotherapy in the future.
Colorectal carcinoma (CRC) is the third most commonly diagnosed cancer and is the second leading cause of cancer-related deaths globally . Although current treatments such as targeted therapy, immunotherapy, and precision treatment have been applied for the treatment of CRC [2, 3], the clinical outcomes are unsatisfactory, and the prognosis of patients with CRC remains bleak. The prognosis of patients with CRC varies widely even among patients with the same stage and therapeutic regimen, largely attributed to the highly heterogeneous nature of CRC . In the current era, immunotherapy using immune checkpoint inhibitors (ICIs) is an encouraging strategy for the treatment of solid tumors, including CRC [5–10]. The heterogeneity of CRC at the genetic and molecular levels has significant ramifications in terms of immunotherapy outcomes . In CRC, microsatellite instability , tumor mutational burden (TMB) , and DNA polymerase epsilon mutations  have emerged as major markers to predict the responses to immunotherapy. However, these indicators do not provide accurate predictions for the current ICIs, thereby facilitating the development of more accurate and reliable biomarkers.
As a common chemical modification of eukaryotic messenger ribonucleic acids (mRNAs), N6-methyladenosine (m6A) can affect various essential biological processes by regulating the expression of target genes [15, 16]. m6A-regulated proteins consist of “writers” (WTAP, METTL3, and METTL14), “erasers” (ALKBH5 and FTO), and “readers” (IGF2BPs and YTHDFs) [17–19]. Currently, numerous studies have offered evidence that m6A regulators perform crucial functions in modulating the maturation, translation, and degradation of RNAs (including mRNAs and noncoding RNAs). Cumulative studies have revealed that m6A regulator dysregulation is correlated with apoptosis, proliferation, self-renewal, developmental defects, malignant tumor progression, and tumor microenvironment (TME) [20–22]. TME refers to the environment surrounding tumor cells, including perivascular cells, immune cells, fibroblasts, molecules, extracellular matrix, and additional stromal components. Currently, in CRC, several reports have suggested that m6A methylation regulators play a crucial role in the tumor immune microenvironment (TIME), which is defined as the immune and the immune-related constituents of the TME. For instance, Cai et al. highlighted that downregulation of the m6A “writer” METTL14 represents an unfavorable patient prognosis and that its low expression level may result in the downregulation of m6A modifications and reduction of the level of immune cell infiltration . Approximately 85% of patients with CRC exhibit mismatch-repair-proficient or microsatellite instability-low (pMMR-MSI-L) tumors, and such patients are poorly responsive to immunotherapy. Wang et al. reported that depletion of the methylation transferases METTL3 and METTL14 could suppress m6A modification and enhance responses of pMMR-MSI-L CRC to anti-PD-1 therapy. Furthermore, in vivo experiments suggested that METTL3- or METTL14-deficient tumors were associated with increased abundance of cytotoxic tumor-infiltrating CD8+ T cells and significantly enhanced the secretion of IFN-γ, CXCL9, and CXCL10 in TME . According to Tsuruta et al., the m6A “erasers” gene-FTO could regulate PD-L1 expression in colon cancer cells in an IFN-γ signaling-independent manner . However, to date, most studies have only focused on the role of single m6A regulators in the regulation of the TIME and have failed to provide a comprehensive perspective on how m6A modification patterns contribute to the TIME.
In this study, we comprehensively investigated m6A modification characteristics and identified three m6A modification patterns (m6Aclusters) with distinct survival benefits, TME immune cell infiltrations, and transcriptome characteristics. Moreover, we constructed an m6A regulator-associated gene signature, m6Ascore, which serves as a risk model for assessing individual m6A modification patterns. Therefore, the m6Ascore is a promising biomarker for predicting prognosis and clinical response to ICIs in patients with CRC.
2. Materials and Methods
2.1. Date Collection
The whole workflow of our present study is displayed in Supplementary Figure S1. RNA sequencing (RNA-seq) data for gene expression analysis, genetic mutations, and clinical information of patients with colon carcinoma () and rectal carcinoma (n =127) were downloaded from The Cancer Genome Atlas (TCGA, http://cancergenome.nih.gov/). Gene expression data and its corresponding clinical information were obtained from the GSE17536 () and GSE78220 () datasets in the Gene Expression Omnibus (GEO) database. The Broad GDAC Firehose (http://gdac.broadinstitute.org/) was retrieved to obtain information on copy number variations (CNVs). In addition, the IMvigor210 cohort, a dataset of advanced urothelial cancer treated with PD-L1 inhibitor (atezolizumab), was collected using the “IMvigor210” R package . GSE78220 (), a dataset of metastatic melanoma treated with anti-PD-1 antibody (pembrolizumab), was collected from the GEO database . The above two immunotherapeutic datasets were collected to assess our findings.
2.2. Analysis of Somatic Mutations
Gene somatic mutation data (MAF files) containing the mutation profile of the CRC cohort were obtained from TCGA database. The visualization procedure was performed with the aid of the “maftools” R package .
2.3. Correlation Analysis between m6A Regulators and Immune Infiltrations
The “CIBERSORT” package in R was employed to assess the infiltration levels of 22 distinct immune cell types across m6A regulators and m6A clusters .
2.4. Unsupervised Clustering for m6A Regulators
A total of 21 m6A modulators were obtained from previous studies  (Supplementary Table S1) and analyzed to identify the characteristics of m6A modification. Unsupervised clustering analysis was performed using the “ConsensusClusterPlus” R package , and tumors were sampled for 1000 bootstrap replications to ensure stability. The cluster number was calculated by the area under the cumulative distribution of the function curve, and the -means method was used with the Euclidean distance metric. Distinct m6A modification profiles were identified according to the expression values of 21 m6A regulators, and patients were classified for subsequent analysis.
2.5. Identification of Differentially Expressed Genes Corresponding to m6A Clusters
The “limma” R package was applied to identify differentially expressed genes (DEGs) across various m6A clusters . DEGs were considered significantly different between groups in the case where and adjusted value < 0.05.
2.6. Gene Set Enrichment Analysis and Gene Enrichment Functional Annotation
Gene set enrichment analysis (GSEA), an approach for interpreting gene expression data, was conducted to thoroughly comprehend the distinct biological functions underlying m6A modification patterns . The Molecular Signature Database (MSigDB) for GSEA (https://www.gsea, http://msigdb.org/gsea/index.jsp) was used to acquire the gene set “h.all.v7.2.symbols.gmt.” The enrichment as well as modeled enrichment scores were calculated using the “GSEA” R package. Values with adjusted were considered statistically significant. The clusterProfiler package in R was used to conduct Gene Ontology (GO) analysis of m6A-associated DEGs. Values with were considered statistically significant.
2.7. Correlation Analysis between m6A Clusters and Immunomodulators
A total of 75 immunomodulators (IMs) were extracted from previous studies . One-way analysis of variance (ANOVA) and Kruskal–Wallis tests were used to investigate the differences in IM potency among different clusters.
2.8. Dimension Reduction and Construction of the m6Ascore Model
The m6Ascore model was constructed based on m6A cluster-associated DEGs. First, a total of 194 m6A cluster-associated DEGs were selected for subsequent analysis. Unsupervised clustering was performed to define DEG clusters. Next, principal component analysis (PCA) was performed for dimension reduction, and the final weight values of PC1 and PC2 were used to construct the m6Ascore model with the formula as follows : , where “” represents the expression of m6A regulator cluster-associated DEGs.
2.9. Statistical Analysis
The correlation coefficient was evaluated by Spearman’s correlation analysis. The optimal cutoff value for each cohort was determined by the “survminer” R package. Statistical significance of comparisons between two groups was computed by Student’s -test or the Wilcoxon signed-rank test. The Kruskal-Wallis test and one-way ANOVA were used for the purpose of examining statistical differences among the three groups. The Kaplan-Meier method was employed to perform prognostic analysis, whereas the log-rank test was used to assess if there were any differences among the groups. A protein-protein m6A regulator interaction (PPI) network was constructed using the STRING database . All statistical values were double-sided. Statistical significance was set at . All statistical analyses were performed using R 3.6.1 (https://www.r-project.org/).
3.1. Characterization and Clinical Value of m6A Regulators in CRC
In the present study, 21 m6A-related genes including “writers” (ZC3H13, VIRMA, WTAP, RBM15, RBM15B, METTL3, METTL14, and METTL16), “erasers” (ALKBH5 and FTO), and “readers” (IGF2BP1, IGF2BP2, IGF2BP3, RBMX, HNRNPC, HNRNPA2B1, YTHDC1, YTHDC2, YTHDF1, YTHDF2, and YTHDF3) were identified from previous studies . In TCGA-CRC dataset, the overall mutation frequency of the above 21 RNA methylation-related genes was high, with 143 of 163 (87.73%) samples having mutations in 21 genes encoding m6A regulators. The genes undergoing the highest mutation frequencies were ZC3H13 (26%), KIAA1429 (17%), and YTHDC2 (17%) (Figure 1(a)). Then, we investigated the CNV of 21 genes encoding m6A regulators and found that CNVs were frequent in TCGA-CRC dataset. The CNV of ALKBH5, HNRNPA2B1, IGF2BP3, KIAA1429, METTL16, YTHDF1, YTHDF3, and ZC3H13 was more than 50%. Among them, ALKBH5 and METTL16 primarily exhibited deletions, and the other genes mostly exhibited amplification (Figure 1(b)). At the transcriptomic level, we compared the differential expression of genes encoding 21 m6A regulators in tumor and healthy samples, and the results showed that 13 out of 21 genes were significantly differentially expressed. The genes IGF2BP1, IGF2BP3, and YTHDF1 were significantly upregulated in tumor samples, whereas the genes YTHDC1, YTHDC2, ALKBH5, METTL14, METTL16, FTO, WTAP, YTHDF2, YTHDF3, and RBM15B were significantly downregulated in tumor tissues (Figure 1(c)). Furthermore, the PPI network constructed using STRING implies extensive protein interactions among these regulators. The interactions and prognostic relevance of the 21 m6A regulators in patients with CRC are illustrated by the regulatory network shown in Figure 1(d), and the results suggested a positive coexpression profile in most genes. Taken together, the above results indicate that the genetic and expression differences of m6A regulators are highly heterogeneous in CRC, suggesting that differences in the expression of m6A regulators might have a key role in the onset and progression of CRC.
3.2. Prognosis and Immune Landscape of m6A Regulators
First, the coexpression and prognostic significance of the 21 m6A regulators in TCGA-CRC cohort were analyzed and further visualized in the network (Figure 2(a)). The results revealed that most of the m6A regulators were significantly positively correlated.
Recent evidence has demonstrated that m6A modifications are involved in tumor immunity . For the purpose of further analyzing the association between m6A modifications and the TME, we investigated the association of 21 m6A gene expression profiles with the infiltration of 22 distinct immune cell types. An extensive correlation was found between the expression of m6A regulators and the infiltration of most immune cell subtypes, with significant differences between the expression of different genes and immune cell infiltration (Figure 2(b)). The m6A eraser YTHDC2 was determined as a tumor suppressor in the above analysis. Results of GSEA demonstrated a predominant enrichment of ADHERENS JUNCTION, ECM RECEPTOR INTERACTION, ERBB SIGNALING PATHWAY, PROSTATE CANCER, and SMALL CELL LUNG CANCER in the YTHDC2 high-expression group and a predominant enrichment of OLFACTORY TRANSDUCTION in the YTHDC2 low-expression group (Figure 2(c)). Results of survival analysis indicated a significantly longer overall survival (OS) duration in patients belonging to the YTHDC2 high expression group compared to those belonging to the YTHDC2 low expression group in both TCGA-CRC and IMvigor210 cohorts (Figures 2(d) and 2(e)). The anatomical location of the primary tumor is also recognized as an important factor for colon cancer prognosis . As illustrated in Figure 2(f), no significant differences were recorded for the expression of the YTHDC2 gene in primary tumors with different anatomical locations. TMB is an important indicator to assess the efficacy of immunotherapy for CRC . The expression of the YTHDC2 gene was significantly higher in the high-TMB group than in the low-TMB group (Figure 2(g)). Given the crucial role played by YTHDC2 in tumor immunity, we thoroughly investigated if YTHDC2 expression could predict immunotherapeutic responses. Unfortunately, in the IMvigor210 dataset, no significant difference was discovered in terms of the expression of YTHDC2 in patients with therapeutic response compared to patients without therapeutic response (Figure 2(h)).
Taken together, these results disclose a significant correlation between the expression of m6A methylation regulators and the TME. The expression of the m6A eraser YTHDC2 is closely associated with immune cell infiltration and might be a favorable biomarker for prognosis. However, the prediction significance of YTHDC2 in terms of immunotherapeutic efficacy still needs further validation.
3.3. Identification of m6A Modification Characteristics Based on 21 m6A Regulators
m6A regulators might be responsible for the CRC heterogeneity, and they are tightly correlated with TME. To further identify novel m6A regulator profiles, we performed unsupervised clustering based on the expression of 21 m6A regulators in TCGA-CRC data. As demonstrated in Supplementary Figure S2, three clusters showed the best clustering performance, and patients were stratified into m6Acluster-1 (), cluster-2 (), and cluster-3 (). Survival analysis revealed that the three m6A clusters had significant differences in the OS. The OS duration of patients in m6Acluster-1 and m6Acluster-2 was significantly longer than that of patients in m6Acluster-3 (Figure 3(a)). The PCA algorithm was employed to visualize the expression profile of 21 m6A regulators, and the results indicated that three clusters could be well distinguished (Figure 3(b)). Therefore, the classification of the m6A clusters was reasonable.
Subsequently, we investigated the differences in immune cell infiltrations among m6Aclusters and found that the infiltration levels of CD8 T cells, activated NK cells, M1 macrophages, M2 macrophages, and neutrophils were significantly higher in m6Acluster-1 and m6Acluster-3 than in the m6A-2 cluster. However, the infiltration levels of plasma cells, M0 macrophages, resting CD4-positive memory T cells, regulatory T cells, and monocytes were higher in m6Acluster-2 than in m6Acluster-1 and m6Acluster-3 (Figure 3(c)). Next, we examined the differences in the role of signaling pathways across the three clusters, and the results indicated that m6Acluster-1 and m6Acluster-3 were mainly enriched in KEGG TERPENOID BACKBONE BIOSYNTHESIS, KEGG ALZHEIMERS DISEASE, KEGG PARKINSONS DISEASE, KEGG OXIDATIVE PHOSPHORYLATION, and KEGG PEROXISOME and that m6Acluster-2 was mainly enriched in KEGG NITROGEN METABOLISM, KEGG PRIMARY BILE ACID BIOSYNTHESIS, and KEGG NEUROACTIVE LIGAND RECEPTOR INTERACTION (Figure 3(d)).
3.4. Clinical and Transcriptomic Characteristics of the Three m6A Clusters
First, we explored the role of m6A regulators in the classification of m6A clusters. The heat map of m6A-associated genes showed that IGF2BP1 and IGF2BP3 played a major role in the classification process (Figure 4(a)). Immune infiltration-associated signatures have a significant role in TME. Therefore, it is necessary to thoroughly probe into the correlation between m6A clusters and various immune infiltration-associated signatures as established by Mariathasan et al. . We first investigated differences in the expression of these immune infiltration-associated signatures in m6A clusters, and the results indicated that m6Acluster-1 and m6Acluster-3 were enriched in immune activation signatures such as CD8 T effector, antigen processing machinery, immune checkpoint, cytolytic activity, type I IFN response, coinhibition T cell, coinhibitory APC, costimulation T cell, MHC-I HLA, and MHC-II HLA. Moreover, m6Acluster-3 was enriched in stroma activation phenotype signatures, including EMT1, EMT3, and pan-F-TBRS (Figure 4(b)). Furthermore, we performed differential expression analysis of IMs among the three clusters and found that multiple genes showed significant differential expression (Figure 4(c)). A comprehensive analysis of the above results and the characteristics of immune cell infiltrations and patient prognosis in the three subtypes revealed that m6Acluster-1 is an immune-inflamed phenotype associated with elevated infiltration levels of immune cells and that m6Acluster-3 is an immune-excluded phenotype associated with the infiltration levels of multiple immune cell types. However, these cells are reserved in the surrounding nests of tumor cells instead of penetrating their stroma. Furthermore, m6Acluster-2 was an immune-desert phenotype characterized by a few immune cell infiltrates and a suppressive immunological landscape .
To further uncover the potential biological profile of distinct m6A patterns, DEGs among m6A clusters were analyzed. A total of 194 genes were identified, out of which 167 genes were significantly upregulated and 27 were significantly downregulated (Figure 4(d)). GO analysis was performed on the upregulated and downregulated DEGs, and the top 10 pathways enriched in the three functional categories (biological process (BP), cellular component (CC), and molecular function (MF)) are shown in bubble plots. As shown in Figures 4(e) and 4(f), the majority of the enriched pathways are linked to BPs such as differentiation of immune cells, ligand and receptor activity of growth factors and cell membranes, and stromal-related and vesicle formation. The above results indicate the different clinical and transcriptomic features of the m6A clusters.
3.5. Construction of the m6Ascore Model
Unsupervised clustering analysis was performed using the obtained expression profiles of 194 tumor m6A cluster-associated DEGs. As a result, the tumor samples of TCGA-CRC could be categorized into three m6Acluster-associated differentially expressed gene clusters (DEG-cluster), named Gene cluster-1, Gene cluster-2, and Gene cluster-3. Results of the prognostic analysis indicated significant survival differences among DEG clusters (Supplementary Figure S3). The above analysis was performed based on m6A methylation modifications, reflecting accurately the expression profiles of m6A regulators in CRC. Furthermore, based on the DEGs among m6A clusters, we could construct a risk score system by PCA algorithm, referred to as the m6Ascore, for the purpose of quantifying m6A modification profiles in individual patients with CRC. The optimal density threshold of the m6Ascore associated with survival was calculated using the “survminer” R package, with a threshold value of 3.23. This threshold value was used to classify TCGA-CRC samples into high- and low-m6Ascore groups (Figure 5(a)). The findings from the Kaplan-Meier survival analysis demonstrated an improved prognosis condition among patients in the low-m6Ascore group compared to those in the high-m6Ascore group (Figure 5(b)).
Then, we sought to explore the correlation between the DEG clusters and clinical characteristics. As shown in Figure 6(a), the patients in Gene cluster-2 and Gene cluster-3 were younger than 65 years, had a better prognosis, and were diagnosed with left-sided colon and rectum cancer. (Figure 6(b)). Additionally, the differential expression of IMs among three clusters was analyzed, and the results indicated that multiple genes (including KIR2DL1, KIR2DL3, GZMA, IFNA2, IL12A, TNFSF9, VEGFA, VEGFB, CD274, VTCN1, and SELP) showed significant differential expression (Figure 6(c)). Then, we compared the m6Ascore in left- and right-sided samples and the results suggested much greater scores among the right-sided patients than those among the left-sided patients (Figure 6(d)). TMB is an important indicator to assess the efficacy of immunotherapy for CRC . The results demonstrated a remarkably elevated m6Ascore among patients belonging to the high-TMB group compared to those belonging to the low-TMB group (Figure 6(e)). For m6A methylation DEG clusters, DEG cluster-1 demonstrated the highest median m6Ascore, while DEG cluster-2 demonstrated the lowest median m6Ascore (Figure 6(f)). These findings might provide novel insights into the mechanisms of m6A modifications and gene mutations among immune checkpoints in CRC. The association among m6A methylation regulators, DEG clusters, anatomical sites, m6Ascore, and OS is presented as a Sankey diagram (Figure 6(g)).
3.6. Validation of the m6Ascore Model in GEO Datasets
The robustness of the m6A risk model for predicting CRC prognosis was further validated on the GSE17536 dataset. The m6Ascore for the GSE17536 cohort was calculated using the previously screened DEGs. Next, the optimal density threshold of the m6Ascore associated with survival was calculated using the “survminer” R package, and the samples of GSE17536 were divided into high and low-m6Ascore groups. Results of the Kaplan-Meier analysis also revealed that the patients in the high-m6Ascore group had a poorer prognosis than those in the low-m6Ascore group (Figure 7(a), ). Moreover, as shown in the heat map (Figure 7(b)), the m6Ascore could serve to assess several clinical features. The group with a high m6Ascore was associated with favorable histological subtypes and earlier staging in the GSE17536 cohort. All these findings showed that the m6Ascore model might be used to predict the prognosis and certain clinical parameters of patients with CRC and to provide potential therapeutic implications in the clinical setting.
3.7. Assessment of the Predictive Power of m6Ascore for ICI Immunotherapy Response
The immunophenogram was developed for predicting checkpoint blockade responses at a pancancer level by Charoentong et al. . To determine the predictive power of the m6Ascore in terms of immunotherapeutic benefit, we obtained the dataset of TCGA-CRC cohort from The Cancer Immunome Atlas (TCIA) database and explored the correlation between immunophenoscore (IPS) and m6Ascore by performing immunophenogram analysis. The results illustrated that CTLA4-negative+PD-1-negative and CTLA4-positive+PD-1-positive cells exhibited significant differences in the immunotherapeutic response between the high- and low-m6Ascore groups (), whereas CTLA4-negative+PD-1-positive and CTLA4-positive+PD-1-negative cells exhibited no differences in the immunotherapeutic response between the high- and low-m6Ascore groups (). Moreover, the high-m6Ascore group tended to have a higher IPS score than the low-m6Ascore group, suggesting that patients in the high-m6Ascore group were more likely to gain benefit from immunotherapy (Figures 8(a)–8(d)).
Next, we investigated whether m6Ascore could predict the immunotherapeutic responses using the GSE78220 and IMvigor210 cohorts. Patients in the IMvigor210 cohort receiving anti-PD-L1 immunotherapy were classified into two groups according to the m6Ascore, namely, the high m6Ascore and the low m6Ascore. The findings demonstrated a significantly higher immunotherapeutic response rate in the high-m6Ascore group than that in the low-m6Ascore group (Figures 8(e) and 8(g)). Moreover, survival analysis demonstrated a significantly longer OS duration among patients in the low-m6Ascore group than that in the high-m6Ascore group (Figure 8(f)). Similarly, the same trend was found in the GSE78220 cohort, in which the patients were treated with anti-PD-1 immunotherapy (Figures 8(h)–8(j)). All these findings confirmed that m6Ascore might serve to predict the efficacy of immunotherapies and the prognosis of patients with CRCs.
Recent studies have linked m6A modifications to the regulation of tumor immunity and response to ICIs . To date, accumulating reports have shown the crucial roles of m6A regulators in carcinogenesis, development, diagnosis, therapy, and prognosis [21, 40]. However, the in-depth analysis of the m6A methylation regulator pattern and its significance on the TIME of CRC needs further elucidation. In this study, we shed light on the association between different m6A modification profiles and the immune landscape in CRC. This study improves our understanding of the TME and facilitates the development of more accurate immune therapy strategies for the treatment of patients with CRC.
In the present study, we first globally assessed somatic mutations, CNVs, and RNA expression of the 21 m6A regulators in TCGA-CRC cohort patients. We found that mutations and CNVs of the 21 m6A regulators were frequent in TCGA-CRC patients. At the transcriptomic level, we found that 2 m6A regulator expression was higher and 8 regulator expression was lower in TCGA-CRC patients. In addition, the correlation between each m6A regulators and TME-infiltrating cells was explored. Our results indicated that significant differences were observed between genes and immune cell infiltrations. Collectively, the results suggested that 21 genes encoding m6A regulators are highly heterogeneous in CRC patients and dysregulation of m6A regulators was associated with the development and antitumor immune response of CRC.
In addition, here, we focused on the regulatory role of the m6A reader YTHDC2 in colorectal cancer. To the best of our knowledge, the m6A readers have been investigated in other cancers [41–43], but YTHDC2 has been poorly studied in CRC. In the present study, the outcomes indicated that YTHDC2 was lowly expressed in colorectal cancer and coupled with prognostic analysis hinted that it might be a tumor suppressor, and this result was in line with Ji et al.’s results . Further, YTHDC2 was significantly correlated with multiple immune cell infiltrations, and the expression of YTHDC2 gene in the high-TMB group was found to be significantly higher than that in the low-TMB group; hence, we inferred that it might be correlated with the efficacy of immunotherapy; unfortunately, it was not found after the validation using IMvigor210 dataset; the reason might be that this external dataset was not from CRC. Thus, in the future, CRC datasets need to be validated and basic experiments need to be performed to further investigate the mechanism of YTHDC2 in the regulation of immune cell infiltration.
Considering the highly individual heterogeneity of m6A regulator modification, then we uncovered the characteristics of m6A methylation modification mediated by 21 genes encoding m6A regulators in CRC. Based on 21 m6A regulators, we identified three distinct models of m6A modification by unsupervised consensus clustering in the CRC cohort. These three patterns termed m6Acluster-1, m6Acluster-2, and m6Acluster-3 are characterized by different prognoses and immune cell infiltrations in the TME. m6Acluster-1 presented a high proportion of effector TME immune cell activation and infiltration level, which are the hallmarks of the immune-inflamed phenotype, also described as a “hot” immune microenvironment . Although m6Acluster-3 was also characterized by both elevated levels of TME immune cell infiltration, it had a high level of EMT1, EMT3, and pan-F-TBRS expression, which could determine the activation of stroma-associated pathways. Moreover, the patients exhibiting the m6Acluster-3 phenotype had worse survival compared to those exhibiting the m6Acluster-1 phenotype. The location and migration of T cells as well as the status of the stroma play crucial roles for antitumor immunosurveillance. In the case of high levels of immune cell infiltration, m6Acluster-3 could be classified as an immune rejection phenotype since T cells are restricted to the stroma, thus failing to penetrate the core tumor nest . In contrast to m6Acluster-3 and m6Acluster-1, m6Acluster-2 was characterized by the lack of immune cell infiltration and immunosuppression, which could be labeled as an immune-desert phenotype. Previous research has suggested that this phenotype was probably correlated with the absence of activated and primed T cells and immune tolerance . The immune rejection phenotype and immune-desert phenotype are also referred to as a “cold” immune microenvironment . Charoentong et al. illustrated that the expression levels of IMs were closely associated with tumor genotypes and that tumor genotypes could determine tumor immunophenotypes and tumor escape mechanisms . We found significant differences in the expression of IMs among the three m6A clusters. Collectively, the above findings confirmed the reliability of m6A modification patterns for identifying and classifying immune phenotypes of CRC.
Next, DEGs among the three m6A clusters were identified as m6A modification pattern-associated genes and were found to be enriched mainly in immune-related and stromal-related biological processes. Similar to the results of m6A patterns, three DEG clusters were recognized by unsupervised consensus clustering. These DEG clusters featured significantly different survival outcomes and IM expression. Moreover, we constructed a risk model, termed m6Ascore, for the purpose of quantifying the pattern of m6A modifications in patients with CRCs. Thus, the m6Ascore risk model could divide TCGA-CRC cohort into two groups, and survival analysis indicated that the patients belonging to the low-m6Ascore group had significantly longer OS duration than those belonging to the high-m6Ascore group. We employed the GSE17536 cohort to verify the reliability and robustness of the risk model in terms of OS prediction.
We next investigated whether m6Ascore is related to the TMB. The key initiator of adaptive immune activation is neoantigen recognition, and mutational profiles are suggested to be promising biomarkers to predict clinical response to ICI-based immunotherapy. However, it is difficult to detect the expression of overall neoantigens. TMB is easily available and commonly employed to estimate neoangiogenic load and is a good marker for predicting immunotherapeutic responses [49, 50]. The TMB was significantly higher among the patients in the high-m6Ascore group than those in the low-m6Ascore group. Subsequently, in both the IMvigor210 and GSE78220 cohorts, patients in the high-m6Ascore group demonstrated better responses to ICI immunotherapy in contrast with those in the low-m6Ascore group, which was consistent with the reports that patients with a high-TMB had a superior clinical response to immunotherapy. Furthermore, the m6Ascore can also be useful for the evaluation of clinical characteristics of patients, such as anatomical location, histological type, and tumor stages. Taken together, the above findings imply that m6Ascore is a reliable and robust biomarker for predicting the efficacy of ICI immunotherapies and is a promising tool to promote individualized immunotherapy for patients with CRC in the future.
There are a few limitations to the present study. First, the findings of this study are based on bioinformatic predictions; the m6A modification patterns and m6Ascore risk model were only established according to the 21 m6A regulator expression in TCGA-CRC cohort. Hence, further large-scale, multicenter, and prospective clinical cohorts are needed to verify our findings. Second, due to the lack of survival data of patients with CRC undergoing immunotherapy, the predictive reliability of m6Ascore in CRC needs further validation. Third, this study is a bioinformatics analysis and can be used as a preliminary reference. Further deep basic experimental researches are required to discover the association between m6A methylation and immune cell regulations, as well as the efficacy of immunotherapy in CRC cohort.
This study categorized CRC into three stable immune subtypes with distinct prognoses and TME based on 21 m6A methylation regulators. Importantly, a risk model-m6Ascore was developed to predict prognosis, clinical features, and the efficacy of immunotherapy in patients with CRC. The predictive power of the m6Ascore was validated on external datasets. Therefore, m6Ascore may be a viable tool to predict the immunotherapy response in patients with CRC, thus paving the way for personalized immunotherapy in the future.
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare that they have no competing interests.
XD designed the study. QY and YZ are the principal investigators and edited the manuscript. FW and JB conducted statistical analysis and data management. FC and FW revised the manuscript. All authors read and approved the final manuscript.
This study was supported by Shaanxi Provincial People’s Hospital Talent Support Program (No. 2021JY-44), Innovation Capability Support Program of Shaanxi (Nos. 2019GHJD-14 and 2021TD-40), and Key Research and Development Program of Shaanxi (No. 2020GXLH-Y-019)..
Figure S1: the workflow of our study design and analytical pipeline. Figure S2: consensus clustering of the expression profiles of 21 m6A methylation regulators in colorectal cancer (CRC). (A–D) Clustering results at classification numbers , 3, 4, and 5, respectively. (E) Distribution of CDF curves for consensus clustering. (F–I) Survival curves at classification numbers , 3, 4, and 5, respectively. (J) Distribution of area under cumulative distribution frequency (CDF) curves for consensus clustering. Figure S3: consensus clustering of differentially expressed genes (DEGs) among tumor m6A clusters. (A–D) Clustering results at classification numbers , 3, 4, and 5, respectively. (E) Distribution of cumulative distribution frequency (CDF) curves for consensus clustering. (F–I) Survival curves at classification numbers , 3, 4, and 5, respectively. (J) Distribution of area under CDF curves for consensus clustering. (Supplementary Materials)
R. L. Siegel, K. D. Miller, and A. Jemal, “Cancer statistics, 2019,” CA: a Cancer Journal for Clinicians, vol. 69, no. 1, pp. 7–34, 2019.View at: Publisher Site | Google Scholar
J. Sandhu, V. Lavingia, and M. Fakih, “Systemic treatment for metastatic colorectal cancer in the era of precision medicine,” Journal of Surgical Oncology, vol. 119, no. 5, pp. 564–582, 2019.View at: Publisher Site | Google Scholar
K. E. Wong, S. C. Ngai, K. G. Chan, L. H. Lee, B. H. Goh, and L. H. Chuah, “Curcumin nanoformulations for colorectal cancer: a review,” Frontiers in Pharmacology, vol. 10, p. 152, 2019.View at: Publisher Site | Google Scholar
C. J. Punt, M. Koopman, and L. Vermeulen, “From tumour heterogeneity to advances in precision treatment of colorectal cancer,” Nature Reviews. Clinical Oncology, vol. 14, no. 4, pp. 235–246, 2017.View at: Publisher Site | Google Scholar
R. L. Ferris, G. Blumenschein Jr., J. Fayette et al., “Nivolumab for recurrent squamous-cell carcinoma of the head and neck,” The New England Journal of Medicine, vol. 375, no. 19, pp. 1856–1867, 2016.View at: Publisher Site | Google Scholar
W. Hugo, J. M. Zaretsky, L. Sun et al., “Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma,” Cell, vol. 165, no. 1, pp. 35–44, 2016.View at: Publisher Site | Google Scholar
K. Shitara, M. Özgüroğlu, Y. J. Bang et al., “Pembrolizumab versus paclitaxel for previously treated, advanced gastric or gastro-oesophageal junction cancer (KEYNOTE-061): a randomised, open-label, controlled, phase 3 trial,” Lancet, vol. 392, no. 10142, pp. 123–133, 2018.View at: Publisher Site | Google Scholar
J. Brahmer, K. L. Reckamp, P. Baas et al., “Nivolumab versus docetaxel in advanced squamous-cell non-small-cell lung cancer,” The New England Journal of Medicine, vol. 373, no. 2, pp. 123–135, 2015.View at: Publisher Site | Google Scholar
T. André, K. K. Shiu, T. W. Kim et al., “Pembrolizumab in microsatellite-instability-high advanced colorectal cancer,” The New England Journal of Medicine, vol. 383, no. 23, pp. 2207–2218, 2020.View at: Publisher Site | Google Scholar
H. J. Lenz, E. Van Cutsem, M. Luisa Limon et al., “First-line nivolumab plus low-dose ipilimumab for microsatellite instability-high/mismatch repair-deficient metastatic colorectal cancer: the phase II CheckMate 142 study,” Journal of Clinical Oncology, vol. 40, no. 2, pp. 161–170, 2022.View at: Publisher Site | Google Scholar
E. Picard, C. P. Verschoor, G. W. Ma, and G. Pawelec, “Relationships between immune landscapes, genetic subtypes and responses to immunotherapy in colorectal cancer,” Frontiers in Immunology, vol. 11, p. 369, 2020.View at: Publisher Site | Google Scholar
D. Y. Lizardo, C. Kuang, S. Hao, J. Yu, Y. Huang, and L. Zhang, “Immunotherapy efficacy on mismatch repair-deficient colorectal cancer: from bench to bedside,” Biochimica Et Biophysica Acta. Reviews on Cancer, vol. 1874, no. 2, article 188447, 2020.View at: Publisher Site | Google Scholar
D. A. Fabrizio, T. J. George Jr., R. F. Dunne et al., “Beyond microsatellite testing: assessment of tumor mutational burden identifies subsets of colorectal cancer who may respond to immune checkpoint inhibition,” Journal of Gastrointestinal Oncology, vol. 9, no. 4, pp. 610–617, 2018.View at: Publisher Site | Google Scholar
E. Domingo, L. Freeman-Mills, E. Rayner et al., “Somatic POLE proofreading domain mutation, immune response, and prognosis in colorectal cancer: a retrospective, pooled biomarker study,” The Lancet Gastroenterology & Hepatology, vol. 1, no. 3, pp. 207–216, 2016.View at: Publisher Site | Google Scholar
K. Chen, Z. Wei, Q. Zhang et al., “WHISTLE: a high-accuracy map of the human N6-methyladenosine (m6A) epitranscriptome predicted using a machine learning approach,” Nucleic Acids Research, vol. 47, no. 7, article e41, 2019.View at: Publisher Site | Google Scholar
N. Liu and T. Pan, “N6-Methyladenosine-encoded epitranscriptomics,” Nature Structural & Molecular Biology, vol. 23, no. 2, pp. 98–102, 2016.View at: Publisher Site | Google Scholar
G. Jia, Y. Fu, X. Zhao et al., “N6-methyladenosine in nuclear RNA is a major substrate of the obesity-associated FTO,” Nature Chemical Biology, vol. 7, no. 12, pp. 885–887, 2011.View at: Publisher Site | Google Scholar
X. L. Ping, B. F. Sun, L. Wang et al., “Mammalian WTAP is a regulatory subunit of the RNA N6-methyladenosine methyltransferase,” Cell Research, vol. 24, no. 2, pp. 177–189, 2014.View at: Publisher Site | Google Scholar
Y. Niu, X. Zhao, Y. S. Wu, M. M. Li, X. J. Wang, and Y. G. Yang, “N6-Methyl-adenosine (m6A) in RNA: an old modification with a novel epigenetic function,” Genomics, Proteomics & Bioinformatics, vol. 11, no. 1, pp. 8–17, 2013.View at: Publisher Site | Google Scholar
X. Liu, L. Liu, Z. Dong et al., “Expression patterns and prognostic value of m<sup>6</sup>A-related genes in colorectal cancer,” American Journal of Translational Research, vol. 11, no. 7, pp. 3972–3991, 2019.View at: Google Scholar
J. Li, L. Liang, Y. Yang, X. Li, and Y. Ma, “N6-methyladenosine as a biological and clinical determinant in colorectal cancer: progression and future direction,” Theranostics, vol. 11, no. 6, pp. 2581–2593, 2021.View at: Publisher Site | Google Scholar
L. Guo, H. Yang, C. Zhou, Y. Shi, L. Huang, and J. Zhang, “N6-Methyladenosine RNA modification in the tumor immune microenvironment: novel implications for immunotherapy,” Frontiers in Immunology, vol. 12, article 773570, 2021.View at: Publisher Site | Google Scholar
C. Cai, J. Long, Q. Huang et al., “M6A "writer" gene METTL14: a favorable prognostic biomarker and correlated with immune infiltrates in rectal cancer,” Frontiers in Oncology, vol. 11, article 615296, 2021.View at: Publisher Site | Google Scholar
L. Wang, H. Hui, K. Agrawal et al., “m<sup>6</sup> a RNA methyltransferases METTL3/14 regulate immune responses to anti-PD-1 therapy,” The EMBO Journal, vol. 39, no. 20, article e104514, 2020.View at: Publisher Site | Google Scholar
N. Tsuruta, K. Tsuchihashi, H. Ohmura et al., “RNA N6-methyladenosine demethylase FTO regulates PD-L1 expression in colon cancer cells,” Biochemical and Biophysical Research Communications, vol. 530, no. 1, pp. 235–239, 2020.View at: Publisher Site | Google Scholar
S. Mariathasan, S. J. Turley, D. Nickles et al., “TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells,” Nature, vol. 554, no. 7693, pp. 544–548, 2018.View at: Publisher Site | Google Scholar
W. Hugo, J. M. Zaretsky, L. Sun et al., “Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma,” Cell, vol. 168, no. 3, p. 542, 2017.View at: Publisher Site | Google Scholar
A. Mayakonda, D. C. Lin, Y. Assenov, C. Plass, and H. P. Koeffler, “Maftools: efficient and comprehensive analysis of somatic variants in cancer,” Genome Research, vol. 28, no. 11, pp. 1747–1756, 2018.View at: Publisher Site | Google Scholar
A. M. Newman, C. B. Steen, C. L. Liu et al., “Determining cell type abundance and expression from bulk tissues with digital cytometry,” Nature Biotechnology, vol. 37, no. 7, pp. 773–782, 2019.View at: Publisher Site | Google Scholar
Z. Tu, L. Wu, P. Wang et al., “N6-Methylandenosine-related lncRNAs are potential biomarkers for predicting the overall survival of lower-grade glioma patients,” Frontiers in Cell and Development Biology, vol. 8, p. 642, 2020.View at: Publisher Site | Google Scholar
M. D. Wilkerson and D. N. Hayes, “ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking,” Bioinformatics, vol. 26, no. 12, pp. 1572-1573, 2010.View at: Publisher Site | Google Scholar
M. E. Ritchie, B. Phipson, D. Wu et al., “limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic Acids Research, vol. 43, no. 7, article e47, 2015.View at: Publisher Site | Google Scholar
A. Subramanian, P. Tamayo, V. K. Mootha et al., “Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles,” Proceedings of the National Academy of Sciences of the United States of America, vol. 102, no. 43, pp. 15545–15550, 2005.View at: Publisher Site | Google Scholar
V. Thorsson, D. L. Gibbs, S. D. Brown et al., “The immune landscape of cancer,” Immunity, vol. 48, no. 4, pp. 812–830 e14, 2018.View at: Google Scholar
D. Zeng, M. Li, R. Zhou et al., “Tumor microenvironment characterization in gastric cancer identifies prognostic and immunotherapeutically relevant gene signatures,” Cancer Immunology Research, vol. 7, no. 5, pp. 737–750, 2019.View at: Publisher Site | Google Scholar
D. Szklarczyk, A. Franceschini, S. Wyder et al., “STRING v10: protein-protein interaction networks, integrated over the tree of life,” Nucleic Acids Research, vol. 43, no. D1, pp. D447–D452, 2015.View at: Publisher Site | Google Scholar
F. Petrelli, G. Tomasello, K. Borgonovo et al., “Prognostic survival associated with left-sided vs right-sided colon cancer: a systematic review and meta-analysis,” JAMA Oncology, vol. 3, no. 2, pp. 211–219, 2017.View at: Publisher Site | Google Scholar
D. S. Chen and I. Mellman, “Elements of cancer immunity and the cancer-immune set point,” Nature, vol. 541, no. 7637, pp. 321–330, 2017.View at: Publisher Site | Google Scholar
P. Charoentong, F. Finotello, M. Angelova et al., “Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade,” Cell Reports, vol. 18, no. 1, pp. 248–262, 2017.View at: Publisher Site | Google Scholar
T. Sun, R. Wu, and L. Ming, “The role of m6A RNA methylation in cancer,” Biomedicine & Pharmacotherapy, vol. 112, article 108613, 2019.View at: Publisher Site | Google Scholar
T. Li, P. S. Hu, Z. Zuo et al., “METTL3 facilitates tumor progression via an m6A-IGF2BP2-dependent mechanism in colorectal carcinoma,” Molecular Cancer, vol. 18, no. 1, p. 112, 2019.View at: Publisher Site | Google Scholar
T. Liu, Q. Wei, J. Jin et al., “The m6A reader YTHDF1 promotes ovarian cancer progression via augmenting EIF3C translation,” Nucleic Acids Research, vol. 48, no. 7, pp. 3816–3831, 2020.View at: Publisher Site | Google Scholar
L. Ma, T. Chen, X. Zhang et al., “The m6A reader YTHDC2 inhibits lung adenocarcinoma tumorigenesis by suppressing SLC7A11-dependent antioxidant function,” Redox Biology, vol. 38, article 101801, 2021.View at: Publisher Site | Google Scholar
L. Ji, S. Chen, L. Gu, and X. Zhang, “Exploration of potential roles of m6A regulators in colorectal cancer prognosis,” Frontiers in Oncology, vol. 10, p. 768, 2020.View at: Publisher Site | Google Scholar
T. F. Gajewski, “The next hurdle in cancer immunotherapy: overcoming the non-T-cell-inflamed tumor microenvironment,” Seminars in Oncology, vol. 42, no. 4, pp. 663–671, 2015.View at: Publisher Site | Google Scholar
H. Salmon, K. Franciszkiewicz, D. Damotte et al., “Matrix architecture defines the preferential localization and migration of T cells into the stroma of human lung tumors,” The Journal of Clinical Investigation, vol. 122, no. 3, pp. 899–910, 2012.View at: Publisher Site | Google Scholar
J. M. Kim and D. S. Chen, “Immune escape to PD-L1/PD-1 blockade: seven steps to success (or failure),” Annals of Oncology, vol. 27, no. 8, pp. 1492–1504, 2016.View at: Publisher Site | Google Scholar
Q. Duan, H. Zhang, J. Zheng, and L. Zhang, “Turning cold into hot: firing up the tumor microenvironment,” Trends Cancer, vol. 6, no. 7, pp. 605–618, 2020.View at: Publisher Site | Google Scholar
R. M. Samstein, C. H. Lee, A. N. Shoushtari et al., “Tumor mutational load predicts survival after immunotherapy across multiple cancer types,” Nature Genetics, vol. 51, no. 2, pp. 202–206, 2019.View at: Publisher Site | Google Scholar
N. A. Rizvi, M. D. Hellmann, A. Snyder et al., “Cancer immunology. Mutational landscape determines sensitivity to PD-1 blockade in non-small cell lung cancer,” Science, vol. 348, no. 6230, pp. 124–128, 2015.View at: Publisher Site | Google Scholar