Analysis of N6-Methyladenosine Methyltransferase Reveals METTL14 and ZC3H13 as Tumor Suppressor Genes in Breast Cancer

Objectives Recently, an increasing number of studies have revealed that N6-methyladenosine (m6A) functions as a significant post-transcriptional modification which plays a critical role in the occurrence and progression of enriched tumors by regulating coding and non-coding RNA biogenesis. However, the biological function of m6A in breast cancer remains largely unclear. Materials and Methods In this study, we used a series of bioinformatic databases and tools to jointly analyze the expression of m6A methylation transferases (METTL3, METTL14, WTAP, RBM15, RBM15B and ZC3H13) and investigate the prognostic value of METTL14 and ZC3H13 in breast cancer. Besides, we analyzed the downstream carcinogenic molecular mechanisms related to METTL14 and ZC3H13 and their relationship with immune infiltration in breast tumor tissues. Results The results showed that METTL14 and ZC3H13 were the down-regulated m6A methylation transferases in breast cancer. Survival outcome analysis suggested that abnormally low expression of METTL14 and ZC3H13 could predict unfavorable prognosis in four breast cancer subtypes. Moreover, their down-regulation was associated with ER-, PR- and triple-negative breast cancer patients, as well as tumor progression (increased Scarff, Bloom and Richardson grade status and Nottingham Prognostic Index classification). Co-expression analysis revealed that METTL14 and ZC3H13 had a strong positive correlation with APC, an antagonist of the Wnt signaling pathway, indicating they might cooperate in regulating proliferation, invasion, and metastasis of tumor cells. METTL14, ZC3H13, and APC expression levels had significant positive correlation with infiltrating levels of CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells, and negative correlation with Treg cells in breast cancer. Conclusions This study demonstrated that down-regulation of METTL14 and ZC3H13 which act as two tumor suppressor genes was found in breast cancer and predicted poor prognosis. Their abnormal expression promoted breast cancer invasion by affecting pathways related to tumor progression and mediating immunosuppression.


INTRODUCTION
The morbidity and mortality of breast cancer rank first and second in global female tumors, respectively. Globally, the incidence of breast cancer has been increasing year by year, and the age of onset has gradually become younger (1). The incidence of breast cancer is relatively low in China, but the number of newly diagnosed breast cancer cases has been increasing continuously in recent years. It has already been the malignant tumor with the highest incidence among women in some large or medium-sized cities, which seriously threatens women's health and life, and brings huge economic and health problems to the society of China. Although the prognosis of breast cancer continues to improve with the progress of treatment, breast cancer is still the main cause of female death at present (2).
The occurrence and development of tumors are driven by the disorder of genetic, epigenetic, and environmental factors, and epigenetic factors play an important role as a bridge between genetic and environmental factors (3,4). RNA modifications are also important types of post-transcriptional epigenetic modification, and N6-methyladenosine (m6A) is the most predominant modification of mRNA. In mammalian cells, there is an average of 1-2 m6A sites per 1,000 nucleotides (5,6). At present, many proteins involved in m6A have been identified, which are classified into m6A methylation transferases ("writer" proteins), m6A demethylases ("eraser" proteins), and m6A reading and binding proteins ("reader" proteins) (7). RNA is methylated at the sixth nitrogen atom of adenylate by the action of " writer" proteins, and methyltransferase-like 3 (METTL3) and methyltransferase-like 14 (METTL14) form a hetero complex, which complete this modification process together with Wilms tumor 1 associated protein (WTAP) and other factors, such as Putative RNAbinding protein 15 (RBM15) and Zinc Finger CCCH-Type Containing 13 (ZC3H13). By contrast, the m6A modified regions can be demethylated under the action of two "eraser" proteins, fat mass, and obesity-associated protein (FTO) and alkB homolog 5 (ALKBH5), making m6A a reversible modification. These m6A modified sites can be recognized by "readers" proteins to regulate RNA metabolism, including translation, splicing, translocation, degradation, and processing (7).
The abnormal expression of m6A "writer" proteins results in abnormal levels of m6A modification in tumor cells, which further abnormalizes the metabolism of tumor-related genes mRNA, thereby promoting the development of a variety of cancers (8,9). The up-regulated METTL3 in gastric cancer was positively correlated with the poor prognosis of patients, and METTL3 promoted epithelial-mesenchymal transition of gastric cancer, leading to tumor invasion and metastasis (10). WTAP expression was up-regulated in hepatocellular carcinoma; WTAP-mediated m6A modification led to epigenetic silencing of the tumor suppressor gene ETS1, and promoted the progression of hepatocellular carcinoma by regulating the cell cycle (11). The expression of METTL14 was down-regulated in colorectal cancer, which indicated a poor prognosis of patients, and promoted the proliferation and invasion of tumor cells by inhibiting the degradation of the oncogene XIST (12). However, the abnormal expression of m6A "writer" proteins in breast cancer still remains largely unknown, and the gene targets and molecular mechanisms involved downstream also need to be further elucidated.
In this study, we used Oncomine and The Cancer Genome Atlas (TCGA) databases to jointly analyze the expression of m6A "writer" in breast cancer including METTL3, METTL14, WTAP, RBM15, RBM15B, and ZC3H13, indicating that the mRNA expression levels of METTL14 and ZC3H13 were lower in tumor samples. We further combined the data of TIMER, cBioPortal, Kaplan-Meier plotter and other databases to analyze the prognostic value of METTL14 and ZC3H13 in breast cancer, their co-expressed genes and related signaling pathways, as well as the relationship between their expression levels and infiltrating immune cells in tumor tissues. In total, we highlighted the important role of METTL14 and ZC3H13 in breast cancer, provided promising markers for predicting prognosis of patients, and explained the potential molecular mechanism of the two as tumor suppressor genes.

Gene Expression Data Mining
ONCOMINE gene expression array dataset (13) (https://www. oncomine.org), an online cancer microarray database, was utilized to analyze the transcription levels of different m6A RNA methylation writer genes: METTL3, METTL14, WTAP, RBM15, RBM15B, and ZC3H13 in different cancers. The cut-off of P value was 0.05 and fold change were 2.0, respectively. Then, METTL3, METTL14, RBM15B, ZC3H13, and APC gene expression levels in various cancer types were further analyzed and verified via The Tumor IMmune Estimation Resource 2.0 (TIMER 2.0) algorithm database (14) (http://timer.cistrome.org/) based on data from The Cancer Genome Atlas (TCGA). The heatmaps of the selected gene mRNA expression were conducted and visualized by the UCSC Xena (15,16) (https://xena.ucsc.edu/public), an interactive online visualization of seminal cancer genomic dataset, based on the data from TCGA.

Survival Outcome Analysis
The prognostic values of METTL14, ZC3H13, APC, and selected hub genes in breast cancer samples were assessed by the online tools and database, Kaplan-Meier plotter (19) (http://kmplot. com/analysis/index.php?p=service&cancer=breast). To analyze the overall survival (OS) and progression-free survival (RFS) of breast cancer patients, patient samples were divided into two groups by the best cut-off value by the tool automatically and calculated via the Kaplan-Meier analysis and Logrank-P test. Besides, the prognostic roles of METTL14, ZC3H13, and APC in enriched datasets from Gene Expression Omnibus database (GEO) were further analyzed and verified by PrognoScan database (20) (http://dna00.bio.kyutech.ac.jp/PrognoScan/ index.html). The prognostic values of METTL14, ZC3H13, APC, and selected hub genes were also estimated with TCGA breast cancer samples by using Kaplan-Meier plotter.

The Associations of Gene Expression and Clinical Parameters
The METTL14 and ZC3H13 mRNA expression data and the clinical data of breast cancer patients in TCGA were extracted from the UCSC Xena (15,16). After deleting incomplete cases, the high and low expression groups were divided by the median value of the mRNA expression level and then for further analysis. Besides, the associations of METTL14 and ZC3H13 gene expression and clinical parameters were also evaluated and verified using the Breast Cancer Gene-Expression Miner v4.0 (bc-GenExMiner v4.0) online dataset (21) (http://bcgenex. centregauducheau.fr/BC-GEM/GEM-Accueil.php?js=1), in accordance with the RNA-seq data and clinical data from TCGA and GEO datasets.

Gene Correlation Analysis
The top 200 co-expressed genes of METTL14 and ZC3H13 were analyzed and obtained from the Breast Invasive Carcinoma (TCGA, PanCancer Atlas) on cBioPortal by Spearman's Correlation method (17,18). Then, the co-expressed genes of the two genes were cross-referenced to obtain a cohort of 76 common co-expressed genes. The correlation between APC and METTL14, ZC3H13 were also analyzed on GEPIA (22) (http:// gepia.cancer-pku.cn/detail.php) and bc-GenExMiner v4.0 (21) by Spearman's and Pearson's correlation methods, respectively.

Pathway and Gene Ontology Enrichment
Functional enrichment analysis including Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment was performed based on the WebGestalt website (23) (http://www.webgestalt.org/option. php) to understand the function of the extracted 76 common co-expressed genes. GO analysis can be divided into three categories: biological processes (BPs), cellular components (CCs), and molecular functions (MFs). P <0.05 was considered significant via Fisher's exact test.

Protein-Protein Interaction Network and Sub-Module Construction
The STRING database (24) was utilized to construct the PPI network of 76 common co-expressed genes with an interaction score above 4.0 and hiding the unconnected genes. Then, ClueGO APP on Cytoscape 3.7.2 (25) was employed to assess the function of the given PPI network, and the most important module and hub genes of the PPI network were selected and constructed by the MCODE APP according to the rules as follows: degree cut-off = 2, node score cut-off = 0.2, Max depth = 100, and k-score = 2.

Infiltrating Immune Cell Analysis
The correlation of mRNA expression level of METTL14, ZC3H13, APC with the abundance of seven types of infiltrating immune cells (CD4+ T cells, CD8+ T cells, Treg cells, B cells, neutrophils, macrophages, and dendritic cells) and the association between immune infiltrates and somatic CNV of METTL14 and ZC3H13 in breast cancer patients were calculated using The Tumor IMmune Estimation Resource 2.0 (TIMER 2.0) algorithm database (26,27) (http://timer.cistrome.org/), noting that the abundance of Treg cells was analyzed by the CIBERSORT method. Tumor purity is an important factor affecting the analysis of immune infiltration in tumor samples by genomic methods.

Tissue Samples and Immunohistochemical Staining
The breast cancer tissue microarrays (F601101) comprising 50 breast tumor tissues and 10 adjacent normal tissues were purchased from Zhong Ke Guang Huang Biotech (http:// bioaitech.com). IHC was performed on the tissue microarrays using the UltraSensitive ™ S-P Methods. Briefly, after the tissue samples were first dewaxed, H 2 O 2 was used to block the activity of endogenous peroxidase. The primary antibodies for METTL14 and ZC3H13 were incubated with a diluted solution, and then a secondary antibody was added and developed by diaminobenzoquinone (DAB). At last, hematoxylin was used to counterstain the microarrays.

Statistical Analysis
Students' t test was employed to compare mRNA expression in Oncomine and TIMER 2.0, also to analyze the staining index. ANOVA was used to identify the significance of METTL14 and ZC3H13 mRNA levels in different putative copy-number alterations, mutations groups. The Kaplan-Meier curve and Logrank P test in Kaplan-Meier (KM) plotter, and Univariable COX in PrognoScan were used to analyze the survival outcomes. The chi-squared test was used to investigate the significance of the correlation of METTL14 and ZC3H13 expression with clinical parameters in different breast cancer groups. Gene expression significant difference between subgroups is assessed by Welch's test on bc-GenExMiner v4.0. The correlation of gene expression or infiltrating immune cells was evaluated using the Spearman's and Pearson's correlation method. Fisher's exact test was adopted to measure the gene enrichment. P <0.05 was considered statistically significant.

Down-Regulation of METTL14 and ZC3H13 mRNA Expression in Breast Cancer
First off, the expression profiles of N6-methyladenosine (m6A) RNA methylation associated methyltransferases were analyzed via Oncomine database. The expression of METTL3, METTL14, RBM15B, and ZC3H13, but not of WTAP and RBM15, was significantly decreased in breast cancer ( Figure 1A). More details, Oncomine analysis of cancer vs normal samples in enriched breast cancer patient datasets showed that their expression was reduced in invasive breast carcinoma stroma, invasive ductal breast carcinoma stroma, invasive mixed breast carcinoma, and ductal breast carcinoma in situ, respectively ( Figures 1B-F). Furthermore, mining of TIMER database based on TCGA breast cancer samples demonstrated that only the mRNA expression levels of METTL14 and ZC3H13 were also lower in tumor samples ( Figure 1G), and there was not any difference of METTL3 and RBM15B ( Figure S2). Considering that WTAP did not carry a consistent gene abnormal expression trend in Oncomine database ( Figures 1A and S1), and that METTL3 and RBM15B were only downregulated in Oncomine database but not in TCGA breast cancer samples ( Figures 1A  and S2), we excluded these three genes and only included METL14 and ZC3H13 for subsequent studies ( Figure 2).

Rare Genomic Alteration of METTL14 and ZC3H13 in Breast Cancer
We evaluated the genomic alteration of METTL14 and ZC3H13 among breast cancer samples by the cBioPortal online database (TCGA, PanCancer Atlas). Overall, about 36 (4%) samples, a very small ratio, had genetic changes found in breast cancer, with 0.9% in METTL14 and 2.7% in ZC3H13, respectively ( Figure  3A). Moreover, there were 0.5% samples with amplification and 0.4% with mutation in METTL14, and 1.4% with deep deletion, 0.2% with amplification, and 1.1% with mutation in ZC3H13 ( Figure 3B). The mutation ratio of ZC3H13 was relatively higher, including 11 missense and six truncating points, whereas only four missense mutation points were found in METTL14 ( Figure 3C). In addition, for the relationship between expression and genomic alteration, the expression of METTL14 and ZC3H13 was much higher in samples with amplification and lower in those with deletion compared with normal samples ( Figure 3D). However, we did not find any differences between mutant and wild type samples ( Figure 3E).

Low METTL14 and ZC3H13 Expression Demonstrate Poor Prognosis in Breast Cancer
In order to explore the roles of METTL14 and ZC3H13 in the survival of breast cancer patients, the Kaplan-Meier plotter was  used to calculate the correlation between the mRNA expression levels and survival. The Kaplan-Meier curve and Log-rank P test confirmed that the low expression of METTL14 (Prob ID: 235552 at) and ZC3H13 (Prob ID: 209851 at) was negatively correlated with the overall survival (OS) and progression-free survival (RFS) (Figures 4A-D). But we only found that reduced expression of METTL14 but not ZC3H13 was associated with poor RFS, and there was no correlation between OS and the expression of the two in TCGA breast cancer samples ( Figure  S3). Generally, breast cancer can be divided into four types according to the molecular characteristics, including luminal type A, luminal type B, HER2 enriched type, and triple negative type, which could influence the adjuvant treatments and the prognosis of patients (28). Survival analysis showed that the decreased expression levels of METTL14 and ZC3H13 were related to RFS in all of four types ( Figures 4E, F). Besides, we also made a validation of their prognostic roles using Prognoscan database and got a consistent trend in various survival types ( Table 1). Thus, the two genes play a role in tumor-suppressing and can be used as potential prognostic markers of breast cancer.

The Association of METTL14 and ZC3H13 Expression With Clinical and Pathological Features
After excluding cases with incomplete clinical data of TCGA breast cancer patients from UCSC Xena, 637 cases were included to analyzed the association of METTL14 and ZC3H13 expression with clinical and pathological features through chi-squared test.
The results showed that the expression of METTL14 and ZC3H13 was significantly associated with ER status, PR status, PAM50 subtype. Besides, the proportion of patients with N2-N3 stages is higher in the low expression group of ZC3H13 ( Table  2). Further, the profiles of 4,712 breast cancer patient cohorts in bc-GenExMiner 4.0 were examined for validation. METTL14 and ZC3H13 expression was significantly decreased in ER-, PRand basal-like, TNBC groups. However, HER2− patients had somewhat increase in METTL14 and ZC3H13 mRNA expression. In addition, compared with p53 wild type samples, the genes' mRNA expression was remarkably decreased in p53 mutated samples. As for the Scarff, Bloom and Richardson (SBR) grade status and Nottingham Prognostic Index (NPI) classification, the expression of the two genes gradually decreased with increasing grade and Index ( Figure 5 and Figure S4). Generally, with higher rate of SBR or NPI, the lower of the survival rate was associated. There was no significant relationship between nodal status.

KEGG and GO Enrichment of Co-Expressed Genes
The cBioPortal database was applied to select the top 200 positively co-expressed genes of METTL14 and ZC3H13 based on the data from Breast invasive carcinoma (TCGA, PanCancer  Atlas). The co-expressed genes obtained from the two genes were used as intersection to build a set of 76 common co-expressed genes ( Table 3 and Figure 6A). Then, we used WebGestalt for functional and pathway enrichment analysis of the 76 coexpressed genes. The results of GO analysis indicated that these genes mainly play a vital role in the regulation of chromatin organization, cellular response to DNA damage stimulus, cell cycle, histone modification, and regulation of microtubule, consistent with their molecular function and cellular component. It has been confirmed that these processes and functions are involved in tumor development and drug sensitivity ( Figures 6B-D , Table S1). Meanwhile, the pathways enriched by KEGG were also related to tumors, including hepatocellular carcinoma, Hippo signaling pathway, MAPK signaling pathway, and pluripotency of stem cells ( Figure 6E, Table S1).

METTL14 and ZC3H13 Associated PPI Network Construction
The PPI network was constructed based on the 76 co-expressed genes by using the STRING database, and the crucial modules were built by Cytoscape ( Figure 7A). The function of this PPI network is mainly involved in transcription factor binding, chromatin organization, and cellular response to abiotic stimulus ( Figure 7B). There were 10 genes in three crucial modules, including GOLGB1, TRIP11, GCC2, ASH1L, WDFY3, KIAA1109, RNF111, KLHL11, HECTD1, and UBR1 ( Figure 7C). The clustered heatmap of these 10 genes and METTL14, ZC3H13 was analyzed and visualized by UCSC Xena, indicating their related expression pattern ( Figure 7D, Figure S6). Furthermore, the prognostic roles of the 10 genes were performed using Kaplan-Meier plotter. Accordingly, WDFY3, KLHL11, GOLGB1, KIAA1109, UBR1, ASH1L, GCC2, and TRIP11 exhibited poor RFS in the lower expression groups ( Figure 7E). In TCGA breast cancer samples, only reduced expression of GCC2 and ASH1L could show poor RFS ( Figure S5).

Co-Expression of APC and METTL14, ZC3H13
Based on the enrichment results of KEGG pathway (Table S1), we noticed a familiar tumor suppressor gene, APC, which is associated with several carcinoma related pathways: Hippo signaling pathway, MAPK signaling pathway, and pluripotency of stem cells. It is confirmed that APC acts as an antagonist of the Wnt signaling pathway in several cancer types, and the low expression of APC promotes tumor progression. Besides, it is also involved in other processes, like cell migration and adhesion, transcriptional activation, and apoptosis (29,30). The data and results from different databases and datasets showed that both METTL14 and ZC3H13 had high correlation coefficients with APC (R > 0.6, P < 0.05) (Figures 8A-F). Then, a thorough inquiry of the expression profile of APC was made using TIMER database, indicating that its mRNA expression was significantly decreased in enriched cancer types, such as breast cancer, colorectal cancer, lung squamous cell carcinoma, lung adenocarcinoma, etc ( Figure 8G). Subsequently, the prognostic value of APC (Prob ID: 215310 at) in breast cancer was also performed using Kaplan-Meier plotter; the results confirmed that the lower expression of APC mRNA was associated with the worse RFS of different breast cancer types ( Figure 8H). However, there was no relationship between the expression of APC and RFS in TCGA breast cancer samples ( Figure S5).

Correlation of METTL14 and ZC3H13 Expression With Seven Types of Immune Cells
Studies have confirmed that the excessive activation of the Wnt signaling pathway in tumors can lead to tumor immunosuppression and immune escape, further promoting tumor invasion and metastasis (31,32). Based on the correlation between METTL14, ZC3H13, and APC expression pattern, we speculated that METTL14 and ZC3H13 would also be indirectly involved in mediating the immune response of tumors. So, we analyzed the correlation between the genes expression and seven types of infiltrating immune cells (CD4+ T cells, CD8+ T cells, Treg cells, B cells, neutrophils, macrophages and dendritic cells) by TIMER database. METTL14, ZC3H13, and APC were all significantly positively correlated with CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells in BRCA, negatively with Treg cells (P < 0.05) ( Figure 9A). In different breast cancer subtypes, the correlations were not all the same. Importantly, APC showed a better correlation with these immune cells than METTL14 and ZC3H13 (Figures 9B-D). Besides, we also calculated the association between immune infiltrates and somatic CNV of METTL14 and ZC3H13 in breast cancer, and these immune cells showed diverse enrichment trends in different CNV types across different types of breast cancer ( Figure S7).

Validation of METTL14 and ZC3H13 Expression in Breast Cancer Tissues
To verify the results of bioinformatics analysis, immunohistochemical (IHC) staining was performed on tissue microarray slides containing 50 breast cancer tissues and 10 adjacent normal tissues. The staining patterns of METTL14 and ZC3H13 in tumor and normal tissues were shown in Figures 10A, B. METTL14 was mainly localized in the   nucleus but not in other parts of cells, and ZCH3H13 was observed not only in the nucleus but also in the cytoplasm in cells. Considering that m6A modification took place in the nucleus, we compared the expression level of the two proteins in the nucleus. The levels of the expression were quantitated by the Staining index based on tumor cell proportion and staining intensity. The result further confirmed that both METTL14 and ZC3H13 were significantly down-regulated in the nucleus in tumor tissues relative to normal control, which was consistent with the results of bioinformatics analysis on RNA levels ( Figures 10C, D).

DISCUSSION
The pathological and biological functions of enriched chemical modifications of proteins and DNA have been well verified and  extensively studied. Recently, m6A modification, as a new layer of RNA epigenetic regulation, has also been increasingly studied in human diseases including cancer (33). With the development of m6A sequence technology, m6A modification sites (RRACH, R = A or G; H = A, C, or U) which are mainly enriched near stop codons, in 3′-untranslated regions (3′-UTRs) and within long internal exons have been identified in a great many of coding or non-coding RNAs (34). With the in-depth study of the process of m6A, researchers have identified a series of proteins related to this modification, including "writers" (METTL3, METTL14, WTAP, RBM15, RBM15B, and ZC3H13), "erasers" (FTO and ALKBH5) and "readers" (IGF2BP1/2/3, YTHDF1/2/3, YTHDC2, HNRNPA2B1, etc) (7). By regulating the level of m6A modification, the abnormal expression of these proteins in tumors affects the fate and metabolism of tumor-related mRNAs, including translation, splicing, translocation, degradation, and processing and plays an important role in the occurrence and development of cancer (35).  leukemia progression (36,37). METTL3 was also over-expressed in pancreatic cancer, bladder cancer, glioma and gastric cancers, and promotes proliferation, invasion, and drug resistance of cancer cells (38)(39)(40)(41). By contrast, METTL14 was down-regulated in colorectal cancer and hepatocellular carcinoma and could inhibit cancer cells' growth and metastasis by regulating m6A-dependent primary microRNA processing (12,42).
In agreement with previous reports (12,42), our study also identified METTL14 as a tumor suppressor gene in breast cancer, as well as ZC3H13, while other m6A "writers" (METTL3, RBM15, and RBM15B) were not dysregulated in tumor tissues compared with normal control based on the data from TCGA and Oncomine databases. We also compared their expressions in tissue samples, and the obtained expression was consistent with the bioinformatics analysis. What is more, survival analysis showed that low METTL14 and ZC3H13 expression could demonstrate poor prognosis in breast cancer for OS and RFS. Low expression of METTL14 and ZC3H13 was related to breast cancer progression (increased SBR grade status and NPI classification) and was significantly decreased in ER-, PR-and basal-like, TNBC patients. The above results suggested that the down-regulated METTL14 and ZC3H13 led to a decrease in m6A modification levels in breast cancer, which promoted cancer cells' invasion and metastasis. Interestingly, it has been demonstrated that m6A "eraser" FTO was significantly upregulated in breast cancer, which could promote breast cancer cell proliferation, colony formation and reduce apoptosis (43). In addition, it is confirmed that ALKBH5 was also up-regulated in breast carcinoma, promoted demethylation of m6A, and then improved the stemness of breast cancer cells (44). In summary, these studies combined with our results have confirmed the reduction of m6A modification levels in breast cancer.
M6A "writers" and "erasers" determine the m6A modification on a specific mRNA. It is also important that m6A "readers" can specifically recognize and bind the modification sites to sort the mRNA and transmit information, thereby establishing an effective m6A associated network (7). To further elucidate the downstream molecular mechanism involved in METTL14 and ZC3H13 by regulating m6A in breast cancer, we selected a set of 76 common positively co-expressed genes of METTL14 and ZC3H13. The pathways enriched by KEGG of these genes were associated with tumors, such as hepatocellular carcinoma, Hippo signaling pathway, MAPK signaling pathway and pluripotency of stem cells. Importantly, we noticed a familiar tumor suppressor gene, APC, among these genes. It was confirmed that APC could inhibit the Wnt signaling pathway, a carcinogenic signaling pathway in several cancers, and the low expression of APC promoted tumor progression (45). Afterwards, we analyzed and confirmed that there were multiple m6A sites on the APC mRNA via the Whistle database (Table S2). In addition, by analyzing RNA binding proteins that can bind to APC mRNA by Starbase database, we found that m6A "reader" Insulin Like Growth Factor 2 MRNA Binding Protein 1/2/3 (IGF2BP1/2/3) could bind to several regions of APC mRNA ( Figure S8). IGF2BP1/2/3 can recognize and bind the conserved GG(m6A)C sequence on mRNA under normal and stress conditions and enhance the stability of mRNA. In HEK293T cells, IGF2BP1/2/3 mainly recognizes and binds to the m6A sites near the stop codon of transcription, leading to the accumulation of carcinogenic products such as MYC (46). In colorectal cancer, IGF2BP2 can recognize m6A sites in the CDS region of SOX2, enhancing its stability and preventing its degradation, leading to the occurrence and development of colorectal cancer (47). Based on these documents and results, we proposed a scientific hypothesis ( Figure 10E): METTL14 and ZC3H13 can promote the m6A modification process and level of APC mRNA, after which IGF2BP1/2/3 recognizes the m6A sites and improves the stability of APC mRNA. However, since the expressions of METTL14 and ZC3H13 in breast cancer were significantly reduced, this process is reversed, resulting in a decrease in the stability of APC mRNA and inhibiting the protein level of APC. Thus, a decrease in APC protein level leads to abnormal activation of Wnt signaling pathway in breast cancer.
Abnormal activation of the Wnt pathway is thought to positively regulate various properties of tumor cells, including proliferation, survival, and metastasis (48). In addition, Wnt signaling also plays an important role in the regulation of tumor immune microenvironment, which can regulate the infiltration and activity of various types of T cells in tumors. Usually excessive activation of Wnt signaling pathway in tumors can lead to tumor immunosuppression and immune escape (49,50). In colorectal cancer, inhibiting the activation of the Wnt signaling pathway can reduce the infiltration of Treg cells in the tumor and up-regulate the content of CD8+ T cells in the tumor, thereby enhancing the response to immunotherapy (51). Based on the correlation of APC and METTL14, ZC3H13 in breast cancer, we speculated that METTL14 and ZC3H13 were indirectly involved in mediating tumor immune responses. Analysis of the correlation between gene expression and seven infiltrating immune cells via the TIMER database showed that METTL14, ZC3H13, and APC were all significantly positively correlated with CD4+ T cells, CD8+ T cells, neutrophils, macrophages, and dendritic cells in BRCA, negatively with Treg cells. This indicated that low expression of METTL14 and ZC3H13 promoted immunosuppression in the breast cancer, which is also an important downstream mechanism for their role in the progression and metastasis of breast cancer.
In conclusion, our study revealed that METTL14 and ZC3H13, as two tumor suppressor genes, were down-regulated in breast cancer, and the low expression of METTL14 and ZC3H13 was negatively correlated with the OS and RFS. Their abnormal expression promoted the development of breast cancer by affecting pathways related to tumor development. Among one of these downstream molecular mechanisms, they could affect the m6A modification level of APC mRNA, then activate the Wnt signaling pathway. However, further experimental studies are necessary to provide solid confirmation of our scientific hypothesis and results.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
P-JG and J-WZ contributed to the conception of the study. P-JG and Y-CS contributed to experimental technology and experimental design. P-JG, YY, XH, and W-JS performed the data analyses. P-JG, S-RH, and Y-FZ wrote the manuscript. LW and J-WZ supervised the study. All authors contributed to the article and approved the submitted version.

FUNDING
This work was supported by funds from Health commission of Hubei Province scientific research project (WJ2019H020, WJ2019H028).

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2020.578963/ full#supplementary-material SUPPLEMENTARY TABLE 1 | Pathway and GO enrichment results of 76 common co-expressed genes of METTL14 and ZC3H13.
FIGURE S1 | The expression of WTAP in breast cancer tissues vs normal control in the datasets of Oncomine databas.