Transcriptome Mapping of the Internal N7-Methylguanosine Methylome in Messenger RNAs in Human Oral Squamous Cell Carcinoma

Background : Internal N7-methylguanosine (m7G) methylation in mammalian messenger RNAs (mRNAs) is essential in disease development. However, the status of internally m7G-modified mRNAs in oral squamous cell carcinoma (OSCC) remains poorly understood. Methods : Methylated RNA immunoprecipitation sequencing (MeRIP-seq) was used to identify the m7G modification level of mRNAs and the expression of mRNAs between OSCC and normal tissues. These differentially methylated and expressed genes were subjected to Gene Ontology (GO) enrichment, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis and construction of protein-protein interaction (PPI) networks. Quantitative real-time PCR (qPCR) assay was performed to detect the expression of Baculoviral IAP Repeat Containing 3 ( BIRC3 ) in vitro . The biological function of BIRC3 in OSCC was clarified using CCK-8, Transwell migration and Western blot assays. Results : The m7G-mRNA profile showed 9514 unique m7G peaks within 7455 genes in OSCC tissues. In addition, the most conserved m7G motif within mRNAs in OSCC was GGARG (R = G/A). The identified m7G peaks were mainly distributed in the coding sequence region within mRNAs in OSCC. GO enrichment and KEGG pathway analyses showed that m7G-modified genes were closely related to cancer progression. m7G-modified hub genes were screened from the constructed PPI networks. Furthermore, BIRC3 with high m7G methylation showed high expression in OSCC cell lines, as confirmed by qPCR assay. Functionally, the knockdown of BIRC3 significantly inhibited the proliferation and migration ability of CAL-27 cells in vitro functional assays. In addition, the relative expression of E-cadherin expression was elevated, while Vimentin and N-cadherin protein expression was decreased in CAL-27 cells transfected with si-BIRC3. This study suggests that BIRC3 could promote OSCC proliferation and migration, which may be associated with involvement in epithelial-mesenchymal transition (EMT) progression. Conclusions : This paper constructed a transcriptome map of internal m7G in mRNAs, which provides potential research value to study the role of m7G methylation in OSCC.


Introduction
Oral squamous cell carcinoma (OSCC) is one of the most common malignant tumors worldwide [1].Recent studies have confirmed that the incidence and mortality rate of OSCC are increasing annually due to the high local recurrence rate and susceptibility to metastasis [2,3].Although technological advances in surgery, radiotherapy, and chemotherapy have improved the local control of OSCC, the 5-year survival rate of OSCC patients has remained at less than 60% [2,4].Therefore, exploring the pathogenesis of OSCC and identifying effective therapeutic targets are crucial to improving OSCC patients' prognoses.
RNA methylation is an essential epigenetic modification in post-transcriptional regulation [5].
N7-methylguanosine (m7G), one of the most common types of RNA methylation, adds a methyl group to the 7th N of RNA guanine in the presence of methyltransferase [6].
Many previous studies have confirmed that m7G methylation is found in ribosomal RNA (rRNA) [7], transfer RNA (tRNA) [8], microRNA (miRNA) [9] and long noncoding RNA (lncRNA) [10] structures.Recently, Zhang et al. [11] discovered internal m7G sites within mammalian mRNAs using advanced m7G-methylated RNA immunoprecipitation sequencing (m7G-MeRIP-seq) technology.Furthermore, m7G methylation within mRNAs can affect disease development [12,13].Specifically, m7Gmodified mRNAs can promote angiogenesis in peripheral arterial disease [12].Internal m7G methylation within mR-NAs is identified in acute myeloid leukemia (AML) using m7G-MeRIP-seq, and hyper m7G-modified genes are associated with multidrug resistance in AML [13].Moreover, it has been confirmed that m7G-modified tRNAs can participate in tumor progression in head and neck squamous cell carcinoma (HNSCC) [14,15].For example, m7G-modified tRNAs can regulate critical signaling pathways (PI3K/AKT/mTOR and WNT/β-catenin), affecting the progression of HNSCC [14,15].Many studies have confirmed that N6-methyladenosine (m6A) modifications are common in mRNAs, and can influence cancer progression through regulation of methylation in OSCC [16].However, the status of internal m7G methylation within mRNAs in OSCC remains poorly understood.Therefore, to study the role of m7G methylation in OSCC, we aimed to construct and analyze the transcriptome mapping of internal m7G in mRNAs.The comprehensive bioinformatics analysis showed the difference in mRNA m7G methylation and the expression of mRNAs between OSCC and normal tissues.Moreover, Gene Ontology (GO) enrichment, Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis and construction of proteinprotein interaction (PPI) networks for m7G-modified genes were performed.Ultimately, the biological functions of the screened m7G-modified Baculoviral IAP Repeat Containing 3 (BIRC3) gene in OSCC were clarified in vitro.In conclusion, our study analyzed the characteristics and potential functions of these m7G-modified genes, which can provide potential research value to study the role of m7G methylation in OSCC.

Tissue Collection, Library Construction and Sequencing
Six tissue samples (C1/N1, C2/N2, and C2/N3) from three OSCC patients were collected from the First Affiliated Hospital of Weifang Medical University.This study was approved by the Ethics Committee of Weifang Medical University [17], and informed consent was obtained from all patients.Library construction and sequencing were performed by CloudSeq Inc. (Shanghai, China), and concrete details about that sequencing procedure have been described in the study [17].In brief, by using m7G-specific antibodies to enrich the m7G-methylated RNA fragments and combined with high-throughput sequencing technology, m7G methylation on RNA was located and quantified.

Sequencing Analysis
The raw reads generated were first subjected to Q30 quality control.Then, high-quality clean reads were obtained using cutadapt software (Version: 1.9.3,Berlin, Germany) [18].The high-quality reads were compared with the reference genome using Hisat2 software (Version: 2.0.4,Baltimore, MD, USA) [19].Raw counts at the gene level were obtained as mRNA expression profiles using HTSeq software (Version: 0.9.1, Barcelona, Spain).Finally, fold change ≥2 and p ≤ 0.05 were the criteria used to screen differentially expressed genes between OSCC and normal tissues using edgeR software (Version: 3.16.5,Melbourne, VCP, Australia) [20].The MACS (Bergisch Gladbach, NRW, Germany) default parameters were used to identify methylated sites, and the p-value of peaks was calculated using the Poisson distribution [21].DiffReps package was used to identify differentially methylated sites, and p-value was calculated using the negative binomial test [22].Differentially methylated genes were screened out with default screening criteria of fold change ≥2 and p ≤ 0.00001.

Bioinformatics Analysis
The RCircos package was used to show the distribution of m7G on chromosomes in OSCC and normal tissues [23].The motif of methylation peaks was identified using Dreme software (Brisbane, QLD, Australia); the motif with the smallest E value was selected to be the most significant [24].
GO and KEGG analyses were performed using the Database for Annotation, Visualization, and Integrated Discovery (DAVID; https://david.ncifcrf.gov/)for differentially expressed genes and differentially m7G-modified genes [25].GO and KEGG terms with p ≤ 0.05 were considered statistically significant.In addition, the top 10 terms of GO enrichment and KEGG pathway were shown according to the enrichment score [-log 10 (P)].The background gene set used for GO and KEGG analysis was the homo sapiens genome (HG19).
PPI networks were constructed for differentially m7Gmodified mRNAs using Search Tool for Retrieval of Interacting Genes/Proteins software (STRING; https://cn.string-db.org/)[26]; an interaction score of 0.7 (median confidence) was the standard.In addition, the CytoNCA plugin of Cytoscape software (Version: 3.8.0,Palo Alto, CA, USA) was used to screen hub genes based on betweenness centrality scores [27].

Cell Culture
Five human OSCC cell lines (SCC-9, SCC-25, HN6, HN30, CAL-27) used in this study were obtained from the Shanghai Jiao Tong University School of Medicine, and the Research Resource Identifiers (RRIDs) were listed in previous studies [31][32][33].The HN6, HN30, CAL-27 cell lines were cultured in Dulbecco's modified Eagle's medium (DMEM, GIBCO, NY, USA) with 10% fetal bovine serum (FBS, GIBCO, NY, USA) and 1% penicillin and streptomycin (GIBCO, NY, USA), while the SCC-9 and SCC-25 cell lines were cultured in DMEM/F12 (GiBCO, NY, USA) with 10% FBS and 1% penicillin and streptomycin.Normal control cells were obtained from primary cultures of oral mucosal epithelial cells.Oral mucosal epithelial cells were cultured in keratinocyte serum-free medium (GIBCO, NY, USA) with 0.2 ng/mL of recombinant epidermal growth factor (Thermo Fisher Scientific, Waltham, MA, USA).All cells were maintained at 37 °C in a 5% CO 2 atmosphere and were mycoplasma-free (Mycoplasma PCR Detection Kit, Beyotime, Shanghai, China).The study was carried out in accordance with the guidelines of the Declaration of Helsinki and approved by the Ethics Committee of the Weifang Medical University [17].All cell lines were authenticated shortly before use by the STR profiling technique, carried out by Shanghai Biowing Applied Biotechnology Co., Ltd.

Total RNA Extraction and Quantitative Real-Time PCR (qPCR)
Total RNA was extracted from cultured cells using TRIzol reagent.Total RNA was reverse transcribed to cDNA using the PrimeScript RT kit (TaKaRa, Tokyo, Japan).Target genes were amplified using a real-time fluorescent quantitative PCR kit (Life Technologies, Carlsbad, CA, USA).Gene expression was normalized for all samples using β-actin, and the gene expression was quantified by the comparative threshold cycle (2 −∆∆CT ) method.The primer sequences were as follows:

Transfection of Small Interfering RNA (siRNA)
Small interfering RNA (siRNA) was designed and synthesized (Sangon Biotech, Shanghai, China).The targeted sequence of si-BIRC3 was 5 ′ -CAGUUCGUACAUUUCUUUCAT-3 ′ .Transfection was performed using Lipo8000 TM (Beyotime, Shanghai, China) according to the manufacturer's instructions, and the knockdown efficiency of BIRC3 was subsequently verified by qPCR.CAL-27 cells transfected with si-BIRC3 were used for the following functional assays.

Transwell Migration Assay
CAL-27 cells (2.5 × 10 4 /well) transfected with si-BIRC3 for 24 h were placed in the upper chamber in 150 µL of serum-free DMEM; 600 µL of DMEM containing 10% FBS was added to the lower chamber and incubated in an incubator at 37 °C and 5% CO 2 for 24 h.After that, the surface cells at the bottom of the chamber were fixed with 4% paraformaldehyde (Beyotime, Shanghai, China), stained with 0.5% crystal violet (Beyotime, Shanghai, China), and photographed in 5 randomly selected fields under an anatomical microscope (Leica S8APO Stereo Microscope, Leica Microsystems, Wetzlar, Hessen, Germany).

Statistical Analysis
SPSS 16.0 (IBM Corp., Armonk, NY, USA) was used for statistical analysis, and GraphPad Prism 7.0 (La Jolla, CA, USA) was used to graph the data.An independent samples t-test (two groups) was used to compare the significant differences between the two groups; one-way ANOVA was used to compare multiple groups.The results (p < 0.05) were considered statistically significant.

Overview of Internal m7G Methylation within mRNAs in OSCC and Normal Tissues
To clarify the overview of internal m7G methylation within mRNAs in OSCC and normal tissues, we took the intersection of C1, C2, and C3 (Supplementary Fig. 1A,C).Meanwhile, we took the intersection of N1, N2, and N3 (Supplementary Fig. 1B,D).The analysis showed 9514 unique m7G peaks within 7455 genes in OSCC tissues (Fig. 1A,B).Internal m7G methylation features within mR-NAs in two groups were further analyzed.The internal m7G peaks within mRNAs in both groups were mainly distributed on chromosome 1 (chr1), followed by chromosome 19 (chr19) (Fig. 1C).Furthermore, the distribution of m7G methylation on all chromosomes in two groups was shown using the RCircos package (Fig. 1D).In addition, the percentage of genes with more than four m7G peaks within mRNAs was found to be 34.54% in OSCC tissues, indicating that internal m7G methylation sites are not unique (Fig. 1E).

Analysis of the Regional Distribution of Internal m7G Peaks in OSCC and Normal Tissues
Identification of the regional distribution of m7G methylation showed that m7G methylation occurred in all regions within mRNAs, including coding sequence (CDS), 3 ′ untranslated region (UTR), 5 ′ UTR, Start C and Stop C in OSCC tissues and normal tissues (Fig. 1F,H).In OSCC tissues, the proportion of m7G methylation sites distributed in the 5 ′ UTR region was 15.58%.Meanwhile, in normal tissues, the proportion of m7G methylation sites distributed in the 5 ′ UTR region was 14.60% (Fig. 1F).By analyzing the data on the distribution of m7G methylation sites in the mRNA region in the samples (C1, C2, C3, N1, N2, N3), the results showed that there was no difference in the distribution of m7G methylation sites in all mRNA region in the OSCC tissues compared with the normal tissues (Fig. 1G).

Analysis of Internal m7G Motifs within mRNAs in OSCC and Normal Tissues
The motif sequences were analyzed using Dreme software in OSCC and normal tissues.Notably, GGARG (R = G/A) sequences were highly clustered at the m7G sites (E value of cancer: 7.0 × 10 −12 ; E value of normal: 3.8 × 10 −6 ) in OSCC and normal tissues (Fig. 1I).

Potential Biological Functions of Differentially Methylated Genes in OSCC
Differentially methylated genes (fold change ≥2, p ≤ 0.00001) were identified in OSCC tissues using MeRIPseq.The top 10 hypermethylated genes were listed according to the fold change (Table 1).
To clarify the biological functions of differentially methylated genes in OSCC, GO analysis based on biological processes (BP), cellular components (CC), molecular functions (MF) and KEGG enrichment analysis for these genes were performed.Hyper m7G-modified genes in OSCC were associated with ion transport, multicellular organismal process and transmembrane transport (GO term: BP); an integral component of the plasma membrane, an intrinsic component of the membrane, and cell projection (GO term: CC); transmembrane transporter activity, and ion transmembrane transporter activity (GO term: MF) (Fig. 2A).In contrast, hypo m7G-modified genes in OSCC were related to localization, regulation of signaling and cell communication (GO term: BP); the cytoplasm, intracellular organelle, and the cytosol (GO term: CC); protein binding, enzyme binding, and cytoskeletal protein binding (GO term: MF) (Fig. 2C).Moreover, KEGG pathway analysis showed that hyper m7G-modified in OSCC were enriched in the calcium signaling pathway and ATP-binding cassette (ABC) transporters (Fig. 2B).In contrast, hypo m7Gmodified genes in OSCC were associated with the Hippo and adenosine 5 monophosphate (AMP)-activated protein kinase (AMPK) signaling pathway (Fig. 2D).

Functional Analysis of Differentially Expressed Genes in OSCC
The heatmap showed gene expression profiles identified by transcriptome sequencing (Supplementary Fig. 2).Volcano and scatter plots showed 2185 upregulated genes and 2022 downregulated genes in OSCC tissues (fold change ≥2 and p ≤ 0.05) (Fig. 3A,B).To investigate the role of differentially expressed genes in OSCC, GO and KEGG pathway enrichment analyses of these genes were performed.Upregulated genes in OSCC were closely associated with response to stress, the cell cycle, and the immune system process (GO term: BP); the cytoplasm, extracellular exosome, and nucleoplasm (GO term: CC); protein, DNA replication origin and ATP binding (GO term: MF) (Fig. 3C).In contrast, downregulated genes in OSCC were related to a muscle system process, muscle contraction, and muscle structure development (GO term: BP); contractile fibers, myofibrils, and sarcomeres (GO term: CC); cytoskeletal protein binding, actin binding, and structural constituent of muscle (GO term: MF) (Fig. 3E).Furthermore, KEGG pathway analysis revealed that upregulated genes in OSCC were involved in the p53 signaling pathway and cell cycle (Fig. 3D).In contrast, downregulated genes in OSCC were associated with CGMP-PKG and the calcium signaling pathway (Fig. 3F).

Combined Analysis of Differentially Methylated Genes and Differentially Expressed Genes in OSCC
To investigate the biological role of differentially expressed genes affected by m7G methylation in OSCC, the combined analysis were completed.The result showed 1221 genes with upregulated mRNA expression with m7G hypermethylation sites (referred to as hyper-upregulated genes) and 1645 genes of downregulated mRNA levels with m7G hypomethylation sites (referred to as hypodownregulated genes; fold change ≥2, p ≤ 0.00001).The top 10 hyper-upregulated genes were shown in Table 2.In addition, GO and KEGG pathway analyses were performed on these genes.GO analysis revealed that the hyper-upregulated genes in OSCC were mostly involved in the mitotic cell cycle, apoptotic process, and cell junction (Fig. 4A); KEGG pathway analysis showed these genes were enriched in the cell cycle, cellular senescence and DNA replication (Fig. 4B).In contrast, GO analysis showed that hypo-downregulated genes in OSCC were involved in muscle contraction, stress fiber, and actin-binding (Fig. 4C); KEGG pathway analysis showed that these genes were enriched in the calcium, apelin, and cGMP-PKG signaling pathways (Fig. 4D).To explore the protein interactions of differentially expressed methylation genes, PPI networks of these genes were generated using the STRING software.Among these 1221 hyper-upregulated genes, 792 nodes and 440 edges (p: 1.0 × 10 −16 ) were found (Fig. 4E).Among these 1645 hypo-downregulated genes, 1070 nodes and 482 edges were found (p: 1.0 × 10 −16 ) (Fig. 4F).In addition, hub genes were filtered based on intergenic centrality scores using the CytoNCA plugin of Cytoscape software.The top three core genes screened were TP73, CASP1, and RUNX3 in the hyper-upregulated genes and ACTN2, PPKACA, and MEF2C in hypo-downregulated genes (Fig. 4E,F).

Bioinformatic Characteristics and Biological Behavior of the m7G-modified Gene BIRC3 in OSCC Tissues
By the co-analysis, the heatmap showed the top 10 upregulated genes with high m7G methylation based on fold change (Fig. 5A).The result showed that the highly expressed BIRC3 gene had the highest m7G methylated peaks (logFC: 3.062) among the top 10 screened genes.The m7G peaks were visualized by Integrative Genomics Viewer software (IGV, Version: 2.14.1, Cambridge, MA, USA), and the BIRC3 transcript has higher m7G methylation peaks at the sites (chr11: 102208521-102208840) in OSCC tissues than in normal tissues (Fig. 5B).The STRING software showed that the BIRC3 protein could interact with TNF receptor-associated factor (Fig. 5C).Furthermore, UALCAN (Fig. 5D), TIMER and GEPIA (Supplementary Fig. 3A,B) showed that BIRC3 was highly expressed in HNSCC tissues.Then, the expression of BIRC3 was upregulated in OSCC cell lines (SCC-9, SCC-25, HN6, HN30, and CAL-27), as confirmed by qPCR assay (Fig. 5E).To further investigate the biological function of BIRC3 in OSCC, CAL-27 cells transfected with si-BIRC3 were constructed.The BIRC3 expression in CAL-27 cells was knocked down by the qPCR assay (Fig. 5F).Functionally, the migration of CAL-27 cells transfected with si-BIRC3 was significantly inhibited in the Transwell migration assay (Fig. 5G).The proliferation of CAL-27 cells transfected with si-BIRC3 was significantly decreased, as determined by the CCK-8 assay (Fig. 5H).The relative expression of (epithelial-mesenchymal transition, EMT)-related marker protein E-cadherin was significantly increased in the CAL-27 cells transfected with si-BIRC3.In contrast, the expression of N-cadherin and Vimentin proteins was significantly decreased in the Western blot assay (Fig. 5I).These results suggests that the knockdown of BIRC3 could inhibit the biological behavior of OSCC in vitro, which may be related to the process of EMT.

Discussion
m7G methylation is a common RNA modification that can affect tumor progression, recurrence, and metastasis  [35].Additionally, previous studies have found that m7G methylation usually acts through methyltransferases, and the methyltransferases that have been identified and studied include methyltransferase-like1 (METTL1), WD repeat domain 4 (WDR4), Williams Beuren syndrome chromosome region 22 (WBSCR22) and RNA guanine-7 methyltransferase (RNMT) [35,36].The METTL1/WDR4 complex, identified as an m7G methyltransferase of tRNAs, may influence the metabolic process of RNA, which affects the progression of tumors [14].A recent study reported that m7G-modified tRNA catalyzed by the METTL1/WDR4 complex could stabilize tRNAs, enhancing the transla-tion of oncogenic genes in HNSCC [14].Importantly, METTL1 has been revealed to be a methyltransferase that installs a subset of m7G into mRNAs [11].Whether the METTL1/WDR4 complex in OSCC can affect the gene expression by mediating their internal m7G methylation levels needs to be further explored.
To explore the role of m7G methylation in diseases progession, researchers developed the MeRIP-seq technology based on specific m7G antibodies for transcriptional profiling of m7G methylation in human cells or tissues [11,13,37].Likewise, to understand the status of m7G methylation in OSCC, the m7G MeRIP-seq method was   used to measure the internal m7G methylated levels of mR-NAs in OSCC and normal tissue in this study.However, m7G MeRIP-seq method still has limitations in detecting the internal mRNAs, which might result in the minor overlap in m7G peaks [11,38].This result revealed the difference in m7G methylation between OSCC and normal tissue, suggesting that m7G methylation might be associated with the progress of OSCC.Importantly, the mechanism by which m7G regulates gene expression was closely related to the distribution region of m7G methylation within mR-NAs [39].It was found that gene expression was closely associated with elevated levels of m7G methylation in the CDS region, which might promote mRNA translation efficiency and protein expression [40].Interestingly, it showed that m7G methylation was mainly enriched in the CDS region in OSCC tissues.Notably, the results also showed that there were more m7G peaks in the 5 ′ UTR region of gene transcripts in OSCC tissues than in normal tissues, which might lead to abnormalities in various cellular interactions.For example, elevated m7G methylation in the 5 ′ UTR region that assesses for complementarity with AUG might affect mRNA translational initiation, affecting critical cancer-related gene expression [41].The sequencing analysis showed the distribution ratio of internal m7G methylation in the region of mRNAs, which provides the theoretical basis for further investigating the pathogenic mechanism of m7G methylation in OSCC.
The Bioinformatics function analysis showed that differentially methylated genes were associated with ABC transporters and the AMPK signaling pathway in OSCC tissues.Significantly, these cellular biological functions are associated with cancer progression.OSCC patients often develop multidrug resistance during clinical treatment, which significantly reduces the efficacy [42].Previous studies have shown that high expression of ABC transporter (ABCC1 and ABCG2) promotes multidrug resistance in OSCC [42].Moreover, the AMPK signaling pathway could participate in the proliferation migration in OSCC [43] and is critical in tumor metabolism [44].We hypothesize that hypo m7G-modified genes in OSCC tissues might lead to abnormalities in tumor metabolism, which would be related to cancer progression.Thus, many cellular biological functions of these differentially methylated genes may provide clues to further clarify further the specific role of m7G methylation in OSCC.
This study showed the correlation between mRNA expression and m7G methylation, and 1221 hyperupregulated genes and 1645 hypo-downregulated genes were identified.The hub genes TP73, CASP1, RUNX3, ACTN2, PPKACA and MEF2C were screened from the PPI network constructed for these genes.Previous studies demonstrated these hub genes could affect tumor development [45][46][47].For example, TP73 is highly expressed in OSCC and can influence cell biological processes of tumorigenesis, such as apoptosis, autophagy and metabolism in tumors [45,46].The downregulated expression of ACTN2 in OSCC is associated with increased survival [47].However, these genes' RNA methylation and regulatory roles have been poorly studied.Therefore, the m7G-modified hub genes provides significant value for studying the specific regulatory mechanisms of m7G methylation in OSCC.Notably, we screened the essential m7G-modified gene BIRC3 by co-analysis.
BIRC3 is a member of the inhibitors of apoptosis proteins (IAPs) family, which can be involved in cellular bioprocesses through the regulation of cysteine asparaginase and the signaling of mitotic kinases [48].Notably, the highly expressed BIRC3 is significantly correlated with poor prognosis and chemotherapy resistance among multiple cancers, including non-small cell lung cancer [49], esophageal squamous cell carcinoma [50], HNSCC [51], and colorectal cancer [52].Previous studies have reported that the BIRC3 gene can act as an oncogenic factor in HN-SCC [51].Our sequencing results showed that the highly expressed BIRC3 gene had high m7G methylation, and it was identified that the expression of BIRC3 was higher in OSCC cell lines than in normal cells, consistent with the sequencing results.Furthermore, high-expressed BIRC3 in OSCC was associated with poor prognosis, lymph node metastasis, and radiotherapy resistance [53,54].The expression of BIRC3 up-regulated by circDOCK1 was found to inhibit OSCC cell apoptosis [55].In this study, the function of BIRC3 with high m7G methylation in OSCC was also explored.BIRC3 could promote the ability of proliferation and migration and the progress of EMT.Previous studies found that the Snail transcription factor could interact with BIRC3, thus promoting ovarian carcinogenesis [56].The mechanisms by which BIRC3 affects EMT progression need further investigation.

Conclusions
The m7G methylation difference in internal mRNAs between OSCC and normal tissues were analyzed.GO and KEGG pathway analyses revealed that m7G-modified genes were closely associated with cancer progression.Furthermore, BIRC3 could promote the biological behavior of OSCC in vitro, which may be related to the process of EMT.This article identified critical m7G-modified genes and investigated their potential functions, providing a foundation for further research on the role of m7G methylation in OSCC.

Fig. 1 .
Fig. 1.Overview of internal m7G methylation within mRNAs in oral squamous cell carcinoma (OSCC).(A) Venn diagram showed the unique m7G peaks in OSCC and normal tissues.(B) Venn diagram showed the specific m7G-modified genes in OSCC and normal tissues.(C) Chromosome distribution characteristics of m7G peaks in OSCC and normal tissues.(D) The distribution of m7G at the chromosome was displayed using the Rcircos package in OSCC and normal tissues.(Red: cancer; Blue: normal).(E) The percentage of genes with different m7G peaks in OSCC and normal tissues.(F) Regional distribution of internal m7G methylation sites within mRNAs in OSCC and normal tissues.(CDS, coding sequence; UTR, untranslated region).(G) The bar plot showed the distribution percentage of m7G methylation sites in the mRNA region in the samples (C1, C2, C3, N1, N2, N3).(H) m7G peaks density analysis of m7G methylation in OSCC and normal tissues.(C, cancer; N, normal).(I) The most conserved m7G motifs in OSCC and normal tissues, respectively.

Fig. 2 .
Fig. 2. The bioinformatic function of differentially methylated genes in OSCC.The top 10 Gene Ontology (GO) terms (A) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (B) of hyper m7G-modified genes genes in OSCC.GO enrichment included biological process (BP), cellular component (CC) and molecular function (MF).The top 10 GO terms (C) and KEGG pathways (D) of hypo m7G-modified genes in OSCC.Enrichment score: -log10 (p); Size: gene count.

Fig. 3 .
Fig. 3.The bioinformatic function was analyzed for differentially expressed genes identified in OSCC.(A) The volcano map showed differentially expressed genes in OSCC.The bar chart showed statistic of differential expressed genes between OSCC and normal tissues.(red: upregulated, grey: not changed, blue: downregulated).(B) Scatter plot showed differentially expressed genes.(Red: upregulated, Blue: not changed, Green: downregulated; C: cancer; N: normal).The top 10 GO terms (C) and KEGG pathways (D) of upregulated genes in OSCC.The top 10 GO terms (E) and KEGG pathway (F) on downregulated genes in OSCC.Enrichment score: -log10 (p); Size: gene count.

Fig. 4 .
Fig. 4. Bioinformatics analysis of differentially expressed methylation genes identified by co-analysis.The top 10 GO terms.(A) and KEGG pathways (B) of hyper-upregulated genes in OSCC.The top 10 GO terms (C) and KEGG pathways (D) of hypo-downregulated genes in OSCC.Enrichment score: -log10 (p); Size: gene count.(E) The protein-protein interaction (PPI) network was constructed among 1221 hyper-upregulated genes, and TP73, CASP1 and RUNX3 were identified as hub genes.(F) The PPI network was constructed among 1645 hypo-downregulated genes, and ACTN2, PPKACA and MEF2C were identified as hub genes.The red circles represent the screened hub genes.

Fig. 5 .
Fig. 5.The Bioinformatic characterization and biological functions of the BIRC3 gene in OSCC.(A) The heatmap showed the top 10 upregulated genes with hypermethylation sites.(B) Visualization of the m7G peak site (chr11:102208521-102208840) in the BIRC3 gene in OSCC tissues.(IP: immunoprecipitation; Red: IP; Blue: Input).(C) The PPI network of BIRC3 was constructed by the STRING software.(D) Gene differential expression analysis from UALCAN showed that the mRNA expression of BIRC3 was significantly upregulated in HNSCC tissues.(E) The high expression of BIRC3 in OSCC cell lines was verified by qPCR assay.(F) The knockdown of BIRC3 at the mRNA level in CAL-27 cells was confirmed by qPCR assay.(G) The Transwell migration assay showed the migration of CAL-27 cells transfected with si-BIRC3 was significantly inhibited.(H) The proliferation ability of CAL-27 transfected with si-BIRC3 was significantly reduced in the CCK-8 assay.OD, optical density.(I) The relative expression of EMT marker protein E-cadherin was significantly increased in the CAL-27 cells transfected with si-BIRC3, while the N-cadherin and Vimentin proteins were significantly decreased using Western blot assay.(*: p < 0.05, **: p < 0.01, ***: p < 0.001, ****: p < 0.0001).