Genetic alterations in m6A regulators were widespread in pancreatic cancer
m6A methylation is dynamically mediated by m6A “writers” and “erasers,” and is recognized by “readers.” We initially used cBioPortal to evaluate genetic alterations in 20 main m6A regulators, including 8 writers, 2 erasers, and 10 readers. This was conducted in 776 cases from 4 pancreatic cancer studies, including Pancreatic Cancer (UTSW, Nat Commun 2015) [5], Pancreatic Adenocarcinoma (TCGA, Firehose Legacy), Pancreatic Adenocarcinoma (QCMG, Nature 2016) [20], and Pancreatic Adenocarcinoma (ICGC, Nature 2012) [21](Figure 1). We found that the total alteration frequency of these 4 studies was 15.21% (118/776). Among them, the alteration frequencies of amplification, mutation, and deep deletion were 6.06% (47/776), 5.41% (42/776) and 2.84% (22/776), respectively. Only 0.90% (7/776) of these cases had two or more alterations (Figure 1A). However, among these 4 studies, Pancreatic Cancer (UTSW, Nat Commun 2015) and Pancreatic Adenocarcinoma (ICGC, Nature 2012) used whole exome sequencing (WES); Pancreatic Adenocarcinoma (TCGA, Firehose Legacy) and Pancreatic Adenocarcinoma (QCMG,Nature 2016) used whole genome sequencing (WGS). So we have divided the cases and identified the alteration frequency in WGS and WES cases respectively. The results showed that the total alteration frequency of 208 cases using WES was 27.88% (58/208). Among them, the alteration frequencies of amplification, mutation, and deep deletion were 5.29% (11/208), 11.06% (23/208) and 9.13% (19/208), respectively. Only 2.40% (5/208) of these cases had two or more alterations (Figure 1B). In 568 cases using WGS the total alteration frequency of was 10.56% (60/568). Among them, the alteration frequencies of amplification, mutation, and deep deletion were 5.46% (31/568), 4.23% (24/568) and 0.51% (3/568), respectively. Only 0.35% (2/568) of these cases had two or more alterations (Figure 1C). Thus, amplification was the most common type of genetic alteration amongst the 20 main m6A regulators in pancreatic cancer. The overall average alteration frequencies of m6A writers, erasers, and readers were 1.03%, 1.42%, and 1.31%, respectively (Figure 1D). The alteration frequencies of m6A writers, erasers, and readers in WGS and WES cases were shown in Figure 1E-F respectively. Moreover, m6A writers VIRMA (2.32%) and WTAP (1.29%), m6A eraser ALKBH5 (1.80%), and m6A readers YTHDF1 (2.32%), IGF2BP1 (1.93%), YTHDF3 (1.55%), and YTHDC1 (1.55%) showed higher alteration frequencies. The total alteration frequencies of Pancreatic Cancer (UTSW, Nat Commun 2015), Pancreatic Adenocarcinoma (TCGA, Firehose Legacy), Pancreatic Adenocarcinoma (QCMG, Nature 2016), and Pancreatic Adenocarcinoma (ICGC, Nature 2012) were 50.46% (55/109), 21.62% (40/185), 5.22% (20/383), and 3.03% (3/99), respectively. Pancreatic Adenocarcinoma (QCMG, Nature 2016), and Pancreatic Adenocarcinoma (ICGC, Nature 2012) only had mutation alterations (Additional file 2).
There were 741 cases with mutation alterations amongst the 4 studies. The top five genes with the highest mutation frequency were KRAS, TP53, SMAD4, TTN, and CDKN2A. Their mutation frequencies were 91.00%, 60.10%, 21.60%, 15.50%, and 13.60%, respectively. SMAD4 and TTN showed higher mutation frequency in m6A regulators altered group, while KRAS, TP53 and CDKN2A showed lower mutation frequency in m6A regulators altered group. Though, there was no statistical significance between them (Figure 2A). Next, we analyzed the mutation frequencies of 20 main m6A regulators in the altered and unaltered groups of these five genes. We found that YTHDF2 did not mutate in any group. Except for YTHDF1, YTHDC1, EIF3A and IGF2BP1, the mutation frequencies of other m6A regulators were higher in the altered KRAS group than unaltered group. The mutation frequencies of most m6A regulators were higher in the altered SMAD4 group than unaltered group, except for VIRMA, RBM15B, YTHDF1, YTHDF3, YTHDC2 and HNRNPA2B1. The mutation frequencies of METTL3, WTAP, RBM15, ZC3H13, FTO, ALKBH5, YTHDC1, YTHDC2, EIF3A and IGF2BP1 were higher in the altered TTN group than unaltered group; while the other m6A regulators showed lower mutation frequencies in the altered TTN group. Except for WTAP, VIRMA, RBM15B, YTHDC1, YTHDC2 and IGF2BP1, the mutation frequencies of other m6A regulators were lower in the altered TP53 group than unaltered group. The mutation frequencies of most m6A regulators were lower in the altered CDKN2A group than unaltered group, except for RBM15 and IGF2BP2. Details were shown in Figure 2B-F and Additional file 3. We then collected the mutation data for 20 main m6A regulators across these 4 studies. We found that there were 68 mutations in these m6A regulators in all cases. The most frequent mutations were missense mutations (60 missense mutations, 3 splice sites, 1 nonsense mutation, 1 in frame deletion, 1 frame shift ins, 1 frame shift deletion, and 1 splice region). The number of cases with missense mutations were as follows: seven in ZC3H13, EIF3A, and IGF2BP1; five in YTHDC1 and YTHDC2; four in WTAP, RBM15, and IGF2BP3; three in METTL14; two in METTL3, FTO, YTHDF1, and IGF2BP2; and 1 in VIRMA, RBMX, RBM15B, ALKBH5, YTHDF3, and HNRNPA2B1 (Figure 1G). The overall average number of mutations in m6A writers, erasers, and readers was 3.5, 1.5, and 3.7, respectively. The readers therefore showed relatively high mutation frequencies. Moreover, IGF2BP1 showed the largest overall number of cases with mutations (7 missense mutations, 1 frame shift deletion, and 1 splice region, Figure 1G).
We investigated the copy number alteration (CNA) data across these 4 studies, and found that there were 293 cases with CNAs. VIRMA, YTHDF1, IGF2BP2, IGF2BP3, HNRNPA2B1, RBMX, EIF3A, and ZC3H13 only demonstrated amplified (AMP) alterations, and the alteration frequencies were 5.80%, 5.50%, 1.70%, 1.40%, 1.40%, 0.70%, 0.70%, and 0.30%, respectively. YTHDC1, YTHDF2, METTL14, WTAP, and YTHDF3 had only homozygously deleted (HOMDEL) alterations, and the alteration frequencies were 3.10%, 2.70%, 1.40%, 1.40%, and 0.30%, respectively. METTL3, RBM15, RBM15B, ALKBH5, FTO, IGF2BP1, and YTHDC2 showed both AMP and HOMDEL alterations, and the alteration frequencies were 1.00% and 0.30%, 0.30% and 0.70%, 0.30% and 0.70%, 2.40% and 2.00%, 1.70% and 0.30%, 2.70% and 0.30%, and 0.30% and 1.00%, respectively (Figure 1H). The specific genetic alterations of the 776 cases were shown in Figure 3.
Transcription levels of m6A regulators in pancreatic cancer
We sought to determine whether genetic alterations affect the expression of m6A regulators in pancreatic cancer. We therefore explored the mRNA expression of m6A regulators between PAAD and normal samples by using the GEPIA. We found that the mRNA expression levels of almost all m6A regulators were higher in PAAD samples except for METTL3. Moreover, the expression differences of METTL14, VIRMA, RBM15, ZC3H13, FTO, ALKBH5, YTHDF1, YTHDF2, YTHDF3, HNRNPA2B1, IGF2BP2 and IGF2BP3 between tumor and normal tissues were statistically significant (Figure 4A). Among them, IGF2BP3 showed the greatest difference in mRNA expression between PAAD patients and healthy individuals (Additional file 4A, Additional file 5). Further, we used the UALCAN to analyze the relationship between IGF2BP3 expression and multiple clinicopathological characteristics of 178 PAAD samples and 4 normal samples in TCGA. The results showed that this difference was not associated with tumor grade, but was closely related to disease stages, lymph node metastasis, race, gender, age, drinking habits, diabetes status, and pancreatitis status. Although only stage 2 PAAD patients showed statistically significant differences in IGF2BP3 expression compared with healthy individuals, the median IGF2BP3 expression was higher in stage 2 to 4 patients than in healthy individuals. This may be due to the other groups having insufficient samples. There was almost no statistically significant difference in IGF2BP3 expression in PAAD patient samples among different tumor grades, disease stages, lymph node metastasis, race, gender, age, drinking habits, diabetes status, and pancreatitis status (Additional file 4B-I, Additional file 6).
To verify the abnormal expression of m6A regulators in pancreatic cancer, we performed qRT-PCR in five human PC cell lines (AsPC-1, BxPC-3, Capan-2, Panc-1, and SW1990) and the control cell line HPDE6-C7. The results showed that the mRNA expression levels of most m6A regulators were higher in PC cell lines except for METTL3, RBMX, RBM15B, FTO and ALKBH5. Moreover, the expression differences of METTL14, VIRMA, RBM15, ZC3H13, FTO, ALKBH5, YTHDF2, YTHDF3, IGF2BP1 and IGF2BP3 between PC cell lines and HPDE6-C7 cell were statistically significant (Figure 4B-G). Among m6A writers and erasers, the expression of METTL14 showed the greatest difference between PC cells and HPDE6-C7 cell. Compared with HPDE6-C7 cell, the mRNA expression of METTL14 in all the five PC cell lines were significantly up-regulated (Figure 4H).
To identify the function of METTL14 in tumor proliferation, we knocked down METTL14 in PC cell lines (BxPC-3 and SW1990) and HPDE6-C7 cell using shRNAs (shMETTL14-09, shMETTL14-10, and shMETTL14-11). METTL14 knockdown was confirmed by RT-qPCR and western blot (Figure 5 A-F). ShMETTL14-11 suppressed the expression of METTL14 most. IncuCyte live cell imaging system results showed that knockdown of METTL14 significantly suppressed cell proliferation in BxPC-3, SW1990 and HPDE6-C7 cells. Similar results of CCK8 (48 hours after treatment) were observed in HPDE6-C7 cell. However, the differences between shMETTL14-09 and negative control group in BxPC-3, SW1990 cells were not significant. Further, we transfected BxPC-3, SW1990 and HPDE6-C7 cells with METTL14 overexpression (METTL14-OE) plasmid. Remarkably, IncuCyte live cell imaging system and CCK8 results both showed that METTL14 overexpression significantly increased the proliferation abilities of BxPC-3, SW1990 and HPDE6-C7 cells (Figure 5G-L).
Alterations in m6A regulators affected DFS/PFS and OS in pancreatic cancer
To determine whether alterations in m6A regulators affected the DFS/PFS and OS of pancreatic cancer patients, we used the cBioPortal to evaluate the survival data of patients with or without genetic m6A regulator alterations from 4 pancreatic cancer studies. There were 30 cases with altered m6A regulators and 110 cases with unaltered m6A regulators in the DFS/PFS data. There were 39 cases with altered m6A regulators and 145 cases with unaltered m6A regulators in the OS data. Kaplan-Meier curve and log-rank test analyses revealed that genetic alterations in m6A regulators were associated with poorer DFS/PFS and OS in pancreatic cancer patients (log-rank test, P = 0.035 and 0.007, respectively, Figures 6A–B). We then analyzed the survival data of 177 pancreatic cancer patients in TCGA. After Kaplan-Meier plot and log-rank test analyses, we found that decreased METTL3 and ALKBH5 mRNA expression levels and increased VIRMA, RBM15B, YTHDF3, YTHDC2, HNRNPA2B1, EIF3A, IGF2BP2, and IGF2BP3 mRNA expression levels were significantly associated with poor DFS/PFS in pancreatic cancer patients (P < 0.05). Furthermore, pancreatic cancer patients with lower METTL3, METTL14, RBMX, FTO, ALKBH5, YTHDF1, and YTHDC1 mRNA expression levels or higher VIRMA, HNRNPA2B1, EIF3A, IGF2BP1, IGF2BP2, and IGF2BP3 mRNA expression levels had significantly poorer OS (P < 0.05, Figures 6C–D). Interestingly, almost all m6A readers had a negative effect on the DFS/PFS and OS in pancreatic cancer patients, while m6A erasers had a positive effect. Most m6A writers simultaneously had a negative effect on the DFS/PFS, but a positive effect on the OS.
Biological interaction network of m6A regulators in pancreatic cancer
Previous research has shown crosstalk amongst the m6A writers, readers, and erasers regulating cancer growth and progression [22]. We investigated whether this coexpression existed amongst m6A regulators in pancreatic cancer. We found that genes in the same functional category had highly correlated expression patterns. Moreover, the expression of most writers, erasers, and readers had a highly positive correlation. There were some exceptions to this observation. YTHDF1 expression was negatively correlated with the expression of most genes, except for a slight positive correlation with IGF2BP1/2/3. The expression of IGF2BP1/2/3, especially IGF2BP2 and IGF2BP3, was negatively correlated with the expression of most writers and erasers (Figure 7A). We used STRING to explore PPI networks between these writers, erasers, and readers in humans. The results showed that EIF3A was expected from these networks, and the PPI enrichment P-value was less than 1.0e-16. Significant GO term analysis showed that m6A regulators were mainly located in the intracellular membrane-bounded organelles, nucleus, nuclear lumen, nucleoplasm, nuclear body, and catalytic complex. Here, the regulators participate primarily in N6-methyladenosine-containing RNA binding, mRNA 3'-UTR and 5'-UTR binding, catalytic activity on RNA, translation regulator activity, mRNA (2'-O-methyladenosine-N6-)-methyltransferase activity, and oxidative RNA demethylase activity. m6A regulators have an important role in the regulation of biological processes, including cellular processe, nitrogen compound metabolic processe, primary metabolic processe, macromolecule metabolic processe, cellular metabolic processe, and gene expression. Reactome pathway enrichment analysis showed that these m6A regulators were enriched in the metabolism of RNA, carrying out processing of capped intron-containing pre-mRNA, insulin-like growth factor-2 mRNA binding proteins (IGF2BPs/IMPs/VICKZs) bind RNA, and reversal of alkylation damage by DNA dioxygenases (Figures 7B–C).
Enrichment of coexpressed genes related to m6A regulators in pancreatic cancer
To investigate the coexpressed genes related to m6A regulators in pancreatic cancer, we used LinkedOmics to analyze the mRNA sequencing data of 20 main m6A regulators from 178 PAAD patients in TCGA. For instance, the volcano plot (Figure 8A) showed that 4102 genes (dark red dots) were significantly positively correlated with ALKBH5, whereas 2970 genes (dark green dots) showed significant negative correlations (false discovery rate [FDR] < 0.01). The top 50 genes significantly positively and negatively correlated with ALKBH5 are shown in the heat map (Figures 8B–C). We used DAVID to analyze significantly enriched GO terms of the top 50 genes that are positively and negatively correlated with m6A regulators. The results indicated that proteins encoded by genes that are significantly positively correlated with m6A regulators are localized mainly to the nucleoplasm, chromosome, nucleolus, microtubule cytoskeleton, ribonucleoprotein complex, spindle, and spliceosomal complex. The molecule functions of these genes were mainly regulation of RNA binding, poly(A) RNA binding, DNA binding, and organic cyclic compound binding. They were primarily involved in biological processes including RNA biosynthetic processe, RNA metabolic processe, transcription, DNA-templated transcription, protein modification processe, and gene expression. KEGG pathways analysis showed enrichment in the cell cycle, spliceosome, RNA transport, mRNA surveillance pathway, ubiquitin mediated proteolysis, viral carcinogenesis, and cancer pathways. However, proteins encoded by genes that are significantly negatively correlated with m6A regulators were located mainly in the membrane-bound vesicles, vacuole, mitochondrial envelope, ubiquitin ligase complex, and Golgi apparatus. They mainly carry out molecular functions such as ubiquitin-protein transferase activity, transferase activity, transferring glycosyl groups, U4atac snRNA binding, snRNP binding, and N6-methyladenosine-containing RNA binding. These affect biological processes such as proteolysis, proteolysis involved in cellular protein catabolic processe, ubiquitin-dependent protein catabolic processe, and modification-dependent protein catabolic processe. KEGG pathways analysis showed enrichment mainly in focal adhesion, lysosome, endocytosis, proteasome, protein digestion and absorption, and metabolic pathways (Figure 9).
We used a Venn diagram to show that genes positively correlated with m6A writers were enriched to 23 pathways, readers were enriched to 32 pathways, and erasers were enriched to 1 pathway. There were 9 intersections between pathways enriched by genes positively correlated with m6A writers and readers, but there were no intersections between pathways enriched by genes positively correlated with erasers and writers, nor erasers and readers. Simultaneously, genes negatively correlated with m6A writers were enriched to 22 pathways, readers were enriched to 17 pathways, and erasers were enriched to 15 pathways. There were 5 intersections between pathways enriched by genes negatively correlated with m6A writers and erasers, 1 intersection between writers and readers, and 1 intersection between erasers and readers. Interestingly, there were 1 intersection between pathways enriched by genes positively correlated with m6A readers and genes negatively correlated with m6A erasers, but there were no intersections between pathways enriched by genes positively correlated with m6A writers and genes negatively correlated with m6A erasers (Figures 10A-C). KEGG pathways enrichment analysis of the top 50 genes significantly positively and negatively correlated with each m6A regulator revealed that m6A writers, erasers, and readers affect a number of signaling pathways in pancreatic cancer (Figure 10D). These also have important roles in the occurrence and development of pancreatic cancer.