The landscape of m6A regulators in healthy and cervical cancer samples
In this study, we investigated the role of 23 m6A RNA methylation regulatory genes in cervical cancer. As shown in Figure 1B, the expression level of 23 m6A regulators varied greatly between two groups. The expression levels of METTL16 and FTO were significantly decreased, while expression levels of RBM15, YTHDF2, and HNRNPA2B1 were significantly increased in tumor tissues (P <0.05). Further, we determined the prevalence of somatic mutations of 23 m6A regulators in cervical cancer. A total of 37 of the 289 (12.8%) samples had genetic alterations of m6A regulators, including 6 types of mutations which are depicted in different colors (Figure 1A). LRPPRC showed the highest mutation frequency, followed by ZC3H13, YTHDC2, and RBM15, while VIRMA, RNPA2B1, and GFBP1/2/3 did not show any mutations in cervical cancer samples (Fig. 1A). Next, the 4 highest mutated genes containing 2 readers (LRPPRC and YTHDC2) and 2 writers (ZC3H13 and RBM15) were divided into a wild-type group and mutation group. The gene expression of other m6A regulators was calculated for these 2 groups. We found that tumors with mutations of YTHDC2 showed a low expression of writer gene RBM15 (Fig. 1E), and reader genes IGFBP2, RBMX, and YTHDF1 were significantly downregulated in ZC3H13-mutant as well as RBM15-mutant tumors compared to wild-type tumors (Fig. 1F-H). Other m6A regulators had no expression differences in the wild-type and mutation groups (P >0.05).
Further analysis of 23 m6A regulators revealed that CNV mutations were common. FMR1, RBMX, HNRNPC, METTL3, YTHDC1, VIRMA, LRPPRC, and HNRNPA2B1 showed widespread CNV amplification. In contrast, YTHDF1/2/3, ALKBH5, METTL14/16, IGFBP1/2/3, WTAP, HNRNPA2B1, YTHDC2, RBM15/15B, ZC3H13, and FTO had frequent CNV deletions (Fig. 1D). Interestingly, YTHDF1 and YTHDF3 had CNV deletions only. The location of CNV alterations of 23 m6A regulators on chromosomes is shown in Figure 1C. These results indicate significant differences between normal and cervical cancer samples in the genomic and transcriptomic landscape of m6A regulators, indicating that expression alterations and genetic variation in m6A regulators played a important role in regulating the malignant progression of cervical cancer.
Unsupervised clustering of 23 m6A regulators in the early cervical cancer cohort
OS data and clinical information from 1 GEO dataset (GSE44001) and the TCGA-CC cohort were included in our study. To further explore the characteristics of the m6A modification phenotypes, we selected m6A regulators to identify subgroups of early cervical cancer samples. We utilized consensus clustering analysis to stratify samples with qualitatively different m6A modification patterns based on the expression of 23 m6A regulators. The coefficient of variation among clusters was calculated according to the number of categories, and K = 3 for early cervical cancer was identified as the optimal choice (Fig. 2B-C). We named these 3 clusters as m6A clusters A, B, and C, respectively, and included 208 cases in pattern cluster A, 195 cases in cluster B, and 201 cases in cluster C (Fig. 2A). However, Kaplan-Meier survival analysis showed no significant difference in prognosis among the 3 clusters (P =0.067, Fig. 2D).
m6A methylation modification regulatory network and enrichment analysis
The comprehensive landscape of the interactions of the 23 m6A regulators, the regulator connections, and their prognostic significance in early CC patients was shown in an m6A regulator network (Fig. 3B). In this network, positive relationships between m6A regulators are indicated with pink lines, and the sole negative correlation between HNRNPA2B1 and METTL16 is shown with a blue line. Writers METTL14 and ZC3H13 and the reader YTHDF3 were prognostic risk-related genes, while no relationship with prognosis was found among the rest of the regulators. The results revealed that cross-talk among the regulators of writers and readers may play a critical role in the formation of different m6A modification patterns.
We further investigated the potential m6A-related transcriptional expression change across 3 m6A modification patterns of early cervical cancer to estimate the underlying genetic alterations and expression perturbations within these phenotypes. A total of 115 DEGs representing the significant distinguishing index of the 3 m6A modification patterns were regarded as an m6A-related signature and depicted in the Venn diagram (Fig. 3A). We then compared the GO analyses and KEGG pathways between each of them to investigate biological responses in the 3 m6A modification patterns. The activation state of biological pathways was evaluated by GSVA enrichment analysis.
Biological functions of the m6A regulator were classified as biological processes, cellular components, or molecular functions (Fig. 3C-D). GO analysis showed that DEGs were mainly enriched in the biological processes associated with the regulation of protein serine / threonine kinase activity, cell leading edge, and SH3 domain binding. Furthermore, the analysis using "clusterProfiler" showed that these genes were significantly enriched in the C-type lectin receptor signaling pathway, an immune pathway associated with antigen presentation (Fig. 3E-F). These findings suggested that m6A modifications had a role in immune regulation in the TME.
m6A phenotype-related DEGs in early cervical cancer
Based on the 115 most representative m6A phenotype-related signature genes, unsupervised consensus clustering analysis was performed to obtain 4 stable transcriptomic phenotypes. Then the patients were divided into 4 distinct m6A gene signature subgroups with different clinicopathologic features and were defined as m6A gene clusters A, B, C, and D (Fig. 4B-C). The heatmap indicated that downregulated m6A regulators were mainly concentrated in m6A gene cluster A while most m6A regulators in cluster D were upregulated. There was no significant difference in the expression of m6A regulators in clusters B or C (Fig. 4A). Further survival analysis showed a significant prognostic difference among the 4 m6A gene signatures in early cervical cancer patients (P <0.001). The signature of m6A gene cluster D was related to the best prognosis, while the m6A gene cluster A was associated with the worst survival outcomes (Fig. 4D).
Expression levels of the 23 m6A regulators were also compared among the 4 gene signature subgroups (Fig. 4E). We observed significant differences in the expression of 17 m6A regulators among the 4 m6A gene signature subgroups (Fig. 4E). Thirteen prognostic-related m6A regulators were extracted using the Kaplan-Meier method (P <0.05, Fig. 5A-M). Low expression of FMR1, HNRNPA2B1, IGFBP2, and METTL16 was related to a worse prognosis. High expression of FTO, HNRNPC, IGFBP1, METTL14, RBM15B, ZC3H13, and YTHDF1/2/3 also showed worse OS outcomes.
Construction of the m6Sig score
Although our study identified the role of m6A modification in prognosis, the findings were based on a patient population and did not accurately predict the patterns of m6A methylation modification in individual tumors. Therefore, we established m6Sig, a scoring system to assess such impacts systematically. Based on an optimal cut-off value of 0.3120783, early cervical cancer patients were separated into a high- or low-m6Sig score group. Considering the complexity of m6A modification quantification, we depicted the workflow of m6Sig score construction with an alluvial diagram (Fig. 6C). Most patients in m6A gene cluster D were in the high-m6Sig score group, whereas most patients in m6A gene cluster A were in the low-m6Sig score group. Survival analysis indicated a significant prognostic difference between the 2 groups (P =0.002, Fig. 6D). The low-m6Sig score group was associated with a worse prognosis, while the high-m6Sig score group had better survival outcomes. We further evaluated the relationship between known biological signatures and the m6Sig score. Significant differences were found among the m6A gene cluster groups in the m6Sig score (P <0.001, Fig. 6B). Notably, m6A gene cluster D showed the highest m6Sig score, followed by m6A gene clusters B, C, and A. There was also a significant difference in m6Sig score between m6A clusters A and B and between A and C (P <0.001), but there was no differences between m6A clusters B and C (P =0.93, Fig. 6A). Figure 6A shows that m6A cluster A had the lowest m6Sig score.
Immune microenvironment characteristics in distinct m6A modification patterns
The immune infiltration of 23 immunocytes was evaluated to detect differences in immune microenvironment characteristics among the distinct m6A modification patterns. We found that 15 types of immune cells presented differential expression profiles among the 3 unsupervised clusters (Fig. 6E). The results showed that compared with m6A clusters B and C, cluster A had a higher level of infiltrated activated CD4 T cells, activated CD8 T cells, activated dendritic cells, gamma delta T cells, immature B cells, MDSC, macrophages, monocytes, natural killer T cells, neutrophils, plasmacytoid dendritic cells, regulatory T cells, T follicular helper cells, and type I T helper cells (P <0.05, Fig. 6E). These results revealed that the activity of most immunocytes increased in cluster A. Further, we analyzed the correlation between the m6Sig score and 23 immunocytes. We found that except for plasmacytoid dendritic cells, the immunocytes all had a negative relationship with the m6Sig score (Fig. 6F). These results strongly suggested that the unsupervised clustering of 23 m6A regulators was significantly correlated with TME.
Significantly mutated genes and tumor mutational signatures
Then, we analyzed the difference in the distribution of somatic mutations between the low- and high-m6Sig score groups in early cervical cancer using the "maftools" package in R. The 20 genes with the highest mutation rates and their mutation classification of early cervical cancer were summarized in a waterfall map. As shown in figures 7A and 7B, the top 5 most frequently mutated genes were TTN (29%), PIK3CA (24%), KMT2C (23%), MUC4 (20%), and MUC16 (17%) in the low-m6Sig score group, and PIK3CA (31%), TTN (27%), MUC16 (14%), MUC4 (13%), and FBXW7 (13%) in the high-m6Sig score group. Interestingly, mutations of the TTN, PIK3CA, MUC4, and MUC16 genes were more common in both subgroups. The missense mutation was the highest mutation in the 7 variant classifications. The box plot in Figure 7D shows that no difference between the low- and high-m6Sig score groups was found. m6Sig score and TMB exhibited a significant negative correlation (Fig. 7E). However, Kaplan-Meier analysis indicated that there was no significant difference in OS between the high- and low-TMB groups (P =0.165, Fig. 7C), suggesting that TMB may not be a risk-independent prognostic factor in early cervical cancer. The results provided a new perspective to explore the mechanisms of m6A methylation modification in tumor somatic mutations and the shaping of the TME landscape.