Abstract LncRNAs and tumor microenvironment (TME) exert an important effect in antitumor immunity. Nonetheless, the role of m^6A-related lncRNA clustering patterns in prognosis, TME and immunotherapy of cervical cancer (CC) remains unknown. Here, based on 7 m^6A-related prognostic lncRNAs obtained from TCGA-CC dataset, two m^6AlncRNA clustering patterns were determined. m^6AlncRNA clusterA was characterized by immune cell infiltrates and immune activation. m^6AlncRNA clusterB was characterized by enrichment of immune evasion and tumorigenic activation pathways as well as survival and clinical stage disadvantage. Then, principal component analysis algorithms were used to construct m^6AlncRNAscore based on prognostic differentially expressed genes between two m^6AlncRNA clusters to quantify m^6AlncRNA clustering patterns. m^6AlncRNAscore was an independent prognostic protective factor. Higher Th2 and Treg cells and enrichment of immunosuppressive pathways were observed in the low-m^6AlncRNAscore group, with poorer survival. High-m^6AlncRNAscore was characterized by increased infiltration of activated CD8 T cell, enrichment of immune activation pathways, lower IL-10 and TGF-beta1 levels, and higher immunophenscore values, indicating inflamed TME and better anti-tumor immunotherapy efficacy. Quantitative Real-Time Polymerase Chain Reaction was used for detection of m^6A-related prognostic lncRNAs. Collectively, we identified two m^6AlncRNA clustering patterns which play a nonnegligible role in the prognosis, TME heterogeneity and immunotherapy of CC patients. Subject terms: Cancer, Genetics, Immunology, Biomarkers, Oncology Introduction Cervical cancer (CC) is the fourth leading cause of cancer deaths in women globally^[32]1. In China, 59,000 women died of CC in 2020, accounting for about 17% of all CC deaths worldwide in the same year^[33]2. High mortality of CC is associated with recurrence and an advanced stage at diagnosis. Patients with recurrent and/or advanced CC have limited treatment options and poor prognosis, with a 5-year survival probability of 17%^[34]3. Immunotherapy represented by immune checkpoint inhibitors (ICIs) provides a promising perspective for cancer treatment. However, the overall response rate of ICIs was only 10–25% in previous clinical trials for CC^[35]4. The mechanisms behind the poor response of ICIs deserve further investigation. Increasing evidence indicates that tumor microenvironment (TME) not only influences tumor cell growth and metastasis, but also strongly affects tumor immune escape and immunotherapy efficacy^[36]5,[37]6. Therefore, the heterogeneity and complexity of TME and novel biomarkers associated with TME should be further explored to predict immunotherapeutic response and provide new therapeutic targets for tumors. N6-methyladenosine (m^6A) modification, the most common epigenetic modification in eukaryotic messenger RNAs and long non-coding RNAs (lncRNAs), plays a crucial role in RNA processing, splicing, degradation, and translation, thereby affecting cell self-renewal, differentiation, tumorigenesis, and tumor progression^[38]7–[39]9. m6A modification is a dynamic reversible process regulated by methyltransferases, binding proteins, and demethylases^[40]10. Methyltransferases are composed of METTL3/14/16, ZC3H13, RBM15, WTAP, VIRMA, and RBM15B, and catalyze the RNA methylation modification process^[41]10–[42]12. Demethylases include FTO^[43]12 and ALKBH3/5^[44]12, and mediate the RNA methylation removal process. Binding proteins consist of YTHDF1/2/3, YTHDC1/2, HNRNPA2B1, LRPPRC, FMR1, TRMT112, ZCCHC4, NUDT21, CPSF6, SETD2, SRSF3, SRSF10, XRN1, NXF1, PRRC2A, IGF2BP1/2/3, IGFBP3, and RBMX, exerting a vital role in carcinogenesis, invasion, and metastasis by combining with m^6A motif^[45]10–[46]12. Although over 85% of the human genome is transcribed, less than 3% of the transcripts encode protein, and the remaining transcripts mainly are non-coding RNAs^[47]13. LncRNAs, more than 200 nucleotides in length, constitute the largest group of ncRNAs and play a key role in transcriptional and post-transcriptional regulation^[48]14. It has been reported that lncRNAs exert an important effect on oncogenesis, metastasis, TME, and tumor immune escape and might be potential therapeutic targets for cancer^[49]14–[50]18. However, the relationship between m^6A-related lncRNA clustering patterns and TME immune infiltration remains unclear. In this study, we established a scoring system, m^6AlncRNAscore, to quantify the m^6AlncRNA clustering patterns in individual patients with CC. Further, we explored the independent prognostic value of m^6AlncRNAscore in the overall survival (OS), progression-free survival (PFS), and disease-specific survival (DSS), and the potential predictive role in immunotherapy efficacy. Additionally, we explored the correlation of m^6AlncRNA clustering patterns with clinicopathologic characteristics, TME immune infiltration, and somatic mutation, as well as the potential mechanisms in CC. Finally, we validated the expression of 4 m^6A-related prognostic lncRNAs in tumor samples and normal tissues. Materials and methods Data acquisition The RNA sequencing and somatic mutation data were downloaded from the Cancer Genome Atlas (TCGA, [51]https://portal.gdc.cancer.gov/) database. The immunophenscore (IPS) data were downloaded from the Cancer Immunome Atlas (TCIA, [52]https://tcia.at/) database. The clinical information was downloaded from UCSC Xena ([53]https://xenabrowser.net/). The sequencing data of 306 CC samples and 3 normal tissues were downloaded. The clinical information was summarized in Supplementary Table [54]1. Patients with OS less than or equal to 30 days were excluded, leaving the remaining 273 CC patients were enrolled into the further survival analysis. All methods were performed in accordance with the relevant guidelines and regulations. Identification of m^6A-related lncRNAs The transcriptome sequencing genes were divided into mRNA genes and lncRNA genes according to the human genome annotation data. Expression levels of 34 m^6A regulators were extracted from the mRNA data. Pearson correlation coefficient was then used to assess the correlation between m^6A regulators and lncRNAs. The lncRNAs with absolute correlation coefficient > 0.40 and P value < 0.001 were regarded as m^6A-related lncRNAs. Next, univariate Cox regression analysis was applied to determine lncRNAs associated with prognosis. The m^6A-related lncRNAs with P value < 0.05 were considered as m^6A-related prognostic lncRNAs. Consensus clustering analysis The “ConsensusClusterPlus” package (1000 iterations) was utilized to divide patients into different clustering patterns, referred to as m^6AlncRNA clusters, based on the expression levels of m^6A-related prognostic lncRNAs. According to the expression of prognosis-associated DEGs between different m^6AlncRNA clusters, patients were again classified into different clustering subtypes, termed as gene clusters. Gene set variation analysis (GSVA) To explore the difference of biological process activity between different subgroups, we conducted GSVA enrichment analysis by using the “GSVA” package. The “c2.cp.kegg.v7.4.symbols” gene sets were downloaded from MSigDB database for running GSVA analysis. Adjusted P value < 0.05 was regarded as statistically significant. Estimation of TME immune cell infiltration The single-sample gene-set enrichment analysis (ssGSEA) algorithm was used to quantify the relative infiltration levels of TME immune cells. The gene set for marking 23 immune cell types was acquired from the published study^[55]19,[56]20. The ssGSEA score was applied to represent the relative abundance of each infiltrating immune cell in each patient. Moreover, ESTIMATE algorithm was utilized to quantify the stromal and immune components for each patient. Identification of DEGs and KEGG pathway enrichment analysis The “limma” package was utilized to determine differentially expressed genes (DEGs) between different m^6AlncRNA clusters. The significance criterion for determining DEGs was set as adjusted P value < 0.001. The “clusterProfiler” package was employed to perform KEGG enrichment analysis for the DEGs to explore their potential biological behaviors. Generation of m^6AlncRNAscore To further investigate the role of m^6AlncRNA clustering patterns in CC, we constructed a scoring system, namely m^6AlncRNAscore, for individual patients based on the prognosis-associated DEGs between different m^6AlncRNA clusters. The procedures for m^6AlncRNAscore establishment were as follows: first, we extracted prognosis-associated DEGs by univariate Cox regression analysis; second, principal component analysis (PCA) was used to construct m^6AlncRNAscore based on each prognostic DEG expression for each sample. The advantage of this method is that the score is focused on the largest well-correlated (or anti-correlated) gene block in the set, while the contribution weight from genes that are not tracked with other set is reduced. Similar to GGI establishment^[57]20,[58]21, the m^6AlncRNAscore formula was as follows: [MATH: m6AlncRNAscore=∑1n(PC1i+PC2i) :MATH] where n is the total number of prognosis-associated DEGs, and i is the expression of the ith prognostic DEG. Somatic mutation analysis The “maftools” package was used to analyze the somatic mutation data of patients. Tumor mutation burden (TMB), mutations per million bases, was calculated for each patient. Then, we compared TMB between different m^6AlncRNAscore groups. Prediction of response to ICIs TCIA database provides immune profiles and antigenomes for 20 solid tumors including CC. IPS ranges from 0 to 10, represents tumor immunogenicity. The larger the IPS value, the stronger the immunogenicity. It has been validated that IPS could predict the response of tumor patients to ICIs^[59]19–[60]23. Drug sensitivity prediction We predicted the chemotherapeutic drug sensitivity based on the Genomics of Drug Sensitivity in Cancer (GDSC) database ([61]https://www.cancerrxgene.org/). oncoPredict package was used to estimate the half-maximal inhibitory concentration (IC50). Construction of ceRNA network We firstly obtained the miRNAs interacting with the m^6A-related prognostic lncRNAs by co-expression method. The miRNAs with absolute correlation coefficient < -0.20 and P value < 0.001 were regarded as related miRNAs. Then, we predicted the miRNA target genes (mRNA) by miRanda, miRDB, miRTarBase and TargetScan software. When all four kinds of software consider this gene as the target gene of miRNA, we regard this gene as the final target gene. The lncRNA-miRNA and miRNA-mRNA regulatory relationships were integrated to construct the competing endogenous RNA (ceRNA) network using Cytoscape software. Sample collection We totally collected 14 cervical tissue specimens, including 6 cervical cancer samples and 6 healthy controls in the Gynecology Department of Cancer Hospital Affiliated to Shanxi Medical University. Ethical approval was obtained from the Science Research Ethics Committee of Cancer Hospital Affiliated to Shanxi Medical University (No: SJJ202105). Informed consent and approval were provided by all participants. 6 patients with cervical cancer were newly diagnosed FIGO stage I/II patients without receiving any treatment. Quantitative real-time polymerase chain reaction (qRT-PCR) Total RNA was isolated from 12 samples using RNA TRIzol reagent (Tiangen Biotech Co., Ltd., Beijing, China, #DP451). cDNA synthesis was conducted with PrimeScriptTM RT Master Mix (Takara Biomedical Technology Co., Ltd., Beijing, China, #RR036Q). Real-time PCR was then performed with TB Green Premix Ex Taq (Takara Biomedical Technology Co., Ltd., Beijing, China, #RR820A). Relative expression of lncRNAs were normalized to GAPDH and calculated by 2-ΔΔCt method. Primers sequences are listed in Supplementary Table [62]2. Statistical analysis All statistical analyses were done in R version 4.0.4. Pearson correlation test was employed for assessing the relationship between m^6A regulators and lncRNAs. Wilcoxon rank sum test was applied to compare the quantitative data such as m^6A-related prognostic lncRNAs, immune cell infiltration, and m^6AlncRNAscore between groups. Kaplan–Meier method was utilized to draw survival curves, and log-rank test was performed to compare the survival difference between groups. The predictive accuracy of m^6AlncRNAscore was evaluated using the receiver operating characteristic (ROC) curve and area under curve (AUC). Multivariable Cox regression model was applied to ascertain the independent prognostic factors of CC. According to the association between m^6AlncRNAscore, TMB, and OS, we used the “survminer” package to find the optimal cutoff values of m^6AlncRNAscore and TMB, respectively. Patients were then divided into different groups according to the optimal cutoff value. Unless otherwise specified, a two-sided P value < 0.05 was considered statistically significant. Results Identification of m^6A-related lncRNAs Pearson correlation analysis was used to assess the relationship between 14,086 lncRNAs and 34 m^6A regulators. Total 112 lncRNAs with absolute correlation coefficient > 0.40 and P value < 0.001 were considered as m^6A-related lncRNAs. Univariate Cox regression analysis was used to explore the prognostic roles of m^6A-related lncRNAs. Of the 112 m^6A-related lncRNAs, 7 were associated with the OS (Table [63]1). These results indicated that the 7 m^6A-related lncRNAs, including [64]AC024270.4, [65]AC008124.1, [66]AL109811.2, [67]AC015922.2, [68]AC099850.4, [69]AC025176.1, and RPP38-DT, might be potential prognostic biomarkers of CC, termed as m^6A-related prognostic lncRNAs. Table 1. m^6A-related lncRNAs associated with prognosis of cervical cancer. m^6A-related prognostic lncRNAs HR (95% CI) P value [70]AC024270.4 0.048 (0.004, 0.581) 0.017 [71]AC099850.4 1.042 (1.007, 1.079) 0.018 [72]AC025176.1 0.806 (0.659, 0.984) 0.034 [73]AC008124.1 0.628 (0.422, 0.935) 0.022 [74]AL109811.2 0.801 (0.672, 0.954) 0.013 [75]AC015922.2 1.088 (1.022, 1.159) 0.009 RPP38-DT 0.068 (0.005, 0.838) 0.036 [76]Open in a new tab Expression profiles of m^6A-related prognostic lncRNAs To explore the potential biological function of m^6A-related lncRNAs in the occurrence of CC, we compared the expression profiles of 7 m^6A-related prognostic lncRNAs between CC samples and normal tissues. Notably, the tumor samples showed significantly lower expression levels of [77]AC024270.4, [78]AC008124.1, [79]AL109811.2, and [80]AC015922.2, but higher levels of [81]AC099850.4, [82]AC025176.1, and RPP38-DT, compared with the normal samples (Supplementary Fig. [83]1). These findings suggested that the 7 m^6A-related prognostic lncRNAs might possess important biological roles in the development of CC. Consensus clustering patterns of m^6A-related prognostic lncRNAs The “ConsensusClusterPlus” package, using the 7 m^6A-related prognostic lncRNAs, was utilized to explore the molecular subtypes of patients. According to the cumulative distribution function (CDF), the area under the CDF curve, the tracking plot from k = 2 to 9 (Supplementary Fig. [84]2a-c), and the number of cases in any cluster cannot be too small, the k = 2 was identified as the cluster number in our study to divide patients into two different m^6A-related lncRNA clustering patterns (Fig. [85]1a), including 209 cases in m^6AlncRNA clusterA and 64 cases in m^6AlncRNA clusterB. m^6AlncRNA clusterA had a notably better outcome compared with clusterB (Fig. [86]1b). In addition, the heatmap revealed that m^6AlncRNA clusterB was preferentially related to a high FIGO stage (Fig. [87]1c). Figure 1. [88]Figure 1 [89]Open in a new tab Differences in OS, clinicopathological and biological features, and TME characteristics between two different m^6AlncRNA clustering patterns constructed based on 7 m^6A-related prognostic lncRNAs. (A) Consensus clustering matrix for k = 2. (B) Kaplan–Meier curves of OS between m^6AlncRNA clusterA and m^6AlncRNA clusterB. (C) Heatmap and clinicopathological features of two m^6AlncRNA clusters (*, P < 0.05). (D) Heatmap and the activation states of biological pathways in two different m^6AlncRNA clustering patterns. (E) The abundance of each TME infiltrating cell in two m^6AlncRNA clustering patterns (*, P < 0.05; **, P < 0.01; ***, P < 0.001). (F) Immune score, (G) ESTIMATE score, and (H) stromal score in two m^6AlncRNA clusters. TME immune infiltration characteristics of different m^6AlncRNA clustering patterns GSVA was applied to explore the biological behaviors between different m^6AlncRNA clustering patterns. m^6AlncRNA clusterA presented enrichment pathways related to oxidative phosphorylation, cardiac muscle contraction, histidine catabolism, and arachidonic acid metabolism (Fig. [90]1d). m^6AlncRNA clusterB was significantly enriched in immune evasion, stromal, and tumorigenic activation pathways such as TGF beta signaling pathway, ubiquitin mediated proteolysis, focal adhesion, and pathways in cancer. Subsequently, we further compared TME cell infiltrates between two m^6AlncRNA clusters. ClusterA showed higher infiltration levels of multiple immune cells such as activated B cell and activated CD8 T cell than clusterB (Fig. [91]1e). The TME cell-infiltrating characteristic of clusterA was consist with its matching survival advantage. As expected, clusterA exhibited higher immune score (Fig. [92]1f) and ESTIMATE score (Fig. [93]1g), suggesting that clusterA had a significantly higher immune cell content and lower tumor purity. However, no significant difference of stromal score was displayed between two clusters (Fig. [94]1h). These results indicated that the two distinct m^6AlncRNA clustering patterns had markedly different TME. Generation of m^6AlncRNA genes and KEGG pathway enrichment analysis To further explore the potential biological behaviors of each m^6AlncRNA clustering pattern, we determined 786 DEGs between two m^6AlncRNA clusters using the limma package, named as m^6AlncRNA genes. Then, we used the clusterProfiler package to perform KEGG enrichment analysis for the DEGs. Figure [95]2a showed the pathways with significant enrichment. To our surprise, these genes presented enrichment of pathways associated with PD-L1 expression and PD-1 checkpoint pathway and infection-related pathways such as viral carcinogenesis and Epstein-Barr virus infection. Afterward, we utilized univariate Cox regression analysis to explore the effect of DEGs on the survival of patients. Among the 786 genes, 140 were positively or negatively related to the OS with P value < 0.05, regarded as m^6AlncRNA prognostic genes (Supplementary Table [96]3). Figure 2. [97]Figure 2 [98]Open in a new tab KEGG pathway analysis of DEGs between two m^6AlncRNA clustering patterns and differences in OS and clinicopathological characteristics between two gene clusters constructed based on m^6AlncRNA prognostic genes. (A) KEGG pathway analysis for DEGs between two m^6AlncRNA clusters. (B) Consensus clustering matrix of m^6AlncRNA prognostic genes for k = 2. (C) Kaplan–Meier curves of OS between gene clusterA and gene clusterB. (D) Heatmap and clinicopathological characteristics of two gene clusters (***, P < 0.001). Consensus clustering of m^6AlncRNA prognostic genes To further assess the regulation mechanism of m^6AlncRNA clustering pattern in CC, we subsequently performed consensus clustering analysis based on the 140 prognostic DEGs so as to divide patients. The consensus clustering of the 140 m^6AlncRNA prognostic genes classified patients into two different genomic subtypes, considered as gene clusterA (n = 224) and gene clusterB (n = 49), respectively (Supplementary Fig. [99]3a-c and Fig. [100]2b). We found that 4 out of the 7 m^6A-related prognostic lncRNAs showed significantly different expression levels in the two gene clusters (Supplementary Fig. [101]4a). Gene clusterB had significantly better prognosis than gene clusterA (Fig. [102]2c). Moreover, the heatmap showed gene clusterA was preferentially associated with m^6AlncRNA clusterB (Fig. [103]2d). Generation of m^6AlncRNAscore and prognostic value To reveal the role of 140 prognostic DEGs in CC, we used PCA to construct a scoring system to quantify the m^6AlncRNA clustering pattern in each patient, termed as m^6AlncRNAscore. We then divided patients into the high-m^6AlncRNAscore group (n = 147) and the low-m^6AlncRNAscore group (n = 126) according to the cutoff value -0.85 determined by the survminer package. Of the 7 m^6A-related prognostic lncRNAs, 6 displayed significantly different levels between two different m^6AlncRNAscore groups (Supplementary Fig. [104]4b). A better prognosis was observed in the high-m^6AlncRNAscore subgroup (Fig. [105]3a). The alluvial diagram showed the corresponding relationship between m^6AlncRNA cluster grouping, gene cluster grouping, m^6AlncRNAscore grouping, and survival outcomes (Fig. [106]3b). The matching rates of m^6AlncRNA clusterA with high-m^6AlncRNAscore and m^6AlncRNA clusterB with low-m^6AlncRNAscore were 67.5% and 90.6%, respectively. To assess the accuracy of m^6AlncRNAscore in predicting the OS, we performed ROC analysis and found that the 3-year AUC value was 0.708, implying that m^6AlncRNAscore had a good prognostic discrimination performance (Fig. [107]3c). Subsequently, we compared the m^6AlncRNAscore between different clustering subtypes. The m^6AlncRNAscore in m^6AlncRNA clusterA, as expected, was dramatically higher than that in m^6AlncRNA clusterB (Fig. [108]3d). Similarly, gene cluster B had significantly higher m^6AlncRNAscore compared with gene cluster A (Fig. [109]3e). Further stratified survival analysis results showed that the OS time in the low-m^6AlncRNAscore group was dramatically shorter compared with the high-m^6AlncRNAscore group, no matter for patients with grade1/2, grade3/4, age ≤ 60 years, age > 60 years, stage I/II, or stage III/IV (Fig. [110]3f-k). Figure 3. [111]Figure 3 [112]Figure 3 [113]Open in a new tab Construction of m^6AlncRNAscore and prognostic value. (A) Kaplan–Meier curves of OS between the high- and low-m^6AlncRNAscore groups. (B) Alluvial diagram showing the changes of m^6AlncRNA clusters, gene clusters, m^6AlncRNAscore, and survival state. (C) The 3-year ROC curve of m^6AlncRNAscore in the OS. (D) Comparison of m^6AlncRNAscore between two m^6AlncRNA clusters. (E) Comparison of m^6AlncRNAscore in two gene clusters. Kaplan–Meier curves of OS between the high- and low-m^6AlncRNAscore groups in (F) grade1/2 patients, (G) grade3/4 patients, (H) patients with age ≤ 60 years, (I) patients with age > 60 years, (J) stage I/II patients, and (K) stage III/IV patients. Independent prognostic value of m^6AlncRNAscore in the prognosis of CC As shown in Fig. [114]4a, the univariate Cox analysis results showed that m^6AlncRNAscore, age, and FIGO stage were significantly related to the OS of CC patients. Subsequent multivariate Cox analysis results displayed that age was an independent risky factor (HR = 2.157, P value = 0.015), but m^6AlncRNAscore was an independent protective factor (HR = 0.918, P value < 0.001) for the OS of CC patients (Fig. [115]4b). Figure 4. [116]Figure 4 [117]Open in a new tab Independent prognostic value of m^6AlncRNAscore in CC patients. (A) Univariate and (B) multivariate Cox regression analyses for the OS. Univariate Cox regression analyses for the (C) PFS and (D) DSS. Multivariate Cox regression analyses for the (E) PFS and (F) DSS. Considering the significance of PFS and DSS in tumor prognosis, we further validated the prognostic value of m^6AlncRNAscore in the PFS and DSS. In the univariate analysis, high-m^6AlncRNAscore was significantly associated with better PFS (Fig. [118]4c) and DSS (Fig. [119]4d). Moreover, further multivariate Cox regression analysis results showed that m^6AlncRNAscore was not only an independent prognostic factor for the PFS (Fig. [120]4e), but also an independent prognostic factor for the DSS (Fig. [121]4f). Our results strongly indicated that m^6AlncRNAscore had good prognostic value in CC. TME immune infiltration characteristics of different m^6AlncRNAscore groups To verify the biological behaviors of m^6AlncRNA clustering patterns in TME, we performed GSVA and ssGSEA analyses in two m^6AlncRNAscore groups. The high-m^6AlncRNAscore group was characterized by enrichment of hallmark pathways such as oxidative phosphorylation and cardiac muscle contraction (Fig. [122]5a) and infiltration of activated CD8 T cell, CD56dim natural killer cell, and monocyte (Fig. [123]5b). The low-m^6AlncRNAscore group was characterized by enrichment of immunosuppressive, stromal, and carcinogenic activation pathways such as wnt signaling pathway, TGF beta signaling pathway, MAPK signaling pathway, ERBB signaling pathway, focal adhesion, extracellular matrix (ECM)-receptor interaction, and pathways in cancer (Fig. [124]5a). Besides, the low-m^6AlncRNAscore group was rich in T helper 2 (Th2) and regulatory T (Treg) cells, two types of tumor immunosuppressive T cells (Fig. [125]5b). We then explored the expression profiles of immunosuppressive factors IL-10 (Fig. [126]5c) and TGF-beta1 (Fig. [127]5d) and found that their levels in the low-m^6AlncRNAscore group were significantly higher. The above results indicated again that m^6A-related lncRNA clustering patterns played a vital role in shaping TME landscape. Figure 5. [128]Figure 5 [129]Open in a new tab TME cell infiltration characteristics in the high- and low-m^6AlncRNAscore groups. (A) Heatmap and the activation states of biological pathways in two m^6AlncRNAscore groups. (B) The abundance of each TME infiltrating cell in two m^6AlncRNAscore groups (*, P < 0.05; **, P < 0.01; ***, P < 0.001). The expression levels of (C) IL-10 and (D) TGF-beta1 in two m^6AlncRNAscore groups. Clinical and somatic mutation characteristics of different m^6AlncRNAscore groups As expected, the m^6AlncRNAscore was significantly higher in patients with stage I/II than those with stage III/IV (Fig. [130]6a). However, no significant m^6AlncRNAscore difference was observed in different age or grade subgroups. TMB quantification analysis results showed that the low-m^6AlncRNAscore group presented no significant TBM difference in relative to the high-m^6AlncRNAscore group (Fig. [131]6b). Next, the survminer package was applied to classify patients with information of somatic mutation and survival into the high-TMB group (n = 28) and the low-TMB group (n = 226) according to the cutoff value 6.32. A better prognostic tendency was observed in the high TMB group, while no significant difference was displayed between the high- and low-TMB groups (Fig. [132]6c). Moreover, we found that patients with low-m^6AlncRNAscore and low-TMB had the worst prognosis, and the prediction power of m^6AlncRNAscore was not disturbed by TMB during the first 5 years of follow-up (Fig. [133]6d). Figure 6. [134]Figure 6 [135]Open in a new tab Clinical and somatic mutation characteristics in the high- and low-m^6AlncRNAscore groups. (A) Comparison of m^6AlncRNAscore between patients with stage I/II and patients with stage III/VI. (B) Comparison of TMB between the high- and low-m^6AlncRNAscore groups. (C) Kaplan–Meier curves of OS in high- and low-TMB groups. (D) Survival analyses for subgroup patients stratified by m^6AlncRNAscore and TMB using Kaplan–Meier curves. Patient’s response to ICIs in different m^6AlncRNAscore groups It has been reported that IPS values could predict the response of patients to ICIs. Whether patients received anti-CTLA-4 (Fig. [136]7a), anti-PD-L1 (Fig. [137]7b) or anti-CTLA-4 and anti-PD-L1 combination treatments (Fig. [138]7c), the IPS values of the high-m^6AlncRNAscore group were dramatically higher compared with the low-m^6AlncRNAscore group, suggesting that the corresponding ICI therapy responses in the high-m^6AlncRNAscore group were significantly better than those of the low-m^6AlncRNAscore group. These results indicated that patients with high-m^6AlncRNAscore were more likely to benefit from ICIs. Figure 7. [139]Figure 7 [140]Open in a new tab Comparison of IPS values between the high- and low-m^6AlncRNAscore groups. Comparison of IPS values in (A) patients receiving anti-CTLA-4 treatment, (B) patients receiving anti-PD-L1 treatment, and (C) patients receiving anti-CTLA-4 and anti-PD-L1 combination treatment. Drug sensitivity prediction in different m^6AlncRNAscore groups Then, we used the GDSC database to predict the valid drugs of high- and low-m^6AlncRNAscore groups. Supplementary Fig. [141]5 showed that the CC patients in the high-m^6AlncRNAscore group sensitively responded to 12 drugs (AZD3759, BI-2536, [142]CDK95038, Dasatinib, ERK2440, Erlotinib, Gefitinib, Ibrutinib, NU7441, Osimertinib, Sapitinib, and UMI-77). The therapy responses to 18 drugs (Afuresertib, Axitinib, AZD6482, AZD8055, Dactolisib, GNE-317, GSK269962A, Ipatasertib, Leflunomide, MK-2206, Navitoclax, Nilotinib, OSI-027, Oxaliplatin, Palbociclib, PF-4708671, Ribociclib, and SB505124) in the low-m^6AlncRNAscore group were significantly better than those of the high-m^6AlncRNAscore group. Construction of the ceRNA network of the 7 m^6A-related prognostic lncRNAs To further identify the mechanism of the 7 m^6A-related prognostic lncRNAs in CC patients, we constructed the lncRNAs-miRNAs-mRNAs ceRNA network. First, we obtained 30 miRNAs by the co-expression method. Then, miRanda, miRDB, miRTarBase and TargetScan software were used to identify 166 mRNA. Furthermore, we constructed and visualized the ceRNA network by incorporating 7 m^6A-related prognostic lncRNAs, 30 miRNAs, and 166 mRNA (Supplementary Fig. [143]6). Validation of the expression levels of four m^6A-related lncRNAs in CC samples qRT-PCR assay was used to detect the expression of [144]AC024270.4, [145]AC008124.1, [146]AC025176.1 and RPP38-DT in 6 tumor tissues and 6 normal samples. As shown in Supplementary Fig. [147]7, compared with normal tissues, cervical cancer tissues had higher [148]AC024270.4 and [149]AC008124.1, but lower [150]AC025176.1. There was no difference in [151]AC024270.4 expression between tumor samples and normal samples. Discussion Increasing evidence indicates that TME plays an indispensable role in tumor immune escape and immunotherapy efficacy^[152]5,[153]6. Therefore, identifying the role and potential regulatory mechanisms of m^6A-related lncRNA clustering patterns in survival prediction and immune infiltration will deepen our understanding of tumor immune escape and enrich the effective population for cancer immunotherapy. Here, we revealed two distinct m^6AlncRNA clustering patterns based on 7 m^6A-related prognostic lncRNAs. m^6AlncRNA clusterA was characterized by low tumor purity and high infiltration level of immune cells, such as activated B cell and activated CD8 T cell, which are key effectors of anti-tumor immunity^[154]24,[155]25. Further, m^6AlncRNA clusterA were mainly involved in immune activation pathways such as oxidative phosphorylation, cardiac muscle contraction, and arachidonic acid metabolism. Inhibition of oxidative phosphorylation alone limits the proliferation of T cells exposed to persistent antigen and promotes T cell exhaustion by upregulating genes associated with T cell exhaustion^[156]26. Therefore, we speculated that the enriched oxidative phosphorylation signaling pathway in m^6AlncRNA clusterA might promote the self-renewal of T cells, thereby enhancing antitumor immunity. Mediators released from arachidonic acid metabolic pathway play vital roles in maintaining the immune system normal function^[157]27,[158]28. Cardiac muscle contraction pathway has been reported to be associated with autoimmune diseases characterized by abnormally activated immune response^[159]29,[160]30. However, m^6AlncRNA clusterB was characterized by enrichment of immune evasion and tumorigenic activation pathways such as TGF beta signaling pathway, ubiquitin mediated proteolysis, and pathways in cancer. Existing studies imply that TGF beta signaling inhibits not only the innate immunity but also the adaptive immune system, leading to tumor immune evasion and poor response to ICIs^[161]31,[162]32. Thus, TGF beta signaling pathway is a potential tumor therapeutic target worthy of in-depth study^[163]33. Ubiquitin mediated proteolysis is involved in multiple biological processes including immune regulation and inflammatory response^[164]34. Melanoma patients with high level HECTD2, the E3 ubiquitin ligase involved in ubiquitin mediated proteolysis, had worse antitumor immunity and worse outcome of ICI treatment than those with low level HECTD2^[165]34. Consistent with TME immune infiltration characterizations, m^6AlncRNA clusterA had better clinical outcomes in relative to m^6AlncRNA clusterB, which was preferentially related to a higher FIGO stage. Further, we explored the transcriptome difference between two m^6AlncRNA clustering patterns. These DEGs were significantly linked to PD-L1 expression and PD-1 checkpoint pathway and immune-related infection pathways. PD-1/PD-L1 axis negatively regulates T cell activation by inhibiting Ras-Raf-MEK-ERK^[166]35. These findings demonstrated again that the m^6AlncRNA clustering pattern dissimilarity was associated with tumor immunity difference. Moreover, we classified patients into two different genomic subtypes and two distinct m^6AlncRNAscore groups based on the prognostic DEGs. Gene clusterA with poor prognosis was preferentially associated with m^6AlncRNA clusterB with poor prognosis. To our surprise, up to 72.9% of patients had m^6AlncRNA cluster grouping consistent with m^6AlncRNAscore grouping. By integrated analyses, we found that the m^6AlncRNAscore was a reliable and independent prognostic protective biomarker for the OS, PFS, and DSS of CC patients. No matter in the overall or in stratified survival analysis, patients with high-m6AlncRNAscore had better prognosis than low-m6AlncRNAscore. Previous studies have reported that high-TMB predicts a better clinical outcome and a higher ICI response rate in some tumors^[167]36,[168]37. In our study, patients with high-TMB had longer OS than patients with low TMB, while no significant difference was observed. Although it could not be considered that the long-term prediction ability of m^6AlncRNAscore was not affected by TMB due to the small sample size, we could determine that its prediction ability within 5 years was unaffected by TMB. Our findings indicated that m^6AlncRNA clustering patterns might affect tumor immune escape by regulating TME, and finally affected the prognosis of patients. Similar to GSVA results of m^6AlncRNA clusters, the high-m^6AlncRNAscore group was enriched in immune-inflamed pathways such as oxidative phosphorylation and cardiac muscle contraction, while the low-m^6AlncRNAscore group was significantly related to immunosuppressive pathways such as wnt signaling pathway, TGF beta signaling pathway, and MAPK signaling pathway. As discussed earlier in this study, oxidative phosphorylation signaling and cardiac muscle contraction signaling are involved in activated immune response in humans^[169]26,[170]38,[171]39. TGF beta signaling^[172]31–[173]34,[174]40, Wnt signaling^[175]41, and MAPK signaling^[176]42 have been reported to promote tumor immune escape and limit antitumor immune response. Besides, the abundance of activated CD8 T cell and CD56dim natural killer cell was higher in the high-m^6AlncRNAscore group. Consistent with the enrichment of immunosuppressive pathways, the fractions of Th2, Treg, IL-10, and TGF-beta1 were higher in patients with low-m^6AlncRNAscore. Th2 and Treg are tumor immunosuppressive cells. Researchers have already found that patients with cervical cancer express higher Th2 and Treg than women with normal cervix^[177]43. Treg, a major barrier to effective anti-tumor immunotherapy, promotes tumor immune escape by production of immunosuppressiv cytokines such as IL-10 and TGF-beta^[178]44,[179]45. TGF-beta in turn promotes the expansion of Treg^[180]44,[181]45. It has been reported recently that m^6A-related lncRNAs are novel prognostic biomarkers in lung cancer and breast cancer and are associated with TME^[182]46,[183]47. Therefore, m^6AlncRNAscore dissimilarity was significantly associated with TME difference. These findings could provide novel insights for cancer immunotherapy, that is, targeting m^6A-related lncRNAs or m^6A-related lncRNA relevant genes to reverse adverse TME, and then developing new immunotherapeutic drugs. IPS comprehensively represents tumor immunogenicity and has been verified to predict the ICI treatment efficacy^[184]19–[185]23. In our study, the high-m^6AlncRNAscore group had higher IPS values than the low-m^6AlncRNAscore group regardless of ICI therapy regimen, indicating that patients with high-m^6AlncRNAscore might be more likely to benefit from ICIs. Additionally, IL-10 and TGF-beta might be potential new targets for anti-tumor immunotherapy in patients with high-m^6AlncRNAscore, deserving further study. Furthermore, drug sensitivity analysis results found that the high-m^6AlncRNAscore group had different sensitive drugs compared with the high-m^6AlncRNAscore group. These findings demonstrated that m6AlncRNA clustering patterns could affect the therapeutic efficacy of ICIs. Previous research on the 7 m^6A-related lncRNAs was few. [186]AC099850.4 was reported to be in relation to the prognosis of ovarian cancer through lncRNA-miRNA-mRNA competing triplets^48. In our study, we also constructed the ceRNA network, which could help us better understand the mechanisms of the 7 m^6A-related lncRNAs in CC. Consistent with the lower expression levels of [187]AC024270.4 and [188]AC008124.1 in cervical cancer tissues in relative to normal tissues in TCGA database, our qRT-PCR results displayed that [189]AC024270.4 and [190]AC008124.1 had significantly lower levels in tumor samples than normal tissues. The level of [191]AC025176.1 in tumor samples was remarkably higher compared with normal tissues both in our experiment and TCGA analysis results. However, the expression of RPP38-DT in cervical cancer tissues was higher in TCGA database, but lower in our experiment in comparison with normal controls. The expression difference of RPP38-DT may be associated with the small sample size in TCGA-cervical cancer database with only three controls and our experiment with 6 pairs of case–control samples. One limitation of our study was that our results were not validated in another database. Another limitation was that the full mechanisms of the 7 m^6A-related lncRNAs in CC remained unclear. More research is needed in the future. In summary, we firstly revealed the important regulation role of m^6A-related lncRNA clustering patterns on prognosis, TME, and ICI therapeutic efficacy in CC. The difference of m^6A-related lncRNA clustering patterns was a vital factor leading to the heterogeneity and complexity of individual TME. The systematic evaluation of m^6A-related lncRNA clustering patterns will help improve our understanding of TME immune infiltration characteristics and might provide novel potential approaches for immunotherapy response prediction and patient prognostic stratification in CC. Supplementary Information [192]Supplementary Information 1.^ (92.9KB, tif) [193]Supplementary Information 2.^ (124KB, tif) [194]Supplementary Information 3.^ (119.4KB, tif) [195]Supplementary Information 4.^ (114.4KB, tif) [196]Supplementary Information 5.^ (14.1MB, tif) [197]Supplementary Information 6.^ (14.1MB, tif) [198]Supplementary Information 7.^ (73.4KB, tif) [199]Supplementary Information 8.^ (16KB, docx) [200]Supplementary Information 9.^ (31KB, docx) Acknowledgements