- 1Hospital of Stomatology, Guanghua School of Stomatology, Guangdong Provincial Key Laboratory of Stomatology, Sun Yat-sen University, Guangzhou, China
- 2Guangdong Key Laboratory of Pharmaceutical Functional Genes, Department of Biochemistry, School of Life Sciences, Sun Yat-sen University, Guangzhou, China
- 3State Key Laboratory for Biocontrol, Department of Biochemistry, School of Life Sciences, Sun Yat-sen University, Guangzhou, China
As the most common post-transcriptional RNA modification, m6A methylation extensively regulates the structure and function of RNA. The dynamic and reversible modification of m6A is coordinated by m6A writers and erasers. m6A reader proteins recognize m6A modification on RNA, mediating different downstream biological functions. mRNA m6A 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 m6A modification and major m6A 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 m6A peaks and 354 differentially expressed mRNAs with differential m6A peaks in DOK cells. Although the change of m6A modification in DOK cells was less than that in SCC-9 cells, mRNAs with differential m6A in both cell lines were enriched into many identical GO terms and KEGG pathways. Among the 20 known m6A 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 m6A modification. Identifying some key m6A 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
N6-methyladenine (m6A) 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 m6A modifications in each mRNA (1, 2). Most m6A 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).
m6A is decorated by m6A methyltransferase complex, and its components include METTL3/14/16, WTAP, RBM15, RBM15B, VIRMA, and ZC3H13. They are collectively referred to as m6A writers. On the contrary, demethylases such as FTO, ALKBH5 and ALKBH3 can remove m6A and act as erasers. The cooperation of m6A writers and erasers regulates the dynamic and reversible modification of m6A. m6A modification of RNAs must be recognized by m6A reader proteins to mediate different downstream biological functions. Several categories of proteins function as m6A 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 m6A modification plays an important role in regulating RNA splicing, translation, stability, translocation, and advanced structure (3, 4).
There is increasing evidence that RNA m6A methylation and its regulators are associated with a variety of human diseases, especially cancer (3–5). m6A regulators may act as writers to catalyze or erasers to remove m6A 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 m6A modification involves many aspects of tumors. For example, demethylase FTO-mediated m6A demethylation of cytidine deaminase APOBEC3B mRNA promotes arsenic-induced mutagenesis (6). FTO promotes growth and metastasis of gastric cancer through m6A 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). m6A writer METTL3 promotes chemo- and radioresistance in pancreatic cancer cells (11). METTL3 stabilizes HK2 and SLC2A1 (GLUT1) expression in colorectal cancer and m6A-dependent glycolysis enhances colorectal cancer progression (12). METTL3 promotes tumor development by decreasing APC gene expression mediated by APC mRNA m6A-dependent YTHDF binding (13).
The significance of m6A modification in head and neck squamous cell carcinoma (HNSCC) has been well reviewed (4). m6A 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 m6A 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). m6A 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 m6A demethylase ALKBH5, thereby reducing m6A 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 m6A RNA methylation regulators can predict the prognosis of patients with HNSCC (23).
Importantly, m6A 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 m6A regulator can exert biological function via different target genes in the same cancer (26, 27). As described above, the same m6A regulator may act different functions in different tumors (6–10). These findings demonstrated that regulatory networks of m6A methylation in cancers are extremely complex and need to be further explored.
Although m6A 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 m6A modification in premalignant stage. When studying the role of m6A 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 m6A level in OSF tissues was increased, suggesting that m6A modification contributes to OSF (29). However, to our knowledge, the change of m6A profile and the expression of m6A regulators in precancerous cells have not been well studied. In this study, we analyzed the transcriptome-wide m6A methylome of oral precancerous DOK cells (30) by MeRIP-seq and RNA-seq, determined the expression of major m6A regulators of these cells by RT-PCR, and compared their m6A 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% CO2 at 37°C. About 5×107 cells were harvested for MeRIP-seq and RNA-seq.
MeRIP sequencing and RNA sequencing
Transcriptome-wide m6A 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 m6A-methylated RNA fragments by immunoprecipitation with m6A-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). m6A 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 |log2fold 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-cm2 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 log2FC ≥ 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 |log2FC| ≥ 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 m6A methylation in different oral cells
In order to explore the characteristics of m6A 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 m6A peaks were identified in HOEC, DOK and SCC-9 cells, respectively. In all three cell lines, m6A 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 m6A peaks have typical conserved motifs (Figures 1B–D).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/939449/fonc-12-939449-HTML/image_m/fonc-12-939449-g001.jpg)
Figure 1 Distribution of m6A peaks across the length of mRNA transcripts and the representative m6A motifs with typical conserved sequence. (A) The m6A peaks were highly enriched in 3`UTR and stop codon regions. (B–D) The m6A motifs with typical conserved sequence in HOEC, DOK and SCC-9 cells, respectively.
m6A profile changed in oral pre-malignant cells
When analyzing the differential m6A 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 m6A peaks in DOK cells and 3,916 hyper-methylated and 1,349 hypo-methylated m6A peaks in SCC-9 cells (|log2FC|≥1.0 and p < 0.05). Compared with DOK cells, 5,120 hyper-methylated and 1,243 hypo-methylated m6A peaks were identified in SCC-9 cells. Obviously, the m6A 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 m6A peak changes in DOK and SCC-9 cells were shown in Table 1 and the m6A 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 m6A in the early phase of carcinogenesis has its own characteristics.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/939449/fonc-12-939449-HTML/image_m/fonc-12-939449-g002.jpg)
Figure 2 The change of m6A modification in DOK and SCC-9 cells compared with HOEC cells. m6A 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 m6A 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 m6A modification.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/939449/fonc-12-939449-HTML/image_m/fonc-12-939449-g003.jpg)
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 (|log2FC|≥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.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/939449/fonc-12-939449-HTML/image_m/fonc-12-939449-g004.jpg)
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. |log2FC|≥1.0 and p < 0.05.
By jointly analyzing the MeRIP-seq and RNA-seq data, we further identified 354 differentially expressed mRNAs with differential m6A 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) m6A peaks and downregulated mRNAs with hypermethylated (hyper-down) or hypomethylated (hypo-down) m6A peaks (Figure 5A). This indicates that the expression changes of some genes are related to the changes of m6A modification in DOK cells. Obviously, there are more such genes in SCC-9 cells (Figure 5B), although the correlation between mRNA m6A 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).
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/939449/fonc-12-939449-HTML/image_m/fonc-12-939449-g005.jpg)
Figure 5 Distribution of differentially expressed genes with differential m6A peaks in DOK and SCC-9 cells, compared with HOEC cells. (A) HOEC VS DOK; (B) HOEC VS SCC-9. Hyper-up, m6A peak upregulated and mRNA expression upregulated; Hyper-down, m6A peak upregulated and mRNA expression downregulated; Hypo-up, m6A peak downregulated and mRNA expression upregulated; Hypo-down, m6A peak downregulated and mRNA expression downregulated.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/939449/fonc-12-939449-HTML/image_m/fonc-12-939449-g006.jpg)
Figure 6 GO function and KEGG pathway enrichment of differentially expressed genes with differential m6A 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 m6A 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 m6A peaks, which can be hypermethylated or hypomethylated. The effects of different m6A modificaton on gene expression may be consistent or inconsistent. For example, TNFSF15, CPM and ENG all have hyper- and hypo- m6A peaks in a single transcript (Table 3), but they are all related to the upregulation of gene expression. ARHGEF26 has two hypo- m6A peaks and RGMB has three hypo- m6A peaks, which are all related to the downregulation of gene expression.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/939449/fonc-12-939449-HTML/image_m/fonc-12-939449-t003.jpg)
Table 3 The first 10 differentially expressed genes with differential m6A peaks in DOK cells or SCC-9 cells compared with HOEC cells.
Expression of m6A regulatory genes
We wanted to know whether m6A methylation changes in DOK and SCC-9 cells were related to the regulatory genes of m6A modification. Firstly, we analyzed expression of key m6A regulators based on our sequencing data. Compared with HOEC cells, there was no significant change in the expression of m6A 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 (|log2FC|≥1.0 and p < 0.05) (Table S3).
We further determined the expression levels of these m6A 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 m6A regulatory genes was consistent with the sequence data (Table S3). Obviously, more m6A regulators in SCC-9 cells changed their expression level compared with DOK cells. These also accounted for the more extensive change in m6A modification in cancer cells compared with oral precancerous cells.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/939449/fonc-12-939449-HTML/image_m/fonc-12-939449-g007.jpg)
Figure 7 Expression of m6A 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
m6A 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 m6A sequencing, we found that oral precancerous DOK cells showed, to some extent, changes of m6A 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 m6A modification changes. A total of 354 differentially expressed mRNAs with differential m6A peaks were identified between DOK and HOEC cells. Among the 20 m6A 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 m6A 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 m6A modification in the precancerous stage, or whether m6A 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 m6A regulatory genes was found to be upregulated (28). However, the status of m6A modification and its relationship with m6A regulator expression have not been studied. Although the overall m6A 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 m6A modification to the initiation of cell transformation. This is the first time to reveal the characteristics of m6A 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 m6A 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 m6A modification state. In addition to analyzing more cell lines, we can also directly analyze m6A 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 m6A 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 m6A target genes identified in tumorigenesis is worth studying. As mentioned above, many m6A regulators are also associated with tumorigenesis and progression. In this study we also found that the expression of several m6A 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 m6A modification, involving more m6A peaks and more genes. On the one hand, this suggests that the epigenetic modification of m6A 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 m6A 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 m6A reader YTHDF1 bound to m6A-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 m6A 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 m6A 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
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
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
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
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
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
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
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
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
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.
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
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
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
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
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
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
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
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
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
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
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
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
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.
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
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
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
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
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
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
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
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
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
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
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
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
36. Wang L, Wang S, Li W. RSeQC: quality control of RNA-seq experiments. Bioinformatics (2012) 28:2184–5. doi: 10.1093/bioinformatics/bts356
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
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
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
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
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
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
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
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
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
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
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
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
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
Keywords: m6A modification, m6A regulatory genes, MeRIP sequencing, precancerous cells, dysplastic oral keratinocyte (DOK), human oral epithelial cell (HOEC), oral squamous cell carcinoma cell
Citation: Chen X, Chen L, Tang Y, He Y, Pan K, Yuan L, Xie W, Chen S, Zhao W and Yu D (2022) Transcriptome-wide m6A methylome analysis uncovered the changes of m6A modification in oral pre-malignant cells compared with normal oral epithelial cells. Front. Oncol. 12:939449. doi: 10.3389/fonc.2022.939449
Received: 09 May 2022; Accepted: 29 August 2022;
Published: 28 September 2022.
Edited by:
Zexian Liu, Sun Yat-sen University Cancer Center (SYSUCC), ChinaReviewed by:
Mojgan Alaeddini, Tehran University of Medical Sciences, IranZhicheng Pan, University of California, Los Angeles, United States
Copyright © 2022 Chen, Chen, Tang, He, Pan, Yuan, Xie, Chen, Zhao and Yu. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Wei Zhao, zhaowei3@mail.sysu.edu.cn; Dongsheng Yu, yudsh@mail.sysu.edu.cn