Skip to main content
Download PDF

Single-cell N6-methyladenosine regulator patterns guide intercellular communication of tumor microenvironment that contribute to colorectal cancer progression and immunotherapy

Journal of Translational Medicine volume 20, Article number: 197 (2022) Cite this article

Abstract

Background

N6-methyladenosine (m6A) RNA methylation plays a critical role in key genetic events for various cancers; yet, how m6A functions within the tumor microenvironment (TME) remains to be elucidated.

Methods

A total of 65,362 single cells from single-cell RNA-seq data derived from 33 CRC tumor samples were analyzed by nonnegative matrix factorization (NMF) for 23 m6A RNA methylation regulators. CRC and Immunotherapy cohorts from public repository were used to determine the prognosis and immune response of TME clusters.

Results

The fibroblasts, macrophages, T and B cells were respectively grouped into 4 to 5 subclusters and then classified according to various biological processes and different marker genes. Furthermore, it revealed that the m6A RNA methylation regulators might be significantly related to the clinical and biological features of CRC, as well as the pseudotime trajectories of main TME cell types. Bulk-seq analysis suggested that these m6A-mediated TME cell subclusters had significant prognostic value for CRC patients and distinguished immune response for patients who underwent ICB therapy, especially for the CAFs and macrophages. Notably, CellChat analysis revealed that RNA m6A methylation-associated cell subtypes of TME cells manifested diverse and extensive interaction with tumor epithelial cells. Further analysis showed that ligand-receptor pairs, including MIF −  (CD74 + CXCR4), MIF −  (CD74 + CD44), MDK–NCL and LGALS9 − CD45, etc. mediated the communication between m6A associated subtypes of TME cells and tumor epithelial cells.

Conclusions

Taken together, our study firstly revealed the m6A methylation mediated intercellular communication of the tumor microenvironment in the regulation of tumor growth and antitumor immunomodulatory processes.

Introduction

Colorectal cancer (CRC) ranks third in incidence and second in mortality. Although tremendous efforts have been made to facilitate screening strategies, the prevalence of CRC has been increasing, and 1.9 million new cases were estimated in 2020. Among them, a large population is diagnosed at an advanced stage; nevertheless, the efficacy of current therapies for late-stage CRC is limited [1]. Thus, a deeper understanding of the molecular mechanisms of CRC might help to bring novel strategies to CRC prevention and therapy.

While DNA-related aberrations are widely studied in the literature, post-transcriptional alterations have been relatively less studied. Yet, these changes play a critical role in regulating the initiation and progression of CRC [2, 3]. Among these post-transcriptional alterations, N6-methyladenosine (m6A) is the most common internal modification in transcripts, comprising more than 50% of eukaryote methylations now identified as a new level of crucial epigenetic regulation for mRNA stability, splicing, and translation [4], as well as the generation of small and long non-coding RNAs [5]. Dynamic and reversible m6A methylation consists of readers (m6A-binding proteins), writers (methyltransferases), and erasers (demethylases) responsible for m6A's functions, methylation, and demethylation, respectively [6]. Importantly, m6A mRNA regulation plays a vital role in tumorigenesis, tumor development, and metastasis, and the dysregulation of m6A is closely correlated with the development and pathogenesis of CRC [7,8,9].

Nowadays, emerging evidence has proven the critical role of tumor microenvironment (TME) in the progress and metastasis of tumor. In addition, single-cell transcriptomics further revealed the complex intercellular communication between diverse subtypes of TME cells and tumor cells [10, 11]. The TME cells consists of multiple cell types in addition to tumor cells, such as cancer-associated fibroblasts (CAFs), tumor-associated macrophages (TAMs), T cells, and B cells. Recently, Huilong et al. discovered that ablation of Mettl3 in myeloid cells promotes tumor growth and metastasis in vivo [12]. Dali et al. reported that loss of YTHDF1 in dendritic cells can lead to enhanced cross-presentation of tumor antigens and the cross-priming of CD8+ T cells in vivo [13]. However, less research was performed to investigate the cell–cell interaction between m6A mRNA modification associated subtypes of TME cells with tumor cells.

Here, we investigated the influence of m6A mRNA methylation on the main TME cells, including stromal cells, myeloid cells, T cells, and B cells, based on 65,362 single-cell sequencing data derived from 33 CRC tumor samples [14]. By nonnegative matrix factorization (NMF) clusters of 23 m6A RNA methylation regulators, as described previously [15, 16], it was observed that different patterns of m6A mRNA methylation in each CRC TME cell type subpopulations manifested extensive and diversity communication with tumor epithelial cells, and associated with different immune characteristics, metabolic pathways, transcription characteristics and prognosis. To the best of our knowledge, the present study reveals, for the first time, that m6A mRNA may guide intercellular communication of TME cells with tumor cells to contribute to colorectal cancer progression based on our comprehensive single-cell analysis.

Materials and methods

Study design and data collection

Single-cell mRNA sequence (scRNA-seq) data from 23 CRC patients with 33 samples in the SMC dataset, 23 tumor tissues, and 10 normal adjacent tissues were collected to analyze the landscape of 23 m6A RNA methylation modification regulators. After preliminary sample integration, we generated a gene expression and phenotype matrix for 65,362 scRNA-seq datasets. Full data were downloaded from GSE132465 in the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo). In addition, 11 public datasets of bulk mRNA sequence or microarray data for 2653 CRC patients were also obtained from The Cancer Genome Atlas (TCGA) and GEO databases (Additional file 1: Table S1). All data generated or analyzed during this study are freely available in previous publications or the public domain.

Visualization of TME cell types and subtypes in CRC

Using the Seurat package in R software, we created Seurat objects for total and individual cell types belonging to the scRNA-seq gene expression matrix based on acquisition. Then, the top 2000 genes, selected as the top variable features, were used as the basis for normalizing the scRNA-seq data for each cell by using the FindVariableFeatures of the Seurat package. Furthermore, we performed ScaleData and RunPCA functions to obtain the number of principal components (PC) based on the Seurat objects. We used “t-SNE (t-distributed stochastic neighbor embedding)” dimensionality reduction to further summarize the top principal components. Finally, with the annotated information for each cell in CRC supported by the previous article, the Idents and DimPlot functions were used to annotate and visualize the cells of the major TME cell types or subtypes.

Pseudotime trajectory analysis of m6A mRNA regulators for TME cells

To investigate the relationship of cell pseudotime trajectories with m6A regulators, we employed the Monocle R package for single-cell RNA data for all cell types in CRC [17]. Highly variable genes were set according to the following filtering criteria: mean expression  ≥  0.1 and dispersion_empirical ≥ 1* dispersion_fit. The DDRTree method was used for dimensionality reduction. Then we used the ‘plot_pseudotime_heatmap’ function to visualize heatmaps showing the dynamic expression of m6A regulators in the pseudotime trajectories of different TME cell types in CRC.

Non-negative matrix factorization of m6A mRNA regulators in TME cells

To best observe the effect of m6A-mediated regulator expression on TME cell types, we carried out a dimension reduction analysis for 23 m6A regulators in all TME cell types, using the non-negative matrix factorization algorithm (NMF R package, version 0.20.6), and identified different cell subtypes for these cell types, depending on the scRNA expression matrix. All these steps were performed in a manner similar to the previous studies [10, 18].

Identification of the marker genes of m6A-related cell subtypes in TME cells

We used the FindAllMarkers function to list the markers of each NMF cluster of each cell type in CRC. The parameters for min.pct and logfc.threshold were set as 0.15, and then we filtered the genes by using the adjusted p value < 0.05 for further research. The Dotplot function was performed to show the top highest gene expressions in each NMF cluster. The AddModuleScore function calculated the signature scores based on differentially expressed genes (DEGs) among these NMF cell clusters. The FeaturePlot function was used to show the distribution of specific signatures of NMF cluster scores in the TME of CRC. The special gene sets used in the comparisons among m6A-related clusters were listed in Additional file 1: Table S2.

Functional Enrichment Analysis for NMF m6A-related subtypes

Based on these marker genes among NMF clusters in different TME cell types, the clusterProfiler R package was used to detect Kyoto Encyclopedia of Genes and Genomes (KEGG) and the Reactome pathway database. Ctyoscape enrichment map function was used to cluster the pathways. Only gene sets with adjusted p-value < 0.05 were considered significantly enriched. In addition, 50 hallmark signatures and 113 metabolic pathways were collected from the molecular signature database (MSigDB) and a previous study [19], followed by performing gene set variation analysis (GSVA) to calculate the enrichment scores for these NMF clusters.

SCENIC analysis for NMF m6A-related subtypes

The pySCENIC package (version 0.9.0), a Python-based implementation of the SCENIC pipeline was used to investigate the gene regulatory network of transcription factors (TFs) in CRC [20]. Two gene-motif rankings (hg19-tss-centered-10 kb and hg19-500 bp-upstream) from the RcisTarget database were used to detect the transcription start site (TSS) and the gene regulatory networks in the scRNA-seq data in CRC. TFs with adjusted by Benjamini–Hochberg false discovery rate (BH-FDR) < 0.05 were considered for further research.

Cell–cell communication analysis for NMF m6A-related subtypes

CellChat, an R package described previously, contains ligand-receptor interaction databases for human and mouse that can analyze the intercellular communication networks from scRNA-seq data annotated as different cell clusters [21]. First, we used CellChat to evaluate the major signaling inputs and outputs among all NMF TME cell clusters using CellChatDB.human. Then, we used the netVisual_circle function to show the strength or weakness of cell–cell communication networks from the target cell cluster to different cell clusters in all NMF clusters. Finally, the netVisual_bubble function shows the bubble plots of significant ligand-receptor interactions between the target cell cluster and other TME NMF clusters.

Survival analyses with m6A-related signatures in public bulk RNA-sequence datasets

Based on the FindAllmarker function in the Seurat R package, we generated m6A-related gene signatures for all NMF cell clusters. Also, the main cell type of CRC TME were also calculated based on the scRNA data. Then the GSVA function was used to calculate these gene signature scores in all 12 public datasets of CRC. The log-rank test and Cox proportional hazard regression were conducted to explore the relationship between m6A-related NMF signatures and patients’ prognosis, including overall survival (OS) rate and recurrence-free survival (RFS) rate. The cutoff values of different NMF cell signatures in the different public datasets were determined by the survminer R package used to plot Kaplan–Meier curves. To obtain the prognosis of NMF m6A-related signatures, we used the RMA function from the metafor R package to pool the Cox-regression results of the same signatures from all available public datasets. Finally, the forestplot R package was used to show meta-analysis results.

Collection of immunotherapy transcriptomic

12 Immune checkpoint blockade immunotherapeutic (ICB) cohorts with FPKM or CPM transcriptomic were collected from the public database, included 8 melanoma datasets (Ulloa et al. (2013, MAGE A3,Melanoma [22]);Gide et al. (2019, anti-PD1 or anti-PD1 + CTLA4, Melanoma [23]); Nathanson (2017 CTLA4, Melanoma [24]); Hugo et al. (2016, anti-PD1, Met Melanoma [25]); Lauss et al. (2017, ACT, Melanoma [26]); Liu et al. (2019, anti-PD1, Met Melanoma [27]); Riaz et al. (2017, anti-PD1, Melanoma [28]); VanAllen (2016, CTLA4, MetMelanoma [29])), and other 4 non-melanoma cohorts (IMvigor210 (2018, anti-PDL1,Urothelial Cancer [30]); Braun et al. (2020, anti-PD1, CCRCC [31]); JaeWon et al. (2020, anti-PD1, NSCLC [32]);and Rose et al. (2021, ICB, Bladder Cancer [33])). All patients had the immune response in these cohorts.

Statistical analysis

Standard tests included Student’s t-test, Wilcoxon rank-sum test, Kruskal–Wallis test, and Chi-square test for the differences of continuous target or category variables in these cell subgroups. Pearson analysis was used to show the correlation of different cell signatures or gene expressions among TME CRC cell types. To compare the biological features of each m6A-related subtype for each TME cell type in CRC, we collected many current CRC-related gene signatures or TME cell function-related gene lists from the previous public publication. Then, we used the ComplexHeatmap or pheatmap packages to visualize the different expressions of the scale data of target variables among the NMF clusters from the TME CRC cell types. Routine statistical analyses of the present study were performed in R 4.0 software, and a two-sided p-value below 0.05 was considered statistically significant.

Results

The landscape of m6A regulators in TME cells in CRC

Overall, the scRNA-seq dataset of CRC, as described previously, was used to explore the landscape of m6A RNA methylation regulators (Fig. 1A). The SMC dataset contained 65,362 TME cells annotated with major cell types, including epithelial cells, mast cells, myeloid cells, stromal cells, T cells, and B cells, in 33 samples from 23 CRC patients (Fig. 1B). Cell-chat analysis showed diverse and distinct interactions among these cell types (Fig. 1C). Here, we then fully evaluated the significant different expression between the mean RNA expression of m6A regulators and common variables of the CRC samples by using the AverageExpression function of Seurat, such as class type (normal vs. tumor), MSI status (MSI-H vs. MSS), age group (older > 60 vs. young <  = 60), AJCC stage (I, IIA, IIIA, IIIB IIIC, and IVA) and gender (female vs. male) (Fig. 1D). Obviously, the expression of m6A regulators was indeed different among six cell types in CRC from the SMC dataset (Fig. 1E). The heatmap also shows the differential expression of m6A regulators in 33 CRC samples with significantly different percentages of major cell types (Additional file 2: Figure S1A). Furthermore, to assess the relationship between m6A regulators and immunologic state in TME cells, we also estimated the ImmuneScore using the “estimate” R package [34] for 65,362 TME cells of CRC. Among six TME cell types, we found different strong associations of writer and eraser regulators using ImmuneScore (Additional file 2: Figure S1B).

Fig. 1
figure 1

Overview of m6A RNA methylation regulators in the single-cell data for colorectal cancer. A The overall design of the present study and the data sourced from the SMC dataset (GSE132465). B Cell type annotations by using the Seurat t-distributed stochastic neighbor embedding (t-SNE) plot of 65,362 cells; C Cell–Cell communications between main six cell types by Cell chat analysis. D Average expression of m6A RNA methylation regulators in 65,362 cells according to different clinical variables in the SMC dataset containing cell types, including class types (normal vs. tumor), MSI (MSI-H vs. MSS), Age (old vs. young), Stage (I, IIA, IIIA, IIIB IIIC, and IVA), and Gender (female vs. male) by using z-score. E Heatmap distribution of m6A RNA methylation regulators in B cells, epithelial cells, mast cells, myeloid cells, stromal cells, and T cells. F NMF cluster by using the 23 m6A regulators RNA expression respectively for the main four types of TME cells (CAFs, Macrophage cells, T cells, and B cells) in the scRNA data

Full size image

Lastly, Fig. 1F showed the proportions of the m6A special NMF clusters for the interested four cell types (Fibroblasts, Macrophages, T cells, and B cells) by using the expression of 23 m6A regulators in scRNA data, respectively. All of their marker genes of these m6A related clusters were listed in the Additional file 1: Table S2 in all special TME cells.

Novel m6A-mediated fibroblasts contributed to the TME of CRC

Stromal cells in the CRC dataset could be grouped into fibroblast and non-fibroblast cells in both tumor and normal tissues of CRC (Additional file 2: Figure S2A, B). By combining TCGA-COAD and TCGA-READ into one dataset, we used the xCell algorithm to calculate tumor infiltration of fibroblast cells, and their high abundance showed a poor prognosis for CRC patients (Additional file 2: Figure S2C, p = 0.014). Also, the pseudotime analysis revealed that the m6A RNA regulators had a critical role in the trajectory process of TME cells including fibroblasts, NK cells, macrophages, CD4 + T cells, and CD8 + T cells, etc. (Fig. 2A, and Additional file 2: Figures S2D, S3). Thus, by the cell-chat analysis, we found that there were different number of ligand-receptor links between these m6A-related fibroblast clusters (named as HNRNPA2B1 + CAF-C1 (n = 1939), WTAP + CAF-C2 (n = 245), HNRNPC + CAF-C3 (n = 1194), and NoneMethy-CAF-C4 (n = 84)) and epithelial cells (Fig. 2B and Additional file 2: Figure S2E). Among them, the WTAP + CAF-C2 population had a higher percentage in tumor samples (n = 1501) than that in normal samples (n = 1961, Chi-square test p < 0.001) (Fig. 2C). Enrichment analysis with KEGG showed that the WTAP + CAF-C2 cell cluster exhibited IL-17 activity, TNF signaling, and neutrophil-related functions based on DEGs (Additional file 1: Table S3). Pan-CAF signatures, originated from a previous study [35], were also calculated, and we found that the WTAP + CAF-C2 score was strongly associated with inflammatory CAF (pan-iCAF) (Fig. 2E).

Fig. 2
figure 2

source class, reveals that the percentage of m6A-fib-C2 in tumor is higher than that in normal mucosa (p < 0.001). D Heatmap showing the activated KEGG pathway in main m6A-fib-clusters by using the DEGs among these groups (p < 0.05). E Different m6A-related fibroblast clusters were correlated with the previous signatures (p < 0.05). F Heatmap showing the significantly different activities of TFs among four m6A fibroblast cell clusters by comparing the average AUC using pySCENIC in Python software (Kruskal–Wallis test, p < 0.001). TF activity is scored using AUCell. G Heatmap showing the different average expression of common signaling pathway genes in the four m6A-fib-clusters, including collagens, ECM, MMPs, TGFb, Neo-Angio, Contractile, RAS and Proinflammatory. H Enrichment cluster analysis for activated signaling ways and functions of m6A-related fibroblast types in the Cytoscape by the REAC database

m6A regulators modified the features of fibroblast cells. A Trajectory Analysis reveals the role of m6A genes in for fibroblast cells (3462 cells). B Cell–Cell communications from m6A-related fibroblast cells to epithelial cells. C Bar plot for four m6A-fib-clusters, along with

Full size image

Next, gene regulatory network analysis showed the expression of 26 TFs was significantly different among the four clusters. It is noteworthy that TFs of REL, CEBPB, FOSL1, POU2F2, NFKB1, IRF1, and ETS1 were upregulated in the WTAP + CAF-C2 cluster (Fig. 2F). We additionally compared the expression of some surface protein genes among the four m6A-mediated CAF NMF clusters and found that CD34, LRRN3, BAMB1, P2RY6, CLDN1, and ICAM1 were higher in the WTAP + CAF-C2 (Additional file 2: Figure S2F). From the pathway heatmap (Fig. 2G), these the WTAP + CAF-C2 and NoneMethy-CAF-C4 had significantly different expressions of these pathway genes. Lastly, the enrichment map revealed that the HNRNPA2B1 + CAF-C1, WTAP + CAF-C2, and NoneMethy-CAF-C4 had different REACTOME pathway features (Fig. 2H).

m6A-mediated macrophages resembled classical features

A total of 5822 macrophages were extracted from myeloid cells and divided into tumor-associated macrophages (5586 cells) and normal macrophages (236) (Additional file 2: Figure S4A). Then, we obtained five main m6A-mac clusters (Additional file 2: Figure S4B, C), including four clusters with the expression of m6A regulators (WTAP + mac-C1, n = 1432; HNRNPC + mac-C2, n = 1538; HNRNPA2B1 + mac-C3, n = 2169; YTHDC1 & YTHDF3 + mac-C4, n = 442) and one cluster without expression of m6A regulators (NoneMethy-mac-C5, n = 241) (Fig. 3A). The average number and cell proportion of m6A-mac NMF clusters showed a significant difference between normal and tumor (Chi-squared test p < 0.001, Additional file 2: Figure S4D). Similar to fibroblasts, we also observed different number of ligand-receptor links between these m6A related macrophage clusters to tumor epithelia cells.

Fig. 3
figure 3

m6A regulators contributed to the production of tumor-associated macrophages5822 macrophages. A Cell–Cell communications between main m6A-related macrophage cells to epithelia cells by Cell chat analysis. B t-SNE plots of methy-mac-C1, methy-mac-C3, proinflammatory, C1q+, proliferating, SPP1 + mac, M1, and M2 macrophage signatures for 5582 macrophage cells. C Correlation plot for the above eight gene signatures in 5582 tumor macrophage cells. Methy-mac-C1 clusters are significantly related to proinflammatory macrophages, and methy-mac-C3 clusters are significantly related to SPP1 + and C1q + macrophage cells (p < 0.001, r > 0.5). D Heatmap showing significantly different TFs among m6A macrophage clusters by using pySCENIC in Python software to compare their average AUCs (Kruskal–Wallis test, p < 0.001). TF activity is scored using AUCell. E Heatmap showing significantly different activity of 41/113 metabolic signaling pathway scores by GSVA for 5582 cells among five methy-mac clusters (Kruskal–Wallis test, p < 0.001). F Enrichment cluster analysis for activated signaling ways and functions of m6A-related macrophage types in the Cytoscape by the REAC database

Full size image

Next, using the AddModulScore function of Seurat for these signatures (Fig. 3B, C, Additional file 2: Figure S4E and F) among TAMs, we found that WTAP + mac-C1 were significantly related to proinflammatory macrophages and HNRNPA2B1 + mac-C3 were significantly related to SPP1 + and C1q + macrophages [14]. Checkpoints expression also showed significantly different expression among five m6A-mac clusters (Additional file 2: Figure S4G). Furthermore, the SCENIC analysis of macrophages showed different activation of potential TFs among WTAP + mac-C1 and HNRNPA2B1 + mac-C3 clusters (Fig. 3D).

To assess the relationship between our m6A-mac clusters and special pathways, we used GSVA and detected that 41 out of 113 metabolic pathways were significantly different among the five m6A-mac clusters (Fig. 3E). Then, the 50 hallmark pathways showed broadly different activity among the five clusters (Additional file 2: Figure S4H). The enrichment map also revealed that the WTAP + mac-C1, HNRNPA2B1 + mac-C3 and TTHDC1 & YTHDF3 + mac-C4 had the different REACTOME pathway features (Fig. 3F).

m6A-mediated T/B cell phenotypes underscored the antitumor immune response in CRC

Among the detected 23,115 T cells, 8 main cell types, including CD4+, CD8+, Treg, NK, T helper 17, T follicular helper, etc., were identified to further analysis (Fig. 4A). A total of 5 m6A-related cell clusters were recognized by the NMF algorithm, and named as methy-T-C1 to methy-T-C5 (Fig. 4B), with the different numbers of ligand-receptor links between these m6A-related T cell clusters and tumor epithelia cells (Fig. 4C). Checkpoints expression analysis also showed significantly different expression among these m6A-related T cell clusters (Additional file 2: Figure S5A). Network regulatory analysis showed significant differential expression of TFs among these m6A clusters of T cells (Fig. 4D). In addition, to assess the overall effects of m6A-related T clusters on T cells, we found many differences in the average expression of immune genes of co-stimulation, co-inhibition, and some function-related markers. We also found many differences in the average expression of signatures among these m6A clusters of CD4 + T, CD8 + T, Treg, and NK T cells, including T exhaustion score, T cytotoxic score, T effector score, and T evasion score (Fig. 4E). Also, according to the DEGs, as listed in Additional file 1: Table S3, the enrichment map releveled that CD8 + T-C1, Treg-C1, CD8 + T-C4 and Treg-C4 had more immune function related terms (Additional file 2: Figure S5B). For 9146 B cells, NMF m6A clusters had similar ligand-receptor links to epithelial cells. No significant relationship was found between m6A-related B cell groups with IgG + plasma B cells, IgA + plasma B cells, and CD19 + CD20 + B cells, respectively (Fig. 4F–H). However, the heatmap still revealed the significantly different TFs among the m6A clusters (Fig. 4I).

Fig. 4
figure 4

NMF clusters of m6A methylation regulators for T cells and B cells. A t-SNE plot for 23,115 T cells by 8 cell types in the SMC dataset, including CD + 4, CD + 8, Treg, NK, T helper 17, T follicular helper, gamma delta T, and Unknown. B Cell–Cell communications from main m6A-related T cells to epithelia cells by Cell chat analysis. (C1, HNRNPA2B1 dominant; C2, IGF2BP1 and HNRNPC dominant; C3, IFG2BP3, WTAP, and FTO dominant; C4, ALKBH5, YTHDF2, and YTHDC1 dominant). C Bar plot showing the number and percentage of methy-T cell clusters, including the cluster without m6A methylation regulator expression (n = 1265), among different cell types of T cells. D Heatmap showing significantly different TFs among m6A clusters in CD + 4, CD + 8, Treg, and NK cells by comparing their average AUC by using pySCENIC in Python (Kruskal–Wallis test, p < 0.001). TF activity is scored using AUCell. E Heatmap showing significantly different features among methy-T clusters in CD + 4, CD + 8, Treg and NK cells, including four T function signatures (T exhaustion score, T cytotoxic score, T effector score, and T evasion score), as well as some immune stimulators, inhibitors and T cell function marker genes (Kruskal–Wallis test, p < 0.001). F Cell–Cell communications between main m6A-related B cells types by Cell chat analysis. G, H Bar plot for the number and percentage of m6A clusters by using the NMF clustering algorithm for 9146 B cells as above. J Heatmap showing significantly different TFs among m6A clusters in total B cells

Full size image

m6A-mediated TME patterns contributed the CRC prognosis and immunotherapy

To obtain the signature of main CRC TME cell types, we re-calculated the DEGs of them in the CRC scRNA data and extracted the top 30 as the cell markers. Then, according to all DEGs of m6A-mediated TME cells (Additional file 1: Table S3), we used the GSVA to calculate the m6A sub-score, and to explore the prognosis of them in the CRC patients and pan-cancer by using the meta-analysis for OS and RFS in the 1892 and 2315 CRC patients from 8 and 11 CRC cohorts, as listed in the Fig. 5, respectively. All score were divided into two groups to conduct the cox regression analysis. Interestingly, along with the changing of the main dominated m6A genes in special m6A-mediated sub-cell types, we found that the recurrence-free survival (Fig. 5A, Additional file 2: Table S6) and overall survival (Fig. 5B, Additional file 2: Table S7) rates of them were significantly different among these sub-clusters, including CAF, macrophages, CD8 + T, Treg and B cells. In addition, we used the logistic regression method to observe the same meaningful phenomena of m6A sub TME cells for predicting the immune response for patients who underwent the immunotherapy in the 13 public cancer cohorts, including clear cell renal cell carcinoma (ccRCC), non-small cell lung cancer (NSCLC), met-melanoma, melanoma, urothelial cancer, and bladder cancer (Fig. 5C, and Additional file 1: Table S8). Finally, we inspected the prognosis of m6A sub clusters in the pan-cancer patients which were listed in Additional file 2: Figure S6A and B, and found that different cell sub clusters contributed to different cancers with the significant prognosis.

Fig. 5
figure 5

Overall of the prognosis and immunotherapy response of m6A-related cells types (GSVA score) in the bulk sequence from public cohorts. The cut-off were calculated by the survival R packages (A). RFS analysis (data from 11 CRC cohorts); B OS analysis (data from 8 CRC cohorts); C immunotherapy response analysis (data from 13 immunotherapy cohorts with response rate)

Full size image

m6A-mediated TME patterns enhanced the intercellular communication

By the cell-chat analysis, we listed all ligand-receptor pairs of intercellular communication, including MIF −  (CD74 + CXCR4), MIF −  (CD74 + CD44), MDK–NCL, LGALS9 − CD45, CLEC2D − KLRB1, CLEC2C − KLRB1, CLEC2B − KLRB1, APP − CD74, CD99 − CD99, and ADGRE5 − CD55, were existed from m6A sub clusters to the tumor epithelial cells (Fig. 6A and Additional file 1: Table S4). Herein, the hypothesis of potential mechanism between them was demonstrated in Fig. 6B. Each m6A subtypes might have different strengths and ligand-receptor pairs with tumor epithelial clusters, which suggests that m6A-mediated TME cells might have more interactions with tumor cells and thus contribute to the progress of CRC.

Fig. 6
figure 6

Cell–Cell communications from main m6A-related TME cells to epithelial cells. A The significantly related ligand–receptor interactions from main m6A-related clusters to epithelial cells. B Hypothesis of the mechanism of m6A clustering TME cells affecting cell communication

Full size image

Discussion

To date, several studies have revealed the correlation between RNA m6A modification and the pathogenesis of colorectal cancer [36,37,38,39,40]. However, only a few have investigated the potential tumorigenic role of m6A-modified single cells. In the present study, we, for the first time, have comprehensively explored m6A modification regulators of main cell types in the TME of colorectal cancer and further identified the diversity cell–cell interaction between m6A associated TME cell subtypes and tumor cells at the 10X Genomic single-cell sequence level. This unique and new perspective allowed us to understand how RNA m6A modification of these diverse cellular components of TME affects the fate of individual CRC patients.

Cancer epithelial cells constitute the majority of tumor tissue and drive tumor development. Meanwhile, the heterogeneity of cancer epithelial cells indicates the responsiveness of patients to treatment and determines prognosis. Besides cancer epithelial cells, TME cells, such as multiple types of stromal cells, vascular endothelial cells and infiltrating immune cells, all support growth and promote immune evasion of solid tumors. In the study, we found that the TME cells, including stromal cells, macrophages, T cells and B cells all manifested the diverse m6A regulatory patterns and the extensive communication with tumor epithelial cells based on the single-cell analysis. Furthermore, Cellphone analysis showed that ligand-receptor pairs, including MIF −  (CD74 + CXCR4), MIF −  (CD74 + CD44), MDK–NCL, LGALS9 − CD45, CLEC2D − KLRB1, CLEC2C − KLRB1, CLEC2B − KLRB1, APP − CD74, CD99 − CD99, and ADGRE5 − CD55 mediated the communication between m6A associated subtypes of TME cells and tumor epithelial cells.

Cancer-associated fibroblasts (CAFs), as one of the critical components of stromal cells, were classified as pan-myCAFs, pan-dCAFs, pan-iCAFs, pan-nCAFs, and pan-pCAFs, according to specific molecular characteristics [35]. To date, few studies have revealed the potential role of RNA m6A modification in CAFs. In our study, we found that m6A-mediated fibroblasts manifested more extensive communications with tumor epithelial cells compared with non-m6A-mediated fibroblasts. Furthermore, WTAP-fibroblasts had a strong relationship with inflammatory-CAFs and the elevated expression of proinflammatory factors, such as CXCL1, CXCL2, CXCL3, CCL2, IL-6 and IL-7. Pathway analysis also revealed the participation of CAFs in IL-17 signaling pathway, TNF signaling pathway and the innate immune associated pathways. CAFs may shape an immunosuppressive microenvironment through the secretion of CXCL1, IL6 and CCL2 [41,42,43,44]. Therefore, we speculated that RNA m6A modification CAFs may form the immunosuppressive interaction with tumor cells to promote the progress and metastasis of tumor.

Nowadays, increasing research has revealed the significant role of RNA m6A methylation in the regulation and reprogramming of immune cells [12, 13, 45,46,47]. In the study of Lihui Dong and colleagues [45], C1q + TAMs, a subtype of TAMs, were found to be regulated by an RNA m6A program and promote CD8 + T cell dysfunction by expressing Ebi3 transcript with decreased m6A level. In addition, Lei Zhang et al. [48] reported that SPP1 + TAMs had a pro-angiogenic signature and weakened tumor immunity. By NMF clusters, we found that m6A-mediated subtypes of macrophages all manifested extensive communications with tumor cells. Correlation analysis showed that HNRNPA2B1 + mac-C3 was significantly related to SPP1 + and C1q + macrophages, and further prognosis analysis revealed the inverse correlation between the expression with survival probability. The metabolic process had a profound influence on TAMs and thus modulated cancer progression and immune responses, including glucose, glutamine and fatty acid metabolism [49]. Our research found that m6A-mediated macrophage, especially for HNRNPA2B1 + mac-C3 subtypes, manifested obvious activation of metabolism-related pathways, such as purine biosynthesis, gluconeogenesis, and cysteine and methionine metabolism et al. Moreover, except for macrophages, we found that m6A-mediated T cells (CD8+, CD4+ and regulatory T cells), NK cells and B cells also showed extensive interaction with tumor cells. Furthermore, m6A-mediated subtypes of four main T cells exhibited variable T cell active and inactive characteristics. These findings all indicate the significant role of RNA m6A methylation in immune escape and the tumor-promoting effect of macrophages and T cells.

To identify cell-specific gene regulatory networks, we performed an analysis of TFs at the single-cell level. In general, each subtype of CAFs, macrophages, B cells and various types of T cells all manifested distinct TFs characteristics. For CAFs, WTAP-fibroblasts exhibited a unique TF gene signature, such as ETS1, CEBPB, IRF1, REL and NFKB1. Previous studies have revealed the relationship between m6A modification and the expression of ETS1, CEBPB, IRF1 and REL, which suggests the role of m6A in the regulation of CAFs [50,51,52,53]. In addition, for macrophages, we observed higher activity of SPI1 and STA1 on HNRNPA2B1-mac. Similarly, the correlation between m6A modification and both SPI1 and STA1 was reported in previous research [52, 54,55,56]. Moreover, for B and T cells, we also found distinct TF characteristics of m6A-mediated cell subtypes. To sum up, m6A-mediated cell subtypes may modulate distinct TF regulatory networks to reshape and reprogram the TME. Finally, cell network analysis revealed that these m6A-mediated TME cells were closely connected and communicated with tumor cells. Notably, either m6A-mediated CAFs or immune cell subtypes had more communication with cancer epithelial cells, indicating that formation of an immunosuppressive tumor microenvironment might partially be determined by RNA m6A methylation.

Considering the complex intrinsic patterns of RNA m6A methylation in TME cells, we comprehensive summarized the relationships of these sub clusters' scores with prognosis and immune response from the public bulk RNA-seq cohorts. Clearly, the patients with different domination of m6A regulators of the TME cells had huge prognosis differences of CRC and exceedingly distinguished the immune response for patients who underwent ICB therapy, especially for the CAFs and macrophages, which revealed that the critical role of TME m6A for CRC patients in further research.

As a preliminary study, the major limitations of our analysis were that the low depth of scRNA-seq and the inadequate samples, and our conclusion need to acquire verification in more patients. Compared with bulk RNA-seq, the scRNA-seq of some m6A regulators in CRC would typically be minor and had more zero observation, which might contribute to the bias of the clustering method in our study. Nonetheless, the scRNA-seq analysis still provides us a novel view to reveal the characteristics of m6A methylation regulators in various TME single cells to reduce the tumor heterogeneity in CRC, which is a key forward step for clinical practice.

Conclusions

We, for the first time, identified specific RNA m6A-modificated cell subtypes of TME cells by using the single-cell sequencing analysis method and revealed the m6A methylation mediated intercellular communication of tumor microenvironment in the regulation of tumor growth and antitumor immunomodulatory processes.

Availability of data and materials

All data are available in a public, open access repository. The datasets used and/or analyzed during the current study are available in the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) and The Cancer Genome Atlas (TCGA) network (https://cancergenome.nih.gov/). R and other custom scripts for analyzing data are available upon reasonable request.

References

  1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, Bray F. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin. 2021;71:209–49.

    PubMed  Article  Google Scholar 

  2. Garcia-Cardenas JM, Guerrero S, Lopez-Cortes A, Armendariz-Castillo I, Guevara-Ramirez P, Perez-Villa A, Yumiceba V, Zambrano AK, Leone PE, Paz YMC. Post-transcriptional regulation of colorectal cancer: a focus on RNA-Binding Proteins. Front Mol Biosci. 2019;6:65.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  3. Gao Y, Wang H, Li H, Ye X, Xia Y, Yuan S, Lu J, Xie X, Wang L, Zhang J. Integrated analyses of m (1)A regulator-mediated modification patterns in tumor microenvironment-infiltrating immune cells in colon cancer. Oncoimmunology. 2021;10:1936758.

    PubMed  PubMed Central  Article  Google Scholar 

  4. Li M, Zha X, Wang S. The role of N6-methyladenosine mRNA in the tumor microenvironment. Biochim Biophys Acta Rev Cancer. 2021;1875: 188522.

    CAS  PubMed  Article  Google Scholar 

  5. Roignant JY, Soller M. m (6)A in mRNA: An Ancient Mechanism for Fine-Tuning Gene Expression. Trends Genet. 2017;33:380–90.

    CAS  PubMed  Article  Google Scholar 

  6. Jiang X, Liu B, Nie Z, Duan L, Xiong Q, Jin Z, Yang C, Chen Y. The role of m6A modification in the biological functions and diseases. Signal Transduct Target Ther. 2021;6:74.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  7. Shen C, Xuan B, Yan T, Ma Y, Xu P, Tian X, Zhang X, Cao Y, Ma D, Zhu X, et al. m (6)A-dependent glycolysis enhances colorectal cancer progression. Mol Cancer. 2020;19:72.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  8. Yang X, Zhang S, He C, Xue P, Zhang L, He Z, Zang L, Feng B, Sun J, Zheng M. METTL14 suppresses proliferation and metastasis of colorectal cancer by down-regulating oncogenic long non-coding RNA XIST. Mol Cancer. 2020;19:46.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  9. Chen X, Xu M, Xu X, Zeng K, Liu X, Sun L, Pan B, He B, Pan Y, Sun H, et al. METTL14 Suppresses CRC progression via regulating n6-methyladenosine-dependent primary miR-375 Processing. Mol Ther. 2020;28:599–612.

    CAS  PubMed  Article  Google Scholar 

  10. Chen YP, Yin JH, Li WF, Li HJ, Chen DP, Zhang CJ, Lv JW, Wang YQ, Li XM, Li JY, et al. Single-cell transcriptomics reveals regulators underlying immune cell diversity and immune subtypes associated with prognosis in nasopharyngeal carcinoma. Cell Res. 2020;30:1024–42.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  11. Obradovic A, Chowdhury N, Haake SM, Ager C, Wang V, Vlahos L, Guo XV, Aggen DH, Rathmell WK, Jonasch E, et al. Single-cell protein activity analysis identifies recurrence-associated renal tumor macrophages. Cell. 2021;184 (2988–3005): e2916.

    Google Scholar 

  12. Yin H, Zhang X, Yang P, Zhang X, Peng Y, Li D, Yu Y, Wu Y, Wang Y, Zhang J, et al. RNA m6A methylation orchestrates cancer growth and metastasis via macrophage reprogramming. Nat Commun. 2021;12:1394.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  13. Han D, Liu J, Chen C, Dong L, Liu Y, Chang R, Huang X, Liu Y, Wang J, Dougherty U, et al. Anti-tumour immunity controlled through mRNA m (6)A methylation and YTHDF1 in dendritic cells. Nature. 2019;566:270–4.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  14. Lee HO, Hong Y, Etlioglu HE, Cho YB, Pomella V, Van den Bosch B, Vanhecke J, Verbandt S, Hong H, Min JW, et al. Lineage-dependent gene expression programs influence the immune landscape of colorectal cancer. Nat Genet. 2020;52:594–603.

    CAS  PubMed  Article  Google Scholar 

  15. Zaccara S, Ries RJ, Jaffrey SR. Reading, writing and erasing mRNA methylation. Nat Rev Mol Cell Biol. 2019;20:608–24.

    CAS  PubMed  Article  Google Scholar 

  16. Chong W, Shang L, Liu J, Fang Z, Du F, Wu H, Liu Y, Wang Z, Chen Y, Jia S, et al. m (6)A regulator-based methylation modification patterns characterized by distinct tumor microenvironment immune profiles in colon cancer. Theranostics. 2021;11:2201–17.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  17. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, Trapnell C. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods. 2017;14:979–82.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  18. Puram SV, Tirosh I, Parikh AS, Patel AP, Yizhak K, Gillespie S, Rodman C, Luo CL, Mroz EA, Emerick KS, et al. Single-cell transcriptomic analysis of primary and metastatic tumor ecosystems in head and neck cancer. Cell. 2017;171:1611-1624.e1624.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  19. Rosario SR, Long MD, Affronti HC, Rowsam AM, Eng KH, Smiraglia DJ. Pan-cancer analysis of transcriptional metabolic dysregulation using The Cancer Genome Atlas. Nat Commun. 2018;9:5330.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  20. Kumar N, Mishra B, Athar M, Mukhtar S. Inference of gene regulatory network from single-cell transcriptomic data Using pySCENIC. Methods Mol Biol. 2021;2328:171–82.

    PubMed  Article  Google Scholar 

  21. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, Myung P, Plikus MV, Nie Q. Inference and analysis of cell-cell communication using Cell Chat. Nat Commun. 2021;12:1088.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  22. Ulloa-Montoya F, Louahed J, Dizier B, Gruselle O, Spiessens B, Lehmann FF, Suciu S, Kruit WH, Eggermont AM, Vansteenkiste J, Brichard VG. Predictive gene signature in MAGE-A3 antigen-specific cancer immunotherapy. J Clin Oncol. 2013;31:2388–95.

    CAS  PubMed  Article  Google Scholar 

  23. Gide TN, Quek C, Menzies AM, Tasker AT, Shang P, Holst J, Madore J, Lim SY, Velickovic R, Wongchenko M, et al. Distinct immune cell populations define response to anti-PD-1 Monotherapy and anti-PD-1/Anti-CTLA-4 combined therapy. Cancer Cell. 2019;35 (238–255): e236.

    Google Scholar 

  24. Nathanson T, Ahuja A, Rubinsteyn A, Aksoy BA, Hellmann MD, Miao D, Van Allen E, Merghoub T, Wolchok JD, Snyder A, Hammerbacher J. Somatic mutations and neoepitope homology in melanomas treated with CTLA-4 blockade. Cancer Immunol Res. 2017;5:84–91.

    CAS  PubMed  Article  Google Scholar 

  25. Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, Berent-Maoz B, Pang J, Chmielowski B, Cherry G, et al. Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma. Cell. 2017;168:542.

    CAS  PubMed  Article  Google Scholar 

  26. Lauss M, Donia M, Harbst K, Andersen R, Mitra S, Rosengren F, Salim M, Vallon-Christersson J, Torngren T, Kvist A, et al. Mutational and putative neoantigen load predict clinical benefit of adoptive T cell therapy in melanoma. Nat Commun. 2017;8:1738.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  27. Liu D, Schilling B, Liu D, Sucker A, Livingstone E, Jerby-Arnon L, Zimmer L, Gutzmer R, Satzger I, Loquai C, et al. Integrative molecular and clinical modeling of clinical outcomes to PD1 blockade in patients with metastatic melanoma. Nat Med. 2019;25:1916–27.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  28. Riaz N, Havel JJ, Makarov V, Desrichard A, Urba WJ, Sims JS, Hodi FS, Martin-Algarra S, Mandal R, Sharfman WH, et al. Tumor and microenvironment evolution during Immunotherapy with Nivolumab. Cell. 2017;171 (934–949): e916.

    Google Scholar 

  29. Van Allen EM, Miao D, Schilling B, Shukla SA, Blank C, Zimmer L, Sucker A, Hillen U, Foppen MHG, Goldinger SM, et al. Genomic correlates of response to CTLA-4 blockade in metastatic melanoma. Science. 2015;350:207–11.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  30. Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, Kadel EE III, Koeppen H, Astarita JL, Cubas R, et al. TGFbeta attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature. 2018;554:544–8.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  31. Braun DA, Hou Y, Bakouny Z, Ficial M, Sant’ Angelo M, Forman J, Ross-Macdonald P, Berger AC, Jegede OA, Elagina L, et al. Interplay of somatic alterations and immune infiltration modulates response to PD-1 blockade in advanced clear cell renal cell carcinoma. Nat Med. 2020;26:909–18.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  32. Cho JW, Hong MH, Ha SJ, Kim YJ, Cho BC, Lee I, Kim HR. Genome-wide identification of differentially methylated promoters and enhancers associated with response to anti-PD-1 therapy in non-small cell lung cancer. Exp Mol Med. 2020;52:1550–63.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  33. Rose TL, Weir WH, Mayhew GM, Shibata Y, Eulitt P, Uronis JM, Zhou M, Nielsen M, Smith AB, Woods M, et al. Fibroblast growth factor receptor 3 alterations and response to immune checkpoint inhibition in metastatic urothelial cancer: a real world experience. Br J Cancer. 2021;125:1251–60.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  34. Yoshihara K, Shahmoradgoli M, Martinez E, Vegesna R, Kim H, Torres-Garcia W, Trevino V, Shen H, Laird PW, Levine DA, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.

    PubMed  Article  CAS  Google Scholar 

  35. Galbo PM, Zang X, Zheng D. Molecular Features of Cancer-associated Fibroblast Subtypes and their Implication on Cancer Pathogenesis, Prognosis, and Immunotherapy Resistance. Clin Cancer Res. 2021;27:2636–47.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  36. Chen RX, Chen X, Xia LP, Zhang JX, Pan ZZ, Ma XD, Han K, Chen JW, Judde JG, Deas O, et al. N (6)-methyladenosine modification of circNSUN2 facilitates cytoplasmic export and stabilizes HMGA2 to promote colorectal liver metastasis. Nat Commun. 2019;10:4695.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  37. Li T, Hu PS, Zuo Z, Lin JF, Li X, Wu QN, Chen ZH, Zeng ZL, Wang F, Zheng J, et al. METTL3 facilitates tumor progression via an m (6)A-IGF2BP2-dependent mechanism in colorectal carcinoma. Mol Cancer. 2019;18:112.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  38. Wu Y, Yang X, Chen Z, Tian L, Jiang G, Chen F, Li J, An P, Lu L, Luo N, et al. m (6)A-induced lncRNA RP11 triggers the dissemination of colorectal cancer cells via upregulation of Zeb1. Mol Cancer. 2019;18:87.

    PubMed  PubMed Central  Article  Google Scholar 

  39. Chen H, Gao S, Liu W, Wong CC, Wu J, Wu J, Liu D, Gou H, Kang W, Zhai J, et al. RNA N (6)-Methyladenosine Methyltransferase METTL3 Facilitates Colorectal Cancer by Activating the m (6)A-GLUT1-mTORC1 Axis and Is a Therapeutic Target. Gastroenterology. 2021;160 (1284–1300): e1216.

    Google Scholar 

  40. Sun L, Wan A, Zhou Z, Chen D, Liang H, Liu C, Yan S, Niu Y, Lin Z, Zhan S, et al. RNA-binding protein RALY reprogrammes mitochondrial metabolism via mediating miRNA processing in colorectal cancer. Gut. 2021;70:1698–712.

    CAS  PubMed  Article  Google Scholar 

  41. Zhao Q, Huang L, Qin G, Qiao Y, Ren F, Shen C, Wang S, Liu S, Lian J, Wang D, et al. Cancer-associated fibroblasts induce monocytic myeloid-derived suppressor cell generation via IL-6/exosomal miR-21-activated STAT3 signaling to promote cisplatin resistance in esophageal squamous cell carcinoma. Cancer Lett. 2021;518:35–48.

    CAS  PubMed  Article  Google Scholar 

  42. Miyake M, Hori S, Morizawa Y, Tatsumi Y, Nakai Y, Anai S, Torimoto K, Aoki K, Tanaka N, Shimada K, et al. CXCL1-mediated interaction of cancer cells with tumor-associated macrophages and cancer-associated fibroblasts promotes tumor progression in human bladder cancer. Neoplasia. 2016;18:636–46.

    CAS  PubMed  PubMed Central  Article  Google Scholar 

  43. Kalluri R. The biology and function of fibroblasts in cancer. Nat Rev Cancer. 2016;16:582–98.

    CAS  PubMed  Article  Google Scholar 

  44. Kobayashi H, Enomoto A, Woods SL, Burt AD, Takahashi M, Worthley DL. Cancer-associated fibroblasts in gastrointestinal cancer. Nat Rev Gastroenterol Hepatol. 2019;16:282–95.

    PubMed  Article  Google Scholar 

  45. Dong L, Chen C, Zhang Y, Guo P, Wang Z, Li J, Liu Y, Liu J, Chang R, Li Y, et al. The loss of RNA N (6)-adenosine methyltransferase Mettl14 in tumor-associated macrophages promotes CD8 (+) T cell dysfunction and tumor growth. Cancer Cell. 2021;39 (945–957): e910.

    Google Scholar 

  46. Wang H, Hu X, Huang M, Liu J, Gu Y, Ma L, Zhou Q, Cao X. Mettl3-mediated mRNA m (6)A methylation promotes dendritic cell activation. Nat Commun. 1898;2019:10.

    Google Scholar 

  47. Tong J, Cao G, Zhang T, Sefik E, Amezcua Vesely MC, Broughton JP, Zhu S, Li H, Li B, Chen L, et al. m (6)A mRNA methylation sustains Treg suppressive functions. Cell Res. 2018;28:253–6.

    PubMed  PubMed Central  Article  Google Scholar 

  48. Zhang L, Li Z, Skrzypczynska KM, Fang Q, Zhang W, O’Brien SA, He Y, Wang L, Zhang Q, Kim A, et al. Single-cell analyses inform mechanisms of myeloid-targeted therapies in colon cancer. Cell. 2020;181 (442–459): e429.

    Google Scholar 

  49. Vitale I, Manic G, Coussens LM, Kroemer G, Galluzzi L. Macrophages and metabolism in the tumor microenvironment. Cell Metab. 2019;30:36–50.

    CAS  PubMed  Article  Google Scholar 

  50. Chen Y, Peng C, Chen J, Chen D, Yang B, He B, Hu W, Zhang Y, Liu H, Dai L, et al. WTAP facilitates progression of hepatocellular carcinoma via m6A-HuR-dependent epigenetic silencing of ETS1. Mol Cancer. 2019;18:127.

    PubMed  PubMed Central  Article  CAS  Google Scholar 

  51. Bechara R, Amatya N, Bailey RD, Li Y, Aggor FEY, Li DD, Jawale CV, Coleman BM, Dai N, Gokhale NS, et al. The m (6)A reader IMP2 directs autoimmune inflammation through an IL-17- and TNFalpha-dependent C/EBP transcription factor axis. Sci Immunol. 2021;6:89.

    Article  CAS  Google Scholar 

  52. Wang L, Hui H, Agrawal K, Kang Y, Li N, Tang R, Yuan J, Rana TM. m (6) A RNA methyltransferases METTL3/14 regulate immune responses to anti-PD-1 therapy. EMBO J. 2020;39: e104514.

    CAS  PubMed  PubMed Central  Google Scholar 

  53. He J, Zhou M, Yin J, Wan J, Chu J, Jia J, Sheng J, Wang C, Yin H, He F. METTL3 restrains papillary thyroid cancer progression via m (6)A/c-Rel/IL-8-mediated neutrophil infiltration. Mol Ther. 2021;29:1821–37.

    CAS  PubMed  Article  Google Scholar 

  54. Weng H, Huang H, Wu H, Qin X, Zhao BS, Dong L, Shi H, Skibbe J, Shen C, Hu C, et al. METTL14 inhibits hematopoietic stem/progenitor differentiation and promotes leukemogenesis via mRNA m (6)A Modification. Cell Stem Cell. 2018;22 (191–205): e199.

    Google Scholar 

  55. Wu G, Suo C, Yang Y, Shen S, Sun L, Li ST, Zhou Y, Yang D, Wang Y, Cai Y, et al. MYC promotes cancer progression by modulating m (6) A modifications to suppress target gene translation. EMBO Rep. 2021;22: e51519.

    CAS  PubMed  PubMed Central  Google Scholar 

  56. Liu Y, Liu Z, Tang H, Shen Y, Gong Z, Xie N, Zhang X, Wang W, Kong W, Zhou Y, Fu Y. The N (6)-methyladenosine (m (6)A)-forming enzyme METTL3 facilitates M1 macrophage polarization through the methylation of STAT1 mRNA. Am J Physiol Cell Physiol. 2019;317:C762–75.

    CAS  PubMed  Article  Google Scholar 

Download references

Acknowledgements

Not available.

Funding

This work was supported by the National Natural Science Foundation of China (No: 81972012), and Health Bureau Foundation of Zhejiang Province (No. 2022RC185).

Author information

Author notes

  1. Yuzhen Gao, Hao Wang and Shipeng Chen contributed equally to this work

Affiliations

  1. Department of Clinical Laboratory, Sir Run Run Shaw Hospital of Zhejiang University School of Medicine, Qingchun East Road, Jianggan District, Hangzhou, 310016, Zhejiang, China

    Yuzhen Gao, Rui An, Yadong Chu, Guoli Li, Yanzhong Wang, Xinyou Xie & Jun Zhang

  2. Department of Gastroenterology, Second Affiliated Hospital of Zhejiang University School of Medicine, Hangzhou, 310009, Zhejiang, China

    Hao Wang

  3. Department of Laboratory Medicine, Shanghai Eastern Hepatobiliary Surgery Hospital, Shanghai, 200438, China

    Shipeng Chen

  4. Institution of Gastroenterology, Zhejiang University, Hangzhou, 310009, Zhejiang, China

    Hao Wang

  5. Key Laboratory of Precision Medicine in Diagnosis and Monitoring Research of Zhejiang Province, Hangzhou, 310013, Zhejiang, China

    Yuzhen Gao, Rui An, Yadong Chu, Guoli Li, Yanzhong Wang, Xinyou Xie & Jun Zhang

Authors
  1. Yuzhen Gao
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Hao Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Shipeng Chen
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Rui An
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Yadong Chu
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Guoli Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Yanzhong Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Xinyou Xie
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Jun Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Conception and Design: JZ, YG, HW, and SC; Financial Support: JZ; Provision of Study Materials: YG, HW, and SC; Collection and Assembly of Data: All authors; Data Analysis and Interpretation: YG, HW, and SC; Manuscript Writing: YG, HW, and SC; Manuscript Supervised: JZ and YX. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Xinyou Xie or Jun Zhang.

Ethics declarations

Ethics approval and consent to participate

Not available.

Consent for publication

All authors agree to publish.

Competing interests

The authors declare that they have no conflicts of interest for this work.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Table S1.

The information of including cohorts in the study. Table S2. The geneset lists using in the study. Table S3. the list of methylation-related genes in each NMF sub cluster of CRC. Table S4. Complex intercellular communication networks between the m6A-mediated immune cells clusters and epithelia cells. Table S5. REAC analysis for Each m6A-mediated subClusters. Table S6. OS of m6A-related Cell types In the Bulk sequences. Table S7. RFS of m6A-related Cell types In the Bulk sequences. Table S8. m6A-related Cell type signatures to response of ICB immunotherapy in the Bulk sequences.

Additional file 2: Figure S1.

The expression of m6A regulators and the association of them with ImmuneScore in CRC (Related to Fig. 1). A The expression of m6A regulators in CRC samples with main TME cell types. B The association of m6A expression with ImmuneScore in different TME cell types in scRNA-seq data. Figure S2. Exploration of cancer-association fibroblast cells in CRC. Related to Fig. 2. A t-SNE plot for stromal cells reveals different distribution of fibroblasts and non-fibroblasts between normal mucosa (2197 cells) and CRC tissue (2736 cells). B Heatmap showing the top four genes in normal mucosa and tumor fibroblasts. C The prognosis of fibroblast cells by xCell in TCGA-COAD and READ. D Trajectory analysis for cancer-association fibroblast cells. E Top genes in four NMF clusters for 3462 fibroblast cells, including m6A-fib-C1 (1939 cells, HNRNPA2B1 dominant), m6A-fib-C2 (245 cells, WTAP dominant), m6A-fib-C3 (1194 cells, HNRNPC dominant), and m6A-fib-C4 (84 cells, no m6A methylation). F Heatmap showing the average expression of cell surface protein genes in four m6A fibroblast cell clusters (Kruskal–Wallis test, p < 0.001). Figure S3. Heatmap of pseudotime Trajectory analysis for TME cell subtypes in TME main cell types of CRC, including macrophage, and B cells, as well as four types of T cells (CD8 + T, CD4 + T, Treg, and NK T cells). Figure S4. Features of m6A-mac clusters in CRC. Related to Fig. 3. A t-SNE plot of 5822 macrophages by their source class in SMC dataset. B Heatmap showing the correlation of four m6A methylation regulator clusters by NMF, named as methy-mac-C1 (n = 1432), methy-mac-C2 (n = 1538), methy-mac-C3 (n = 2169), and methy-mac-C4 (n = 442) for 5822 macrophage cells. C Heatmap showing the different expressions for m6A methylation regulators in macrophages cells. D Bar plot showing the number and percentage of methy-mac clusters, including the cluster without m6A methylation regulator expression between tumor and normal samples. E and F The m6A-mac clusters were related to proinflammatory, proliferating and SPP1 + macrophage cells. G The distribution of checkpoints gene expression among five m6A-mac clusters. H The activity of hallmark pathways among five m6A-mac clusters. Figure S5. The immune features and prognosis of m6A-T cell clusters in CRC. Related to Fig. 4. A The distribution of checkpoints gene expression among m6A-related T sub cluster cells, including CD4+, Cd8 + NK, and Regulatory T cells. B Enrichment cluster analysis for activated signaling ways and functions of m6A-related macrophage types in the Cytoscape by the REAC database. Figure S6. The dynamic effects to the prognosis of m6A-related subtype TME cells in the Pan Cancer patients. Related to Fig. 5. A The significant prognosis of the m6A-related cells, including CAF subtype cells (HNRNPA2B1 + CAF-C1; WTAP + CAF-C2; HNRNPC + CAF-C3; NoneMethy-CAF-C4), macrophage subtype (WTAP + mac-C1;HNRNPA2B1 + mac-C3;YTHDC1 &YTHDF3 + mac_C4), T subtype cells (HNRNPA2B1 + CD8 + T cells-C1;ALKBH5 + CD8 + T cells-C4; HNRNPA2B1 + Reg T cells-C1; YTHDF2 + Reg T cells-C4), and B subtype cells (HNRNPA2B1 + B-C1; WTAP + B-C2; NoneMethy + B-C4). B Tables showed dynamic changes of the number of the prognosis among different subtypes in Pan-Cancer patients.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Gao, Y., Wang, H., Chen, S. et al. Single-cell N6-methyladenosine regulator patterns guide intercellular communication of tumor microenvironment that contribute to colorectal cancer progression and immunotherapy. J Transl Med 20, 197 (2022). https://doi.org/10.1186/s12967-022-03395-7

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12967-022-03395-7

Keywords

  • Single-cell
  • m6A
  • Tumor microenvironment
  • Colorectal cancer
  • Prognosis
  • Immunotherapy

Journal of Translational Medicine

ISSN: 1479-5876

Contact us