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 analyze via Oncomine database. The expression of METTL3, METTL14, RBM15B and ZC3H13, but not of WTAP and RBM15, were 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 (Figure 1B-F). Furthermore, mining of TIMER database based on TCGA BRCA samples demonstrated that only the mRNA expression levels of METTL14 and ZC3H13 were also lower in tumor samples (Figure 1G, H), and there was not any difference of METTL3 and RBM15B (Figure S1).
Rare genomic alteration of METTL14 and ZC3H13 in Breast cancer
We evaluated the genomic alteration of METTL14 and ZC3H13 by the cBioPortal online database for BRCA (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 2A). Moreover, there were 0.5% samples with amplification and 0.2% with mutation in METTL14, and 1.4% with deep deletion, 0.2% with amplification and 1.1% with mutation in ZC3H13 (Figure 2B). The mutation ratio of ZC3H13 was relatively higher, including 11 missense and 6 truncating points, whereas only 4 missense mutation points were found in METTL14 (Figure 2C). 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 2D). However, we did not find any differences between mutant and wild type samples (Figure 2E).
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 the survival. The Kaplan-Meier curve and Log-rank P test confirmed that the low expression of METTL14 and ZC3H13 was negatively correlated with the overall survival (OS) and progression-free survival (RFS) (Figure 3A-D). 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 (Figure 3E, 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 BRCA from UCSC Xena, 637 cases were included to analyzed the association of METTL14 and ZC3H13 expression with clinical and pathological features through chi-square test. The results showed that the expression of METTL14 and ZC3H13 was significantly associated with ER statues, PR statues, PAM50 subtype. Besides, the proportion of patients with N2-N3 stage is higher in the low expression group of ZC3H13 (Table 2). Further, the profiles of 4712 breast cancer patients cohorts in bc-GenExMiner 4.0 were examined for validation. METTL14 and ZC3H13 mRNA expression was significantly decreased in ER-, PR- and 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 mRNA expression of the two genes gradually decreased with increasing grade and Index (Figure 4 and Figure S2). Generally, with higher rate of SBR or NPI, the lower of the survival rate was associated. There was no significant relationship between nodal statues.
KEGG and GO enrichment of co-expressed genes
The cBioPortal database was applied to select 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 gene were used intersection to build a set of 76 common co-expressed genes (Table 3 and Figure 5A). Then, we used WebGestalt for functional and pathway enrichment analysis of the 76 co-expressed 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 (Figure 5B, C, D). 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 5E).
METTL14 and ZC3H13 associated PPI network construction
The PPI network was construction based on the 76 co-expressed genes using the STRING database, and the crucial modules were built by Cytoscape (Figure 6A). The function of this PPI network mainly involved in transcription factor binding, chromatin organization and cellular response to abiotic stimulus (Figure 6B). There were 10 genes in 3 crucial modules, including GOLGB1, TRIP11, GCC2, ASH1L, WDFYS, KIAA1109, RNF111, KLHL11, HECTD1 and UBR1 (Figure 6C). The clustered heatmap of these 10 genes and METTL14, ZC3H13 was analyzed and visualize by UCSC Xena, indicating their related expression pattern (Figure 6D, Figure S3). Furthermore, the prognostic roles of the 10 genes were performed using KM plotter. Accordingly, WDFYS, KLHL11, GOLGB1, KIAA1109, UBR1, ASH1L, GCC2 and TRIP11 exhibited poor RFS in lower expression groups (Figure 6E).
Co-expression of APC and METTL14, ZC3H13
Based on the enrichment results of KEGG pathway, we noticed a familiar tumor suppressor gene, APC, which associated in 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) (Figure 7A-F). Then, the expression profile of APC was made a thorough inquiry 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 7G). Subsequently, the prognostic value of APC in breast cancer was also performed using KM plotter, the results confirmed that the low expression of APC mRNA was associated with the worse RFS of different breast cancer types (Figure 7H).
Correlation of METTL14 and ZC3H13 expression with 7 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 7 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 8A). 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 (Figure 8C, D, E). 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 S4).
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 Figure 8A & B, METTL14 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 among tumor tissues relative to normal control, which was consistent with the results of bioinformatics analysis (Figure 8C-F).