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, 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 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 , 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 ( ). Some enriched m6A peaks have typical conserved motifs ( ).

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 and the m6A peak changes of a representative gene were shown in . 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.

Table 1

The top 20 altered m6A peaks in DOK and SCC-9 cells compared with HOEC cells.

DOK				SCC-9
Hyper-methylated		Hypo-methylated		Hyper-methylated		Hypo-methylated
Genes	Peak region	Genes	Peak region	Genes	Peak region	Genes	Peak region
EIF3C	exonic	SULT1A4	UTR3	SLC9A3	exonic	SLX1B	UTR3
NOMO1	exonic	NOMO2	UTR3	SRCIN1	UTR5	AC148477	ncRNA_exonic
TMLHE	UTR3	FAM72D	UTR3	SERF1B	UTR3	PHTF1	intronic
TBC1D3D	exonic	F8A1	UTR5	AL513218	ncRNA_exonic	HSPA1B	exonic
ETS1	UTR3	HES2	UTR3	SCAMP1	UTR5	SMG1P5	ncRNA_exonic
VGF	exonic	ZNF10	UTR3	PWP2	exonic	SP140L	intronic
ROR1	UTR3	SLC30A4	UTR3	NCK1-DT	ncRNA_exonic	ZNF467	UTR5
MAGEA9B	UTR5	LINC02341	ncRNA_exonic	FAM169B	ncRNA_exonic	BRSK1	UTR3
CD177	UTR3	C2orf42	UTR5	EIF3C	UTR3	STEAP2	UTR3
TTC28-AS1	ncRNA_exonic	CYTH3	exonic	MIR137HG	ncRNA_exonic	ITPR1	UTR3

UTR3, 3’ untranslated region; UTR5, 5’ untranslated region.

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 ( ). 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 ( ). The above data may suggest that precancerous DOK cells and SCC-9 cancer cells shared some common features in the changes of m6A modification.

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; ), suggesting that SCC-9 cells have more obvious transcriptional changes than DOK cells. The top 20 ( ) or top 100 ( ) 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

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.

Table 2

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

DOK				SCC-9
Upregulated	log2FC	Downregulated	log2FC	Upregulated	log2FC	Downregulated	log2FC
AC004922	7.30	AC006064	-4.60	HDGFL3	11.19	CBLC	-13.12
IGFL2-AS1	5.18	RPPH1	-3.33	INA	11.05	FAR2	-12.31
CXCL8	3.87	SCARNA5	-3.28	SPINK13	10.43	CLDN3	-11.60
NMRAL2P	3.73	SCARNA10	-3.07	STC1	10.38	MAGED1	-11.40
SLC7A11	3.50	UPK3BL1	-2.78	COL4A5	10.33	HNF1A	-11.10
DHRS9	3.38	PRSS2	-2.52	LAMA4	10.27	RUBCNL	-11.04
AL137800	3.37	RNU1-3	-2.30	NNMT	10.27	FOXA3	-10.80
ZBED2	3.36	VWA5B2	-2.29	CPM	10.21	SMIM22	-10.78
TNFSF9	3.21	RN7SL3	-2.19	ANXA8	10.21	GALNT5	-10.73
DDIT4	3.03	MYL9	-2.04	XIST	10.21	PHGR1	-10.71

FPKM>2 in up-regulated genes.

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 ( ). 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 ( ), 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 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.

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 . 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 ( ), 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.

Table 3

The first 10 differentially expressed genes with differential m6A peaks in DOK cells or SCC-9 cells compared with HOEC cells.

DOK
Hypo-up	log2FC	Hyper-up	log2FC	Hypo-down	log2FC	Hyper-down	log2FC
CXCL8	3.87	TNFSF15	2.57	EPHB3	-1.94	WFDC21P	-1.73
SLC7A11	3.50	UCA1	2.22	EDN1	-1.67	SNN	-1.59
ZBED2	3.36	NELL2	2.05	SNN	-1.59	GAL3ST2	-1.54
DDIT4	3.03	CALB2	2.04	ARHGEF26	-1.55	PBXIP1	-1.38
ANKRD37	2.87	ANGPTL4	1.95	RGMB	-1.36	PARS2	-1.35
TNFSF15	2.57	CD177	1.77	AC159540	-1.24	RHOV	-1.31
CXCL1	2.51	TM4SF1	1.70	CDHR1	-1.21	CBX6	-1.25
HMGCS2	2.36	SLC38A2	1.62	DIPK1A	-1.20	ZNF696	-1.25
TNFAIP3	2.27	BIRC3	1.57	IFITM10	-1.20	B3GNT9	-1.22
HMGCS1	2.22	IDI1	1.51	MECOM	-1.18	MYPOP	-1.22
SCC-9
Hypo-up	log2FC	Hyper-up	log2FC	Hypo-down	log2FC	Hyper-down	log2FC
CPM	10.21	CPM	10.21	UTS2	-10.14	FAR2	-12.31
XIST	10.21	PDZD2	9.29	VIL1	-9.64	MAGED1	-11.40
SALL4	9.95	SMC1B	9.27	BICDL2	-9.46	RUBCNL	-11.04
FN1	9.71	WNT5A	9.08	KRT20	-9.25	SMIM22	-10.78
ARHGEF10	9.69	PSCA	9.04	PLEKHG6	-8.97	FUT2	-9.92
ENG	8.53	FSTL1	8.99	CYP3A5	-8.92	VIL1	-9.64
LOX	8.37	SHISAL1	8.92	MUC13	-8.81	MACC1	-9.38
GNB4	8.33	ENG	8.53	P2RY1	-8.39	TSTD1	-8.93
FERMT2	8.22	ACKR3	8.46	B3GNT3	-8.01	CYP3A5	-8.92
IGFBP7	8.17	LINC02593	8.44	GPX2	-7.88	PTPRH	-8.39

FPKM>2 in up-regulated genes.

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) ( ).

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, ). 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, ). 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 ( ). 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.

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 ( ) 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 ( ). 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 ( ). 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.