Global RNA m6A modification patterns differ in WT and NC samples
To determine the global effect on m6A RNA modification in Wilms’ tumor, we collected 23 tumor (WT) and adjacent non-tumorous control (NC) tissue from patients undergoing tumor resection. Histopathologic analysis confirmed that the cancer cells are densely arranged in tubular shape with small size, deep staining of the nuclei, and sparse cytoplasm. Representative images were shown in Fig. 1A.
To determine the global m6A levels in NC and WT samples, we performed LC-MS/MS. The total m6A methylation levels in total RNA of WT samples and NC samples were (0.21 ± 0.01, n = 23)% and (0.22 ± 0.01, n = 23)%, respectively. The WT samples exhibited relatively lower total m6A levels than the NC samples (Fig. 1B). Hower, these was no statistically significant difference between the two groups. To evaluate transcript-specific m6A changes in the WT group, we performed microarray hybridization. Immunoprecipitated m6A-methylated RNAs were isolated from WT and NC group, and the methylated RNAs were profiled using microarray with probes for 41,246 mRNAs. The results showed that 59 transcripts were differentially m6A-methylated between the WT and NC groups (fold change > 1.5, p < 0.05). Nine of the transcripts were significantly hypermethylated and 50 of the transcripts were significantly hypomethylated in WT samples. The hierarchical clustering and the scatter plot of the microarray hybridization data are shown in Fig. 2A and B.
We further evaluated potential differences in the expression levels among all m6A-methylated genes in nephroblastoma. A total of 1335 m6A-methylated transcripts were differentially expressed between the WT and NC groups (fold change > 1.2, p < 0.05). These results suggest that WT samples tend to have lower levels of m6A methylation and higher expression of m6A-hypomethylated genes.
Abnormal m6A-modified genes are enriched in cancer related pathways
To determine the function of the 1335 genes that are differentially m6A-methylated in WT, we performed protein-protein interaction (PPI) network analysis using Cytoscape software. Network analysis of 117 key differentially expressed genes was carried out using STRING software (https://string-db.org/), with the confidence value set at 0.9, and the degree ≥ 5 (Fig. S1). These results suggest that the m6A-methylated genes that are highly expressed in WT may be functionally related.
To explore the pathological and physiological roles of m6A-modified genes in the WT group, we analyzed the functions of the 1335 differentially expressed m6A-methylated transcripts (fold change > 1.2, p < 0.05) by gene ontology (GO) and KEGG pathway analyses. The GO analysis identified “GTPase activity” and “identical protein binding” as the main biological process. “Cytosol” and “extracellular exosome” were significantly enriched GO terms in the cellular component category, and “positive regulation of gene expression” and “protein transport” were major functions (Fig. 3A). The KEGG pathway analysis revealed that among the enriched pathways, the differently methylated genes in the WT group were enriched in the “Wnt signaling pathway”, “mTOR signaling pathway”, and “Proteoglycans in cancer”. These results suggest that differentially m6A-modified RNAs and/or m6A-modified genes with differential expression may participate in the occurrence and development of cancer developmental processes (Fig. 3B).
Among the altered m6A-modified KEGG pathways enriched in WT, the mTOR signaling pathway has been established to be associated with carcinogenesis signaling[9, 10]. Closer examination suggested that STK11, CAB39L, AKT1S1, PDPK1, STRADA, BRAF, PIK3CB, EIF4E2, and PTEN were altered m6A-modified genes within the mTOR signaling pathway. Among the 9 genes, STK11, CAB39L, AKT1S1, STRADA, PIK3CB, EIF4E2 and PTEN were hyper-methylated, PDPK1and BRAF were hypo-methylated. To validate the result of microarray data for these genes, we conducted gene-specific m6A qPCR assays for the 9 genes. The results suggest that these 9 mRNAs were uniformly hypomethylated in WT, though statistical differences were observed for PDPK1, STRADA, PIK3CB and PTEN (Fig. 4A). Then we analyzed the total mRNA expression levels of these 9 genes by qPCR (Fig. 4B). The results showed that PDPK1 and PIK3CB were significantly lowly expressed (Fold Change > 1.5, P < 0.05), and EIF4E2 was significantly highly expressed (Fold Change > 1.5, P < 0.05) in WT. Overall, these results suggest that m6A modification may function in the process of WT development by changing the regulation of the mTOR signaling pathway and other cancer-related signaling pathways.
Conjoint analysis of MeRIP and RNA in WT and NC samples
For a more comprehensive analysis of the mRNA expression profile in the WT group, we also profiled mRNA expression using microarray with probes for 41,246 mRNAs. There were 2091 mRNAs with differential expression levels between the WT and NC groups, of which 615 mRNAs were upregulated and 1167 mRNAs were downregulated (Fold Change > 1.5, P < 0.05). The hierarchical clustering and scatter plot of the microarray hybridization data are shown in Fig. 5A and 5B. The abnormally regulated genes were selected for gene ontology and KEGG pathway analyses. The results suggest that the dysregulated genes in WT were significantly enriched in “small molecular metabolic process”, “regulation of biology quality”, and “transport”. Moreover, pathway analysis showed that “proteoglycans in cancer”, “PPAR signaling pathway” and “PI3K-Akt signaling pathway” were significantly altered in WT (Fig. 5C and 5D).
To further clarify the correlation between altered mRNAs and m6A-modified mRNAs, we performed conjoint analysis of the differently expressed methylated RNAs and differently expressed RNAs (Fold Change > 1.5, P < 0.05 for mRNA and fold change > 1.2, p < 0.05 for m6A). A total of 120 m6A hyper-methylated mRNAs were identified, among which 51 transcripts were significantly downregulated (hyper-down), and 59 transcripts were significantly upregulated (hyper-up). We also detected a total of 668 hypo-methylated m6A mRNAs, with 212 transcripts obviously upregulated (hypo-up) and 456 transcripts obviously downregulated (hypo-down, Fig. 5E). The ranking of the top three genes in each quadrant of Fig. 5E is listed in Table S1.
The m6A “readers” YTHDF1, YTHDF2 and IGF2BP3 are upregulated in WT tissues
We speculated that the altered m6A abundance in the WT group could potentially be the result of altered m6A modification-related genes. To evaluate this possibility, we detected the expression of m6A methylase complex subunits and m6A demethylase by mRNA qPCR. The “writer” METTL3, METTL14 and WTAP were only slightly upregulated while the “eraser” FTO and ALKBH5 were downregulated; however, these differences were not significant (P > 0.05, Fig. 6). More importantly, expression of the m6A “readers” YTHDF1, YTHDF2 and IGF2BP3 was statistically greater in the WT group than in the NC group (1.5 fold, 1.2 fold and > 2.0 fold, respectively, p < 0.05. Figure 6). These results are consistent with the possibility that the overall levels of hypomethylation in the WT tissues are associated with higher expression of m6A “readers”.