Genetic alterations in m6A regulators are 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) [22], and Pancreatic Adenocarcinoma (ICGC, Nature 2012) [23] (Fig. 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 (Fig. 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 (Fig. 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 (Fig. 1C). Thus, amplification was the most common type of genetic alteration among 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 (Fig. 1D). The alteration frequencies of m6A writers, erasers, and readers in WGS and WES cases were shown in Fig. 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 among 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 statistically significant difference between them (Fig. 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 Fig. 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 (Fig. 1G). The overall average number of mutations in m6A writers, erasers, and readers was 3.5, 1.5, and 3.7, respectively. Among these three, the readers showed the highest mutation frequencies. Moreover, IGF2BP1 showed the largest overall number of cases with mutations (7 missense mutations, 1 frame shift deletion, and 1 splice region, Fig. 1G).
We investigated the copy number alteration 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 (Fig. 1H). The specific genetic alterations of the 776 cases were shown in Fig. 3.
Transcription levels of m6A regulators in pancreatic cancer
We sought to determine whether or not genetic alterations affect the expression of m6A regulators in pancreatic cancer. We 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 (Fig. 4A). Among them, IGF2BP3 showed the greatest difference in mRNA expression between PAAD patients and healthy individuals (Additional files 4A, 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 files 4B-I, 6).
To verify the abnormal expression of m6A regulators in pancreatic cancer, we performed RT-qPCR 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 (Fig. 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 (Fig. 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 blotting (Fig. 5A–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 h after treatment) were observed in these cells, except for the differences between shMETTL14-09 and NC groups in BxPC-3 and SW1990 cells. Further, we transfected BxPC-3, SW1990 and HPDE6-C7 cells with 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 (Fig. 5G–L).
Alterations in m6A regulators affect DFS/PFS and OS of pancreatic cancer patients
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, Fig. 6A, B, Additional file 7). 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, Fig. 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 among the m6A writers, readers, and erasers regulating cancer growth and progression [24]. We investigated whether this coexpression existed among 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 (Fig. 7A, Additional file 8). We used STRING to explore PPI networks between these writers, erasers, and readers in humans. The results show that EIF3A is expected from these networks, and the PPI enrichment P-value is less than 1.0e−16. Significant GO term analysis shows that m6A regulators are 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 process, nitrogen compound metabolic process, primary metabolic process, macromolecule metabolic process, cellular metabolic process, and gene expression. Reactome pathway enrichment analysis shows that these m6A regulators are 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 (Fig. 7B, C, Additional file 9). In addition, we used Cytoscape to perform the crosstalk among m6A writers, readers, and erasers in humans, and visualized GO enrichment analysis of these 20 main m6A regulators in Fig. 7D. The results show that m6A readers participate in almost all the enriched pathways.
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 (Additional file 10). For instance, the volcano plot (Fig. 8A) shows that 4102 genes (dark red dots) are significantly positively correlated with ALKBH5, whereas 2970 genes (dark green dots) show 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 (Fig. 8B, C). We used DAVID to analyze significantly enriched GO terms of the top 50 genes positively or negatively correlated with m6A regulators. The results indicate that proteins encoded by genes significantly positively correlated with m6A regulators are localized mainly to the nucleoplasm, chromosome, nucleolus, microtubule cytoskeleton, ribonucleoprotein complex, spindle, and spliceosomal complex. The molecular functions of these genes are mainly regulation of RNA binding, poly (A) RNA binding, DNA binding, and organic cyclic compound binding. They are primarily involved in biological processes including RNA biosynthetic process, RNA metabolic process, transcription, DNA-templated transcription, protein modification process, and gene expression. KEGG pathways analysis shows enrichment in the cell cycle, spliceosome, RNA transport, mRNA surveillance pathway, ubiquitin mediated proteolysis, viral carcinogenesis, and cancer pathways. However, proteins encoded by genes significantly negatively correlated with m6A regulators are 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 process, ubiquitin-dependent protein catabolic process, and modification-dependent protein catabolic process. KEGG pathways analysis shows enrichment mainly in focal adhesion, lysosome, endocytosis, proteasome, protein digestion and absorption, and metabolic pathways (Fig. 9, Additional file 11, 12).
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 are 9 intersections between pathways enriched by genes positively correlated with m6A writers and readers, but there is no intersection between pathways enriched by genes positively correlated with erasers and writers, nor erasers and readers. Simultaneously, genes negatively correlated with m6A writers are enriched to 22 pathways, readers are enriched to 17 pathways, and erasers are enriched to 15 pathways. There are 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 is 1 intersection between pathways enriched by genes positively correlated with m6A readers and genes negatively correlated with m6A erasers, but there is no intersection between pathways enriched by genes positively correlated with m6A writers and genes negatively correlated with m6A erasers (Fig. 10A–C, Additional file 13, 14, 15). KEGG pathways enrichment analysis of the top 50 genes significantly positively and negatively correlated with each m6A regulator reveals that m6A writers, erasers, and readers affect a number of signaling pathways in pancreatic cancer (Fig. 10D). These also have important roles in the occurrence and development of pancreatic cancer.