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.

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 $\ge$2 and p $\le$ 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 $\ge$2 and p $\le$ 0.00001.

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 $\le$ 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 m7G-modified 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].

Tumor Immune Estimation Resource (TIMER; https://cistrome.shinyapps.io/timer/), the University of Alabama at Birmingham Cancer Data Analysis Portal (UALCAN; https://ualcan.path.uab.edu/) and Gene Expression Profiling Interactive Analysis (GEPIA; http://gepia.cancer-pku.cn/) were used to analyze the expression of the BIRC3 in tumor tissues [28, 29, 30].

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 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 $\beta$-actin, and the gene expression was quantified by the comparative threshold cycle (2${}^{-\mathrm{\Delta }\mathrm{\Delta }\text{CT}}$) method. The primer sequences were as follows: BIRC3: F: 5${}^{\prime }$-TCGCTTGAAAAGACTGGGCT-3${}^{\prime }$, R: 5${}^{\prime }$-CCCGAGATTAGACTAAGTCCCTT-3${}^{\prime }$; $\beta$-actin: F: 5${}^{\prime }$-CATGTACGTTGCTATCCAGGC-3${}^{\prime }$, R: 5${}^{\prime }$-CTCCTTAATGTCACGCACGAT-3${}^{\prime }$ (F: forwards, R: rewards).

Small interfering RNA (siRNA) was designed and synthesized (Sangon Biotech, Shanghai, China). The targeted sequence of si-BIRC3 was 5${}^{\prime }$-CAGUUCGUACAUUUCUUUCAT-3${}^{\prime }$. Transfection was performed using Lipo8000${}^{\text{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.

CAL-27 cells (1.5 $\times$ 10${}^3$/well) transfected with si-BIRC3 for 24 h were added to 96-well plates in triplicate; subsequently, cells were incubated in a 37 °C, 5% CO${}_2$ incubator for 2 h with 10 µL of CCK-8 solution (Dojindo, Kumamoto, Japan). The optical density (OD) at 450 nm was measured daily for 5 days using a microplate reader (Thermo Fisher Scientific, Waltham, MA, USA).

CAL-27 cells (2.5 $\times$ 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).

Western blot analysis was performed as previously described [31]. The primary antibodies included anti-E cadherin (Cat# GB11868; 1:1000; Servicebio, Wuhan , China), anti-N cadherin (Cat# GB11135; 1:1000; Servicebio, Wuhan, China), anti-Vimentin (Cat# GB11192; 1:1000; Servicebio, Wuhan, China) and anti-$\beta$-actin (Cat# GB11001; 1:1000; Servicebio, Wuhan, China); HRP Goat anti-rabbit IgG (Cat# BL003A; 1:10,000; Biosharp, Guangzhou, China) was used as the secondary antibody. Immunostrips were visualized with ECLUltra (New Cell and Molecular Biotech, Suzhou, China), and grayscale values were analyzed with ImageJ software (Version 1.46R, National Institutes of Health, Bethesda, MD, USA) after imaging [34].

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.

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

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.

Identification of the regional distribution of m7G methylation showed that m7G methylation occurred in all regions within mRNAs, including coding sequence (CDS), 3${}^{\prime }$untranslated region (UTR), 5${}^{\prime }$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${}^{\prime }$UTR region was 15.58%. Meanwhile, in normal tissues, the proportion of m7G methylation sites distributed in the 5${}^{\prime }$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).

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 $\times$ 10${}^{-12}$; E value of normal: 3.8 $\times$ 10${}^{-6}$) in OSCC and normal tissues (Fig. 1I).

Differentially methylated genes (fold change $\ge$2, p $\le$ 0.00001) were identified in OSCC tissues using MeRIP-seq. 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 m7G-modified genes in OSCC were associated with the Hippo and adenosine 5 monophosphate (AMP)-activated protein kinase (AMPK) signaling pathway (Fig. 2D).

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: –log${}_{10}$ (p); Size: gene count.

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 $\ge$2 and p $\le$ 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).

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: –log${}_{10}$ (p); Size: gene count.

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 hypo-downregulated genes; fold change $\ge$2, p $\le$ 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).

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: –log${}_{10}$ (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.

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 $\times$ 10${}^{-16}$) were found (Fig. 4E). Among these 1645 hypo-downregulated genes, 1070 nodes and 482 edges were found (p: 1.0 $\times$ 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).

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

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

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.