m6A RNA Methylation Regulator YTHDF1 Correlated With Immune Microenvironment Predicts Clinical Outcomes and Therapeutic Efficacy in Breast Cancer

Objective: Increasing evidence highlights the roles of N⁶-methyladenosine (m⁶A) and its regulators in oncogenesis. Herein, this study observed the associations of m⁶A regulators with breast cancer.

Methods: RNA-seq profiles of breast cancer were retrieved from the Cancer Genome Atlas (TCGA) database. The expression of m⁶A regulators was analyzed in tumor and normal tissues. Their expression correlations were analyzed by Spearson test. Overall survival (OS) analysis of these regulators was then presented. Gene set enrichment analysis (GSEA) was performed in high and low YTHDF1 expression groups. The correlations of YTHDF1 expression with immune cells and tumor mutation burden (TMB) were calculated in breast cancer samples. Somatic variation was assessed in high and low YTHDF1 expression groups.

Results: Most of m⁶A regulators were abnormally expressed in breast cancer compared to normal tissues. At the mRNA levels, there were closely relationships between them. Among them, YTHDF1 up-regulation was significantly related to undesirable prognosis (p = 0.025). GSEA results showed that high YTHDF1 expression was associated with cancer-related pathways. Furthermore, YTHDF1 expression was significantly correlated with T cells CD4 memory activated, NK cells activated, monocytes, and macrophages. There were higher TMB scores in YTHDF1 up-regulation group than its down-regulation group. Missense mutation and non-sense mutation were the most frequent mutation types.

Conclusion: Our findings suggested that dysregulated m⁶A regulator YTHDF1 was predictive of survival outcomes as well as response to immunotherapy of breast cancer, and were closely related to immune microenvironment.

Introduction

The incidence of breast cancer continues to rise globally (1). It represents the most mortal among the female population (2). This malignancy is heterogeneous on clinical, molecular behaviors as well as response to therapies (3). This management is multidisciplinary, including locoregional (surgery or radiotherapy) as well as systemic therapies (4). At present, advanced breast cancer with distant metastasis is incurable. Bone, lung, liver, and brain are common metastatic sites (5). Individualized therapy is future therapeutic goal for breast cancer. Thus, it is vital to elucidate the mechanisms of breast cancer initiation as well as progression.

m⁶A is the most abundant modification in eukaryotic mRNAs, occupying 0.1–0.4% of the total adenine residues (6). The m⁶A modification takes on varied biological functions (7) like RNA splicing, RNA stabilities, nuclear export as well as translation (8). This process is involved in three kinds of m⁶A regulators, called as “writers,” “erasers,” and “readers,” containing “writers” (METTL3, METTL14, WTAP, RBM15, KIAA1429, ZC3H13), “readers” (YTHDC1, YTHDC2, YTHDF1, YTHDF2, HNRNP), as well as “erasers” (FTO, ALKBH5). “Writers” are responsible for catalyzing the formation of m⁶A. “Readers” are charge of decoding m⁶A methylation as well as producing functional signals. “Erasers” can remove the methyl code from targeted mRNAs. The m⁶A modification participates in carcinogenesis through regulating RNA production as well as metabolism. For instance, Niu et al. (9) demonstrated that FTO promoted breast tumor progress by inhibition of BNIP3. As in the study of Cai et al. (10), METTL3 up-regulation accelerates breast cancer cellular proliferation. Recent studies have emphasized the key implications of immune microenvironment in breast cancer progression as well as response to immunotherapy (11). Zhang et al. (12) reported that m⁶A modification was involved in tumor immune microenvironment formation. ALKBH5 m⁶A reader could regulate the response to anti-PD-1 immunotherapy through inhibiting the accumulation of tumor-infiltrating immune cells (13). Hence, evaluation of m⁶A regulators in individual breast cancer may enhance our understanding about characteristics of tumor immune microenvironment as well as improve the therapeutic efficacy of immunotherapies. In this study, we evaluated the expression patterns of m⁶A regulators and their correlations with survival outcomes, immune microenvironment, response to immunotherapy as well as somatic mutation in breast cancer.

Materials and Methods

Data Sourcse

From TCGA database (http://cancergenome.nih.gov/), this study collected RNA-seq profiles of breast cancer subjects on January 15, 2021. Meanwhile, the matched clinical information was also retrieved. After removing samples with survival time of 0, 902 samples including Basal, Her2, LumA, and LumB subtypes were retained for further analysis (Table 1). The microarray expression profiling of 17 normal breast tissues and 104 breast cancer tissues were retrieved from the GSE42568 dataset of the Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/gds/) repository. Also, expression profiles and follow-up information of 266 breast cancer patients were retrieved from the GSE21653 dataset.

TABLE 1

Table 1. Clinicopathological characteristics of breast cancer patients in this study.

m⁶A Regulators

Totally, this study gathered 13 m⁶A regulators including five writers (METTL3, METTL14, WTAP, RBM15, ZC3H13, KIAA1429), five readers (YTHDC1, YTHDC2, YTHDF1, YTHDF2, HNRNPC), two erasers (FTO, ALKBH5). The expressions of these m⁶A regulators were assessed in breast cancer and normal tissue samples by Wilcoxon test. P-value was corrected with Bonferroni method. Spearson correlation analysis was presented between these regulators at the mRNA levels. |r| > 0.5 indicated a significant correlation.

Survival Analysis

Samples with survival status of zero were removed. Overall survival (OS) is defined as the time from the date of diagnosis to death due to any cause. The patients were classified into high and low expression of each m⁶A regulator groups based on the median value of its expression. OS analyses were carried out between groups through univariate cox regression analysis. For each regulator, this study calculated p-value, hazard ratio (HR) as well as 95% confidence interval (CI). Regulators with p-value <0.05 and HR >1 were risk factors of breast cancer prognosis. Meanwhile, those with p-value <0.05 and HR <1 were protective factors.

GSEA

Enrichment analysis was presented with GSEA software (14). Here, breast cancer subjects were separated into high and low YTHDF1 expression groups based on the media value of YTHDF1. Thousand gene set permutations were presented for each analysis. The expression level of YTHDF1 was considered a phenotype. According to the normalized enrichment score (NES), normalized p-value and false discovery rate (FDR), the enrichment pathways of each phenotype were classified. The absolute value of NES ≥1.0, normalized p-value ≤ 0.05 and FDR ≤ 0.25 were confirmed as meaningful gene sets.

Cell Culture and Transfection

Human breast cancer cells MCF-7 (ATCC, USA) were grown in DMEM (Thermo, USA) containing 10% FBS in an incubator with 5% CO₂ at 37°C. MCF-7 cells were seeded in 6-well plates. The culture medium was changed 1 day before transfection, and when the cells reached 70–90% of the growth density, the cells were transfected with synthetic siRNAs (si-YTHDF1; GenePharma, Suzhou, China) through the Lipofectamine^™2000 Transfection Kit (Invitrogen, USA). Simultaneously, non-interfering siRNA was transfected as a negative control. The transfection process strictly followed the instructions of the kit, and the cells transfected for 48 h were collected for next research.

Western Blot

Transfected cells were lysed on ice with cell lysate to extract total protein. The BCA method was used for protein quantification. The proteins were separated by SDS-PAGE gel electrophoresis, and transferred to PVDF membranes at a constant voltage of 80 V. After blocking with 5% skimmed milk powder TBST at room temperature for 2 h, the membranes were incubated with YTHDF1 (1/500; ab230330; Abcam), Axin2 (1:1000; ab32197; Abcam), c-myc (1:5000; ab152146; Abcam), β-catenin (1:10000; ab81305; Abcam), cyclin D1 (1:10000; ab134175; Abcam), and β-actin (1:5000; ab179467; Abcam) overnight at 4°C. After washing, the membranes were incubated with secondary antibodies at room temperature for 1 h. The color was developed by chemiluminescence, and the gel imaging system was used to analyze images. The gray value of the bands was measured.

Somatic Variation

Somatic variant data of breast cancer that were stored in the mutation annotation format (MAF) were obtained from TCGA database. According to VarScan2 variant aggregation as well as masking data, somatic variation analyses were presented through “maftools” package (15). TMB score was determined for each patient as follows: (total mutation/total covered bases).

CIBERSORT

The CIBERSORT algorithm was employed for estimating the fractions of 22 phenotypes of immune cells in each specimen based on gene expression profiles (16). LM22 leukocyte gene signature matrix was utilized in conjunction. Specimens with p-value < 0.05 were retained. For each specimen, the sum of the estimated fractions of immune cells was equal to 1. The enriched scores of each immune cell were compared in high and low YTHDF1 expression groups by Wilcoxon test.

TMB

TMB is defined as the total number of somatic mutations per Mb base in the coding region of an exon, which is an emerging biomarker for judging the efficacy of tumor immunotherapy (17). The greater the TMB score, the better the therapy efficacy. TMB is calculated as the total number of somatic mutations/the size of the target area, and the unit is mutations/Mb. The somatic mutation data were in MAF format. In this study, the somatic mutation data processed by vanscan software were downloaded from TCGA. The “maftools” package was applied to calculate the TMB score of each sample. The difference in TMB scores between high and low YTHDF1 expression groups was assessed by Wilcoxon test.

Differential Expression Analysis

Differentially expressed genes (DEGs) were screened between high and low YTHDF1 expression groups via the limma package. The screening thresholds were as follows: |log2 fold change (FC)| > 1 and adjusted p-value <0.05.

Functional Enrichment Analysis

The “clusterProfiler” package was applied to present Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enrichment analysis based on YTHDF1-related DEGs (18). Terms with adjusted p-value < 0.05 were significantly enriched.

Connectivity Map

CMap database (https://clue.io/) was employed for exploring candidate chemical compounds against breast cancer (19). Based on a list of DEGs, this study searched for the compounds through this database. The CMap connectivity score (range: −1 to 1) was indicative of the specificity degree related to the DEGs. The connectivity score of a compound tended to −1, suggesting that it negatively correlated with the DEGs. On the contrary, the connectivity score of a compound closer to 1 implied that it exhibited positive correlations with the DEGs. Here, the compounds with |connectivity score| ≥90 were candidate chemical agents.

Scratch Test

The cells were seeded in a 6-well plate and placed in a CO₂ cell incubator. After the cells were overgrown, a sterile toothpick was used to make a vertical mark of uniform thickness in each hole. After taking pictures of the initial state of the scratches at 0 h under an inverted microscope, the samples were put back into the incubator. After 24 h, the scratch state was photographed again. ImageJ software was used to calculate the scratch area. Relative migration rate = (0 h scratch area-24 h scratch area)/0 h scratch area.

Transwell

The cells were added to the Matrigel-coated Transwell chamber (1 × 10⁵ cells/well). A medium containing 10% FBS was added to the well in the lower layer of each chamber. After culturing for 24 h, the chamber was removed. After staining with crystal violet, a cotton swab was used to gently wipe away the non-invasive cells in the upper chamber. The invasive cells were observed under an inverted microscope. Then, ImageJ software was applied to calculate the number of invasive cells.

Statistical Analysis

All statistical analysis was carried out through the packages of R software (http://www.r-project.org/) or Graphpad Prism. Student's t-test, Wilcoxon test or ANOVA test was applied for comparisons between groups. P < 0.05 indicated statistical significance.

Results

Expression Patterns and Correlations of 13 m⁶A Regulators in Breast Cancer

Genomic locations for 13 m⁶A regulatory genes were displayed in Figure 1A, as follows: METTL3 (chr14: 21498133-21511342), METTL14 (chr4: 118685392-118715433), WTAP (chr6: 159725585-159756319), KIAA1429 (chr8: 94487689-94553529), RBM15 (chr1: 110338506-110346681), ZC3H13 (chr13: 45954465-46052759), YTHDC1 (chr4: 68310387-68350090), YTHDC2 (chr5: 113513694-113595285), YTHDF1 (chr20: 63195429-63216139), YTHDF2 (chr1: 28736621-28769775), HNRNPC (chr14: 21209136-21269494), FTO (chr16: 53701692-54158512) and ALKBH5 (chr17: 18183078-18209954). The expression of these regulators was compared in breast cancer and normal tissues. We found that METTL14 (adjusted p-value = 9.10e-08), WTAP (adjusted p-value = 2.86e-07), KIAA1429 (adjusted p-value = 2.33e-06), RBM15 (adjusted p-value = 1.52e-03), ZC3H13 (adjusted p-value = 5.56e-12), YTHDC1 (adjusted p-value = 1.68e-03), YTHDF1 (adjusted p-value = 4.72e-27), YTHDF2 (adjusted p-value = 4.72e-02), HNRNPC (adjusted p-value = 1.60e-24), and FTO (adjusted p-value = 2.62e-33) were significantly abnormally expressed in breast cancer compared to normal tissues (Figure 1B). However, no significant differences in METTL3, YTHDC2, and ALKBH5 expression were found between tumor and normal samples. The correlations between 13 m⁶A regulators were assessed at the mRNA levels among breast cancer specimens. Our data showed that METTL14 had a strong correlation to YTHDC1 (r = 0.7) and YTHDC2 (r = 0.65), as displayed in Figure 1C. Meanwhile, YTHDC1 was strongly associated with YTHDF1 (r = 0.64), and moderately correlated to HNRNPC (r = 0.58).

FIGURE 1

Figure 1. Expression patterns and correlations of 13 m⁶A regulators in breast cancer from TCGA database. (A) Genomic locations of the regulators. (B) Bar diagram for the expression of regulators in breast cancer and normal tissues. (C) Assessment of correlations between 13 m⁶A regulators at the mRNA levels. Ns, not significant; ^*p < 0.05; ^**p < 0.01; ^****p < 0.0001.

YTHDF1 Expression Is Correlated to Breast Cancer Patients' Survival

To determine the clinical implications of 13 m⁶A regulatory mRNAs in breast cancer, this study observed the correlations between the expression of 13 m⁶A regulatory mRNAs and patients' clinical outcomes. Our results showed that the expression of ALKBH5, FTO, HNRNPC, METTL3, METTL14, RBM15, WTAP, YTHDC1, YTHDF2, ZC3H13, KIAA1429, and YTHDC2 were all not associated with patients' survival (Figures 2A–L). Only YTHDF1 expression exhibited a significant correlation to subjects' prognosis (Figure 2M). Its expression was a risk factor for breast cancer (p = 0.025; HR: 1.5; 95% CI: 1.03–2.19). Subjects with high YTHDF1 expression often experienced undesirable survival outcomes than those with its low expression.

FIGURE 2

Figure 2. Correlation between 13 m⁶A regulators and breast cancer patients' survival. OS of (A) ALKBH5; (B) FTO; (C) HNRNPC; (D) METTL3; (E) METTL14; (F) RBM15; (G) WTAP; (H) YTHDC1; (I) YTHDF2; (J) ZC3H13; (K) KIAA1429; (L) YTHDC2; (M) YTHDF1 expression for breast cancer patients.

Verification of Expression and Prognosis of YTHDF1 in Breast Cancer

YTHDF1 up-regulation was confirmed in breast cancer in the GSE42568 dataset (p = 5e-09; Figure 3A). Also, its up-regulation was in relation to unfavorable survival outcomes of subjects in the GSE21653 dataset (p = 0.031; Figure 3B).

FIGURE 3

Figure 3. Verification of expression and prognosis of YTHDF1 in breast cancer. (A) Box plots for YTHDF1 expression in the GSE42568 dataset. (B) Survival analysis of breast cancer with high and low YTHDF1 expression in the GSE21653 dataset.

Enriched Pathways in High and Low YTHDF1 Expression Groups

Then, we evaluated the enriched signaling pathways in high and low YTHDF1 expression groups. We found that cell cycle, ERBB signaling pathway, oocyte meiosis, pathways in cancer, spliceosome, ubiquitin mediated proteolysis, and WNT signaling pathway were distinctly enriched in high YTHDF1 expression group (Figure 4A). Meanwhile, ribosome was significantly enriched in low YTHDF1 expression group (Figure 4B). To verify whether YTHDF1 altered WNT pathway activation, YTHDF1 was successfully silenced by siRNA in MCF-7 cells (Figures 4C,D). Our data showed that YTHDF1 knockdown distinctly lowered the expression of Axin2, c-myc, β-catenin, and cyclin D1 in MCF-7 cells, confirming that YTHDF1 participated in the activation of WNT pathway in breast cancer (Figures 4E–H).

FIGURE 4

Figure 4. Enriched pathways in high and low YTHDF1 expression groups. (A) Enriched pathways in high expression group. (B) An enriched pathway in low expression group. (C) Western blot for detecting the expressions of (D) YTHDF1, (E) Axin2, (F) c-myc, (G) β-catenin and (H) cyclin D1 in MCF-7 cells transfected with si-YTHDF1. Ns, not significant; ^**p < 0.01; ^***p < 0.001.

Correlation Between YTHDF1 Expression and Tumor Immune Microenvironment and Response to Immunotherapy

m⁶A modification is closely related to tumor microenvironment cell infiltration in individual tumors (12). Here, we assessed the correlation between YTHDF1 expression and tumor immune microenvironment in breast cancer. Samples with high YTHDF1 expression distinctly exhibited higher infiltration scores of T cells CD4 memory activated and macrophages M1 those with its low expression (Figure 5A; Table 2). Lower infiltration levels of NK cells activated and monocytes were found in subjects with high YTHDF1 expression compared to those with its low expression. Despite the revolutionization of immune checkpoint blockade (ICB) therapy, most patients cannot benefit from ICB therapy (13). We found that high YTHDF1 expression was significantly correlated to higher TMB score, indicating that patients with its up-regulation had a better effect on immunotherapy (Figure 5B). Thus, YTHDF1 expression might be used for predicting the response to immunotherapy.

FIGURE 5

Figure 5. Correlation between YTHDF1 expression and tumor immune microenvironment and response to immunotherapy. (A) Box plots for the association between YTHDF1 expression and infiltration levels of immune cells in breast cancer. (B) Box plots for the association between YTHDF1 expression and TMB score.

TABLE 2

Table 2. The correlations between YTHDF1 expression and tumor-infiltrating immune cells in breast cancer.

Somatic Mutations in High and Low YTHDF1 Breast Cancer

Somatic mutations were evaluated in high and low YTHDF1 breast cancer samples from TCGA database. Among 986 samples, 260 (26.37%) occurred somatic mutations in high YTHDF1 expression group. Here, we displayed the top 20 genes according to mutation frequency. As shown in Figure 6A, TP53 (12%) had the most frequently mutated gene, followed by PIK3CA (10%), TTN (5%), CDH1 (3%), and GATA3 (3%). Four hundred and forty samples occurred somatic mutations in low YTHDF1 expression group (Figure 6B). Consistent with high expression group, TP53 (16%), PIK3CA (17%), TTN (9%), CDH1 (8%), and GATA3 (7%) were the top five mutated genes. Both in high and low YTHDF1 expression groups, missense mutation was the most common mutation type.

FIGURE 6

Figure 6. Somatic mutations in high and low YTHDF1 breast cancer. (A) The top 20 genes according to mutation frequency in high YTHDF1 expression groups. (B) The top 20 genes according to mutation frequency in low expression groups. Each mutation type is identified by a unique color.

Expression Patterns of YTHDF1 in Pan-Cancer

We comprehensively analyzed the expression of YTHDF1 in pan-cancer and corresponding normal tissues. Up-regulation of YTHDF1 was found in bladder cancer, breast cancer, cholangiocarcinoma, colon adenocarcinoma, esophageal carcinoma, glioblastoma multiforme, head and neck squamous cell carcinoma, kidney renal papillary cell carcinoma, liver hepatocellular carcinoma, lung adenocarcinoma, lung squamous cell carcinoma, pheochromocytoma and paraganglioma, prostate adenocarcinoma, rectum adenocarcinoma, stomach adenocarcinoma, and uterine corpus endometrial carcinoma compared to corresponding normal tissues (Figure 7). However, YTHDF1 was lowly expressed in thyroid carcinoma than normal samples.

FIGURE 7

Figure 7. Expression patterns of YTHDF1 in pan-cancer and matched normal tissues. Red indicates tumor specimens and blue indicates normal tissues. Each dot represents one sample. ^*P < 0.05; ^**p < 0.01; ^***p < 0.001.

Enrichment Analysis of YTHDF1-Related DEGs in Breast Cancer

Sixteen up-regulated and down-regulated genes were identified between high and low YTHDF1 expression breast cancer samples (Table 3). Their biological functions were then discovered through GO and KEGG enrichment analyses. We found that up-regulated genes were only enriched in columnar/cuboidal epithelial cell differentiation (Figure 8A). Down-regulated genes were significantly involved in RNA metabolism-related biological processes, such as cAMP-mediated signaling, cyclic-nucleotide-mediated signaling, and adenylate cyclase-modulating G protein-coupled receptor signaling pathway (Figure 8B). Furthermore, they had the molecular functions of G protein-coupled receptor binding and neuropeptide receptor binding. KEGG pathway enrichment analysis revealed that up-regulated genes significantly participated in GABAergic synapse (Figure 8C). In Figure 8D, neuroactive ligand-receptor interaction and renin secretion were distinctly enriched by down-regulated genes.

TABLE 3

Table 3. DEGs between high and low YTHDF1 expression in breast cancer.

FIGURE 8

Figure 8. Enrichment analysis of YTHDF1-related DEGs in breast cancer. GO enrichment results by (A) up-regulated and (B) down-regulated genes. KEGG pathway enrichment results by (C) up-regulated, and (D) down-regulated genes.

Candidate Therapeutic Agents Against Breast Cancer Based on YTHDF1-Related DEGs

Candidate therapeutic agents were discovered based on YTHDF1-related DEGs (Table 4). We found that indoprofen, nabumetone, nimesulide, and phenacetin shared the MoA of Cyclooxygenase inhibitor. Digitoxigenin, helveticoside, ouabain shared the MoA of ATPase inhibitor. Alclometasone, mometasone, and piretanide shared the MoA of Glucocorticoid receptor agonist (Figure 9). These compounds could become candidate therapeutic agents against breast cancer.

TABLE 4

Table 4. Candidate therapeutic agents against breast cancer.

FIGURE 9

Figure 9. CMap identifies the candidate therapeutic agents related to the DEGs between high and low YTHDF1 expression breast cancer samples. Heatmap shows each inhibitor (perturbagen) and its shared mechanism of action (row).

YTHDF1 Knockdown Restrains Migrated and Invasive Abilities of Breast Cancer

Following YTHDF1 knockdown, migrated as well as invasive abilities of breast cancer were investigated in breast cancer. The scratch test demonstrated that silencing YTHDF1 markedly lowered the migrated levels of MCF-7 cells (Figures 10A,B). Moreover, the number of invasive cells was decreased by YTHDF1 knockdown (Figures 10C,D).

FIGURE 10

Figure 10. YTHDF1 knockdown restrains migrated and invasive abilities of in MCF-7 cells transfected with si-YTHDF1. (A,B) Scratch test for the relative migrated levels of transfected MCF-7 cells. (C,D) Transwell for the number of invasive MCF-7 cells. ^*P < 0.05; ^**p < 0.01; ^***p < 0.001.

Discussion

In this study, we characterized the expression patterns of m⁶A regulators and their implications on survival outcomes, immune microenvironment, response to immunotherapy as well as somatic mutations in breast cancer. Our data highlighted the key roles of m⁶A modification and their regulators in tumor progression and prognosis.

We found that most of m⁶A regulators including METTL14, WTAP, KIAA1429, RBM15, ZC3H13, YTHDC1, YTHDF1, YTHDF2, HNRNPC, and FTO were dysregulated in breast cancer than normal tissues. METTL14 displayed a strong correlation to YTHDC1 and YTHDC2 while YTHDC1 was strongly correlated to YTHDF1 and moderately correlated to HNRNPC. Among them, m⁶A reader YTHDF1 aberration is correlated to undesirable survival outcomes in breast cancer subjects, which exhibited the consistency with the research from Anita et al. (20). The carcinogenic roles of YTHDF1 have been confirmed in previous research. For example, Liu et al. (21) reported that YTHDF1 facilitated ovarian carcinoma progress through controlling EIF3C translation. Bai et al. (22) reported that YTHDF1 accelerated tumorigenicity in colorectal carcinoma. Shi et al. (23) found that YTHDF1 correlated to hypoxia adaptation may contribute to non-small cell lung cancer development. Pi et al. (24) also demonstrated that YTHDF1 promoted gastric carcinogenesis through elevating translation of FZD7. Our pan-cancer analysis revealed that YTHDF1 was up-regulated in most types of cancer. Thus, YTHDF1 could be an oncogene. Our experiments confirmed that silencing YTHDF1 suppressed migrated and invasive capacities of breast cancer cells. The GSEA results demonstrated that high YTHDF1 expression was distinctly correlated to cell cycle, ERBB signaling pathway, oocyte meiosis, pathways in cancer, spliceosome, ubiquitin mediated proteolysis, and WNT signaling pathway. Meanwhile, ribosome was significantly enriched in low YTHDF1 expression group. It has been reported that YTHDF1 could regulate cell cycle progression in hepatocellular carcinoma (25). Recent research has reported that YTHDF1 mediates Wnt pathway activation in intestinal stemness (26). These findings were indicative that YTHDF1 up-regulation could participate in carcinogenesis.

Tumor immune microenvironment affects initiation as well as progress in breast cancer (27). Tumor-infiltrating immune cells are correlated to survival outcomes. Here, we found that YTHDF1 expression was distinctly related to T cells CD4 memory activated, macrophages M1, NK cells activated, and monocytes in breast cancer tissues. Vaccines against dendritic cells have exhibited prolonged survival time in breast cancer patients (28). Tumor-associated macrophages may modulate the efficacy of anti-PD-1/PD-L1 therapy in breast cancer (29). Our data were indicative that YTHDF1 might modulate the immune microenvironment of breast cancer, thereby, affecting tumor progression as well as immunotherapy efficacy. For immunotherapy, the higher the TMB of cancer cells, the increased new antigens may be produced. The higher the immunogenicity of the antigen, the stronger the T cell response, and anti-tumor response, which is more suitable for immunotherapy. Immune-checkpoint inhibitors (ICIs) are novel therapeutic strategies against breast cancer. Nevertheless, only some subjects respond to PD-1 or PD-L1 therapy. As widely accepted, breast cancer patients with high TMB can benefit from immunotherapy (30). Here, we found that TMB score was significantly higher in high YTHDF1 expression group compared to its low expression group, indicating that subjects with high YTHDF1 expression were more likely to benefit from immunotherapy.

The occurrence of tumors is the result of the accumulation of somatic mutations (31). In fact, there are basically non-synonymous mutations in the development of tumors. Because the mutation will increase immunogenicity, but in order to avoid being detected and eliminated by the immune system, tumors often increase immune checkpoints (32). The driver gene mutations can lead to tumors, so that a large number of somatic mutations can produce new antigens, which can activate CD8⁺ cytotoxic T cells, thereby, exerting T cell-mediated anti-tumor effects (33). Therefore, when the number of gene mutations accumulates, more new antigens will be produced, which will be more likely to be recognized by the immune system. Among the 986 breast cancer patients, 26.37% samples in high YTHDF1 expression group and 44.62% samples in low expression group occurred genetic mutations, indicating that the frequency of mutations in breast cancer patients was very high. Among them, missense mutation and non-sense mutation were most frequently. Both in high and low expression of YTHDF1 groups, the five most common mutant genes were TP53, PIK3CA, TTN, CDH1, and GATA3, indicating that these mutated genes contributed to the progression of breast cancer.

Except for breast cancer, up-regulation of YTHDF1 was found in various cancers, indicating that YTHDF1 could be an oncogene. To explore underlying molecular mechanisms of YTHDF1 in breast cancer, we screened DEGs between high and low YTHDF1 expression groups. Our results showed that DEGs were mainly involved in RNA metabolism processes, such as cAMP-mediated signaling, cyclic-nucleotide-mediated signaling, adenylate cyclase-modulating G protein-coupled receptor signaling pathway, second-messenger-mediated signaling as well as adenylate cyclase-activating G protein-coupled receptor signaling pathway. These findings indicated that YTHDF1 promoted tumor progression mainly by m⁶A modification. This study also screened several small molecular inhibitors such as cyclooxygenase inhibitor (indoprofen, nabumetone, nimesulide, and phenacetin), ATPase inhibitor (digitoxigenin, helveticoside, ouabain), glucocorticoid receptor agonist (alclometasone, mometasone, and piretanide), which might be candidate therapeutic agents against breast cancer. More experiments should be carried out to investigate the therapeutic effects of these small molecular inhibitors in breast cancer cells.

Conclusion

Collectively, this study characterized the dysregulated expression patterns of m⁶A regulators in breast cancer. Among them, YTHDF1 overexpression was distinctly indicative of undesirable survival outcomes. Moreover, YTHDF1 up-regulation exhibited a significant association with cancer-related pathways such as cell cycle, pathways in cancer and Wnt signaling pathway. YTHDF1 expression was significantly correlated to tumor-infiltrating immune cells, indicating that it might contribute to the complexity and diversity of immune microenvironment. Furthermore, subjects with YTHDF1 up-regulation were more likely to benefit from immunotherapy. Several underlying small molecular compounds against breast cancer were discovered based on YTHDF1-related DEGs. In conclusion, our data suggested the implications of m⁶A regulators in survival outcomes, immune microenvironment as well as response to immunotherapy in breast cancer.

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

JJ and PT conceived and designed the study. YH, QP, and MW conducted most of the experiments and data analysis and wrote the manuscript. XA, YY, YT, and YJ participated in collecting data and helped to draft the manuscript. All authors reviewed and approved the manuscript.

Funding

This work was funded by Chongqing Basic Research and Frontier Exploration Project (cstc2018jcyjAX0379).

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.

Abbreviations

m⁶A, N⁶-methyladenosine; TCGA, The Cancer Genome Atlas; OS, overall survival; GSEA, gene set enrichment analysis; TMB, tumor mutation burden; HR, hazard ratio; CI, confidence interval; NES, normalized enrichment score; FDR, false discovery rate; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; CMap, Connectivity Map.

References

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

PubMed Abstract | CrossRef Full Text | Google Scholar

2. de Bessa Garcia SA, Araújo M, Pereira T, Mouta J, Freitas R. HOX genes function in breast cancer development. Biochim Biophys Acta Rev Cancer. (2020) 1873:188358. doi: 10.1016/j.bbcan.2020.188358

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Wagner J, Rapsomaniki MA, Chevrier S, Anzeneder T, Langwieder C, Dykgers A, et al. A single-cell atlas of the tumor and immune ecosystem of human breast cancer. Cell. (2019) 177:1330–45.e18. doi: 10.1016/j.cell.2019.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Li N, Deng Y, Zhou L, Tian T, Yang S, Wu Y, et al. Global burden of breast cancer and attributable risk factors in 195 countries and territories, from 1990 to 2017: results from the Global Burden of Disease Study 2017. J Hematol Oncol. (2019) 12:140. doi: 10.1186/s13045-019-0828-0

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Liang Y, Zhang H, Song X, Yang Q. Metastatic heterogeneity of breast cancer: molecular mechanism and potential therapeutic targets. Semin Cancer Biol. (2020) 60:14–27. doi: 10.1016/j.semcancer.2019.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Chen XY, Zhang J, Zhu JS. The role of m(6)A RNA methylation in human cancer. Mol Cancer. (2019) 18:103. doi: 10.1186/s12943-019-1033-z

CrossRef Full Text | Google Scholar

7. Li Y, Gu J, Xu F, Zhu Q, Chen Y, Ge D, et al. Molecular characterization, biological function, tumor microenvironment association and clinical significance of m6A regulators in lung adenocarcinoma. Brief Bioinform. (2020). doi: 10.1093/bib/bbaa225

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Niu X, Xu J, Liu J, Chen L, Qiao X, Zhong M. Landscape of N(6)-methyladenosine modification patterns in human ameloblastoma. Front Oncol. (2020) 10:556497. doi: 10.3389/fonc.2020.556497

CrossRef Full Text | Google Scholar

9. Niu Y, Lin Z, Wan A, Chen H, Liang H, Sun L, et al. RNA N6-methyladenosine demethylase FTO promotes breast tumor progression through inhibiting BNIP3. Mol Cancer. (2019) 18:46. doi: 10.1186/s12943-019-1004-4

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Cai X, Wang X, Cao C, Gao Y, Zhang S, Yang Z, et al. HBXIP-elevated methyltransferase METTL3 promotes the progression of breast cancer via inhibiting tumor suppressor let-7g. Cancer Lett. (2018) 415:11–9. doi: 10.1016/j.canlet.2017.11.018

PubMed Abstract | CrossRef Full Text | Google Scholar

11. 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:1293–308.e36. doi: 10.1016/j.cell.2018.05.060

PubMed Abstract | CrossRef Full Text | Google Scholar

12. 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:53. doi: 10.1186/s12943-020-01170-0

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Li N, Kang Y, Wang L, Huff S, Tang R, Hui H, et al. ALKBH5 regulates anti-PD-1 therapy response by modulating lactate and suppressive immune cell accumulation in tumor microenvironment. Proc Natl Acad Sci U S A. (2020) 117:20159–70. doi: 10.1073/pnas.1918986117

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. (2005) 102:15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Mayakonda A, Lin DC, Assenov Y, Plass C, Koeffler HP. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res. (2018) 28:1747–56. doi: 10.1101/gr.239244.118

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Newman AM, Liu CL, Green MR, Gentles AJ, Feng W, Xu Y, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. (2015) 12:453–7. doi: 10.1038/nmeth.3337

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Ravaioli S, Limarzi F, Tumedei MM, Palleschi M, Maltoni R, Bravaccini S. Are we ready to use TMB in breast cancer clinical practice? Cancer Immunol Immunother. (2020) 69:1943–5. doi: 10.1007/s00262-020-02682-w

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Subramanian A, Narayan R, Corsello SM, Peck DD, Natoli TE, Lu X, et al. A next generation connectivity map: L1000 platform and the first 1,000,000 profiles. Cell. (2017) 171:1437–52.e17. doi: 10.1016/j.cell.2017.10.049

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Anita R, Paramasivam A, Priyadharsini JV, Chitra S. The m6A readers YTHDF1 and YTHDF3 aberrations associated with metastasis and predict poor prognosis in breast cancer patients. Am J Cancer Res. (2020) 10:2546–54.

PubMed Abstract | Google Scholar

21. Liu T, Wei Q, Jin J, Luo Q, Liu Y, Yang Y, et al. The m6A reader YTHDF1 promotes ovarian cancer progression via augmenting EIF3C translation. Nucleic Acids Res. (2020) 48:3816–31. doi: 10.1093/nar/gkaa048

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Bai Y, Yang C, Wu R, Huang L, Song S, Li W, et al. YTHDF1 regulates tumorigenicity and cancer stem cell-like activity in human colorectal carcinoma. Front Oncol. (2019) 9:332. doi: 10.3389/fonc.2019.00332

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Shi Y, Fan S, Wu M, Zuo Z, Li X, Jiang L, et al. YTHDF1 links hypoxia adaptation and non-small cell lung cancer progression. Nat Commun. (2019) 10:4892. doi: 10.1038/s41467-019-12801-6

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Pi J, Wang W, Ji M, Wang X, Wei X, Jin J, et al. YTHDF1 promotes gastric carcinogenesis by controlling translation of FZD7. Cancer Res. (2020) 81:2651–65. doi: 10.1158/0008-5472.CAN-20-0066

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Zhao X, Chen Y, Mao Q, Jiang X, Jiang W, Chen J, et al. Overexpression of YTHDF1 is associated with poor prognosis in patients with hepatocellular carcinoma. Cancer Biomark. (2018) 21:859–68. doi: 10.3233/CBM-170791

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Han B, Yan S, Wei S, Xiang J, Liu K, Chen Z, et al. YTHDF1-mediated translation amplifies Wnt-driven intestinal stemness. EMBO Rep. (2020) 21:e49229. doi: 10.15252/embr.201949229

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Tekpli X, Lien T, Røssevold AH, Nebdal D, Borgen E, Ohnstad HO, et al. An independent poor-prognosis subtype of breast cancer defined by a distinct tumor immune microenvironment. Nat Commun. (2019) 10:5499. doi: 10.1038/s41467-019-13329-5

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Basu A, Ramamoorthi G, Jia Y, Faughn J, Wiener D, Awshah S, et al. Immunotherapy in breast cancer: current status and future directions. Adv Cancer Res. (2019) 143:295–349. doi: 10.1016/bs.acr.2019.03.006

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Santoni M, Romagnoli E, Saladino T, Foghini L, Guarino S, Capponi M, et al. Triple negative breast cancer: key role of tumor-associated macrophages in regulating the activity of anti-PD-1/PD-L1 agents. Biochim Biophys Acta Rev Cancer. (2018) 1869:78–84. doi: 10.1016/j.bbcan.2017.10.007

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Barroso-Sousa R, Jain E, Cohen O, Kim D, Buendia-Buendia J, Winer E, et al. Prevalence and mutational determinants of high tumor mutation burden in breast cancer. Ann Oncol. (2020) 31:387–94. doi: 10.1016/j.annonc.2019.11.010

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Zhang Y, Xiong S, Liu B, Pant V, Celii F, Chau G, et al. Somatic Trp53 mutations differentially drive breast cancer and evolution of metastases. Nat Commun. (2018) 9:3953. doi: 10.1038/s41467-018-06146-9

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Zacharakis N, Chinnasamy H, Black M, Xu H, Lu YC, Zheng Z, et al. Immune recognition of somatic mutations leading to complete durable regression in metastatic breast cancer. Nat Med. (2018) 24:724–30. doi: 10.1038/s41591-018-0040-8

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Kraya AA, Maxwell KN, Wubbenhorst B, Wenz BM, Pluta J, Rech AJ, et al. Genomic signatures predict the immunogenicity of BRCA-deficient breast cancer. Clin Cancer Res. (2019) 25:4363–74. doi: 10.1158/1078-0432.CCR-18-0468

PubMed Abstract | CrossRef Full Text | Google Scholar