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

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. 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. 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. 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 (m 6 A) 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 m 6 A methylation consists of readers (m 6 A-binding proteins), writers (methyltransferases), and erasers (demethylases) responsible for m 6 A's functions, methylation, and demethylation, respectively [6]. Importantly, m 6 A mRNA regulation plays a vital role in tumorigenesis, tumor development, and metastasis, and the dysregulation of m 6 A 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 m 6 A mRNA modification associated subtypes of TME cells with tumor cells.
Here, we investigated the influence of m 6 A 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 m 6 A RNA methylation regulators, as described previously [15,16], it was observed that different patterns of m 6 A 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 m 6 A mRNA may guide intercellular communication of TME cells with tumor cells to contribute to colorectal cancer progression based on our comprehensive single-cell analysis.

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 m 6 A 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 Find-VariableFeatures 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 m 6 A mRNA regulators for TME cells
To investigate the relationship of cell pseudotime trajectories with m 6 A 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 DDR-Tree method was used for dimensionality reduction. Then we used the 'plot_pseudotime_heatmap' function to visualize heatmaps showing the dynamic expression of

Non-negative matrix factorization of m 6 A mRNA regulators in TME cells
To best observe the effect of m 6 A-mediated regulator expression on TME cell types, we carried out a dimension reduction analysis for 23 m 6 A 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 m 6 A-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 m 6 A-related clusters were listed in Additional file 1: Table S2.

Functional Enrichment Analysis for NMF m 6 A-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 m 6 A-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 m 6 A-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 m 6 A-related signatures in public bulk RNA-sequence datasets
Based on the FindAllmarker function in the Seurat R package, we generated m 6 A-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 m 6 A-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 m 6 A-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.

Statistical analysis
Standard tests included Student's t-test, Wilcoxon ranksum 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 m 6 A-related subtype for each TME cell type in CRC, we collected many current CRC-related gene signatures or TME cell functionrelated 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.

The landscape of m 6 A regulators in TME cells in CRC
Overall, the scRNA-seq dataset of CRC, as described previously, was used to explore the landscape of m 6 A 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 m 6 A 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 m 6 A regulators was indeed different among six cell types in CRC from the SMC dataset (Fig. 1E). The heatmap also shows the differential expression of m 6 A regulators in 33 CRC samples with significantly different percentages of major cell types (Additional file 2: Figure S1A). Furthermore, to assess the relationship between m 6 A 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). Lastly, Fig. 1F showed the proportions of the m 6 A special NMF clusters for the interested four cell types (Fibroblasts, Macrophages, T cells, and B cells) by using the expression of 23 m 6 A regulators in scRNA data, respectively. All of their marker genes of these m 6 A related clusters were listed in the Additional file 1: Table S2 in all special TME cells.

Novel m 6 A-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 m 6 A 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 ligandreceptor links between these m 6 A-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).

m 6 A-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 m 6 A-mac clusters (Additional file 2: Figure S4B, C), including four clusters with the expression of m 6 A 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 m 6 A regulators (NoneMethy-mac-C5, n = 241) (Fig. 3A). The average number and cell proportion of m 6 A-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 ligandreceptor links between these m 6 A related macrophage clusters to tumor epithelia cells. 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 m 6 A-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 m 6 A-mac clusters and special pathways, we used GSVA and detected that 41 out of 113 metabolic pathways were significantly different among the five m 6 A-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).  (Fig. 4A). A total of 5 m 6 A-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 m 6 A-related T cell clusters and tumor epithelia cells (Fig. 4C). Checkpoints expression analysis also showed significantly different expression among these m 6 A-related T cell clusters (Additional file 2: Figure S5A). Network regulatory analysis showed significant differential expression of TFs among these m 6 A clusters of T cells (Fig. 4D). In addition, to assess the overall effects of m 6 A-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 m 6 A 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 m 6 A clusters had similar ligand-receptor links to epithelial cells. No significant relationship was found between m 6 A-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 m 6 A clusters (Fig. 4I).

m 6 A-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 m 6 A-mediated TME cells (Additional file 1: Table S3), we used the GSVA to calculate the m 6 A 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 m 6 A genes in special m 6 A-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  20:197 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  NT5E  TMIGD2  CD226  CD86  TNFSF4  TNFSF9  ENTPD1  CD58  CD28  SLAMF1  ICOS  TNFRSF14  TNFRSF4  TNFRSF13B  TNFRSF9  CD70  CD27  TNFRSF8  CD80  CD40  TNFSF15  HAVCR1  BTNL8  EGFR  TNFSF8  TNFRSF13C  TNFSF13B  ICOSLG   BTN3A1  BTN3A2  BTN2A2  PDCD1LG2  LGALS9  TIGIT  CD274  IL10RB  BTLA  CTLA4  IL10  C10orf54  CD276  CSF1R  KDR  LAIR1  CD247  KIR2DL1  TGFB1  KIR2DL3  SLAMF7  CD160  CD244  HAVCR2  LAG3  CD96  ADORA2A  PDCD1  PVRL2  IDO1  CD48   GZMB  TBX21  IKZF2  IL2RA −4 C4 C5 CD4+ T CD8+ T NK Regulatory T   Table S4). Herein, the hypothesis of potential mechanism between them was demonstrated in Fig. 6B. Each m 6 A subtypes might have different strengths and ligandreceptor pairs with tumor epithelial clusters, which suggests that m 6 A-mediated TME cells might have more interactions with tumor cells and thus contribute to the progress of CRC.

Discussion
To date, several studies have revealed the correlation between RNA m 6 A modification and the pathogenesis of colorectal cancer [36][37][38][39][40]. However, only a few have investigated the potential tumorigenic role of m 6 A-modified single cells. In the present study, we, for the first time, have comprehensively explored m 6 A modification regulators of main cell types in the TME of colorectal cancer and further identified the diversity cell-cell interaction between m 6 A 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 m 6 A 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 m 6 20:197 ADGRE5 − CD55 mediated the communication between m 6 A 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 m 6 A modification in CAFs. In our study, we found that m 6 A-mediated fibroblasts manifested more extensive communications with tumor epithelial cells compared with non-m 6 A-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 m 6 A 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 m 6 A 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 m 6 A program and promote CD8 + T cell dysfunction by expressing Ebi3 transcript with decreased m 6 A 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 m 6 A-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 m 6 A-mediated macrophage, especially  20:197 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 m 6 A-mediated T cells (CD8+, CD4+ and regulatory T cells), NK cells and B cells also showed extensive interaction with tumor cells. Furthermore, m 6 A-mediated subtypes of four main T cells exhibited variable T cell active and inactive characteristics. These findings all indicate the significant role of RNA m 6 A 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 m 6 A modification and the expression of ETS1, CEBPB, IRF1 and REL, which suggests the role of m 6 A 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 m 6 A 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 m 6 A-mediated cell subtypes. To sum up, m 6 A-mediated cell subtypes may modulate distinct TF regulatory networks to reshape and reprogram the TME. Finally, cell network analysis revealed that these m 6 A-mediated TME cells were closely connected and communicated with tumor cells. Notably, either m 6 A-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 m 6 A methylation.
Considering the complex intrinsic patterns of RNA m 6 A 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 m 6 A 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 RNAseq, the scRNA-seq of some m 6 A 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 m 6 A 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 m 6 A-modificated cell subtypes of TME cells by using the single-cell sequencing analysis method and revealed the m 6 A methylation mediated intercellular communication of tumor microenvironment in the regulation of tumor growth and antitumor immunomodulatory processes.
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 m 6 A regulators and the association of them with ImmuneScore in CRC (Related to Fig. 1). A The expression of m 6 A regulators in CRC samples with main TME cell types. B The association of m 6 A expression with ImmuneScore in different TME cell types in scRNA-seq data. Figure S2. Exploration of cancer-association fibroblast cells in CRC. Related to  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 m 6 A-mac clusters in CRC. Related to Fig. 3.