Transcriptome-wide m6A methylome analysis uncovered the changes of m6A modification in oral pre-malignant cells compared with normal oral epithelial cells

As the most common post-transcriptional RNA modification, m⁶A methylation extensively regulates the structure and function of RNA. The dynamic and reversible modification of m⁶A is coordinated by m⁶A writers and erasers. m⁶A reader proteins recognize m⁶A modification on RNA, mediating different downstream biological functions. mRNA m⁶A modification and its corresponding regulators play an important role in cancers, but its characteristics in the precancerous stage are still unclear. In this study, we used oral precancerous DOK cells as a model to explore the characteristics of transcriptome-wide m⁶A modification and major m⁶A regulator expression in the precancerous stage compared with normal oral epithelial cell HOEC and oral cancer cell SCC-9 through MeRIP-seq and RT-PCR. Compared with HOEC cells, we found 1180 hyper-methylated and 1606 hypo-methylated m⁶A peaks and 354 differentially expressed mRNAs with differential m⁶A peaks in DOK cells. Although the change of m⁶A modification in DOK cells was less than that in SCC-9 cells, mRNAs with differential m⁶A in both cell lines were enriched into many identical GO terms and KEGG pathways. Among the 20 known m⁶A regulatory genes, FTO, ALKBH5, METTL3 and VIRMA were upregulated or downregulated in DOK cells, and the expression levels of 10 genes such as METTL14/16, FTO and IGF2BP2/3 were significantly changed in SCC-9 cells. Our data suggest that precancerous cells showed, to some extent, changes of m⁶A modification. Identifying some key m⁶A targets and corresponding regulators in precancerous stage may provide potential intervention targets for the prevention of cancer development through epigenetic modification in the future.

Introduction

N⁶-methyladenine (m⁶A) methylation is a dynamic and reversible modification regulated by methyltransferases and demethylases in different types of RNAs including messenger RNA (mRNA). This is the most prevalent post-transcriptional regulatory markers on eukaryotic RNAs and one third of the total mammalian mRNA has 3-5 m⁶A modifications in each mRNA (1, 2). Most m⁶A sites are located in the conserved motif DRACH (D=G/A/U, R=G/A, H=A/U/C), which are usually found in 3′ UTRs and near stop codons in mRNAs (1, 2).

m⁶A is decorated by m⁶A methyltransferase complex, and its components include METTL3/14/16, WTAP, RBM15, RBM15B, VIRMA, and ZC3H13. They are collectively referred to as m⁶A writers. On the contrary, demethylases such as FTO, ALKBH5 and ALKBH3 can remove m⁶A and act as erasers. The cooperation of m⁶A writers and erasers regulates the dynamic and reversible modification of m⁶A. m⁶A modification of RNAs must be recognized by m⁶A reader proteins to mediate different downstream biological functions. Several categories of proteins function as m⁶A reader, including YT521-B homology (YTH) domain-containing proteins YTHDC1/2 and YTHDF1/2/3, IGF2 mRNA binding proteins IGF2BP1/2/3, heterogeneous nuclear ribonucleoproteins HNRNPA2B1 and HNRNPC, and eukaryotic initiation factor 3 (eIF3). RNA m⁶A modification plays an important role in regulating RNA splicing, translation, stability, translocation, and advanced structure (3, 4).

There is increasing evidence that RNA m⁶A methylation and its regulators are associated with a variety of human diseases, especially cancer (3–5). m⁶A regulators may act as writers to catalyze or erasers to remove m⁶A modifications in the mRNAs of oncogenes and tumor suppressor genes, which are then recognized by readers to regulate the expression of these genes and their downstream biological functions. The role of m⁶A modification involves many aspects of tumors. For example, demethylase FTO-mediated m⁶A demethylation of cytidine deaminase APOBEC3B mRNA promotes arsenic-induced mutagenesis (6). FTO promotes growth and metastasis of gastric cancer through m⁶A demethylation of caveolin-1 and metabolic regulation of mitochondrial dynamics (7). FTO-mediated downregulation of DACT1 mRNA stability promotes Wnt signaling to facilitate osteosarcoma progression (8). Downregulation of FTO promotes EMT-mediated progression of epithelial tumors and sensitivity to Wnt inhibitors (9). FTO promotes hepatocellular carcinoma tumorigenesis through mediating PKM2 demethylation (10). m⁶A writer METTL3 promotes chemo- and radioresistance in pancreatic cancer cells (11). METTL3 stabilizes HK2 and SLC2A1 (GLUT1) expression in colorectal cancer and m⁶A-dependent glycolysis enhances colorectal cancer progression (12). METTL3 promotes tumor development by decreasing APC gene expression mediated by APC mRNA m⁶A-dependent YTHDF binding (13).

The significance of m⁶A modification in head and neck squamous cell carcinoma (HNSCC) has been well reviewed (4). m⁶A plays an important role in the tumorigenesis, drug resistance and prognosis of oral squamous cell carcinoma (OSCC). Methyltransferase-like 3 (METTL3) promotes OSCC by regulating m⁶A of p38 (14), BMI1 (15), c-Myc (16), PRMT5 and PD-L1 (17). Bioinformatics analysis reveals that HNRNPC and HNRNPA2B1 facilitate progression of OSCC via EMT (18, 19). m⁶A demethylase fat mass and obesity-associated protein (FTO) plays an oncogenic role in arecoline-induced OSCC progression (20) and regulates autophagy and tumorigenesis by targeting eukaryotic translation initiation factor gamma 1 (eIF4G1) in OSCC (21). DEAD-box helicase 3 X-linked (DDX3) is a human RNA helicase that directly regulates m⁶A demethylase ALKBH5, thereby reducing m⁶A methylation of cancer stem cell transcription factor fork head box protein M1 (FOXM1) and Nanog, leading to chemoresistance (22). Analysis based on the cancer genome atlas (TCGA) data indicates that m⁶A RNA methylation regulators can predict the prognosis of patients with HNSCC (23).

Importantly, m⁶A in cancer seems to function as a double-edged sword. The methylation of some genes is related to tumorigenesis, while the demethylation of others will promote tumorigenesis (24, 25). A single m⁶A regulator can exert biological function via different target genes in the same cancer (26, 27). As described above, the same m⁶A regulator may act different functions in different tumors (6–10). These findings demonstrated that regulatory networks of m⁶A methylation in cancers are extremely complex and need to be further explored.

Although m⁶A modification together with corresponding regulators was involved in the occurrence and progression of a variety of cancers, the characteristics of this epigenetic modification in the precancerous stage are still unclear. A few of studies have attempted to explore the m⁶A modification in premalignant stage. When studying the role of m⁶A regulatory genes in colorectal carcinogenesis, it was found that the expression of YTHDF1, IGF2BP1, IGF2BP3, and EIF3B in adenoma was up-regulated at the protein level (28). Oral submucosal fibrosis (OSF) is a precancerous condition. Compared with normal tissues, the overall m⁶A level in OSF tissues was increased, suggesting that m⁶A modification contributes to OSF (29). However, to our knowledge, the change of m⁶A profile and the expression of m⁶A regulators in precancerous cells have not been well studied. In this study, we analyzed the transcriptome-wide m⁶A methylome of oral precancerous DOK cells (30) by MeRIP-seq and RNA-seq, determined the expression of major m⁶A regulators of these cells by RT-PCR, and compared their m⁶A modification with that of normal oral epithelial cells and oral cancer cells.

Materials and methods

Cells and cell culture

This study involved three cell lines, including immortalized normal human oral epithelial cell (HOEC) (BNCC340217) provided by BeNa Culture Collection (Bnbio, Beijing, China), human dysplastic oral keratinocyte (DOK) and oral squamous cell carcinoma cell SCC-9 provided by Shanghai Guandao Biological Engineering Co., Ltd (Sgdbio, Shanghai, China) (30, 31). DOK cells were established from human dysplastic oral mucosa and were considered to be precancerous cells (30). DOK cells were maintained in RPMI-1640 medium, HOEC and SCC-9 cells were maintained in DMEM medium, supplemented with 10% FBS (GIBCO, Australia), 100 U/mL penicillin G and 100 μg/mL streptomycin, in a humidified atmosphere of 5% CO₂ at 37°C. About 5×10⁷ cells were harvested for MeRIP-seq and RNA-seq.

MeRIP sequencing and RNA sequencing

Transcriptome-wide m⁶A sequencing was performed as described previously (32, 33). Briefly, total RNA was isolated using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and quantified. RNA integrity was evaluated by Bioanalyzer 2100 (Agilent, CA, USA). Poly (A) RNA was purified using Dynabeads Oligo (dT)25-61005 (Thermo Fisher, CA, USA) and fragmented into about 100nt using Magnesium RNA Fragmentation Module (NEB, USA). The RNA fragments were divided into two portions, one was kept as input and the other was used to enrich m⁶A-methylated RNA fragments by immunoprecipitation with m⁶A-specific antibody (Synaptic Systems, Germany). RNA-seq library preparation included double stranded cDNA synthesis, addition of A-tailing, adapter ligation, amplification of ligated products, and library purification. The paired-end sequencing (PE150) of libraries was performed on Illumina Novaseq™ 6000 platform (LC-Bio Technology CO., Ltd., Hangzhou, China), following the manufacturer’s protocol.

Bioinformatics analysis

The fastp tool (https://github.com/OpenGene/fastp) was used to remove the low quality reads and trim adaptors (34). Sequence quality of IP and input samples was verified using FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and RseQC (http://rseqc.sourceforge.net/) (35, 36). Clean reads were then mapped to the reference genome Homo sapiens (Version: v101) using HISAT2 (http://daehwankimlab.github.io/hisat2) (37). m⁶A peak calling and analysis of differentially methylated peaks were performed by R package exomePeak2 (https://bioconductor.org/packages/release/bioc/html/exomePeak2.html) (38), and peaks were annotated by intersection with gene architecture using R package ANNOVAR (http://www.openbioinformatics.org/annovar/) (39). The MEME (http://meme-suite.org) (40) and HOMER (http://homer.ucsd.edu/homer/motif) were used for de novo and known motif finding, followed by localization of the motif with respect to peak summit. The expression level of all transcripts and genes from input libraries was analyzed through calculating FPKM (total exon fragments/mapped reads (millions) × exon length (kb)) using StringTie (https://ccb.jhu.edu/software/stringtie) (41). The R package edgeR (https://bioconductor.org/packages/edgeR) (42) was used to identify differentially expressed transcripts and genes, and the threshold was set to |log₂fold change (FC)|≥1 and p value < 0.05. The differentially expressed genes and differentially methylated coding genes were subjected to Gene Oncology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis (43).

RNA extraction and quantitative real-time RT-PCR

Total RNA was isolated by incubating the cells in a 25-cm² culture flask with 1 mL TRIzol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA) in accordance with the manufacturer’s instructions. After genomic DNA was removed by treatment with gDNA Eraser, RNAs were reverse-transcribed with random hexamer primers using PrimeScript RT Enzyme Mix I (PrimeScriptTM RT reagent Kit with gDNA Eraser, TaKaRa, Shiga, Japan). Real-time PCR were processed in triplicate and finished in 20 μL reaction volumes with SYBR^® Premix Ex Taq II (Tli RNaseH Plus, TaKaRa, Shiga, Japan) and the ABI PRISM^®7900 system (ABI). The threshold cycles and relative fold differences were calculated with 2^-ΔΔCt. The primers used in the study are listed in Table S1 and proved to be effective in previous studies (44, 45).

Statistical analysis

For MeRIP-seq, aligned reads were used for peak calling of the MeRIP regions, and significantly enriched regions (peaks) were determined at a threshold of log₂FC ≥ 1 and p value < 0.05. For MeRIP sequencing and RNA sequencing data, differentially methylated peaks were identified through exomePeak2 (38) and differentially expressed genes were identified by the edgeR in R package (42) according to the criteria |log₂FC| ≥ 1 and p < 0.05. The statistical analysis of RT-PCR data were done using Graphpad Prism 8, and p value were calculated using two tailed unpaired student’s t-test. * represent p < 0.05 and ** represent p < 0.01.

Results

Overall features of m⁶A methylation in different oral cells

In order to explore the characteristics of m⁶A methylation in precancerous cells, we did MeRIP-seq and RNA-seq of three oral cell lines, including normal human oral epithelial cell HOEC, precancerous dysplastic oral keratinocyte DOK, and oral squamous cell carcinoma cell SCC-9. Sequencing data were summarized in Table S2, and the reads containing adaptor, low quality bases and undetermined bases were removed from raw data to generate clean data. More than 90% valid data from IP and input samples can be mapped to exons of genes in the reference genome.

Based on clean data, a total of 43,296, 41,946 and 41,817 m⁶A peaks were identified in HOEC, DOK and SCC-9 cells, respectively. In all three cell lines, m⁶A peaks were highly enriched in 3`UTR and stop codon regions, and their distribution and density across the length of mRNA transcripts were similar (Figure 1A). Some enriched m⁶A peaks have typical conserved motifs (Figures 1B–D).

FIGURE 1

Figure 1 Distribution of m⁶A peaks across the length of mRNA transcripts and the representative m⁶A motifs with typical conserved sequence. (A) The m⁶A peaks were highly enriched in 3`UTR and stop codon regions. (B–D) The m⁶A motifs with typical conserved sequence in HOEC, DOK and SCC-9 cells, respectively.

m⁶A profile changed in oral pre-malignant cells

When analyzing the differential m⁶A peaks between each two cell lines, it was found that compared with HOEC cells, there were 1,180 hyper-methylated and 1,606 hypo-methylated m⁶A peaks in DOK cells and 3,916 hyper-methylated and 1,349 hypo-methylated m⁶A peaks in SCC-9 cells (|log₂FC|≥1.0 and p < 0.05). Compared with DOK cells, 5,120 hyper-methylated and 1,243 hypo-methylated m⁶A peaks were identified in SCC-9 cells. Obviously, the m⁶A modification in DOK cells has changed in general compared with HOEC cells, but these changes seem to be insufficient compared with those in SCC-9 cells. Compared with HOEC cells, the top 20 genes with the most significant m⁶A peak changes in DOK and SCC-9 cells were shown in Table 1 and the m⁶A peak changes of a representative gene were shown in Figure 2. These genes in DOK cells were also inconsistent with those in SCC-9 cells, further indicating that the modification of m⁶A in the early phase of carcinogenesis has its own characteristics.

TABLE 1

Table 1 The top 20 altered m⁶A peaks in DOK and SCC-9 cells compared with HOEC cells.

FIGURE 2

Figure 2 The change of m⁶A modification in DOK and SCC-9 cells compared with HOEC cells. m⁶A peak clusters from three cells are shown along the length of the EIF3C transcript. Compared with HOEC cells, hyper-methylated peaks were observed in exonic region in DOK cells and in 3`UTR in SCC-9 cells. EIF3C, eukaryotic translation initiation factor 3 subunit C.

The GO function and KEGG pathway enrichment of differentially methylated mRNAs were analyzed to explore the biological significance of m⁶A modification. Among the top 20 GO terms enriched with differentially methylated mRNAs between DOK and HOEC cells, 13 terms are consistent with those between SCC-9 and HOEC cells (Figures 3A, B). Those differentially methylated genes were significantly enriched in some important biological processes, such as transferase activity, protein phosphorylation, and protein binding and so on. In the first 20 enriched KEGG pathways of differentially methylated mRNAs between DOK and HOEC cells, 5 pathways were also enriched between SCC-9 and HOEC cells, including phosphatidylinositol signaling system, pancreatic cancer, insulin signaling pathway, colorectal cancer and adherens junction (Figures 3C, D). The above data may suggest that precancerous DOK cells and SCC-9 cancer cells shared some common features in the changes of m⁶A modification.

FIGURE 3

Figure 3 GO function and KEGG pathway enrichment of differentially methylated mRNA. (A, B) Compared with HOEC cells, top 20 significantly enriched GO terms in DOK and SCC-9 cells, respectively. (C, D) Compared with HOEC cells, top 20 significantly enriched KEGG pathways in DOK and SCC-9 cells, respectively.

Conjoint analysis of MeRIP-seq and RNA-seq data

We used RNA-seq (MeRIP-seq input library) data to analyze the difference of gene expression between any two cells. Compared with HOEC cells, DOK and SCC-9 cells had 328 and 3,531 significantly upregulated genes and 563 and 3,550 significantly downregulated genes, respectively (|log₂FC|≥1.0 and p < 0.05; Figure 4), suggesting that SCC-9 cells have more obvious transcriptional changes than DOK cells. The top 20 (Table 2) or top 100 (Figure S1) genes differentially expressed between DOK and HOEC cells or between SCC-9 and HOEC cells also showed great differences, although many genes are related to the development and progression of tumors.

FIGURE 4

Figure 4 Volcano plots of differentially expressed genes in DOK and SCC-9 cells compared with HOEC cells. (A) HOEC VS DOK; (B) HOEC VS SCC-9. |log₂FC|≥1.0 and p < 0.05.

TABLE 2

Table 2 The top 20 differentially expressed genes in DOK and SCC-9 cells compared with HOEC cells.

By jointly analyzing the MeRIP-seq and RNA-seq data, we further identified 354 differentially expressed mRNAs with differential m⁶A peaks in DOK cells, including 161 upregulated genes and 193 downregulated genes, compared with HOEC cells. These genes can be divided into four groups, including upregulated mRNAs with hypermethylated (hyper-up) or hypomethylated (hypo-up) m⁶A peaks and downregulated mRNAs with hypermethylated (hyper-down) or hypomethylated (hypo-down) m⁶A peaks (Figure 5A). This indicates that the expression changes of some genes are related to the changes of m⁶A modification in DOK cells. Obviously, there are more such genes in SCC-9 cells (Figure 5B), although the correlation between mRNA m⁶A methylation and its expression level needs to be further evaluated. The GO function and KEGG pathway enrichment of these genes was also analyzed (Figure 6).

FIGURE 5

Figure 5 Distribution of differentially expressed genes with differential m⁶A peaks in DOK and SCC-9 cells, compared with HOEC cells. (A) HOEC VS DOK; (B) HOEC VS SCC-9. Hyper-up, m⁶A peak upregulated and mRNA expression upregulated; Hyper-down, m⁶A peak upregulated and mRNA expression downregulated; Hypo-up, m⁶A peak downregulated and mRNA expression upregulated; Hypo-down, m⁶A peak downregulated and mRNA expression downregulated.

FIGURE 6

Figure 6 GO function and KEGG pathway enrichment of differentially expressed genes with differential m⁶A peaks. (A, B) Top 20 significantly enriched GO terms in DOK VS HOEC and SCC-9 VS HOEC cells, respectively. (C, D) Top 20 significant KEGG pathways in DOK VS HOEC and SCC-9 VS HOEC cells, respectively.

Compared with HOEC cells, top 10 differentially expressed genes with differential m⁶A peaks in DOK cells or SCC-9 cells were summarized in Table 3. Although many of these genes are associated with tumors, they are not consistent in DOK and SCC-9 cells. Compared with DOK cells, the expression changes of differential genes in SCC-9 cells were also more significant. We noted that a transcript can have multiple differential m⁶A peaks, which can be hypermethylated or hypomethylated. The effects of different m⁶A modificaton on gene expression may be consistent or inconsistent. For example, TNFSF15, CPM and ENG all have hyper- and hypo- m⁶A peaks in a single transcript (Table 3), but they are all related to the upregulation of gene expression. ARHGEF26 has two hypo- m⁶A peaks and RGMB has three hypo- m⁶A peaks, which are all related to the downregulation of gene expression.

TABLE 3

Table 3 The first 10 differentially expressed genes with differential m⁶A peaks in DOK cells or SCC-9 cells compared with HOEC cells.

Expression of m⁶A regulatory genes

We wanted to know whether m⁶A methylation changes in DOK and SCC-9 cells were related to the regulatory genes of m⁶A modification. Firstly, we analyzed expression of key m⁶A regulators based on our sequencing data. Compared with HOEC cells, there was no significant change in the expression of m⁶A regulatory genes in DOK cells, but there was significant upregulation of FTO and IGF2BP1 and downregulation of METTL14 and IGF2BP2/3 in SCC-9 cells (|log₂FC|≥1.0 and p < 0.05) (Table S3).

We further determined the expression levels of these m⁶A regulatory genes in different cells by RT-PCR. Among the 20 genes detected, FTO and ALKBH5 were upregulated and METTL3 and VIRMA were downregulated in DOK cells compared with HOEC cells (0.01< p < 0.05, Figure 7). However, METTL16, FTO and ALKBH5 were significantly upregulated and METTL14, RBM15, VIRMA, ZC3H13, IGF2BP2/3 and HNRNPC were significantly downregulated in SCC-9 cells (p < 0.01, Figure 7). The expression of RBM15B, YTHDF2/3, HNRNPA2B1 and METTL3 in SCC-9 cells also decreased or increased to some extent (0.01< p < 0.05). The expression of most m⁶A regulatory genes was consistent with the sequence data (Table S3). Obviously, more m⁶A regulators in SCC-9 cells changed their expression level compared with DOK cells. These also accounted for the more extensive change in m⁶A modification in cancer cells compared with oral precancerous cells.

FIGURE 7

Figure 7 Expression of m⁶A regulatory genes in HOEC, DOK and SCC-9 cells detected by real-time RT-PCR. All bars show mean ± SD of three independent experiments. *p < 0.05, **p < 0.01.

Discussion

m⁶A methylation is the most common RNA post-transcriptional modification, which widely regulates a variety of cellular functions. Its significance in the development and progression of tumor has attracted extensive attention. In this study, through whole transcriptome m⁶A sequencing, we found that oral precancerous DOK cells showed, to some extent, changes of m⁶A modification compared with normal oral epithelial cells. Although the magnitude of these changes in DOK cells was less than that in SCC-9 cells, the two cell lines shared many GO terms and KEGG pathways, which were enriched by the mRNAs with m⁶A modification changes. A total of 354 differentially expressed mRNAs with differential m⁶A peaks were identified between DOK and HOEC cells. Among the 20 m⁶A regulatory genes detected by RT-PCR, FTO, ALKBH5, METTL3 and VIRMA in DOK cells were upregulated or downregulated to a certain extent, but in SCC-9 cells, the expression changes of 10 genes such as METTL14/16, FTO and IGF2BP2/3 were more significant. Our data suggest that precancerous cells do change their m⁶A modification status, but these changes may be more obvious in cancer cells than in precancerous cells.

At present, it is unclear whether there are changes of m⁶A modification in the precancerous stage, or whether m⁶A modification contributes to the initiation of cell transformation. Relevant research is very limited. In colorectal adenomas, a precancerous condition of colorectal cancer, expression of several m⁶A regulatory genes was found to be upregulated (28). However, the status of m⁶A modification and its relationship with m⁶A regulator expression have not been studied. Although the overall m⁶A modification was enhanced in precancerous oral submucosal fibrosis, the target genes and corresponding regulatory factors were unclear (29). In this study, we used a well-established precancerous cell line (30) as a model to investigate the potential contribution of m⁶A modification to the initiation of cell transformation. This is the first time to reveal the characteristics of m⁶A methylation in premalignant cells. It should be mentioned here that only one precancerous cell line was used in this study. This may only reflect the m⁶A modification characteristics of a single precancerous tissue. The genetic background of a single cell line is homogeneous, which is comparable with other cell lines. However, as immortalized cells, immortalization may change some characteristics of the original cells, resulting in false reflection of m⁶A modification state. In addition to analyzing more cell lines, we can also directly analyze m⁶A modification of some oral precancerous tissues, such as lichen planus, leukoplakia and erythroplakia (31, 46). Since it is usually difficult to obtain enough oral precancerous tissue samples, other similar tissues with precancerous characteristics, such as cervical squamous intraepithelial lesions (31) and colorectal adenomas, can also be selected.

We noted that many differentially expressed genes with differential m⁶A peaks in DOK cells (Table 3) are closely related to tumorigenesis. For example, E2F1 transcriptional activation of neural epidermal growth factor-like 2 (NELL2) accelerates the progression of non-small cell lung cancer (47). Calbindin 2 (CALB2) promotes metastasis of hepatocellular carcinoma (48). The upregulation of transmembrane 4 L six family member 1 (TM4SF1) in papillary thyroid carcinoma patients is associated with lymph node metastases (49). Therefore, the role of the m⁶A target genes identified in tumorigenesis is worth studying. As mentioned above, many m⁶A regulators are also associated with tumorigenesis and progression. In this study we also found that the expression of several m⁶A regulators changed in oral precancerous cells. Their exact function in tumors remains to be evaluated.

Compare with DOK cells, oral squamous cell carcinoma cells have more extensive changes in m⁶A modification, involving more m⁶A peaks and more genes. On the one hand, this suggests that the epigenetic modification of m⁶A may play an important role in carcinogenesis. On the other hand, this may indicate that there is still a long way to go from precancerous lesions to cancer. We even found that the mRNA of the same gene can have different m⁶A modifications in precancerous cells and cancer cells (Table 1). For example, eukaryotic translation initiation factor 3 subunit C (EIF3C), a subunit of the protein translation initiation factor EIF3, was hyper-methylated in exonic region in DOK cells and 3`UTR in SCC-9 cells (Figure 2). It was reported that m⁶A reader YTHDF1 bound to m⁶A-modified EIF3C mRNA and promoted the translation of EIF3C and the overall translational output, consequently facilitating tumorigenesis and metastasis of ovarian cancer (50). Compared with cancer cells, the change of m⁶A regulatory gene expression in precancerous cells also seems to be less obvious. It can be supposed that the regulatory genes in precancerous cells are mainly regulated by changing enzyme activity rather than inducing gene expression. This needs further research to evaluate. Therefore, we must explore the function and significance of the key m⁶A target genes and regulators identified in this study in tumorigenesis. This may provide potential intervention targets for the prevention of cancer development through epigenetic modification in the future.

Data availability statement

The data presented in the study are deposited in the GEO repository, accession number GSE213714.

Author contributions

DY, WZ and SC contributed conception and design of the study. XC, LC, YT, YH, KP, LY, WX collected samples and performed experiments. All authors participated in the writing of the manuscript and confirmed the final review of the manuscript.

Funding

This study was funded by the National Natural Science Foundation of China (No. 81873711 and No. 31670788) and by the Open Fund of Guangdong Key Laboratory of Pharmaceutical Functional Genes (No. 2020B1212060031 and No. 2017B030314021).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor ZL declared a shared parent affiliation with the authors at the time of review.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fonc.2022.939449/full#supplementary-material

Supplementary Figure 1 | Hierarchical clustering heatmap of the top 100 differentially expressed genes in DOK and SCC-9 cells compared with HOEC cells. (A) HOEC VS DOK; (B) HOEC VS SCC-9.

References

1. Dominissini D, Moshitch-Moshkovitz S, Schwartz S, Salmon-Divon M, Ungar L, Osenberg S, et al. Topology of the human and mouse m6A RNA methylomes revealed by m6A-seq. Nature (2012) 485:201–6. doi: 10.1038/nature11112

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Meyer KD, Saletore Y, Zumbo P, Elemento O, Mason CE, Jaffrey SR. Comprehensive analysis of mRNA methylation reveals enrichment in 3’ UTRs and near stop codons. Cell (2012) 22:149:1635–46. doi: 10.1016/j.cell.2012.05.003

CrossRef Full Text | Google Scholar

3. Jiang X, Liu B, Nie Z, Duan L, Xiong Q, Jin Z, et al. The role of m6A modification in the biological functions and diseases. Signal Transduct Target Ther (2021) 6:74. doi: 10.1038/s41392-020-00450-x

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Jing FY, Zhou LM, Ning YJ, Wang XJ, Zhu YM. The biological function, mechanism, and clinical significance of m6A RNA modifications in head and neck carcinoma: A systematic review. Front Cell Dev Biol (2021) 9:683254. doi: 10.3389/fcell.2021.683254

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Wang T, Kong S, Tao M, Ju S. The potential role of RNA N6-methyladenosine in cancer progression. Mol Cancer (2020) 19:88. doi: 10.1186/s12943-020-01204-7

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Gao M, Qi Z, Feng W, Huang H, Xu Z, Dong Z, et al. m6A demethylation of cytidine deaminase APOBEC3B mRNA orchestrates arsenic-induced mutagenesis. J Biol Chem (2022) 298:101563. doi: 10.1016/j.jbc.2022.101563

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Zhou Y, Wang Q, Deng H, Xu B, Zhou Y, Liu J, et al. N6-methyladenosine demethylase FTO promotes growth and metastasis of gastric cancer via m(6)A modification of caveolin-1 and metabolic regulation of mitochondrial dynamics. Cell Death Dis (2022) 13:72. doi: 10.1038/s41419-022-04503-7

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Lv D, Ding S, Zhong L, Tu J, Li H, Yao H, et al. M(6)A demethylase FTO-mediated downregulation of DACT1 mRNA stability promotes wnt signaling to facilitate osteosarcoma progression. Oncogene (2022) 41:1727–41. doi: 10.1038/s41388-022-02214-z

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Jeschke J, Collignon E, Al Wardi C, Krayem M, Bizet M, Jia Y, et al. Downregulation of the FTO m(6)A RNA demethylase promotes EMT-mediated progression of epithelial tumors and sensitivity to wnt inhibitors. Nat Cancer (2021) 2:611–28. doi: 10.1038/s43018-021-00223-7

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Li J, Zhu L, Shi Y, Liu J, Lin L, Chen X. m6A demethylase FTO promotes hepatocellular carcinoma tumorigenesis via mediating PKM2 demethylation. Am J Transl Res (2019) 11:6084–92.

PubMed Abstract | Google Scholar

11. Taketo K, Konno M, Asai A, Koseki J, Toratani M, Satoh T, et al. The epitranscriptome m6A writer METTL3 promotes chemo- and radioresistance in pancreatic cancer cells. Int J Oncol (2018) 52:621–29. doi: 10.3892/ijo.2017.4219

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Shen C, Xuan B, Yan T, Ma Y, Xu P, Tian X, et al. m(6)A-dependent glycolysis enhances colorectal cancer progression. Mol Cancer (2020) 19:72. doi: 10.1186/s12943-020-01190-w

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Wang W, Shao F, Yang X, Wang J, Zhu R, Yang Y, et al. METTL3 promotes tumour development by decreasing APC expression mediated by APC mRNA N(6)-methyladenosine-dependent YTHDF binding. Nat Commun (2021) 12:3803. doi: 10.1038/s41467-021-23501-5

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Xu T, Zhang W, Chai L, Liu C, Zhang S, Xu T. Methyltransferase-like 3-induced N6-methyladenosine upregulation promotes oral squamous cell carcinoma by through p38. Oral Dis (2021). doi: 10.1111/odi.14016

CrossRef Full Text | Google Scholar

15. Liu L, Wu Y, Li Q, Liang J, He Q, Zhao L, et al. METTL3 promotes tumorigenesis and metastasis through BMI1 m(6)A methylation in oral squamous cell carcinoma. Mol Ther (2020) 28:2177–90. doi: 10.1016/j.ymthe.2020.06.024

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Zhao W, Cui Y, Liu L, Ma X, Qi X, Wang Y, et al. METTL3 facilitates oral squamous cell carcinoma tumorigenesis by enhancing c-myc stability via YTHDF1-mediated m(6)A modification. Mol Ther Nucleic Acids (2020) 20:1–12. doi: 10.1016/j.omtn.2020.01.033

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Ai Y, Liu S, Luo H, Wu S, Wei H, Tang Z, et al. METTL3 intensifies the progress of oral squamous cell carcinoma via modulating the m6A amount of PRMT5 and PD-L1. J Immunol Res (2021) 2021:6149558. doi: 10.1155/2021/6149558

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Huang GZ, Wu QQ, Zheng ZN, Shao TR, Chen YC, Zeng WS, et al. M6A-related bioinformatics analysis reveals that HNRNPC facilitates progression of OSCC via EMT. Aging (Albany NY) (2020) 12:11667–84. doi: 10.18632/aging.103333

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Zhu F, Yang T, Yao M, Shen T, Fang C. HNRNPA2B1, as a m(6)A reader, promotes tumorigenesis and metastasis of oral squamous cell carcinoma. Front Oncol (2021) 11:716921. doi: 10.3389/fonc.2021.716921

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Li X, Xie X, Gu Y, Zhang J, Song J, Cheng X, et al. Fat mass and obesity-associated protein regulates tumorigenesis of arecoline-promoted human oral carcinoma. Cancer Med (2021) 10:6402–15. doi: 10.1002/cam4.4188

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Wang F, Liao Y, Zhang M, Zhu Y, Wang W, Cai H, et al. N6-methyladenosine demethyltransferase FTO-mediated autophagy in malignant development of oral squamous cell carcinoma. Oncogene (2021) 40:3885–98. doi: 10.1038/s41388-021-01820-7

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Shriwas O, Priyadarshini M, Samal SK, Rath R, Panda S, Das Majumdar SK, et al. DDX3 modulates cisplatin resistance in OSCC through ALKBH5-mediated m(6)A-demethylation of FOXM1 and NANOG. Apoptosis (2020) 25:233–46. doi: 10.1007/s10495-020-01591-8

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Zhao X, Cui L. Development and validation of a m(6)A RNA methylation regulators-based signature for predicting the prognosis of head and neck squamous cell carcinoma. Am J Cancer Res (2019) 9:2156–69.

PubMed Abstract | Google Scholar

24. Li T, Hu PS, Zuo Z, Lin JF, Li X, Wu QN, et al. METTL3 facilitates tumor progression via an m(6)A-IGF2BP2-dependent mechanism in colorectal carcinoma. Mol Cancer (2019) 18:112. doi: 10.1186/s12943-019-1038-7

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Niu Y, Lin Z, Wan A, Chen H, Liang H, Sun L, et al. RNA N6-methyladenosine demethylase FTO promotes breast tumor progression through inhibiting BNIP3. Mol Cancer (2019) 18:46. doi: 10.1186/s12943-019-1004-4

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Yue B, Song C, Yang L, Cui R, Cheng X, Zhang Z, et al. METTL3-mediated N6-methyladenosine modification is critical for epithelial-mesenchymal transition and metastasis of gastric cancer. Mol Cancer (2019) 18:142. doi: 10.1186/s12943-019-1065-4

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Wang Q, Chen C, Ding Q, Zhao Y, Wang Z, Chen J, et al. METTL3-mediated m(6)A modification of HDGF mRNA promotes gastric cancer progression and has prognostic significance. Gut (2020) 69:1193–205. doi: 10.1136/gutjnl-2019-319639

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Kuai D, Zhu S, Shi H, Yang R, Liu T, Liu H, et al. Aberrant expression of m(6)A mRNA methylation regulators in colorectal adenoma and adenocarcinoma. Life Sci (2021) 273:119258. doi: 10.1016/j.lfs.2021.119258

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Li X, Gao Y, Chen W, Gu Y, Song J, Zhang J, et al. N6-methyladenosine modification contributes to arecoline-mediated oral submucosal fibrosis. J Oral Pathol Med (2022) 51:474–82. doi: 10.1111/jop.13292

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Chang SE, Foster S, Betts D, Marnock WE. DOK, a cell line established from human dysplastic oral mucosa, shows a partially transformed non-malignant phenotype. Int J Cancer (1992) 52:896–902. doi: 10.1002/ijc.2910520612

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Chen X, Yi C, Yang MJ, Sun X, Liu X, Ma H, et al. Metabolomics study reveals the potential evidence of metabolic reprogramming towards the warburg effect in precancerous lesions. J Cancer (2021) 12:1563–74. doi: 10.7150/jca.54252

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Zhang Z, Wang Q, Zhang M, Zhang W, Zhao L, Yang C, et al. Comprehensive analysis of the transcriptome-wide m6A methylome in colorectal cancer by MeRIP sequencing. Epigenetics (2021) 16:425–35. doi: 10.1080/15592294.2020.1805684

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Han Z, Yang B, Wang Q, Hu Y, Wu Y, Tian Z. Comprehensive analysis of the transcriptome-wide m(6)A methylome in invasive malignant pleomorphic adenoma. Cancer Cell Int (2021) 21:142. doi: 10.1186/s12935-021-01839-6

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Chen S, Zhou Y, Chen Y, Gu J. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics (2018) 34:i884–90. doi: 10.1093/bioinformatics/bty560

PubMed Abstract | CrossRef Full Text | Google Scholar

35. de Sena Brandine G, Smith AD. Falco: high-speed FastQC emulation for quality control of sequencing data. F1000Res (2019) 8:1874. doi: 10.12688/f1000research.21142.2

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Wang L, Wang S, Li W. RSeQC: quality control of RNA-seq experiments. Bioinformatics (2012) 28:2184–5. doi: 10.1093/bioinformatics/bts356

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods (2015) 12:357–60. doi: 10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Meng J, Lu Z, Liu H, Zhang L, Zhang S, Chen Y, et al. A protocol for RNA methylation differential analysis with MeRIP-seq data and exomePeak R/Bioconductor package. Methods (2014) 69:274–81. doi: 10.1016/j.ymeth.2014.06.008

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res (2010) 38:e164. doi: 10.1093/nar/gkq603

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Bailey TL, Boden M, Buske FA, Frith M, Grant CE, Clementi L, et al. MEME SUITE: tools for motif discovery and searching. Nucleic Acids Res (2009) 37(Web Server issue):W202–8. doi: 10.1093/nar/gkp335

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Pertea M, Pertea GM, Antonescu CM, Chang TC, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol (2015) 33(3):290–5. doi: 10.1038/nbt.3122

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics (2010) 26:139–40. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc (2009) 4:44–57. doi: 10.1038/nprot.2008.211

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Shen M, Li Y, Wang Y, Shao J, Zhang F, Yin G, et al. N(6)-methyladenosine modification regulates ferroptosis through autophagy signaling pathway in hepatic stellate cells. Redox Biol (2021) 47:102151. doi: 10.1016/j.redox.2021.102151

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Hou J, Wang Z, Li H, Zhang H, Luo L. Gene signature and identification of clinical trait-related m(6) a regulators in pancreatic cancer. Front Genet (2020) 11:522. doi: 10.3389/fgene.2020.00522

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Chen X, Zhao Y. Human papillomavirus infection in oral potentially malignant disorders and cancer. Arch Oral Biol (2017) 83:334–9. doi: 10.1016/j.archoralbio.2017.08.011

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Wang Y, Li M, Zhang L, Chen Y, Zhang S. m6A demethylase FTO induces NELL2 expression by inhibiting E2F1 m6A modification leading to metastasis of non-small cell lung cancer. Mol Ther Oncolytics (2021) 21:367–76. doi: 10.1016/j.omto.2021.04.011

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Chu H, Zhao Q, Shan Y, Zhang S, Sui Z, Li X, et al. All-ion monitoring-directed low-abundance protein quantification reveals CALB2 as a key promoter in hepatocellular carcinoma metastasis. Anal Chem (2022) 94:6102–11. doi: 10.1021/acs.analchem.1c03562

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Wang K, Zheng J, Wu Z, Fang JG, Zhao JY, Yao JM, et al. Hyper-expression and hypomethylation of TM4SF1 are associated with lymph node metastases in papillary thyroid carcinoma patients. Neoplasma (2022) 69:341–51. doi: 10.4149/neo_2021_210819N1191

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Liu T, Wei Q, Jin J, Luo Q, Liu Y, Yang Y, et al. The m6A reader YTHDF1 promotes ovarian cancer progression via augmenting EIF3C translation. Nucleic Acids Res (2020) 48:3816– 31. doi: 10.1093/nar/gkaa048

PubMed Abstract | CrossRef Full Text | Google Scholar