Skip to main content

Comprehensive analysis of the transcriptome‐wide m6A methylome in invasive malignant pleomorphic adenoma



Invasive malignant pleomorphic adenoma (IMPA) is a highly invasive parotid gland tumor and lacks effective therapy. N6-Methyladenosine (m6A) is the most prevalent post-transcriptional modification of mRNAs in eukaryotes and plays an important role in the pathogenesis of multiple tumors. However, the significance of m6A-modified mRNAs in IMPA has not been elucidated to date. Hence, in this study, we attempted to profile the effect of IMPA in terms of m6A methylation in mRNA.


Methylated RNA immunoprecipitation with next-generation sequencing (MeRIP-seq) and RNA sequencing (RNA-seq) were utilized to acquire the first transcriptome-wide profiling of the m6A methylome map in IMPA followed by bioinformatics analysis.


In this study, we obtained m6A methylation maps of IMPA samples and normal adjacent tissues through MeRIP-seq. In total, 25,490 m6A peaks associated with 13,735 genes were detected in the IMPA group, whereas 33,930 m6A peaks associated with 18,063 genes were detected in the control group. Peaks were primarily enriched within coding regions and near stop codons with AAACC and GGAC motifs. Moreover, functional enrichment analysis demonstrated that m6A-containing genes were significantly enriched in cancer and metabolism relevant pathways. Furthermore, we identified a relationship between the m6A methylome and the RNA transcriptome, indicating a mechanism by which m6A modulates gene expression.


Our study is the first to provide comprehensive and transcriptome-wide profiles to determine the potential roles played by m6A methylation in IMPA. These results may open new avenues for in-depth research elucidating the m6A topology of IMPA and the molecular mechanisms governing the formation and progression of IMPA.


Pleomorphic adenoma (PA), one of the most common types of salivary neoplasms, typically occurs in the superficial lobe of the parotid gland and accounts for 60–70% of all parotid tumors [1]. Although initially benign, long-standing PA has the potential to transform into malignant pleomorphic adenoma (MPA) and subsequently infiltrate surrounding facial structures, which occurs in 6% of patients [2, 3]. MPA can be further classified into three subtypes according to the depth of tumor infiltration beyond the capsule of the previous PA, namely, non-invasive, micro-invasive, and invasive malignant pleomorphic adenoma (IMPA) [4]. For IMPA, it is well accepted that surgical excision is the first-line treatment modality, but its high incidence of local recurrence and distant metastasis after surgery is a thorny problem that remains unsolved to date [5]. Besides, traditional molecular-targeted drugs against estrogen receptors (ER) or human epidermal growth factor receptor 2 (HER2) are of limited efficiency in treating IMPA, since a variety of patients are negative for both ER and HER2. To circumvent these obstacles, further research is necessary to identify a new molecular basis of IMPA.

Reversible RNA modifications, including N1-Methyladenosine (m1A), 5-Methylcytosine (m5C), and N6-Methyladenosine (m6A), have garnered substantial attention from oncologists in recent years [6]. Among these modifications, m6A is recognized as the most prevalent internal modification of poly-adenylated mRNAs and long non-coding RNAs in eukaryotic cells. Since it was first identified by Wei et al. in the 1970s, m6A has been deciphered to modulate various biological processes, such as adipogenesis, circadian rhythm and embryonic stem cell self-renewal and differentiation [7]. Mechanistically, the m6A RNA modification is enzymatically induced by m6A “writers”, a group of methyltransferases. METTL3, METTL14 and WTAP form the heterodimer core complex and mediate m6A deposition on mammalian nuclear RNAs [8, 9]. Next, m6A “readers” selectively recognize m6A-modified RNAs in the cytoplasm and exert regulatory functions by affecting multiple aspects of RNA metabolism, such as RNA stability, translation, transport and splicing. YTHDF1, YTHDC2, IGF2BP1/2/3, and HNRNPC are the primary “readers” in eukaryotic cells [10,11,12]. Notably, mammalian m6A modification is dynamic and reversible. Various demethylases, collectively known as m6A “erasers”, have been identified to participate in the process by removing N6-Methyladenosine from RNAs in the nucleus; these demethylases, include FTO and ALKBH5 [13].

Recently, m6A has been demonstrated to play a significant role in the development and progression of many tumors. Chen et al. reported the highly diverse m6A modification patterns between clear cell renal cell carcinoma (ccRCC) and tumor-adjacent normal tissue. Newly characterized m6A peaks in ccRCC have been proven to modulate gene expression and cancer-related pathways [14]. Similarly, thousands of dysregulated m6A peaks, which can be detected in colorectal cancer (CRC) samples, are associated with the prognosis of CRC patients [15]. In addition, the well-studied m6A “writer” METTL3 has been positively associated with leukemogenesis of acute myeloid leukemia and represents a therapeutic target in the disease [16]. Depletion of METTL3 in human myeloid leukemia cell lines induces differentiation and apoptosis and even delays leukemia progression in vivo [16]. Furthermore, the m6A levels of mRNAs regulated by METTL3 may also affect leukemogenesis [16]. Mutations of the two m6A “erasers” FTO and ALKBH5, as demonstrated by Cui et al., can facilitate malignant transformation of glioblastomas by promoting glioblastoma stem cell growth, self-renewal and tumor progression [17]. Another study conducted by Zhao et al. confirmed the oncogenic role of the m6A “reader” protein YTHDF1 in hepatocellular carcinoma, but further research is necessary [18].

Given that m6A modification is a dynamic and reversible process, reversing abnormal m6A patterns in tumors represents a promising strategy for targeted therapy. To date, pharmacologically targeting m6A regulators has been proven to be efficacious in both in vitro and in vivo studies [19]. Meclofenamic acid 2 (MA2) is an FDA-approved NSAID, that binds to the surface area of the FTO active site and selectively inhibits FTO. Treating glioblastoma stem cells with MA2 causes an increase in m6A levels, leading to a significant reduction in cell growth and self-renewal [17, 20]. CHTB, N-CDPCB and R-2HG are other FTO inhibitors with potential for clinical application [21]. No inhibitors of other m6A regulators are currently available, but targeting aberrant m6A modification in tumors is a feasible approach and warrants further research.

Limited knowledge is available regarding RNA m6A modification profiles in IMPA. Given the indispensable role of RNA m6A modification in the onset and progression of a wide variety of tumors, it is reasonable to surmise that the dysregulation of m6A modification might be associated with IMPA. To investigate the relationship between m6A and IMPA, we delineated the first-ever m6A transcriptome-wide landscape of IMPA. The function of the new layer of differentially methylated peaks in IMPA was also determined in our work.

Materials and methods

Specimens collections and ethics statement

Five pairs of IMPA tumors and the corresponding normal tissues were acquired at the time of surgery from Shanghai Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine between 5 Jan 2019 and 5 Oct 2019. After surgery, the candidate tissues were immediately separated into RNase-free centrifuge tubes, snap-frozen in liquid nitrogen and stored at − 80 °C until RNA isolation. Furthermore, all sample sections (4 µm) were stained with hematoxylin and eosin (H&E) and histologically reviewed by two pathologists. This study was approved by the Ethics Committee of Shanghai Ninth People’s Hospital.

MeRIP-seq and RNA-seq

Total RNA from each sample was isolated and purified using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) following the instructions of the manufacturer. Next, the amount and quality of the RNA from each sample were quantified using a NanoDrop ND-1000 (NanoDrop, Wilmington, DE, USA) instrument. RNA integrity was evaluated by a Bioanalyzer 2100 (Agilent, CA, USA) based on the threshold of RIN number > 7.0 and confirmed by electrophoresis with a denaturing agarose gel. Next, poly (A) RNAs were purified from 50 µg of total RNA by applying Dynabeads Oligo (dT) 25-61005 (Thermo Fisher, CA, USA) using two rounds of purification. Then, by using Magnesium RNA Fragmentation Module (NEB, cat.e6150, USA) at 86 ℃ for 7 min, the poly(A) RNA was fragmented into small pieces. After that step, the fragmented RNAs were incubated for 2 h at 4 ℃ in the presence of an m6A-specific antibody (No. 202003, Synaptic Systems, Germany) in IP buffer (50 mM Tris-HCl, 750 mM NaCl and 0.5 % Igepal CA-630). A portion of the initial fragmented RNA was used as the input RNA-seq library for MeRIP-seq. Next, the IP RNA was reverse-transcribed to synthesize cDNA by SuperScript™ II Reverse Transcriptase (Invitrogen, cat. 1896649 USA), which were subsequently used to create U-labeled second-stranded DNAs with E. coli DNA polymerase I (NEB, cat.m0209, USA), RNase H (NEB, cat.m0297, USA) and dUTP Solution (Thermo Fisher, cat.R0133, USA). Furthermore, an A-base was added to the blunt ends of each strand for the purpose of ligation to the indexed adapters. Each adapter contains a T-base overhang for ligating the adapter to the specific A-tailed fragmented DNA. Single- or dual-index adapters were ligated to the piece of fragments, and size selection and purity were performed with protein-A beads. After heat-labile UDG enzyme (NEB, cat.m0280, USA) treatment of the U-labeled second-stranded DNAs, the ligated products were amplified via PCR technology under the following conditions: (1) initial denaturation at 95 ℃ for 3 min; (2) 8 cycles of denaturation at 98 ℃ for 15 s, annealing at 60 ℃ for 15 s, and extension at 72 ℃ for 30 s; (3) and then final extension at 72 ℃ for 5 min. The average insert size for the final cDNA library was 300 ± 50 bp. Finally, both the input tissues without immunoprecipitation and the IP tissues were subjected to the 2 × 150-bp paired-end sequencing (PE150) on an Illumina NovaSeq™ 6000 (LC-Bio Technology Co., Ltd., Hangzhou, China) following the manufacturer’s instructions.

Bioinformatics analysis process

After the paired-end reads were obtained from Illumina NovaSeq™ 6000, the clean data were harvested using the Fastq ( package to perform quality control, including adapter-trimming, and removal of low-quality reads in fastq format. After that step, we utilized HISAT2 ( software to align the clean reads to the reference genome (version: hg38). Next, the mapped reads of IP and input libraries were provided for the R package exomePeak (, which served as m6A peak calling in bed or bigwig format that can be adapted for visualization via the Integrative Genomics Viewer (IGV) software ( Differentially methylated peaks were also identified via exomePeak according to the criteria |log2Fold Change (FC) | > 1 and p-value < 0.05. Furthermore, the MEME ( and HOMER ( packages were used for de novo and known motif finding followed by localization of the motif concerning peak summit. Peaks were annotated by the intersection with gene architecture utilizing the R package ChIPseeker ( Next, StringTie ( was applied to perform expression level for all mRNAs from input libraries by calculating FPKM (total exon fragments/mapped reads (millions) × exon length (kB)). Finally, the differentially expressed genes were identified with |log2FC| > 1 and p-value < 0.05 by the edgeR ( package in R.

Public databases and analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID) database ( was used in Gene Oncology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis for genes associated with m6A-modified peaks. Regarding the functional enrichment, the GO covered three domains: biological process (BP), cellular component (CC), and molecular function (MF). The RMBase v2.0 ( database was applied to predict the m6A modification site in BACE2. All the functional terms were regarded as significant as the threshold of p-value < 0.05.


Overall features of m6A methylation in IMPA

After the raw data were filtered and corrected, more than 35,000,000 high-quality reads acquired from each sample were mapped to the hg38 genome (Table 1). Of these high-value data, 5671 m6A methylated peaks were detected based on the pairwise comparison of read distribution between MeRIP-seq and RNA-seq samples in both tumor and corresponding normal tissues (Fig. 1a). Furthermore, 25,490 m6A modified peaks were identified in the IMPA group, which were related to 13,735 genes. Similarly, 33,930 m6A modified peaks were identified among 18,063 transcripts of genes in the control group (Fig. 1a, b). The differences and overlap in m6A modified genes were demonstrated in Fig. 1b (Additional file 1: Table S1). Among these, we discovered that only 10,046 m6A modified genes were generally expressed in both groups. The top 20 differentially modified peaks are listed in Table 2. Taken together, these results indicate an apparent discrepancy in global m6A modification patterns between IMPA and corresponding normal tissues.

Table 1 Summary of sequencing data and reads alignment statistics from MeRIP-seq and RNA-seq of IMPA tumors and corresponding normal tissues
Fig. 1
figure 1

Overview of m6A methylation modification in IMPA and adjacent normal control tissues. a, b Venn diagrams showing 5671 overlaps of m6A peaks (a) and 10,046 overlaps of m6A-modified genes (b) between IMPA and normal control tissues, respectively. c Number of peaks per transcript in IMPA and normal control tissues. d Identification of 183 hyper-methylated and 150 hypo-methylated m6A peaks in five pairs of IMPA and corresponding normal control, which presented a significant alteration, respectively (|log2 FC | > 1, p < 0.05)

Table 2 The top 20 differently methylated m6A peaks (IMPA tumor/normal)

To characterize the variation in the number of m6A-modified peaks in the transcript, we analyzed MeRIP-seq results and observed that approximately 67% of all modified transcripts had one m6A peak. To go one step further, a relatively large number of genes (93.7%) harbored 1 to 3 m6A-modified sites (Fig. 1c, Additional file 2: Table S2). Next, the abundance of m6A peaks between IMPA and adjacent normal control samples was analyzed with the norm of |log2FC| > 1 and p < 0.05. A total of 183 hyper-methylated peaks and 150 hypo-methylated peaks were depicted in IMPA samples compared with the control groups (Fig. 1d, Additional file 3: Table S3).

Topological patterns of IMPA methylation

The topological patterns of m6A modification in IMPA was investigated, the results presented in Fig. 2a demonstrate that the distribution of m6A peaks was principally enriched in coding sequence (CDS) regions close to the stop codon regardless of the IMPA tumor or normal tissue origin. For a more precise evaluation of m6A peaks, both m6A-modified peaks in the general and unique groups were analyzed according to the transcript’s location. Overall, the m6A sites were primarily deposited in the CDS, 3′-UTR, and stop codon (Fig. 2b i). However, the IMPA groups exhibited a distinct pattern from the control group with a relative increase in m6A enrichment harbored in the vicinity of the CDS region but less m6A enrichment near the stop codon (Fig. 2b ii, 2b iii). Besides, the differentially methylated m6A sites in the two groups most often appeared in the CDS (50.9%), followed by the 3′-UTR (20.6 %) (Fig. 2b iv). The altered m6A peaks were further observed to be mapped to all human chromosomes. The top three chromosomes sharing the most hyper-methylated and hypo-methylated m6A peaks within mRNA were chromosomes 1, 16, and 19, and hyper-methylated m6A peaks were observed in all chromosome but chromosome 4 (Fig. 2c, Additional file 4: Table S4). Additionally, the HOMER package was applied to search for consensus motifs deposited in the region surrounding the m6A peak and to identify the classical AAACH and GGAC motif structures (Fig. 2d, e), which was in accord with the previous findings. The identification of a strong consensus strengthens the authenticity of the identified m6A peaks and suggest the existence of a dominant methylation mechanism.

Fig. 2
figure 2

Distribution and topology pattern of m6A peaks along transcripts and chromatin in IMPA and the corresponding normal tissues. a Difference in the density of m6A peaks along the transcript. Each transcript was divided into three parts: 5’-UTR, CDS, and 3’-UTR. b Pie charts showing the distributions of m6A peaks in the whole transcriptome of IMPAs and controls (i), m6A peaks in IMPAs (ii), m6A peaks in controls (iii), and differentially methylated m6A peaks between IMPAs and controls (iv). c The distribution of differentially methylated m6A peaks with significance along the chromosome of IMPA. d, e Sequencing logo displaying AAACC and GGAC conserved sequence motifs for the m6A peak regions, with p-values of 1*e− 301 and 1*e− 203, respectively

m6A-containing genes are involved in important biological processes and pathways

In this study, transcriptome features of the abnormal m6A-modified genes in IMPA were determined by MeRIP-seq. To investigate the biological significance of m6A-containing genes in IMPA, the hyper-methylated transcripts were subjected the GO and KEGG analyses. With GO analysis, genes were divided into three functional domains: biological process (BP), cellular component (CC), and molecular function (MF). Figure 3a (Additional file 5: Table S5) lists the top 25 BP terms, top 15 CC terms, and top 10 MF terms of the hyper-methylated genes which indicated that those genes were significantly enriched in ATP transmembrane transporter activity, antiporter activity, and nucleic acid binding (Fig. 3b, Additional file 6: Table S6). With respected to the KEGG pathways of genes with up-methylated m6A peaks, we found that the hyper-methylated mRNAs were enriched in biotin metabolism, morphine addiction, nicotine addiction, aldosterone synthesis and secretion, the PPAR signaling pathway, and the mTOR signaling pathway (Fig. 3c, Additional file 7: Table S7). To sum up, the enrichment results may serve to increase the scientific interest in RNA m6A modification in IMPA.

Fig. 3
figure 3

GO function and KEGG pathway enrichment of hyper-methylated m6A genes. a Major enrichment and meaningful GO terms of hyper-methylated m6A genes in IMPA. b The top twenty significant GO enrichment terms. c The top twenty significant KEGG pathways

Overview of transcriptome profiles and conjoint analysis of m6A methylation and gene expression

To detect the transcriptome profiles of IMPA samples versus the corresponding normal tissues from five pairs of IMPA patients, we employed RNA-seq (MeRIP-seq input library) to quantify the differential gene expression in the 2 groups. Compared with normal control samples, 3028 genes were differentially expressed in IMPA (|log2FC| > 1, and p < 0.05), including 2178 up-regulated genes and 850 down-regulated genes, respectively (Fig. 4a, Additional file 8: Table S8). The heatmap of the differentially expressed genes between IMPA and control is shown in Fig. 4b (Additional file 9: Table S9). Next, according to the conjoint analysis of RNA-seq and MeRIP-seq data, all genes were classified into 4 parts, including 24 up-regulated genes that were significantly hyper-methylated (10; up-hyper) or hypo-methylated (14; up-hypo) and 7 down-regulated genes that were significantly hyper-methylated (2; down-hyper) or hypo-methylated (5; down-hypo, |log2FC| > 1, and p < 0.05, Fig. 4c, Additional file 10: Table S10). These findings necessitate further exploration of the correlation of differentially methylated genes and differentially expressed gene levels in IMPA and normal samples. As illustrated in Fig. 4d (Additional file 11: Table S11), the genes were classed evenly into 10 groups according to the degree of expression, from low to high, the m6A density of each category was determined separately. An approximately positive correlation between m6A methylation and gene expression level before the 0.6 thresholds was determined, while the converse trend was observed for genes whose expression is over 0.6.

Fig. 4
figure 4

Comprehensive analysis of the m6A methylome and RNA transcriptome in IMPA. a, b Scatter plot and hierarchical clustering presenting the differentially expressed genes in IMPA. c Four quadrant graph exhibiting the distribution of mRNAs with a significant alteration in both m6A and mRNA levels between IMPA and adjacent normal tissues (|log2 FC| > 1, p < 0.05). d Estimation of m6A peaks density in IMPA. e Expression levels of genes containing m6A-specific peaks in different five regions of transcript between IMPA and normal tissues

Furthermore, we also analyzed the expression of genes harboring differentially methylated m6A peaks in the five transcriptive regions between IMPA and control samples. In comparison with the control groups, the higher expression of mRNA specific to IMPA tissues in all five regions is also depicted in Fig. 4e (Additional file 12: Table S12). These results lay a solid groundwork for an in-depth research on m6A modifications in IMPA.


As an emerging field of post-transcriptional gene regulation, m6A modification has attracted extensive attention and interest in the research community. Studies have identified transcriptome-wide maps of m6A in multiple biological and pathological conditions, including maternal-to-zygotic transition, heart failure, and neurodegenerative disease [22,23,24]. It is noteworthy that m6A modification also plays an essential role in the development and drug responses of various cancers [25]. Through decades of research, global mapping of m6A in most kinds of cancers (e.g., hepatocellular carcinoma, glioblastoma, and breast cancer) has been achieved [25]. Nonetheless, to date, few studies on head and neck tumors have been conducted, and the distribution of m6A modification in IMPA at the transcriptome-wide level has not been fully elucidated.

In this study, we established an antibody-based methodology that presented transcriptome-wide m6A modification patterns in IMPA samples versus tumor-adjacent normal tissues, and we subsequently analyzed the alterations to gene expression and cancer-related pathways caused by aberrant m6A modifications. To the best of our knowledge, this report describes the first-ever high-throughput sequencing to comprehensively investigate RNA methylation in IMPA. Our data identified broad mRNA methylation and demethylation features over one of the most guarded prognosis subtypes, IMPA.

We observed abundant m6A methylation sites in the IMPA transcriptome by using the MeRIP-seq technique. The density was determined to be approximately 1.57 sites per transcript, which was significantly higher than that of sheep and pig, which exhibited approximately < 0.5 and < 0.6 m6A sites per transcript (Fig. 1c), respectively [26, 27]. Taken together, the m6A distribution density results imply that epigenetic modification may vary between different organisms or may be related to species differentiation. It is interesting to note that in all IMPA samples of this study, the m6A sites were also enriched around the CDS and in the vicinity of stop codons, close to the topological pattern of modification in colorectal cancer (Fig. 2b ii) [15]. As mentioned in the literature review, internal exon methylation would take effect on the alternative splicing. At the same time, translational regulation can be influenced by methylation near the stop codon [28]. Therefore, it can be assumed that m6A could play a vital role in controlling protein synthesis in IMPA. These findings, in combination with the results presented in Fig. 4d, offer some conceptual hypotheses in which m6A could regulate gene expression to a certain extent.

The most prominent finding to emerge from our data is that the difference in m6A modification between the tumor tissues and corresponding normal tissues was significant, which was also observed in breast cancer, non-small cell lung cancer, and pancreatic adenocarcinoma with different methylation levels [29,30,31]. This result suggests that the m6A modification could contribute to IMPA progression. A possible explanation for this global alteration of m6A modification patterns might be the unique expression of the key m6A regulators. Surprisingly and contrary to expectations, our current study did not detect significant differences among the key m6A regulators (writers, readers, and erasers) based on RNA-sEq. This result might be due to the relatively small sample size of our study or post-translational modifications [32]. Subsequent RNA-seq studies in IMPA are warranted to determine the actual reason.

RRACH and GGAC are proved to be over-represented motifs of m6A that are strongly preserved in humans [33,34,35]. Accordingly, in our current study, we identified multiple similar m6A consensus motifs in IMPA (Fig. 2d, e). These analyzed results were in keeping with the findings of a large number of previous studies in oral carcinomas, which suggested that we called specific m6A sites successfully [36].

As the most prevalent epigenetic modification in RNA, m6A has been proven to function as an essential mediator in tumorigenesis [37, 38]. Recently, Chen et al. reported that SOCS2 m6A methylation could contribute to its mRNA expression and liver carcinogenesis [39]. Additionally, the methylation of FOXMs serves to maintain the tumorigenicity of glioblastoma [40]. Indeed, in our study, we successfully identified a significant hyper-methylated characteristic in BACE2 (Table 2). Interestingly, in glioblastoma, patients with higher BACE2 expression exhibited an unfavorable outcome [41]. The GSEA database indicated that BACE2 could be involved in cell invasion and cell migration in gliomas [41]. Besides, compared with normal tissue, BACE2 was also observed to be over-expressed in breast and colon tumors [42, 43]. Furthermore, according to the post-transcriptional RNA modifications database, RMBase v2.0, BACE2 also shares an m6A methylation motif (Additional file 13: Fig. S1). Therefore, it may suggest that BACE2 affects IMPA progression via m6A modification manners. Further study is necessary to elucidate whether BACE2 methylation triggers the IMPA oncogenic process.

Prior studies have demonstrated that signaling pathways are significant targets to modulate tumor progression [44,45,46]. In our study, many critical biological pathways related to differentially methylated genes were enriched via GO, such as ATP transmembrane transporter activity, antiporter activity, and nucleic acid binding. Studies have revealed that these biological functions are closely implicated in the onset and development of tumors, including tumor growth, metastasis and chemo-resistance. For instance, high expression of ATP-binding cassette transmembrane transporters results in the efflux of cytotoxic agents from cancer cells, which leads to drug resistance. Qiu et al. demonstrated that the elevated expression of cellular nucleic acid-binding protein in fibrosarcoma cells down-regulated the transcription and expression of hnRNP K, which enhanced tumor cell death and suppressed tumor cell metastasis [47]. Based on KEGG analysis, genes with upregulated m6A modification were associated with mTOR signaling pathways in IMPA. It is intriguing to note that PI3K/AKT/mTOR pathways would contribute to cancer cell proliferation in head and neck squamous cell carcinoma (HNSCC), possibly serving a pro-tumorigenic function [48]. In accordance with the above HNSCC study, the most common malignant salivary gland neoplasm, adenoid cystic carcinoma, exhibits similar characteristics in the PI3K/AKT/mTOR axis [49]. One therapeutic method for tumors was to target the PI3K/AKT/mTOR pathways with an inhibitor [50]. Hence, it could be hypothesized that modulating the m6A modification of the transcripts of genes in the mTOR signaling pathways might provide a novel means of targeting IMPA for treatment.

Additionally, the carbohydrate and lipid metabolism have always been closely related to cancer cell proliferation and migration [51]. The KEGG pathway enrichment of coding genes harboring the hyper-methylated peaks in IMPA samples was related to the PPAR signaling pathways (Fig. 3c). The PPAR pathway may be involved in transcriptional regulation in many pathophysiological processes, which not only participate in fat and sugar metabolism, but also play roles in tumor development [52,53,54,55]. Moreover, researches have indicated a close association between the alteration of m6A modification and lipid metabolism, thereby, providing a link between m6A modification and IMPA tumorigenesis through metabolic ways [56,57,58]. Further studies, seeking to elucidate the potential mechanisms underlying PPAR signaling pathways in IMPA, are warranted.

In summary, we analyzed the differential m6A methylome in IMPA relative to the corresponding normal controls, and demonstrated a strong association between m6A modification and IMPA progression. Our study paves the way for future investigations aimed to explore the mechanism in the pathogenesis of IMPA.

Availability of data and materials

The datasets generated and/or analyzed during the current study are from the DAVID ( and RMBase v2.0 ( MeRIP-seq and RNA-seq were uploaded to the GEO database (GSE161879), and can be accessed with the secure token key ejsdoiaujpehxyb.



Pleomorphic adenoma


Malignant pleomorphic adenoma


Invasive malignant pleomorphic adenoma


Methylated RNA immunoprecipitation with next-generation sequencing




RNA sequencing






Head and neck squamous cell carcinoma


Gene Oncology


Kyoto Encyclopedia of Genes and Genomes


Biological process


Cellular component


Molecular function


The Database for Annotation Visualization and Integrated Discovery


Integrative Genomics Viewer


  1. Erickson LA. Pleomorphic adenoma in the parotid gland. Mayo Clin Proc. 2017;92(3):e55–6.

    Article  PubMed  Google Scholar 

  2. Xia L, Wang Y, Hu Y, Zhang C, Gu T, Wang L, Li J, Yu W, Tian Z. Human epithelial growth factor receptor 2 in human salivary carcinoma ex pleomorphic adenoma: a potential therapeutic target. Cancer Manag Res. 2018;10:6571–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Hashimoto K, Yamamoto H, Shiratsuchi H, Nakashima T, Tamiya S, Nishiyama K, Higaki Y, Komune S, Tsuneyoshi M, Oda Y. HER-2/neu gene amplification in carcinoma ex pleomorphic adenoma in relation to progression and prognosis: a chromogenic in-situ hybridization study. Histopathology. 2012;60(6b):E131–42.

    Article  PubMed  Google Scholar 

  4. Hu Y, Xia L, Zhang C, Xia R, Tian Z, Li J. Clinicopathologic features and prognostic factors of widely invasive carcinoma ex pleomorphic adenoma of parotid gland: a clinicopathologic analysis of 126 cases in a Chinese population. J Oral Maxillofac Surg. 2020;78(12):2247–57.

    Article  PubMed  Google Scholar 

  5. Abu-Ghanem Y, Mizrachi A, Popovtzer A, Abu-Ghanem N, Feinmesser R. Recurrent pleomorphic adenoma of the parotid gland: Institutional experience and review of the literature. J Surg Oncol. 2016;114(6):714–8.

    Article  PubMed  Google Scholar 

  6. Huang AZ, Delaidelli A, Sorensen PH. RNA modifications in brain tumorigenesis. Acta Neuropathol Commun. 2020;8(1):64.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Karthiya R, Khandelia P. m6A RNA methylation: ramifications for gene expression and human health. Mol Biotechnol. 2020;62(10):467–84.

    Article  CAS  PubMed  Google Scholar 

  8. Liu J, Yue Y, Han D, Wang X, Fu Y, Zhang L, Jia G, Yu M, Lu Z, Deng X, et al. A METTL3-METTL14 complex mediates mammalian nuclear RNA N6-adenosine methylation. Nat Chem Biol. 2014;10(2):93–5.

    Article  CAS  PubMed  Google Scholar 

  9. Ping XL, Sun BF, Wang L, Xiao W, Yang X, Wang WJ, Adhikari S, Shi Y, Lv Y, Chen YS, et al. Mammalian WTAP is a regulatory subunit of the RNA N6-methyladenosine methyltransferase. Cell Res. 2014;24(2):177–89.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  10. Wang X, Zhao BS, Roundtree IA, Lu Z, Han D, Ma H, Weng X, Chen K, Shi H, He C. N(6)-methyladenosine modulates messenger RNA translation efficiency. Cell. 2015;161(6):1388–99.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Hsu PJ, Zhu Y, Ma H, Guo Y, Shi X, Liu Y, Qi M, Lu Z, Shi H, Wang J, et al. Ythdc2 is an N(6)-methyladenosine binding protein that regulates mammalian spermatogenesis. Cell Res. 2017;27(9):1115–27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Alarcón CR, Goodarzi H, Lee H, Liu X, Tavazoie S, Tavazoie SF. HNRNPA2B1 is a mediator of m(6)A-dependent nuclear RNA processing events. Cell. 2015;162(6):1299–308.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  13. Jia G, Fu Y, Zhao X, Dai Q, Zheng G, Yang Y, Yi C, Lindahl T, Pan T, Yang YG, et al. N6-methyladenosine in nuclear RNA is a major substrate of the obesity-associated FTO. Nat Chem Biol. 2011;7(12):885–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Chen Y, Zhou C, Sun Y, He X, Xue D. m(6)A RNA modification modulates gene expression and cancer-related pathways in clear cell renal cell carcinoma. Epigenomics. 2020;12(2):87–99.

    Article  PubMed  CAS  Google Scholar 

  15. Zhang Z, Wang Q, Zhang M, Zhang W, Zhao L, Yang C, Wang B, Jiang K, Ye Y, Shen Z, et al. Comprehensive analysis of the transcriptome-wide m6A methylome in colorectal cancer by MeRIP sequencing. Epigenetics. 2020.

    Article  PubMed  PubMed Central  Google Scholar 

  16. Vu LP, Pickering BF, Cheng Y, Zaccara S, Nguyen D, Minuesa G, Chou T, Chow A, Saletore Y, MacKay M, et al. The N(6)-methyladenosine (m(6)A)-forming enzyme METTL3 controls myeloid differentiation of normal hematopoietic and leukemia cells. Nat Med. 2017;23(11):1369–76.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Cui Q, Shi H, Ye P, Li L, Qu Q, Sun G, Sun G, Lu Z, Huang Y, Yang CG, et al. m(6)A RNA methylation regulates the self-renewal and tumorigenesis of glioblastoma stem cells. Cell Rep. 2017;18(11):2622–34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Zhao X, Chen Y, Mao Q, Jiang X, Jiang W, Chen J, Xu W, Zhong L, Sun X. Overexpression of YTHDF1 is associated with poor prognosis in patients with hepatocellular carcinoma. Cancer Biomark. 2018;21(4):859–68.

    Article  CAS  PubMed  Google Scholar 

  19. Lan Q, Liu PY, Haase J, Bell JL, Hüttelmaier S, Liu T. The critical role of RNA m(6)A methylation in cancer. Cancer Res. 2019;79(7):1285–92.

    Article  CAS  PubMed  Google Scholar 

  20. Huang Y, Yan J, Li Q, Li J, Gong S, Zhou H, Gan J, Jiang H, Jia GF, Luo C, et al. Meclofenamic acid selectively inhibits FTO demethylation of m6A over ALKBH5. Nucleic Acids Res. 2015;43(1):373–84.

    Article  CAS  PubMed  Google Scholar 

  21. He W, Zhou B, Liu W, Zhang M, Shen Z, Han Z, Jiang Q, Yang Q, Song C, Wang R, et al. Identification of a novel small-molecule binding site of the fat mass and obesity associated protein (FTO). J Med Chem. 2015;58(18):7341–8.

    Article  CAS  PubMed  Google Scholar 

  22. Tang J, Zheng C, Zheng F, Li Y, Wang YL, Aschner M, Guo Z, Yu G, Wu S, Li H. Global N6-methyladenosine profiling of cobalt-exposed cortex and human neuroblastoma H4 cells presents epitranscriptomics alterations in neurodegenerative disease-associated genes. Environ Pollut. 2020;266(Pt 2):115326.

    Article  CAS  PubMed  Google Scholar 

  23. Wu S, Zhang S, Wu X, Zhou X. m(6)A RNA methylation in cardiovascular diseases. Mol Ther. 2020;28(10):2111–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Zhao BS, Wang X, Beadell AV, Lu Z, Shi H, Kuuspalu A, Ho RK, He C. m(6)A-dependent maternal mRNA clearance facilitates zebrafish maternal-to-zygotic transition. Nature. 2017;542(7642):475–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Deng X, Su R, Weng H, Huang H, Li Z, Chen J. RNA N(6)-methyladenosine modification in cancers: current status and perspectives. Cell Res. 2018;28(5):507–17.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Lu Z, Liu J, Yuan C, Jin M, Quan K, Chu M, Wei C. m(6)A mRNA methylation analysis provides novel insights into heat stress responses in the liver tissue of sheep. Genomics. 2020;113(1):484–92.

    Article  PubMed  CAS  Google Scholar 

  27. Tao X, Chen J, Jiang Y, Wei Y, Chen Y, Xu H, Zhu L, Tang G, Li M, Jiang A, et al. Transcriptome-wide N (6) -methyladenosine methylome profiling of porcine muscle and adipose tissues reveals a potential mechanism for transcriptional regulation and differential methylation pattern. BMC Genom. 2017;18(1):336.

    Article  CAS  Google Scholar 

  28. Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, Cesarkas K, Jacob-Hirsch J, Amariglio N, Kupiec M, et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-sEq. Nature. 2012;485(7397):201–6.

    Article  CAS  PubMed  Google Scholar 

  29. Chang G, Shi L, Ye Y, Shi H, Zeng L, Tiwary S, Huse JT, Huo L, Ma L, Ma Y, et al. YTHDF3 induces the translation of m(6)A-enriched gene transcripts to promote breast cancer brain metastasis. Cancer Cell. 2020;38(6):857–71.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Sun S, Han Q, Liang M, Zhang Q, Zhang J, Cao J. Downregulation of m(6) A reader YTHDC2 promotes tumor progression and predicts poor prognosis in non-small cell lung cancer. Thorac Cancer. 2020;11(11):3269–79.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Tang R, Zhang Y, Liang C, Xu J, Meng Q, Hua J, Liu J, Zhang B, Yu X, Shi S. The role of m6A-related genes in the prognosis and immune microenvironment of pancreatic adenocarcinoma. PeerJ. 2020;8:e9602.

    Article  PubMed  PubMed Central  Google Scholar 

  32. Sheng Q, Wang C, Li X, Qin H, Ye M, Xiong Y, Wang X, Li X, Lan M, Li J, et al. Highly efficient separation of methylated peptides utilizing selective complexation between lysine and 18-Crown-6. Anal Chem. 2020;92(23):15663.

    Article  CAS  PubMed  Google Scholar 

  33. Fu Y, Dominissini D, Rechavi G, He C. Gene expression regulation mediated through reversible m6 RNA methylation. Nat Rev Genet. 2014;15(5):293–306.

    Article  CAS  PubMed  Google Scholar 

  34. Wen K, Zhang Y, Li Y, Wang Q, Sun J. Comprehensive analysis of transcriptome-wide m(6)A methylome in the anterior capsule of the lens of high myopia patients. Epigenetics 2020:1–14.

  35. Li Z, Weng H, Su R, Weng X, Zuo Z, Li C, Huang H, Nachtergaele S, Dong L, Hu C, et al. FTO plays an oncogenic role in acute myeloid leukemia as a N(6)-Methyladenosine RNA demethylase. Cancer Cell. 2017;31(1):127–41.

    Article  PubMed  CAS  Google Scholar 

  36. Liu L, Wu Y, Li Q, Liang J, He Q, Zhao L, Chen J, Cheng M, Huang Z, Ren H, et al. METTL3 Promotes tumorigenesis and metastasis through BMI1 m(6)A methylation in oral squamous cell carcinoma. Mol Ther. 2020;28(10):2177–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Zhao Y, Feng F, Guo QH, Wang YP, Zhao R. Role of succinate dehydrogenase deficiency and oncometabolites in gastrointestinal stromal tumors. World J Gastroenterol. 2020;26(34):5074–89.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  38. Wu KJ. The role of miRNA biogenesis and DDX17 in tumorigenesis and cancer stemness. Biomed J. 2020;43(2):107–14.

    Article  PubMed  PubMed Central  Google Scholar 

  39. Chen M, Wei L, Law CT, Tsang FH, Shen J, Cheng CL, Tsang LH, Ho DW, Chiu DK, Lee JM, et al. RNA N6-methyladenosine methyltransferase-like 3 promotes liver cancer progression through YTHDF2-dependent posttranscriptional silencing of SOCS2. Hepatology. 2018;67(6):2254–70.

    Article  CAS  PubMed  Google Scholar 

  40. Zhang S, Zhao BS, Zhou A, Lin K, Zheng S, Lu Z, Chen Y, Sulman EP, Xie K, Bögler O, et al. m(6)A Demethylase ALKBH5 maintains tumorigenicity of glioblastoma stem-like cells by sustaining FOXM1 Eexpression and cell proliferation program. Cancer Cell. 2017;31(4):591–606.e596.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Wang H, Chen Z, Wang S, Gao X, Qian M, Qiu W, Zhang Z, Zhang S, Qi Y, Sun X, et al. TGFβ1-induced beta-site APP-cleaving enzyme 2 upregulation promotes tumorigenesis through the NF-κB signalling pathway in human gliomas. Mol Oncol. 2020;14(2):407–25.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Tsuji N, Kondoh K, Furuya M, Kobayashi D, Yagihashi A, Inoue Y, Meguro T, Horita S, Takahashi H, Watanabe N. A novel aspartate protease gene, ALP56, is related to morphological features of colorectal adenomas. Int J Colorectal Dis. 2004;19(1):43–8.

    Article  PubMed  Google Scholar 

  43. Kondoh K, Tsuji N, Kamagata C, Sasaki M, Kobayashi D, Yagihashi A, Watanabe N. A novel aspartic protease gene, ALP56, is up-regulated in human breast cancer independently from the cathepsin D gene. Breast Cancer Res Treat. 2003;78(1):37–44.

    Article  CAS  PubMed  Google Scholar 

  44. Sun SY. mTOR-targeted cancer therapy: great target but disappointing clinical outcomes, why? Front Med. 2020.

    Article  PubMed  PubMed Central  Google Scholar 

  45. Khongthong P, Roseweir AK, Edwards J. The NF-KB pathway and endocrine therapy resistance in breast cancer. Endocr Relat Cancer. 2019;26(6):R369-r380.

    Article  PubMed  Google Scholar 

  46. Sha YL, Liu S, Yan WW, Dong B. Wnt/β-catenin signaling as a useful therapeutic target in hepatoblastoma. Biosci Rep 2019, 39(9).

  47. Qiu J, Chen S, Su L, Liu J, Xiao N, Ou TM, Tan JH, Gu LQ, Huang ZS, Li D. Cellular nucleic acid binding protein suppresses tumor cell metastasis and induces tumor cell death by downregulating heterogeneous ribonucleoprotein K in fibrosarcoma cells. Biochim Biophys Acta. 2014;1840(7):2244–52.

    Article  CAS  PubMed  Google Scholar 

  48. Bernard M, Cardin GB, Cahuzac M, Ayad T, Bissada E, Guertin L, Bahig H, Nguyen-Tan PF, Filion E, Ballivy O, et al. Dual inhibition of autophagy and PI3K/AKT/MTOR pathway as a therapeutic strategy in head and neck squamous cell carcinoma. Cancers. 2020;12(9):2371.

    Article  CAS  PubMed Central  Google Scholar 

  49. Zboray K, Mohrherr J, Stiedl P, Pranz K, Wandruszka L, Grabner B, Eferl R, Moriggl R, Stoiber D, Sakamoto K, et al. AKT3 drives adenoid cystic carcinoma development in salivary glands. Cancer Med. 2018;7(2):445–53.

    Article  CAS  PubMed  Google Scholar 

  50. Gibbons JJ, Abraham RT, Yu K. Mammalian target of rapamycin: discovery of rapamycin reveals a signaling pathway important for normal and cancer cell growth. Semin Oncol. 2009;36(Suppl 3):3-s17.

    Google Scholar 

  51. Vykoukal J, Fahrmann JF, Gregg JR, Tang Z, Basourakos S, Irajizad E, Park S, Yang G, Creighton CJ, Fleury A, et al. Caveolin-1-mediated sphingolipid oncometabolism underlies a metabolic vulnerability of prostate cancer. Nat Commun. 2020;11(1):4279.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Lefebvre AM, Chen I, Desreumaux P, Najib J, Fruchart JC, Geboes K, Briggs M, Heyman R, Auwerx J. Activation of the peroxisome proliferator-activated receptor gamma promotes the development of colon tumors in C57BL/6J-APCMin/+ mice. Nat Med. 1998;4(9):1053–7.

    Article  CAS  PubMed  Google Scholar 

  53. Rosen ED, Sarraf P, Troy AE, Bradwin G, Moore K, Milstone DS, Spiegelman BM, Mortensen RM. PPAR gamma is required for the differentiation of adipose tissue in vivo and in vitro. Mol Cell. 1999;4(4):611–7.

    Article  CAS  PubMed  Google Scholar 

  54. Martin G, Schoonjans K, Staels B, Auwerx J. PPARgamma activators improve glucose homeostasis by stimulating fatty acid uptake in the adipocytes. Atherosclerosis. 1998;137(Suppl):75–80.

    Article  Google Scholar 

  55. Panasyuk G, Espeillac C, Chauvin C, Pradelli LA, Horie Y, Suzuki A, Annicotte JS, Fajas L, Foretz M, Verdeguer F, et al. PPARγ contributes to PKM2 and HK2 expression in fatty liver. Nat Commun. 2012;3:672.

    Article  PubMed  CAS  Google Scholar 

  56. Wu R, Yao Y, Jiang Q, Cai M, Liu Q, Wang Y, Wang X. Epigallocatechin gallate targets FTO and inhibits adipogenesis in an mRNA m(6)A-YTHDF2-dependent manner. Int J Obes. 2018;42(7):1378–88.

    Article  CAS  Google Scholar 

  57. Mo X, Lei S, Zhang Y, Zhang H. Genome-wide enrichment of m(6)A-associated single-nucleotide polymorphisms in the lipid loci. Pharmacogenom J. 2019;19(4):347–57.

    Article  CAS  Google Scholar 

  58. Li X, Yang J, Zhu Y, Liu Y, Shi X, Yang G. Mouse maternal high-fat intake dynamically programmed mRNA m(6)A modifications in adipose and skeletal muscle tissues in offspring. Int J Mol Sci. 2016;17(8):1336.

    Article  PubMed Central  CAS  Google Scholar 

Download references


We thank lc-bio Ltd. Co. (Hangzhou, China) for the MeRIP-seq service and the subsequent bioinformatics analysis.


The work was supported by the Natural Science Foundation of Shanghai (Grant No. 18ZR1422200).

Author information

Authors and Affiliations



ZT and YQW conceived and developed the outline of the study, and reviewed the draft of the paper. ZH, BY, and QW performed the experiments, provided the statistical analysis, and wrote the manuscript. YHH collected the IMPA samples. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Yuqiong Wu or Zhen Tian.

Ethics declarations

Ethics approval and consent to participate

This study was approved by the Ethics Committee of Shanghai Ninth People’s Hospital.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

Table S1. The differences and overlap of m6A modified genes in IMPA and corresponding normal tissues.

Additional file 2.

Table S2. Number of peaks per transcript in IMPA and normal control tissues.

Additional file 3.

Table S3. The estimation of hyper- and hypo-methylated peaks in IMPA samples compared with the control groups.

Additional file 4.

Table S4. Number of altered m6A peaks in human chromosomes.

Additional file 5.

Table S5. The enriched GO terms of the hyper-methylated genes.

Additional file 6.

Table S6. The enriched GO terms of the hyper-methylated genes.

Additional file 7.

Table S7. The enriched KEGG pathways of the hyper-methylated genes.

Additional file 8.

Table S8. The distribution of genes with significant changes in mRNA levels.

Additional file 9.

Table S9. The differentially expressed genes in IMPA and normal control tissues.

Additional file 10.

Table S10. The distribution of genes with significant changes in both m6A modification and mRNA levels.

Additional file 11.

Table S11. The estimation of m6A peaks density in IMPA.

Additional file 12.

Table S12. Expression levels of genes containing m6A-specific peaks in different five regions of transcript between IMPA and normal tissues.

Additional file 13: Fig. S1.

The RMBase v2.0 database shows the m6A modification motif in BACE2 based on GSE37003. motif sequence: GGACU.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Han, Z., Yang, B., Wang, Q. et al. Comprehensive analysis of the transcriptome‐wide m6A methylome in invasive malignant pleomorphic adenoma. Cancer Cell Int 21, 142 (2021).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • Malignant pleomorphic adenoma
  • N6-Methyladenosine
  • MeRIP-seq