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 N6-methyladenosine (m6A) and its regulators in oncogenesis. Herein, this study observed the associations of m6A regulators with breast cancer. Methods: RNA-seq profiles of breast cancer were retrieved from the Cancer Genome Atlas (TCGA) database. The expression of m6A 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 m6A 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 m6A 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 6 A is the most abundant modification in eukaryotic mRNAs, occupying 0.1-0.4% of the total adenine residues (6). The m 6 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 6 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 6 A. "Readers" are charge of decoding m 6 A methylation as well as producing functional signals. "Erasers" can remove the methyl code from targeted mRNAs. The m 6 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 6 A modification was involved in tumor immune microenvironment formation. ALKBH5 m 6 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 6 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 6 A regulators and their correlations with survival outcomes, immune microenvironment, response to immunotherapy as well as somatic mutation in breast cancer.

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 6 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.
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 2 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 TM 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 pvalue < 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 YTHDF1related 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 2 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 5 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.

YTHDF1 Expression Is Correlated to Breast Cancer Patients' Survival
To determine the clinical implications of 13 m 6 A regulatory mRNAs in breast cancer, this study observed the correlations between the expression of 13 m 6 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)  experienced undesirable survival outcomes than those with its low expression.

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).

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).
Correlation Between YTHDF1 Expression and Tumor Immune Microenvironment and Response to Immunotherapy m 6 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.

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   (Figure 7). However, YTHDF1 was lowly expressed in thyroid carcinoma than normal samples.

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-nucleotidemediated signaling, and adenylate cyclase-modulating G proteincoupled 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.

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.

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).

DISCUSSION
In this study, we characterized the expression patterns of m 6 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 6 A modification and their regulators in tumor progression and prognosis. We found that most of m 6 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 6 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 upregulated 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 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). 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 nonsynonymous 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 6 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 6 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 tumorinfiltrating 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 6 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).