Crosstalk of RNA methylation writers defines tumor microenvironment and alisertib resistance in breast cancer

Background: The five major RNA methylation modifications (m6A, m1A, m6Am, m5C, and m7G) exert biological roles in tumorigenicity and immune response, mediated mainly by “writer” enzymes. Here, the prognostic values of the “writer” enzymes and the TCP1 role in drug resistance in breast cancer (BC) were explored for further therapeutic strategies.

Methods: We comprehensively characterized clinical, molecular, and genetic features of subtypes by consensus clustering. RNA methylation modification “Writers” and related genes_risk (RMW_risk) model for BC was constructed via a machine learning approach. Moreover, we performed a systematical analysis for characteristics of the tumor microenvironment (TME), alisertib sensitivity, and immunotherapy response. A series of experiments in vitro were carried out to assess the association of TCP1 with drug resistance.

Results: One “writer” (RBM15B) and two related genes (TCP1 and ANKRD36) were identified for prognostic model construction, validated by GSE1456, GSE7390, and GSE20685 cohorts and our follow-up data. Based on the patterns of the genes related to prognosis, patients were classified into RMW_risk-high and RMW_risk-low subtypes. Lower RMW_Score was associated with better overall survival and the infiltration of immune cells such as memory B cells. Further analysis revealed that RMW_Score presented potential values in predicting drug sensitivity and response for chemo- and immunotherapy. In addition, TCP1 was confirmed to promote BC alisertib-resistant cell proliferation and migration in vitro.

Conclusion: RMW_Score could function as a robust biomarker for predicting BC patient survival and therapeutic benefits. This research revealed a potential TCP1 role regarding alisertib resistance in BC, providing new sights into more effective therapeutic plans.

1 Introduction

Breast cancer (BC) is a complex disease affected by multiple risk factors, including genetic and environmental factors, recognized as the most common cause of cancer-related deaths among women worldwide (1–3). Recent studies have focused on investigating the underlying association of the tumor microenvironment (TME) with BC. The TME, composed of tumor cells and multiple non-malignant cell types (immune cells, fibroblasts, and endothelial cells), exerts influences at all disease stages, including tumor initiation, metastatic progression, and response to therapies (4–6). Evidence indicates that significant differences could be derived from various cellular compositions of immune infiltration in BC, thus offering prognostic and predictive values (7). In D. Hammerl and T. Karn’s investigations, tumor-infiltrating leukocytes (TILs) could serve as a sensitive biomarker, a high level of which indicated better survival, especially in patients with triple-negative breast cancer (TNBC) (8, 9). However, apart from the apparent progress achieved in immunotherapies, clinical benefits from these applications are restricted to a minority of patients (10). TILs and other components of the TME in BC were demonstrated to predict response to anticancer immunotherapies (11). Thus, a comprehensive understanding of the TME in BC will assist in enhancing clinical prediction and developing effective individual strategies.

Recent years have witnessed the identification of over 100 types of RNA modifications including methylations, cytosine modifications, isomerization of uridine, and ribose modification (12, 13). Among them, RNA methylation modifications are involved in the pathogenesis of various diseases including cancer through the dysregulation of epigenetic pathways (14, 15). Here, most attention has been focused on five major RNA methylation modifications, i.e., N⁶-methyladenosine (m⁶A), N¹-methyladenosine (m¹A), N⁶,2′-O-dimethyladenosine (m⁶Am), 5-methylcytosine (m⁵C), and 7-methylguanosine (m⁷G), reckoned to be activated by “writers” enzymes (16–23). The valid associations of the five major RNA modifications and related enzymes with the TME are reported in recent studies (24–26). However, given only one or two RNA modification “writers” contained in these studies, the antitumor effect of RNA modifications through multiple regulators lacks attention, especially in BC. Thus, our study sheds light on the correlations between numerous “writers” and the TME in BC for a holistic understanding.

Alisertib (MLN8237), an Aurora A inhibitor, has been shown to prevent Aurora A from being phosphorylated. This drug has been reported to prevent N-myc signaling and tumor growth by inhibiting the interaction between N-myc and the factor AURKA (27). By modifying the expression of the BCL-2 family (Bcl-xL), the drug might cause apoptosis (28). In a randomized clinical trial, the safety and effectiveness of weekly paclitaxel therapy with alisertib addition in metastatic breast cancer patients were assessed (29). This phase 2 study found that compared with paclitaxel, weekly paclitaxel plus oral alisertib presented improved progression-free survival (PFS) significantly, and the toxic impacts were manageable. Nevertheless, studies concerning the potential mechanism of BC cells resistant to alisertib remain few.

In this work, the prognostic roles of RNA methylation modification “writers” were explored using bioinformatics and statistical analysis, based on cases from The Cancer Genome Atlas (TCGA), the Gene Expression Omnibus (GEO) dataset, and our follow-up data. One “writer” and two related genes were identified to construct the RNA methylation modification “writers” and related genes_risk score (RMW_Score) model. It was observed that the risk signature was related to immune infiltration, therapeutic response, and genetic and molecular features. Interestingly, we also proposed the TCP1 role in resistance to alisertib, which may offer a novel target for overcoming drug resistance.

2 Materials and methods

2.1 Data collection and processing

Complete clinical information and gene expression data of BC patients were extracted from TCGA and GEO databases. With the use of the “caret” R package, the 1,089 BC patients from TCGA database (https://portal.gdc.cancer.gov/) were randomly assigned to training and validation cohorts in a 1:1 ratio. GSE1456, GSE7390, and GSE20685 cohorts were retrieved from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). A total of 194 BC samples from our follow-up data were used for further validation. The IMvigor210 cohort was included to evaluate the role of RMW_Score in predicting immunotherapy benefits. The expression of data and clinical details of the IMvigor210 cohort were obtained from http://research-pub.gene.com/IMvigor210CoreBiologies.

2.2 Clustering pattern

Unsupervised clustering algorithm was performed to identify the robust clustering of BC using the “ConsensusClusterPlus” package. According to previously published literature, 20 acknowledged RNA methylation modification “writers” were collated to explore potential prognostic biomarkers. These “writers” consist of seven m⁶A enzymes (METTL14, METTL16, WTAP, RBM15, RBM15B, ZC3H13, and VIRMA), five m¹A enzymes (TRMT61A, TRMT61B, TRMT10C, TRMT6, and TRPM5), one m⁶Am enzyme (PCIF1), five m⁵C enzymes (NSUN2, NSUN3, NSUN4, NSUN5, and NSUN6), and two m⁷G enzymes (METTL1 and WDR4).

2.3 Enrichment analysis

We identified differentially expressed genes (DEGs) via Bioconductor packages “DESeq2” and “limma”. “GSEA” and “GSVA” R packages were utilized respectively for gene set enrichment analysis (GSEA) and gene set variation analysis (GSVA) to study the variation of RNA modification patterns in biological processes. The gene sets were available from the MSigDB database.

2.4 Immune infiltration

The infiltrating level of 28 immune cell types was quantified using single-sample gene set enrichment analysis (ssGSEA). ssGSEA employed an enrichment score to represent the relative abundance of each immune cell type by package “gsva”. The normalized distribution was from 0 to 1. Via TIMER, CIBERSORT, and xCell algorithms, the levels of immune infiltrating cells were calculated. The related data were downloaded from the TIMER2.0 website (http://timer.comp-genomics.org/).

2.5 Construction of the RMW_Score and drug sensitivity

The RMW_Score of each BC patient was calculated by the formula Risk score = ${\displaystyle {\sum }_{i=1}^n(Coe{f}_i\ast X_i)}$, where i represents the RNA modification phenotype-related genes and X_i represents the expression value of each gene. Coef_i is the coefficient of each gene in the RMW_Score model. The Coef_i for RBM15B, TCP1, and ANKRD36 were −0.101981953, 0.017033828, and −1.224417758, respectively. Genomics of Drug Sensitivity in Cancer (GDSC; http://www.cancerrxgene.org/downloads) provided the antitumor drugs in cancer cell lines and relative drug targets/pathways for investigation. Spearman’s correlation analysis was applied to calculate the correlation of RMW_Score with drug sensitivity (|Rs| > 0.2 and false discovery rate (FDR)<0.05).

2.6 Cell culture and cell transfection

The MDA-MB-231, BT-549, alisertib-resistant MDA-MB-231, and BT-549 cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) (Gibco, Grand Island, NY, USA) medium containing 10% fetal bovine serum (FBS; Gibco) and 1% penicillin/streptomycin (Gibco, USA) at 37°C in a 5% CO₂ incubator. MDA-MB-231/alisertib and BT-549/alisertib cells were also cultured in 2 μg/ml of alisertib to maintain drug resistance. Small interfering RNA (siRNA) (GenePharma, Shanghai, China) was utilized to silence transiently the expression of TCP1. Based on the producer’s instructions, we conducted transfection with Lipofectamine 2000 (Invitrogen, Carlsbad, CA, USA).

2.7 Alisertib-sensitive assay

Cell Counting Kit-8 (CCK-8) was used to determine the half-maximal inhibitory concentration (IC50) for the drug-sensitivity test to alisertib. Fresh media were used to cultivate 1.5 * 10⁴ MDA-MB-231 and BT-549 cells in 96-well plates. The corresponding concentration of alisertib (0.1, 0.5, 1, 5, 10, 20, and 40 μM) was administrated to cells. Cell Counting Kit-8 (Dojindo, Japan) was used at the specified time to measure the in vitro drug sensitivity using a fluorometric microplate reader (Thermo Fisher Scientific, Waltham, MA, USA) set at 450 nm. With GraphPad 8.0, the IC50 was visually determined.

2.8 Quantitative reverse transcription polymerase reaction

TRIzol reagent (Invitrogen, USA) was used to extract total RNA from cells in accordance with the producer’s protocol. RNA was reverse transcribed into cDNA through a reverse transcription kit (Takara, Maebashi, Japan). The primers sequences were listed as follows: GAPDH-F: 5′-GACAGTCAGCCGCATCTTCT-3′; GAPDH-R: 5′-TTAAAAGCAGCCCTGGTGAC-3′; TCP1-F: 5′-CGGGATCCATGGCGGTGAAGGCCCTT-3′; TCP1-R: 5′-GCTTCTAGATCAGCCTTTAAGAGATGAC-3′.

2.9 EdU assay

BeyoClick™ EdU Cell Proliferation Kit with Alexa Fluor 594 (Beyotime, Shanghai, China) was utilized. Cells were washed with phosphate-buffered saline (PBS). DMEM (Gibco, USA) and 10 μM of EdU were then added to the plate. After 2-h incubation at 37°C/5% CO₂, the medium containing DMEM and EdU was removed by washing the cells with PBS. A solution of 4% paraformaldehyde was used as a fixative solution at room temperature for 30 min. The cells were then stained with DAPI for 3 min. At last, the cells were washed with PBS and observed with an inverted microscope.

2.10 Cell proliferation assay

MDA-MB-231/alisertib cells receiving si-NC or si-TCP1 treatment were seeded in 96-well plates. CCK-8 solution (RiboBio, Guangzhou, China) measuring 10 μl was added per well at a specified time (0, 24, 48 h, etc.) based on the protocol. We used a microplate reading element (Synergy4, USA) to measure the cell absorbance at 450 nm.

2.11 Transwell assay

The methods of the migration assay were similar to those of the invasion assay except for the membrane type in the upper transwell. The membrane for migration detection was normal, while that for invasion detection was a Matrigel-coated membrane (BD Biosciences, San Jose, CA, USA). MDA-MB-231/alisertib cells receiving si-NC or si-TCP1 were cultured with serum-free DMEM medium in the upper chamber. A complete medium including 10% FBS was placed into the bottom chamber as a chemotactic agent. After 48-h incubation, the cells were fixed with 4% paraformaldehyde and then stained with crystal violet. At last, we used a microscope to observe the difference.

2.12 Colony formation assay

The single-layer culture cell suspension during the logarithmic growth phase was diluted due to multiple gradients, and the culture plate was inoculated with the appropriate cell density. After the supernatant was discarded, PBS was used to wash the cells twice. Pure methanol measuring 5 mL was added at room temperature for 15 min to fix cells. After removal of the fixative solution, Giemsa stain was added to establish a staining solution for 10–30 min. Finally, the plate was turned upside down with a grid of transparencies overlaying, and the clones were directly counted by the naked eye.

2.13 Western blotting

radioimmunoprecipitation assay (RIPA) buffer (Beyotime, China) with protease inhibitors (Sigma-Aldrich, St. Louis, MO, USA) was used to lyse the cells. The protein lysates were separated by the 10% sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE) and then transferred by polyvinylidene fluoride (PVDF) membranes (Millipore, Billerica, MA, USA). Membranes were exposed under primary (anti-TCP1 antibody, Proteintech, Wuhan, China; anti-α-tubulin antibody, Proteintech, China) and secondary antibody incubation. After washing, the chemiluminescence system (Bio-Rad, Hercules, CA, USA) and Image Lab Software (NIH) were respectively utilized for obtaining signals and processing them.

2.14 Cell cycle arrest and apoptosis measurement by flow cytometry

MDA-MB-231/alisertib cells treated differently were washed with PBS. The cells were then fixed with 75% ethanol and kept at −20°C overnight. The cell cycle detection kit (MultiSciences, Hangzhou, China) and Annexin V-APC/PI Apoptosis Detection Kit (MultiSciences, China) were used respectively for cell cycle distribution and apoptotic cell staining. The different cells were ultimately analyzed by flow cytometry.

2.15 Statistical analysis

Spearman’s correlation analysis was used to estimate the correlation coefficient of RNA modification “writers”. Moreover, univariate Cox regression assisted in identifying statistically significant genes related to prognostic value. The association between the RMW_Score signature and prognosis was analyzed via the Kaplan–Meier survival analysis with “survival” and “survminer” packages. The cutoff point of survival information was used for dichotomy RMW_Score. Then, potential points were tested to obtain the maximum rank statistic. In addition, the classification performance of the RMW_Score signature in prognosis was evaluated through the receiver operating characteristic (ROC) curve and the area under the curve (AUC) values via the “ROC” R package. Univariate and multivariate Cox regression were employed to demonstrate the role of RMW_Score in predicting prognosis independently. A nomogram model was constructed by the “rms” R package. All statistical analyses were performed via R software (Version 4.2.1), GraphPad Prism 8.0 (GraphPad Software, San Diego, CA, USA), and SPSS 26.0 (IBM, Armonk, NY, USA). Statistical significance was defined as p< 0.05.

3 Results

3.1 Consensus clustering for BC patients in TCGA cohort

The flowchart of our work is shown in Figure 1. With consensus clustering performed based on RNA methylation modification (k = 2), 1,089 BC patients from TCGA dataset were stratified into two subgroups: cluster 1 (n = 321) and cluster 2 (n = 768, Figure 2A). Noticeably, the Kaplan–Meier curves revealed that patients in cluster 2 had better overall survival (OS; p = 0.013). To explain the significant difference in OS, Gene Ontology (GO) enrichment analysis was performed, showing that DEGs between the clusters were identified enriched in organelle fission, nuclear division, and amide binding (Figure 2B). We also observed the proportions of 28 immune cell types between clusters 1 and 2 by leveraging the ssGSEA algorithm. The distribution of 28 immune cell types is illustrated in Figure 2C, supporting that cluster 1 was associated with higher proportions of activated CD4 T cell, gamma delta T cell, memory B cell, and type 2 T helper cell. However, the proportions of most immune cell types were higher in cluster 2 than those in cluster 1, such as central memory CD4 T cell, central memory CD8 T cell, effector memory CD8 T cell, T follicular helper cell, type 1 T helper cell, type 17 T helper cell, CD56bright natural killer cell, CD56dim natural killer cell, eosinophil, macrophage, mast cell, monocyte, natural killer cell, natural killer T cell, neutrophil, and plasmacytoid dendritic cell. In addition, Figure 2D displays the distribution of immune cells in different subtypes, including clusters, tumor purity, estimate score, immune score, and stromal score.

FIGURE 1

Figure 1 Flowchart of the study design.

FIGURE 2

Figure 2 (A) Consensus clustering matrix for k = 2 and Kaplan–Meier curve revealed OS of patients in cluster 1 (red) and cluster 2 (blue). (B) Gene Ontology (GO) based on differentially expressed genes (DEGs) between two clusters. (C) The infiltrating levels of 28 immune cell types were compared between cluster 1 (red) and cluster 2 (blue). ns, no significance,*p < 0.05, **p < 0.01 and ***p < 0.001. (D) Heatmap for the landscape of immune cell types (rows) in patients (columns). Column annotations represent cluster, Tumor Purity, ESTIMATE Score, Immune Score, and Stromal Score. (E) The box plot compared the expression of different immune targets between cluster 1 (red) and cluster 2 (blue). ns,no significance,*p < 0.05, **p < 0.01 and ***p < 0.001. Waterfall maps revealed genetic alterations of BC patients in cluster 1 (F) and cluster 2 (G). (H) Relative distribution of tumor mutational burden (TMB) in cluster 1 versus cluster 2 exhibited in the boxplot. (I) Venn diagram indicating 69 survival-related pathways. (J–L) Kaplan–Meier curves for three pathways highly associated with survival, including regulation of autophagy, JAK/STAT signaling pathway, and drug metabolism cytochrome p450. OS, overall survival; BC, breast cancer.

Subsequently, more molecular and genetic insights were focused on clustering. Related gene expressions such as CD274 and CTLA4 were upregulated in cluster 1 over cluster 2 except TBX2 (Figure 2E). The results from Figures 2F, G revealed that PIK3CA accounted for 36%, the highest incidence in cluster 2, whereas TP53 had the highest incidence (51%) of mutations in cluster 1. Of note, compared with cluster 2, tumor mutational burden (TMB) was significantly enriched in cluster 1 (Figure 2H). Moreover, 69 differentially expressed pathways in two clusters were identified as related to survival (Figure 2I). GSVA was performed to investigate these differentially expressed pathways between the cluster subtypes. In that matter, BC patients with low expression of regulation of autophagy pathway had favorable prognosis, whereas high expressions of JAK/STAT signaling and drug metabolism cytochrome p450 pathways were significantly associated with worse outcomes (Figures 2J–L). Other major pathways were exhibited in the heatmap (Figure 3A). Furthermore, as displayed in Figure 3B, PCIF1 and NSUN2 had a widespread frequency of copy number variation (CNV) gain in both clusters.

FIGURE 3

Figure 3 (A) A heatmap showing KEGG pathways (rows) for each sample (columns) in both clusters via gene set variation analysis (GSVA). (B) The CNV gain and loss frequency of RNA modification “writers” in cluster 1 and cluster 2. (C) Box plots compared the expression distribution of 20 “writers” between paired tumor (red) and normal (blue) tissues. ns, no significance, *p < 0.05, **p < 0.01, ***p < 0.001. (D) Heatmap displaying the correlations among 20 “writers” of five types of RNA methylation modification in BC by Spearman’s correlation analysis. Red, negative correlation; blue, positive correlation. ^*p< 0.05, ^**p< 0.01, ^***p< 0.001. (E, F) Identification of prognostic “writers” and related genes in forest plots using univariate and multivariate Cox regression. (G) Kaplan–Meier curves for BC patients in the training cohort and validation cohort from TCGA database. (H) DEGs in the volcano plot. Red and blue dots represent upregulated and downregulated genes, respectively. (I) Gene set enrichment analysis (GSEA) for the association of hallmarks, biological process, and molecular function with DEGs. (J) Analyses for the expression of 20 “writers” in different RMW_Score groups. ns, no significance, *p < 0.05, **p < 0.01, ***p < 0.001. (K) Heatmap for the landscape of immune cell types (rows) in patients (columns). Column annotations represent risk group, risk score, Tumor Purity, ESTIMATE Score, Immune Score, and Stromal Score. (L) Correlations of RMW_Score with six immune cell types. KEGG, Kyoto Encyclopedia of Genes and Genomes; CNV, copy number variation; BC, breast cancer; TCGA, The Cancer Genome Atlas; DEGs, differentially expressed genes.

3.2 Investigation of five types of RNA methylation modification “writers”

The findings based on the cluster subtype supported the roles of RNA methylation modification in prognosis and immune infiltration. To quantify the RNA methylation modification patterns in individual BC patients, we attempted to develop a prognostic score model. The roles of 20 RNA modification “writers” were explored in this work. According to current investigations (19, 23, 30, 31), seven m⁶A modification “writers”, five m¹A modification “writers”, one m⁶Am modification “writer”, five m⁵C modification “writers”, and two m⁷G modification “writers” were included. First, we observed that compared with normal samples, the expression of “writers” such as TRMT61A, TRPM5, TRMT10C, and NSUN2 was significantly upregulated in BC samples, whereas “writers” like TRMT61B and ZC3H13 were apparently decreased in BC samples (Figure 3C). Additionally, pairwise correlations among the expression of the 20 writers were calculated to estimate their relationships (Figure 3D). Obviously, a remarkable correlation existed in the same categories. Significant correlations were also observed among different types of RNA modification writers. The expressions of METTL1, TRMT61A, NSUN5, and NSUN3 were positively correlated with other writers. On the contrary, negative correlations were present in the expression of NSUN6, RBM15, NSUN2, TRMT6, TRMT10C, and WTAP with other writers. Thus, negative correlations occurred more frequently than positive ones.

3.3 Construction and validation of a prognostic signature for BC

We then performed univariate Cox regression analysis to evaluate the relationship between these writers and the prognosis of BC patients (Figure 3E). Further analysis of these writers, 20 writers, and 82 related genes were curated (p< 0.05, $Cor\gg 0.65$. Then, machine learning was utilized to identify three genes (RBM15B, TCP1, and ANKRD36) linked with BC prognosis (Figure 3F). In that matter, TCP1 was correlated with m⁶A modification “writer” WTAP. According to the results, a scoring model termed as RMW_Score (RNA methylation modification “writers” and related genes_risk score) was constructed to quantify RNA modification patterns of individual BC patients. The RMW_Score of each patient could be calculated, and 1,089 BC cases from TCGA database were divided into RMW_Score-high and RMW_Score-low groups. Then, the Kaplan–Meier survival analysis was performed, suggesting that patients in the RMW_Score-high group had a worse outcome than the RMW_Score-low group in both the training cohort (n = 545) and validation cohort (n = 544) (p< 0.001, Figure 3G).

3.4 Immune landscapes of RMW_Score signature

To explore the potential mechanisms underlying the different prognoses of RMW_Score-high and RMW_Score-low groups, we first investigated the DEGs in the two groups. As indicated in the volcano plot (Figure 3H), a total of 1,272 DEGs consisting of 1,053 upregulated and 219 downregulated genes were identified. Then, we observed that these genes were enriched in biological processes, such as the interleukin-1-mediated signaling pathway and nucleoside triphosphate metabolic process by GSEA (Figure 3I). These DEGs were also found related to hallmarks like mTORC1 signaling and reactive oxygen species pathways. We also compared the expression of 20 “writers” in the two groups (Figure 3J) and found that in the RMW_Score-low group, NSUN6, TRMT61B, PCIF1, ZC3H13, METTL14, NSUN4, RBM15B, METTL16, NSUN3, and VIRMA were upregulated while TRMT10C, METTL1, NSUN5, and WTAP were downregulated.

Next, the ssGSEA approach was applied to estimate the infiltrating levels of immune cells in BC. The association of immune cell distribution with different subtypes (RMW_Score, tumor purity, estimate Score, immune score, and stromal score) is depicted in the heatmap (Figure 3K). In that matter, the infiltrating levels of major immune cells were positively linked with RMW_Score, including monocyte, gamma delta T cell, and activated CD4 T cell. The correlations of the six types of immune cells (activated CD4 T cell, gamma delta T cell, monocyte, activated CD8 T cell, activated dendritic cell, and CD56bright natural killer cell) with RMW_Score could be viewed in Figure 3L. Furthermore, the significant association of RMW_Score with nine immune checkpoints including CTLA4 was observed (Figure 4A). We also set sight on the distribution of TMB and found that patients with high RMW_Score had higher TMB (p< 0.001, Figure 4B). In addition, a positive correlation was observed between RMW_Score and TMB (Figure 4C). The prognosis of BC patients with high immune scores was better than that of patients with low immune scores (p = 0.019, Figure 4D). We also investigated the distribution of immune infiltrating cells among RMW_Score-high and RMW_Score-low groups via TIMER, CIBERSORT, and xCell (Figure 4E). These results suggested the involvement of the RMW_Score signature in the BC immune microenvironment.

FIGURE 4

Figure 4 (A) The difference in the expression of different immune targets between RMW_Score-high (red) and RMW_Score-low (blue) groups. ns, no significance, *p < 0.05, **p < 0.01, ***p < 0.001. (B) The TMB of two RMW_Score groups is compared and plotted. (C) The correlation of TMB with RMW_Score is shown. (D) Kaplan–Meier curves analyses revealed that high immune scores were associated with more favorable OS (p = 0.019). (E) A heatmap designed for the normalized scores of immune and stromal cell infiltrations. Blue, lower infiltration; red, higher infiltration. Wilcoxon test was used to compare the difference between the groups. ^*p< 0.05, ^**p< 0.01, ^***p< 0.001. Meanwhile, gene mutation patterns and clinical characteristics are exhibited as an annotation. TMB, tumor mutational burden; OS, overall survival.

3.5 Mutation status associated with RMW_Score signature

Moreover, 69 differentially expressed pathways linked with OS were selected (Figure 5A). A heatmap with GSVA was established to visualize and evaluate the major relative pathways, including amino sugar and nucleotide sugar metabolism, oxidative phosphorylation, and galactose metabolism pathways (Figure 5B). Of note, CNV mutations existed prevalently. PCIF1, NSUN2, and TRMT6 exhibited an extensive frequency of CNV gain no matter in the RMW_Score-high or RMW_Score-low group, as shown in Figure 5C.

FIGURE 5

Figure 5 (A) Venn diagram for 69 distinct pathways related to survival in two subgroups. (B) The survival-related pathways (rows) for each patient (columns) are shown in the heatmap. Red represents high expression, and blue represents low expression. (C) The CNV gain and loss frequency in RMW_Score-high and RMW_Score-low groups. (D) ROC curve analyses were performed to evaluate the predictive value of RMW_Score in patients (AUC of 0.659, 0.619, and 0.615 for 5-, 7.5-, and 10-year overall survival, respectively). (E, F) Forest plots for the independence of RMW_Score as a prognostic predictor by univariate and multivariate Cox regression. (G) ROC curves for RMW_Score, age, gender, stage, T, M, N, and PAM50. (H) Heatmap depicting the relationships between three patterns (one RNA modification “writer” and two related genes) and RMW_Score and clinical characteristics. (I) RMW_Score, age, stage, T, N and PAM50 were used to construct a nomogram (*, p< 0.05; ***, p< 0.001). (J) Comparison of RMW_Score expression in different PAM50 subgroups. (K) Alluvial diagram of survival status in groups with different clusters and RMW_Score. (L) The Kaplan–Meier curves for 194 patients from our follow-up data. (M) Kaplan–Meier curves suggested that the RMW_Score-low group all had better OS in the GSE1456 cohort (n = 159, p< 0.001), GSE7390 cohort (n = 198, p = 0.005), and GSE20685 cohort (n = 327, p = 0.003). (N) The Kaplan–Meier curves of OS for predicting the immunotherapeutic and chemotherapeutic benefits respectively in the IMvigor210 cohort and TCGA cohort. CNV, copy number variation; ROC, receiver operating characteristic; AUC, area under the curve; OS, overall survival.

3.6 RMW_Score signature for nomogram construction

The AUC values of the signature for predicting 5-, 7.5-, and 10-year OS were respectively 0.659, 0.619, and 0.615 (Figure 5D). To further confirm the predictive capacity of the RMW_Score signature, we employed univariate Cox regression analysis based on clinicopathological features including age, gender, stage, T, M, N, and PAM50 (Figure 5E). It implied that RMW_Score could serve as a prognostic biomarker (HR = 1.210, 95%CI 1.102–1.329, p< 0.001). Moreover, with multivariate Cox regression analysis performed, RMW_Score was an independent and robust predictive factor (Figure 5F, HR = 1.109, 95%CI 1.003–1.227, p = 0.043). The AUC of the RMW_Score signature was 0.632, validating its predictive advantage compared with other factors (Figure 5G). Additionally, a heatmap in Figure 5H showed the association of RNA modification patterns with clinical features and RMW_Score. As expected, a nomogram was constructed according to RMW_Score, PAM50, age, stage, N, and T stage (Figure 5I). As Figure 5J shows, the distribution of RMW_Score in PAM50 subtypes was different significantly. The alluvial diagram in Figure 5K reveals associations between two clusters and two RMW_Score groups in TCGA patients.

3.7 External validation and potential therapeutic value for RMW_Score signature

A total of 194 BC samples from our follow-up data were further analyzed (Figure 5L). As expected, patients with high RMW_Score had worse outcomes, verifying the reliability of the prognostic model. Consistent with these results, patients in the RMW_Score-low group from the GSE1456 cohort (n = 159, p< 0.001), GSE7390 cohort (n = 198, p = 0.005), and GSE20685 cohort (n = 327, p = 0.003) also showed an obvious survival advantage (Figure 5M). We also confirmed the predictive role of RMW_Score in immunotherapeutic benefits for BC patients. Considering the lack of published datasets of BC patients receiving immunotherapy, the urothelial cancer patients (IMvigor210) receiving anti-PD-L1 therapy were used in our study. It was noticed that patients with high RMW_Score in the IMvigor210 cohort had inferior prognoses than those with low RMW_Score (p = 0.030, Figure 5N). Moreover, 236 BC patients receiving chemotherapy selected from TCGA dataset were divided into RMW_Score-high and RMW_Score-low groups. The Kaplan–Meier curves showed that patients in the RMW_Score-low group had more survival benefits (p = 0.011, Figure 5N). Together, RMW_Score could function as an effective biomarker for managing suitable and satisfactory treatment strategies.

We evaluated the relationship between the RMW_Score and the response to medications in BC cell lines in order to better understand the implications of the RMW_Score on drug response. In the Genomics of Drug Sensitivity in Cancer database, we found 43 substantially linked relationships between RMW_Score and drug sensitivity from 345 drugs by using Spearman’s correlation analysis (Figure 6A) (32). In that matter, 41 pairs displayed drug resistance associated with the RMW_Score, including alisertib (Rs = 0.526, p< 0.001), Ara-G (Rs = 0.51, p< 0.001), and LDN-193189 (Rs = 0.492, p< 0.001). Drug sensitivity related to the RMW_Score could be noticed in two pairs, including AZD6482 (Rs = −0.212, p< 0.05) and pictilisib (Rs = −0.273, p< 0.01). As illustrated in Figure 6B, drug sensitivity correlated with RMW_Score-high was mainly targeting PI3K/mTOR signaling pathway, whereas drug sensitivity correlated with RMW_Score-low was targeting apoptosis regulation and RTK signaling pathway.

FIGURE 6

Figure 6 (A) Spearman’s analysis estimating the correlation between RMW_Score and drug sensitivity. Each column indicates a drug. Rs > 0, drug resistance; Rs< 0, drug sensitive. (B) Drug-targeted signaling pathways based on the RMW_Score. Red represents resistance; blue represents sensitivity. The significance of the correlation was indicated by the size of the point. (C) The Kaplan–Meier curves for subgroups with different TCP1 expressions. (D) The cell viability assay was used to compare drug resistance between BC cells and their respective drug-resistant cells. (E) The expression levels of TCP1 were detected by qRT-PCR. (F) The effects of expression of TCP1 on the inhibiting rate were calculated by CCK-8 assay after 48-h treatment of alisertib, and the IC50 was calculated by SPSS. (G) The efficiency of TCP1 knockdown was confirmed by Western blotting analysis. (H) EdU assays indicated the effects of TCP1 knockdown on cell proliferation in MDA-MB-231/alisertib cells. Scale bars = 100 μm. (I) Colony formation assay of siRNA-treated MDA-MB-231/alisertib cells. (J) Apoptosis analysis of the effects of TCP1 knockdown and statistical analysis are shown. (K) Cell cycle analysis of the effects of TCP1 knockdown and the cell cycle distribution is shown. (L) CCK-8 assays showed the effects of TCP1 knockdown on MDA-MB-231/alisertib cell proliferation. (M, N) Effects of TCP1 knockdown on cell migration and invasion in MDA-MB-231/alisertib cells by transwell assay. Scale bars = 100 μm. Data were analyzed with the mean ± SD of three replicates. *p< 0.05, **p< 0.01, and ***p< 0.001. BC, breast cancer; CCK-8, Cell Counting Kit-8.

3.8 TCP1 expression positively correlated with alisertib resistance

The RMW_Score model of this work was constructed based on RBM15B, TCP1, and ANKRD36, in which TCP1 was the most significant element for its lowest p-value (p< 0.001). The Kaplan–Meier analysis of TCP1 supported that TCP1 was a risk factor for BC patients, and its high expression was apparently associated with poor OS (Figure 6C, p< 0.001). Moreover, the potential role of TCP1 in regulating drug resistance of various tumors including acute myeloid leukemia (AML) and lung adenocarcinoma (LUAD) has been investigated in recent studies (33, 34). In our work, the cell viability assay and quantitative reverse transcription polymerase reaction (qRT-PCR) results respectively confirmed that MDA-MB-231/alisertib and BT-549/alisertib cells were more resistant to alisertib and had higher TCP1 expression than their parental cells significantly (Figures 6D, E). These findings suggested that TCP1 may contribute to BC cells’ resistance to alisertib. Therefore, we transiently transfected MDA-MB-231/alisertib and BT-549/alisertib cells, which have considerably greater TCP1 expression, with TCP1-specific siRNA vector to reduce TCP1 expression to investigate whether the cellular TCP1 level is connected to the resistance of BC cells to alisertib. According to evidence from cell viability assays, knocking down TCP1 in MDA-MB-231/alisertib and BT-549/alisertib cells increased alisertib’s capacity to suppress cell viability and reduced the drug’s half-life (IC50) (Figure 6F). We constructed TCP1 siRNA and verified the knockdown efficiency in MDA-MB-231/alisertib cells through Western blotting (Figure 6G).

3.9 The effects of TCP1 knockdown in MDA-MB-231/alisertib

For evaluating the effects of TCP1 knockdown in MDA-MB-231/alisertib, a series of experiments in vitro were conducted. It was observed that TCP1 knockdown significantly weakened MDA-MB-231/alisertib cell proliferation through EdU (Figure 6H). The amount of MDA-MB-231/alisertib colonies was lower in the si-TCP1 group than that in the negative control (NC) group, suggesting that the colony formation capacity could be affected by TCP1 knockdown (Figure 6I). Additionally, we discovered that the quantity of apoptotic MDA-MB-231/alisertib cells was dramatically increased in the si-TCP1 group (Figure 6J). As shown in Figure 6K, TCP1 knockdown increased the percentage of MDA-MB-231/alisertib in the S and G2 phases while decreasing the G1 phase population by flow cytometry analyses. Furthermore, the result in Figure 6L indicated the effect of TCP1 on promoting the MDA-MB-231/alisertib cell proliferation through CCK-8 assays. We also found that compared with the NC group, the si-TCP1 group was less invasive and migrated (Figures 6M, N). Collectively, these data revealed the effects of TCP1 on promoting the proliferation and migration of MDA-MB-231/alisertib and drug resistance.

4 Discussion

Increasing evidence has revealed the essential roles of RNA methylation modifications in the TME through interaction with numerous “writers”. The single type of RNA methylation modification “writer” was the main concern in most reported works of literature. However, a clear view of the functions of multiple types of “writers” in BC has not been established. Herein, the RMW_Score model (containing m⁶A, m¹A, m⁶Am, m⁵C, and m⁷G) was developed via machine learning approaches as an integrative predictor of BC patients’ survival and therapeutic responses in BC. Then, using data from 20 RNA modification enzymes, we determined two unique RNA modification patterns and identified two BC subtypes associated with RNA modification. The further landscape of immune cell, genetic, and molecular characteristics in the two clusters was explored, and cluster 2 was closely correlated with higher infiltrating levels of immune cells compared with cluster 1. We also found that 69 pathways differently expressing in the two clusters were related to survival, most of which like the JAK–STAT signaling pathway appeared to upregulate in BC patients with favorable clinical outcomes. These findings revealed the essential roles of RNA methylation modifications in the TME and BC prognosis.

Tumor cell metastasis and drug resistance are impacted by epithelial–mesenchymal transition (EMT) (35, 36), and naïve CD4⁺ T-cell recruitment might inhibit TI Tregs, restore immune tumor killing, and inhibit cancer growth and metastases (37). According to Bach et al., the JAK–STAT pathway, which regulates the expression of PDL1 and MHC class I, is activated by IFN-γ released from effector T cells in cancer cells. This signaling cascade can cause tumor cell death in a variety of ways (38, 39). Therefore, we developed the RMW_Score scoring model (consisting of one distinct RNA methylation modification pattern (RBM15B) and two related genes (TCP1 and ANKRD36)) to quantify RNA methylation modification patterns in individual tumors, predicting the prognosis of BC patients and facilitating effective therapeutic strategies. RBM15B, the paralog of RBM15, was found involved in prognostic models of several types of tumors, such as alcohol-related hepatocellular carcinoma (A-HCC), small cell lung cancer (SCLC), and melanoma (40–42). The predictive values of TCP1 and ANKRD36 were also elucidated in many investigations (43, 44). For example, it was reported that TCP1 may contribute to HCC cell proliferation and metastasis by regulating the Wnt7b/β-catenin pathway and is a molecular marker for the prognosis (44). Of note, we explored the involvement of the RMW_Score model. The RMW_Score-high subtype indicated poor survival significantly, which was verified in several cohorts including the GSE1456 cohort, GSE7390 cohort, GSE20685 cohort, and 194 samples from follow-up data.

Subsequently, we observed the association of the RMW_Score model with the tumor microenvironment and found that the amounts of immune infiltrating cells were apparently different between the two groups. The RMW_Score-low group was associated with high infiltrating levels of memory B cell, eosinophil, mast cell, natural killer cell, and plasmacytoid dendritic cell. Recent studies reported that tumor-infiltrating B cells (TIL-B) in BC were correlated with improved outcomes (45, 46). An analysis of HER2+ and triple-negative breast cancer patients from the BIG 02-98 clinical trial confirmed the correlation between positive outcomes and higher TIL-B densities (46). Similar to these studies, the results in our work implied that high amounts of memory B cells were associated with better clinical outcomes.

Additionally, the results from enrichment analyses revealed that DEGs were remarkably enriched in the interleukin-1-mediated signaling pathway, nucleoside triphosphate metabolic process, mTORC1 signaling, and reactive oxygen species pathway. In a further analysis of these pathways, we found 69 pathways related to survival, such as Notch signaling and oxidative phosphorylation (OXPHOS) pathways. Recent studies indicated that OXPHOS could be upregulated in many cancers, including lymphomas, pancreatic ductal adenocarcinoma, and breast cancer, based on which OXPHOS inhibitors exerted positive effects on treating cancers (47). For instance, marizomib (Mzb) inhibited TNBC metastasis via OXPHOS inhibition (48). Moreover, we found that BC cases from either the RMW_Score-low group or the RMW_Score-high group had high mutations of TP53 and PIK3CA and an enormous frequency of CNV gain of PCIF1, NSUN2, and TRMT6. Nevertheless, patients with high RMW_Score were dramatically associated with TMB.

Moreover, the RMW_Score model could independently predict the prognosis of BC patients revealed by univariate and multivariate Cox regression analyses, thus showing predominant advantages compared with several clinicopathological parameters. Based on RMW_Score, PAM50, age, stage, T, and N stage, a visualized nomogram was designed to offer risk score prediction for each patient. Furthermore, we also found the remarkable effects of the RMW_Score model on therapeutic response. RMW_Score was correlated with sensitivity to drugs (AZD6482 and Pictilisib) targeting PI3K/mTOR signaling pathway and resistance to drugs targeting apoptosis regulation and RTK signaling pathway. It was discovered that the drug metabolism enzyme P450 plays a significant role in affecting the development of BC. The results demonstrated that patients in the RMW_Score-high group may benefit from drugs targeting PI3K/mTOR signaling pathway. In addition to drug response, we analyzed the effects of RMW_Score as a predictor of chemotherapy and immunotherapy responses. As expected, patients achieving chemotherapy with low RMW_Score presented prolonged survival outcomes. Urothelial cancer patients (IMvigor210) in the RMW_Score-low group achieved apparent benefits from anti-PD-L1 therapy. Despite the veracity of RMW_Score in anti-PD-L1 therapy, further work is needed to assess the response in patients with BC. Together, our model could work as a robust and precise biomarker for improving therapeutic strategies.

T-complex protein 1 (TCP1), also termed CCT1, is one of eight CCT (chaperonin-containing T complex) subunits (CCT1–8). According to recent investigations, TCP1 could promote the tumorigenesis and progression of various tumors (49, 50). For instance, TCP1 has been found overexpressed in ovarian cancer, which is correlated with poor prognosis (49). Activation of the PI3K signaling pathway may contribute to the process of TCP1 upregulation in promoting cancer progression. TCP1 is also found crucial to BC patients’ survival, connected to driving oncogenes (50). Furthermore, increasing studies have reported the effects of TCP1 involved in drug resistance (34). TCP1 presented high expression in acute myeloid leukemia cells, responsible for drug resistance. In the current work, we found that TCP1 was elevated in drug-resistant MDA-MB-231/alisertib and BT-549/alisertib cells. While TCP1 overexpression imparts resistance to the therapeutic benefits of MDA-MB-231/alisertib and BT-549/alisertib, TCP1 downregulation can increase their sensitivity to alisertib. In light of these findings, it appears that alisertib drug resistance was linked to high TCP1 levels.

Although the RMW_Score model presented reliable and efficient prediction in prognosis and therapeutic response of BC patients, there still existed some limitations. Tumor heterogeneity and retrospective datasets were the main reasons. Moreover, anti-PD-L1 therapy response was analyzed in this work based on patients with urothelial cancer (IMvigor210) due to insufficient data on BC patients undergoing anti-PD-L1 therapy. Thereby, further investigation is needed to evaluate immunotherapeutic response in BC patients.

5 Conclusions

In summary, RBM15B, TCP1, and ANKRD36 were identified by machine learning approach to construct RNA methylation modification “writers” and related genes_risk score (RMW_Score) model. We found that RMW_Score could predict BC patient survival and therapeutic benefits as a sensitive marker. Additionally, TCP1 was confirmed to promote BC alisertib-resistant cell proliferation and migration in vitro.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/supplementary material.

Author contributions

XZ, JT and LS conceived and designed the study. XZ, YZ and LS performed bioinformatics analyses. YZ, CZ and HZ analyzed and interpreted data and/or supervised parts of the study. XZ and LS wrote the paper. JT and YZ revised the paper. All authors contributed to the article and approved the submitted version.

Funding

The current study was supported by the National Natural Science Foundation of China (Grant Nos. 82072931 and 82002805).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Arthur RS, Wang T, Xue X, Kamensky V, Rohan TE. Genetic factors, adherence to healthy lifestyle behavior, and risk of invasive breast cancer among women in the uk biobank. J Natl Cancer Inst (2020) 112(9):893–901. doi: 10.1093/jnci/djz241

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global cancer statistics 2020: Globocan estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin (2021) 71(3):209–49. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Britt KL, Cuzick J, Phillips KA. Key steps for effective breast cancer prevention. Nat Rev Cancer (2020) 20(8):417–36. doi: 10.1038/s41568-020-0266-x

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Flister MJ, Bergom C. Genetic modifiers of the breast tumor microenvironment. Trends Cancer (2018) 4(6):429–44. doi: 10.1016/j.trecan.2018.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Quail DF, Joyce JA. Microenvironmental regulation of tumor progression and metastasis. Nat Med (2013) 19(11):1423–37. doi: 10.1038/nm.3394

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Azizi E, Carr AJ, Plitas G, Cornish AE, Konopacki C, Prabhakaran S, et al. Single-cell map of diverse immune phenotypes in the breast tumor microenvironment. Cell (2018) 174(5):1293–308 e36. doi: 10.1016/j.cell.2018.05.060

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Ali HR, Chlon L, Pharoah PD, Markowetz F, Caldas C. Patterns of immune infiltration in breast cancer and their clinical implications: a gene-expression-based retrospective study. PloS Med (2016) 13(12):e1002194. doi: 10.1371/journal.pmed.1002194

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Karn T, Jiang T, Hatzis C, Sanger N, El-Balat A, Rody A, et al. Association between genomic metrics and immune infiltration in triple-negative breast cancer. JAMA Oncol (2017) 3(12):1707–11. doi: 10.1001/jamaoncol.2017.2140

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Hammerl D, Smid M, Timmermans AM, Sleijfer S, Martens JWM, Debets R. Breast cancer genomics and immuno-oncological markers to guide immune therapies. Semin Cancer Biol (2018) 52(Pt 2):178–88. doi: 10.1016/j.semcancer.2017.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Polk A, Svane IM, Andersson M, Nielsen D. Checkpoint inhibitors in breast cancer - current status. Cancer Treat Rev (2018) 63:122–34. doi: 10.1016/j.ctrv.2017.12.008

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Baxevanis CN, Fortis SP, Perez SA. The balance between breast cancer and the immune system: challenges for prognosis and clinical benefit from immunotherapies. Semin Cancer Biol (2021) 72:76–89. doi: 10.1016/j.semcancer.2019.12.018

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Roundtree IA, Evans ME, Pan T, He C. Dynamic rna modifications in gene expression regulation. Cell (2017) 169(7):1187–200. doi: 10.1016/j.cell.2017.05.045

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Barbieri I, Kouzarides T. Role of RNA modifications in cancer. Nat Rev Cancer (2020) 20(6):303–22. doi: 10.1038/s41568-020-0253-2

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Zhou Y, Kong Y, Fan W, Tao T, Xiao Q, Li N, et al. Principles of rna methylation and their implications for biology and medicine. BioMed Pharmacother (2020) 131:110731. doi: 10.1016/j.biopha.2020.110731

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Shi H, Wei J, He C. Where, when, and how: Context-dependent functions of rna methylation writers, readers, and erasers. Mol Cell (2019) 74(4):640–50. doi: 10.1016/j.molcel.2019.04.025

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Yang X, Yang Y, Sun BF, Chen YS, Xu JW, Lai WY, et al. 5-methylcytosine promotes mrna export - nsun2 as the methyltransferase and alyref as an m(5)c reader. Cell Res (2017) 27(5):606–25. doi: 10.1038/cr.2017.55

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Dominissini D, Nachtergaele S, Moshitch-Moshkovitz S, Peer E, Kol N, Ben-Haim MS, et al. The dynamic n(1)-methyladenosine methylome in eukaryotic messenger rna. Nature (2016) 530(7591):441–6. doi: 10.1038/nature16998

PubMed Abstract | CrossRef Full Text | Google Scholar

18. He L, Li H, Wu A, Peng Y, Shu G, Yin G. Functions of n6-methyladenosine and its role in cancer. Mol Cancer (2019) 18(1):176. doi: 10.1186/s12943-019-1109-9

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Lin S, Liu Q, Lelyveld VS, Choe J, Szostak JW, Gregory RI. Mettl1/Wdr4-mediated M(7)G trna methylome is required for normal mrna translation and embryonic stem cell self-renewal and differentiation. Mol Cell (2018) 71(2):244–55 e5. doi: 10.1016/j.molcel.2018.06.001

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Sendinc E, Valle-Garcia D, Dhall A, Chen H, Henriques T, Navarrete-Perea J, et al. Pcif1 catalyzes m6am mrna methylation to regulate gene expression. Mol Cell (2019) 75(3):620–30 e9. doi: 10.1016/j.molcel.2019.05.030

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Mauer J, Luo X, Blanjoie A, Jiao X, Grozhik AV, Patil DP, et al. Reversible methylation of m(6)am in the 5' cap controls mrna stability. Nature (2017) 541(7637):371–5. doi: 10.1038/nature21022

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Nombela P, Miguel-Lopez B, Blanco S. The role of m(6)a, m(5)c and psi rna modifications in cancer: Novel therapeutic opportunities. Mol Cancer (2021) 20(1):18. doi: 10.1186/s12943-020-01263-w

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Sun HL, Zhu AC, Gao Y, Terajima H, Fei Q, Liu S, et al. Stabilization of erk-phosphorylated mettl3 by usp5 increases m(6)a methylation. Mol Cell (2020) 80(4):633–47 e7. doi: 10.1016/j.molcel.2020.10.026

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Gu Y, Wu X, Zhang J, Fang Y, Pan Y, Shu Y, et al. The evolving landscape of n(6)-methyladenosine modification in the tumor microenvironment. Mol Ther (2021) 29(5):1703–15. doi: 10.1016/j.ymthe.2021.04.009

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Zhang B, Wu Q, Li B, Wang D, Wang L, Zhou YL. M(6)a regulator-mediated methylation modification patterns and tumor microenvironment infiltration characterization in gastric cancer. Mol Cancer (2020) 19(1):53. doi: 10.1186/s12943-020-01170-0

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Pan J, Huang Z, Xu Y. M5c rna methylation regulators predict prognosis and regulate the immune microenvironment in lung squamous cell carcinoma. Front Oncol (2021) 11:657466. doi: 10.3389/fonc.2021.657466

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Calvisi DF, Pinna F, Ladu S, Pellegrino R, Simile MM, Frau M, et al. Forkhead box m1b is a determinant of rat susceptibility to hepatocarcinogenesis and sustains erk activity in human hcc. Gut (2009) 58(5):679–87. doi: 10.1136/gut.2008.152652

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Ren BJ, Zhou ZW, Zhu DJ, Ju YL, Wu JH, Ouyang MZ, et al. Alisertib induces cell cycle arrest, apoptosis, autophagy and suppresses emt in ht29 and caco-2 cells. Int J Mol Sci (2015) 17(1). doi: 10.3390/ijms17010041

PubMed Abstract | CrossRef Full Text | Google Scholar

29. O'Shaughnessy J, McIntyre K, Wilks S, Ma L, Block M, Andorsky D, et al. Efficacy and safety of weekly paclitaxel with or without oral alisertib in patients with metastatic breast cancer: a randomized clinical trial. JAMA Netw Open (2021) 4(4):e214103. doi: 10.1001/jamanetworkopen.2021.4103

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Chen H, Yao J, Bao R, Dong Y, Zhang T, Du Y, et al. Cross-talk of four types of rna modification writers defines tumor microenvironment and pharmacogenomic landscape in colorectal cancer. Mol Cancer (2021) 20(1):29. doi: 10.1186/s12943-021-01322-w

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zhang C, Jia G. Reversible rna modification n(1)-methyladenosine (m(1)a) in mrna and trna. Genomics Proteomics Bioinf (2018) 16(3):155–61. doi: 10.1016/j.gpb.2018.03.003

CrossRef Full Text | Google Scholar

32. Yang W, Soares J, Greninger P, Edelman EJ, Lightfoot H, Forbes S, et al. Genomics of drug sensitivity in cancer (gdsc): a resource for therapeutic biomarker discovery in cancer cells. Nucleic Acids Res (2013) 41(Database issue):D955–61. doi: 10.1093/nar/gks1111

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Danni X, Jiangzheng Z, Huamao S, Yinglian P, Changcheng Y, Yanda L. Chaperonin containing tcp1 subunit 3 (cct3) promotes cisplatin resistance of lung adenocarcinoma cells through targeting the janus kinase 2/signal transducers and activators of transcription 3 (jak2/stat3) pathway. Bioengineered (2021) 12(1):7335–47. doi: 10.1080/21655979.2021.1971030

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Chen X, Chen X, Huang Y, Lin J, Wu Y, Chen Y. Tcp1 increases drug resistance in acute myeloid leukemia by suppressing autophagy via activating akt/mtor signaling. Cell Death Dis (2021) 12(11):1058. doi: 10.1038/s41419-021-04336-w

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Nieto MA. Epithelial-mesenchymal transitions in development and disease: old views and new perspectives. Int J Dev Biol (2009) 53(8-10):1541–7. doi: 10.1387/ijdb.072410mn

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Xiao P, Long X, Zhang L, Ye Y, Guo J, Liu P, et al. Neurotensin/Il-8 pathway orchestrates local inflammatory response and tumor invasion by inducing m2 polarization of tumor-associated macrophages and epithelial-mesenchymal transition of hepatocellular carcinoma cells. OncoImmunology (2018) 7(7). doi: 10.1080/2162402x.2018.1440166

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Su S, Liao J, Liu J, Huang D, He C, Chen F, et al. Blocking the recruitment of naive cd4(+) t cells reverses immunosuppression in breast cancer. Cell Res (2017) 27(4):461–82. doi: 10.1038/cr.2017.34

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Schoenfeld AJ, Hellmann MD. Acquired resistance to immune checkpoint inhibitors. Cancer Cell (2020) 37(4):443–55. doi: 10.1016/j.ccell.2020.03.017

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Bach EA, Aguet M, Schreiber RD. The ifn gamma receptor: a paradigm for cytokine receptor signaling. Annu Rev Immunol (1997) 15:563–91. doi: 10.1146/annurev.immunol.15.1.563

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Liu J, Zhou Z, Ma L, Li C, Lin Y, Yu T, et al. Effects of rna methylation n6-methyladenosine regulators on malignant progression and prognosis of melanoma. Cancer Cell Int (2021) 21(1):453. doi: 10.1186/s12935-021-02163-9

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Zhang Y, Zeng F, Zeng M, Han X, Cai L, Zhang J, et al. Identification and characterization of alcohol-related hepatocellular carcinoma prognostic subtypes based on an integrative n6-methyladenosine methylation model. Int J Biol Sci (2021) 17(13):3554–72. doi: 10.7150/ijbs.62168

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Zhang Z, Zhang C, Luo Y, Wu P, Zhang G, Zeng Q, et al. M(6)a regulator expression profile predicts the prognosis, benefit of adjuvant chemotherapy, and response to anti-pd-1 immunotherapy in patients with small-cell lung cancer. BMC Med (2021) 19(1):284. doi: 10.1186/s12916-021-02148-5

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Iqbal Z, Absar M, Akhtar T, Aleem A, Jameel A, Basit S, et al. Integrated genomic analysis identifies ankrd36 gene as a novel and common biomarker of disease progression in chronic myeloid leukemia. Biol (Basel) (2021) 10(11). doi: 10.3390/biology10111182

CrossRef Full Text | Google Scholar

44. Tang N, Cai X, Peng L, Liu H, Chen Y. Tcp1 regulates wnt7b/beta-catenin pathway through p53 to influence the proliferation and migration of hepatocellular carcinoma cells. Signal Transduct Target Ther (2020) 5(1):169. doi: 10.1038/s41392-020-00278-5

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Hu Q, Hong Y, Qi P, Lu G, Mai X, Xu S, et al. Atlas of breast cancer infiltrated b-lymphocytes revealed by paired single-cell rna-sequencing and antigen receptor profiling. Nat Commun (2021) 12(1):2186. doi: 10.1038/s41467-021-22300-2

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Garaud S, Buisseret L, Solinas C, Gu-Trantien C, de Wind A, Van den Eynden G, et al. Tumor infiltrating b-cells signal functional humoral immune responses in breast cancer. JCI Insight (2019) 5. doi: 10.1172/jci.insight.129641

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Ashton TM, McKenna WG, Kunz-Schughart LA, Higgins GS. Oxidative phosphorylation as an emerging target in cancer therapy. Clin Cancer Res (2018) 24(11):2482–90. doi: 10.1158/1078-0432.CCR-17-3070

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Raninga PV, Lee A, Sinha D, Dong LF, Datta KK, Lu X, et al. Marizomib suppresses triple-negative breast cancer via proteasome and oxidative phosphorylation inhibition. Theranostics (2020) 10(12):5259–75. doi: 10.7150/thno.42705

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Weng H, Feng X, Lan Y, Zheng Z. Tcp1 regulates pi3k/akt/mtor signaling pathway to promote proliferation of ovarian cancer cells. J Ovarian Res (2021) 14(1):82. doi: 10.1186/s13048-021-00832-x

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Guest ST, Kratche ZR, Bollig-Fischer A, Haddad R, Ethier SP. Two members of the tric chaperonin complex, cct2 and tcp1 are essential for survival of breast cancer cells and are linked to driving oncogenes. Exp Cell Res (2015) 332(2):223–35. doi: 10.1016/j.yexcr.2015.02.005

PubMed Abstract | CrossRef Full Text | Google Scholar