Skip to main content
Download PDF

High-throughput microarray reveals the epitranscriptome-wide landscape of m6A-modified circRNA in oral squamous cell carcinoma

BMC Genomics volume 23, Article number: 611 (2022) Cite this article

Abstract

Background

Emerging transcriptome-wide high-throughput screenings reveal the landscape and functions of RNAs, such as circular RNAs (circRNAs), in human cancer. In addition, the post-transcriptional RNA internal modifications, especially N6-methyladenosine (m6A), greatly enrich the variety of RNAs metabolism. However, the m6A modification on circRNAs has yet to be addressed.

Results

Here, we report an epitranscriptome-wide mapping of m6A-modified circRNAs (m6A-circRNA) in oral squamous cell carcinoma (OSCC). Utilizing the data of m6A methylated RNA immunoprecipitation sequencing (MeRIP-seq) and m6A-circRNAs microarray, we found that m6A-circRNAs exhibited particular modification styles in OSCC, which was independent of m6A-mRNA. Besides, m6A modification on circRNAs frequently occurred on the long exons in the front part of the coding sequence (CDS), which was distinct from m6A-mRNA that in 3’-UTR or stop codon.

Conclusion

In conclusion, our work preliminarily demonstrates the traits of m6A-circRNAs, which may bring enlighten for the roles of m6A-circRNAs in OSCC.

Highlights

1. m6A-circRNAs exhibited their particular modification style in OSCC, which was independent of m6A-mRNA.

2. m6A on circRNAs frequently occurred on the long exons in the front part of CDS, which was distinct from m6A-mRNA that in 3’-UTR or stop codon.

Peer Review reports

Background

Oral squamous cell carcinoma (OSCC) acts as the most common cancers in head and neck, accounting for 90% of oral malignant tumor [1]. In addition to local hyperplasia, tissues erosion and dysfunction, OSCC frequently results in lymphatic metastasis to the neck area. Although comprehensive therapeutic schedules have made remarkable progress, including surgical excision, chemotherapy and radiotherapy, the survival and prognosis of OSCC sufferer remain poor [2]. The latest researches show that the tumorigenesis of OSCC is a complicated process involved in diverse alterations of genetic and epigenetic. Therefore, unremitting exploration targeting the initiation and oncogenesis of OSCC could provide valuable directions for precise treatment.

Circular RNAs (circRNAs) are groups of covalently closed loop transcripts in noncoding transcriptome [3]. CircRNAs are characterized by a continuous loop generated from the back-splicing of linear RNA, which is different from linear RNA. Benefited from covalently-bonded RNA molecule, circRNAs could resist the digestion of RNA enzyme and stably exist in the cellular microenvironment, which is responsible for their high abundancde in different species. In OSCC, the functions of circRNAs have been preliminarily uncovered. Certain circRNAs regulate the OSCC cellular glycolysis metabolism, e.g. circ_0000140 [4] and hsa_circRNA_100290 [5]. Certain circRNAs regulate the epithelial-mesenchymal transition (EMT) of OSCC, e.g. hsa_circ_0009128 [6] and circEPSTI1 [7]. Moreover, some circRNAs serve as potential diagnostic predictors or biomarkers for OSCC, e.g. hsa_circ_0008309 [8] and hsa_circ_0003829 [9]. Therefore, the evidence illustrates the critical roles participating in OSCC initiation and progression.

N6-methyladenosine (m6A) is one of the most abundant internal modifications in eukaryotic messenger RNA (mRNA) that exerts essential roles in mRNA fate, including mRNA stability, splicing and translation [10]. RNA m6A modification is a well-known chemical modification occurred in the sixth nitrogen. The specific transcriptome modifications have been proved to influence the cellular pathophysiology, thereby regulating tumorigenesis. For example, RNA modification enzyme methyltransferase-like 3 (METTL3) is upregulated in OSCC cohorts and the high expression of METTL3 is associated with poor prognosis. Functionally, METTL3 promotes cell proliferation, migration, invasion and self-renewal through promoting BMI1 translation under the cooperation with IGF2BP1 in OSCC [11]. Besides, m6A demethylase fat mass and obesity-associated protein (FTO) knockdown induces the downregulation of m6A-contained eIF4G1, which is captured by YTHDF2, and enhances the autophagic flux, thus inhibiting OSCC tumorigenesis. With the release of new research work about m6A, there are reasons to believe that the mechanism by which m6A regulates OSCC could be identified.

Given that m6A modification regulates the fate of RNA, while circRNAs are a group of covalent closed-loop RNA, there could be a pivotal connection between m6A and circRNAs [12]. In fact, current researches prove this hypothesis in human cancers. However, the role of m6A on the function of circRNAs has yet to be addressed. Here, our team utilized the data of MeRIP-seq and m6A-circRNAs epitranscriptomic microarray analysis to investigate the epitranscriptome-wide mapping of m6A-modified circRNA in OSCC. Our work focused on the traits of m6A-circRNAs and probed into the association of m6A-circRNAs and m6A-mRNA, which may bring enlighten for the roles of m6A-circRNAs in OSCC.

Results

m6A-circRNA epitranscriptomic microarray analysis revealed the m6A-circRNAs profile in OSCC

m6A-circRNA epitranscriptomic microarray analysis was performed according to the workflow using SCC25 cells and HOK cells. As regarding to the microarray raw data, the differentially m6A-methylated circRNAs were expressed as ‘m6A methylation level’ and ‘m6A quantity’, which was different from conventional circRNA microarray analysis. Thus, in our m6A-circRNA epitranscriptomic microarray analysis, the differentially expressed m6A-circRNAs were presented by two aspects, including ‘m6A-circRNA Methylation Level’ and ‘m6A-circRNA Quantity’. The screening threshold was set to p-value < 0.05 and fold change > 1.5 fold. Based on this threshold, there were 104 up-regulated m6A-circRNA and 145 down-regulated m6A-circRNA on the basis of ‘m6A-circRNA Methylation Level’, meanwhile, there were 2586 up-regulated m6A-circRNA and 472 down-regulated m6A-circRNA on the basis of ‘m6A-circRNA Quantity’. The raw data of ‘m6A-circRNA Methylation Level’ and ‘m6A-circRNA Quantity’ were presented by Heat Map (Fig. 1A, B), Volcano Plot (Fig. 1C, D), and Scatter Plot t (Fig. 1E F). Taken together, m6A-circRNA epitranscriptomic microarray analysis revealed the m6A-circRNAs profile in OSCC.

Fig. 1
figure 1

m6A-circRNA epitranscriptomic microarray analysis revealed the m6A- circRNAs profile in OSCC. A, B Heat Map showed the differentially expressed m6A-circRNAs, including ‘m6A-circRNA Methylation Level’ and ‘m6A-circRNA Quantity’. OSCC cells were SCC25 cells and normal cells were HOK cells. C, D Volcano Plot showed the differentially expressed m6A-circRNAs. Based on this threshold (p-value < 0.05, fold change > 1.5 fold), there were 104 up-regulated m6A-circRNA and 145 down-regulated m6A-circRNA on the basis of ‘m6A-circRNA Methylation Level’ (C). And, there were 2586 up-regulated m6A-circRNA and 472 down-regulated m6A-circRNA on the basis of ‘m6A-circRNA Quantity’ (D). E, F Scatter Plot displayed the correlation distribution of m6A-circRNA in SCC25 cells and HOK cells

Full size image

The significantly expressed m6A-circRNAs

Given that m6A-circRNA epitranscriptomic microarray analysis discovered hundreds or thousands of circRNAs with different ‘m6A methylation level’ or ‘m6A quantity’, further analysis focused on their intersection to determine the potential m6A-methylated circRNAs with significant expression. Circos plot showed the locations of m6A-cricRNAs on human chromosomes, including ‘m6A methylation level’ or ‘m6A quantity’(Fig. 2A and B). In the schematic, the outermost layer was chromosome map of human genome. The inner green layer indicated all the m6A-cricRNAs detected by m6A-circRNA epitranscriptomic microarray. The inner blue layer indicated the m6A-cricRNAs from OSCC cells (SCC25), and the innermost red layer indicated the m6A-cricRNAs from normal control cells (HOK) cells.

Fig. 2
figure 2

The significantly expressed m6A-circRNAs. A, B Circos plot showed the locations of m6A-cricRNAs on human chromosomes, including ‘m6A methylation level’ or ‘m6A quantity’. The outermost layer was chromosome map of human genome. The inner green layer indicated all the m6A-cricRNAs detected by m6A-circRNA epitranscriptomic microarray. The inner blue layer indicated the m6A-cricRNAs from OSCC cells (SCC25). The innermost red layer indicated the m6A-cricRNAs from normal control cells (HOK) cells. C, D Venn diagram showed the up-regulated 100 m6A-circRNAs and 90 down-regulated m6A-circRNAs. Several up-regulated or down-regulated m6A-circRNAs were randomly chosen and listed in the diagram

Full size image

In consideration of the fact that m6A-modified circRNAs were categorized and detected based on ‘m6A methylation level’ and ‘m6A quantity’, our team took the intersection of ‘m6A methylation level’ and ‘m6A quantity’ to zoom out the scope, thereby ascertaining the up-regulated or down-regulated m6A-circRNAs. Venn diagram showed the up-regulated m6A-circRNAs (100) (Fig. 2C) and down-regulated m6A-circRNAs (90) (Fig. 2D). Besides, several up-regulated or down-regulated m6A-circRNAs were randomly chosen and listed in the diagram. Moreover, among these screened m6A-circRNAs, the significantly up-regulated (Top 10) m6A-circRNAs were shown in Table 1 and down-regulated m6A-circRNAs were shown in Table 2.

Table 1 Top 10 up-regulated m6A-circRNAs in OSCC
Full size table
Table 2 Top 10 down-regulated m6A-circRNAs in OSCC
Full size table

The genomic characteristic of m6A-circRNAs in OSCC

An essential literature by Dr. Alan C Mullen [13] (2017) has raised an enlightening point that m6A-modified circRNAs exhibit characteristic rules, thus, we investigate the genomic characteristic of m6A-circRNAs in OSCC with the help of MeRIP-Seq and m6A-circRNA epitranscriptomic microarray. The differentially expressed m6A-modified circRNAs, including up-regulated (100 m6A-circRNAs) and down-regulated (90 m6A-circRNAs), were incorporated into the research. Firstly, exons quantity of up-regulated (Fig. 3A) and down-regulated (Fig. 3B) m6A-modified circRNAs were counted. Results illustrated that the exons of m6A-circRNAs mainly concentrated on 1–3 exons. Moreover, the length (spliced length, not genomic length) was counted by stages (100 bp). Results indicated that the length of m6A-circRNAs primarily distributed at 200–500 bp (Fig. 3C, D). In other words, the shorter circRNAs with a small number of exons were more easily methylated. Our previous research performed MeRIP-seq to detect the m6A profile in OSCC cells [14]. Taking MeRIP-seq as reference, we could find the difference and connection between m6A-mRNA and m6A-circRNAs. Metagene profiles showed the enrichment of m6A modification on mRNA and circRNA across RNA transcriptome. In up-regulated (Fig. 3E) and down-regulated (Fig. 3F) m6A-circRNAs, the m6A modification sites were primarily located in the front part of CDS, which was distinct from m6A-mRNA that in 3’-UTR or stop codon. Taken together, the data of m6A-circRNA epitranscriptomic microarray revealed the genomic characteristic of m6A-circRNAs in OSCC.

Fig. 3
figure 3

The genomic characteristic of m6A-circRNAs in OSCC. A, B The exons quantity of up-regulated and down-regulated m6A-modified circRNAs were counted in differentially expressed m6A-modified circRNAs, including up-regulated (100 m6A-circRNAs) and down-regulated (90 m6A-circRNAs). C, D The length (spliced length, not genomic length) was counted by stages count (~ 100 bp). E, F Metagene profiles of enrichment of m6A modification sites on mRNA and circRNA of gene locus across RNA transcriptome. 5’-untranslated regions (5’-UTR), coding sequences (CDS), 3’-untranslated regions (3’-UTR)

Full size image

The association of m6A-circRNAs and m6A-mRNA in OSCC

To date, the major field of m6A epitranscriptome research focused on mRNA m6A modification. In previous published literature, our team had performed the MeRIP-seq on OSCC prior to present m6A-circRNAs epitranscriptomic microarray analysis [14]. Utilizing the MeRIP-seq data and m6A-circRNAs epitranscriptomic microarray data, we found several traits of m6A-circRNAs and probed into the association of m6A-circRNAs and m6A-mRNA in OSCC. Combined with the existing literature and our findings, we confirmed that exons of certain gene could generate both pre-mRNA transcripts with m6A modification and exon-derived circRNAs with m6A modification (Fig. 4A). Moreover, the m6A modification was installed by identical m6A methyltransferase complex. Indeed, there might be other types of m6A-modified circRNAs, e.g. intron-derived circRNAs or intron/exon-derived circRNAs, which were elusive and will be investigate in our further research.

Fig. 4
figure 4

The association of m6A-circRNAs and m6A-mRNA in OSCC. A Schematic diagram illustrated the biogenesis of m6A-circRNAs and m6A-mRNA in OSCC cells. The m6A modification was installed by identical m6A methyltransferase complex (m6A writers). Red A indicated the m6A modification site. B Schematic diagram presented the genomic structure of m6A-mRNA, unmethylated circRNAs and m6A-circRNAs. C c-Myc mRNA acted as a representative for m6A-mRNA, which was methylated at 3’-UTR and without circRNA derived. D circHIPK3 (hsa_circ_0000284) represented the group of circRNAs without m6A modification in neither circular transcripts nor host gene. E CircFOXM1 (hsa_circ_0025039) represented the group of circRNAs without m6A modification in circRNA transcript, however, their host genes were m6A-modified at 3’-UTR. F CircKRT5 (hsa_circ_0026457) represented the group of circRNAs only m6A-modified in circRNA transcript, instead of host genes. G CircIPO9 (hsa_circ_0015936) represented the group of circRNAs m6A-modified both in circRNA transcript and their host genes

Full size image

All the profound analysis was performed based on our MeRIP-seq data and m6A-circRNAs epitranscriptomic microarray data. According to the difference of methylated location, we concluded several subgroups of m6A-modified circRNAs (Fig. 4B). Firstly, a representative non-circRNA-derived mRNA (c-Myc mRNA) with m6A modification was selected as negative control (Fig. 4C). The m6A modification site was located in the 3’-UTR of c-Myc mRNA. Besides, a circRNA-derived gene (HIPK3) was selected as positive control, whose exon-2 was cyclized to be circHIPK3 (hsa_circ_0000284) (Fig. 4D). Our MeRIP-seq data revealed that there wasn’t any remarkable m6A modified site in the HIPK3 mRNA. The subgroup of circRNA, such as circHIPK3, is a more common type of transcripts without m6A modification neither in circRNA nor host gene.

Furthermore, we investigated the interaction of m6A modification with circRNAs or their host genes. We noticed a feature that whether circRNAs were m6A-modified wasn’t related to their host genes methylated or not. For example, circFOXM1 (hsa_circ_0025039) was a non-m6A-modified circRNA, however, its host gene (FOXM1) was remarkably m6A-modified at 3’-UTR (Fig. 4E). In addition to the circRNAs derived from m6A-modified host genes, another group of m6A-modified circRNAs was cyclized from unmethylated host genes. For example, crcKRT5 (hsa_circ_0026457) was a m6A-modified circRNA, however, there wasn’t any m6A site in its host gene (KRT5) (Fig. 4F). Moreover, the m6A modification could be installed both circRNA transcript and host genes. For instance, circIPO9 (hsa_circ_0015936) was a m6A-modified circRNA, besides, its host gene (IPO9) was m6A-modified at 3’-UTR (Fig. 4G).

Overall, based on MeRIP-seq and m6A-circRNAs epitranscriptomic microarray data, our findings illustrated that m6A-circRNAs exhibited their particular modification style in OSCC, which was independent of m6A-mRNA.

GO and KEGG pathway analysis

To investigate the mechanisms correlated to m6A-circRNAs in OSCC, Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of the host genes of differentially expressed circRNAs were performed. GO terms included the biological process (BP), cellular component (CC), and molecular function (MF) categories. The top 10 enriched GO terms in BP, CC and MF were ‘keratinization’, ‘Keratin filament’ and ‘glutathione binding’ respectively (Fig. 5A). The top 10 KEGG pathways were shown as following (Fig. 5B) and the host genes of differentially expressed circRNAs were mainly associated with necroptosis. Given that circRNAs could act as miRNA sponge to harbor their downstream miRNAs, thereby releasing the fettered target mRNAs to modulate the cancer progression. The regulation pattern was regarded as competing endogenous RNA (ceRNA) (Supplementary Info: supplementary Figure S2).

Fig. 5
figure 5

GO and KEGG Pathway Analysis. A Gene Ontology (GO) enrichment analysis of biological process (BP), cellular component (CC), and molecular function (MF). The vertical axis stands for the enrichment score of GO terms. B Bulb map displayed the top 10 KEGG pathways for host genes of differentially expressed circRNAs. The GO and KEGG analysis got permission of KEGG Database [15,16,17]. The size of the dot indicates the gene counts enriched in the pathway. The of the dot indicates the significance (p value) of the enriched pathway

Full size image

Discussion

As an indispensable part of epigenetic regulation, m6A modification is suggested as a crucial regulator participating in human cancer tumorigenesis. Chemical modification of m6A on RNAs is an efficient way of regulating molecular function, which influences the downstream pathways. Increased knowledge of m6A modification of noncoding RNA highlighted their effect on gene expression. Modifications on RNA contribute to the post-transcriptional regulation of RNA fate. Here, our research focused on the m6A-modified circRNA (m6A-circRNA) epitranscriptome-wide mapping in OSCC, which may provide valuable references for tumor epigenetics research.

CircRNAs have been identified to play important roles in tumorigenesis by virtue of their stability and enzyme-resistance. High-throughput sequencing for circRNAs reveals that many genes, previously considered as protein-coding genes, can produce circRNAs by back-splicing. In a variety of tumor subgroup, circRNAs participate in the cell differentiation, proliferation, energy metabolism and chemoradiotherapy resistance. As the research of circRNA progressed, the functions of circRNAs have been gradually investigated. The achievement is acquired not only in solid cancers, but also in OSCC. Our previous research found that, in OSCC, circUHRF1 up-regulated in OSCC and accelerated the tumorigenesis [18]. Up to now, the major regulatory type of circRNAs is competing endogenous RNA (ceRNA), by which circRNAs sponge micro RNAs to relieve the activity restrain of target proteins, e.g. circHIPK3 [19], circFOXO3 [20] and circIGHG [21]. In addition to ceRNA, more regulation manners of circRNAs have been reported. For instance, circ-ZNF609 is associated with heavy polysomes and could be translated into protein via splicing-dependent and cap-independent manner [22]. Thus, with the further research, the modes of circRNAs will be reported more and more.

Alan C. Mullen et al. (2017) firstly reported the profile of m6A modifications on circRNAs through a computational pipeline (AutoCirc) tool [13]. In this study, researchers defined thousands of m6A-circRNAs with cell-type-specific expression. Besides, m6A-circRNAs share identical m6A readers and writers with mRNAs, however, m6A-circRNAs are frequently derived from exons that unmethylated in mRNAs. Moreover, this study proposed another critical trait that the same exons methylated both on mRNA and m6A-circRNAs exhibit less stability. Overall, this groundbreaking work makes a valuable contribution to the field of m6A-modified circRNAs.

Here, our study performed the newely-developed m6A-circRNA epitranscriptomic microarray to detect the profile of m6A-modified circRNAs in OSCC. In our results, the differentially expressed m6A-circRNAs were presented by two aspects, including ‘m6A-circRNA Methylation Level’ and ‘m6A-circRNA Quantity’ based on threshold of p-value < 0.05 and fold change > 1.5 fold. Finally, there were 104 up-regulated m6A-circRNA and 145 down-regulated m6A-circRNA on the basis of ‘m6A-circRNA Methylation Level’, meanwhile, there were 2586 up-regulated m6A-circRNA and 472 down-regulated m6A-circRNA on the basis of ‘m6A-circRNA Quantity’. In previous study, our team had performed the MeRIP-seq on OSCC and thus we analyzed the genomic characteristic of m6A-circRNAs in OSCC utilizing the data of MeRIP-Seq and m6A-circRNA epitranscriptomic microarray. The differentially expressed m6A-modified circRNAs were selected into the research, including up-regulated (100 m6A-circRNAs) and down-regulated (90 m6A-circRNAs). The role of these up-regulated/down-regulated m6A-circRNAs is a promising research direction in OSCC progression. Combining with existing literature, the circRNAs with m6A modification display critical roles in tumor pression, which need more investigation.

As regarding to the genomic characteristic of m6A-circRNA in OSCC, we found that the exons of m6A-circRNAs mainly concentrated on 1–3 exons, and the length primarily distributed at 200–500 bp. Moreover, shorter circRNAs with a small number of exons were more easily methylated. Then, we analyzed the enrichment of m6A modification on mRNA and circRNA across RNA transcriptome. Results indicated that the m6A modification sites of m6A-circRNAs were primarily located in the front part of CDS, which was distinct from m6A-mRNA that in 3’-UTR or stop codon. Taken together, the genomic characteristic of m6A-circRNAs in OSCC was preliminarily investigated with the help of MeRIP-Seq and m6A-circRNA epitranscriptomic microarray.

To date, emerging literatures report the cross-talking of m6A and circRNA in human cancer [23]. Within the scope of epigenetics, RNA modifications and circRNAs are two rapidly expanding fields, and increasing number of researchers are beginning to turn their attention in this direction. For example, the m6A modification of circNSUN2 increases its export to the cytoplasm and circNSUN2 enhances the stability of HMGA2 mRNA to promote colorectal carcinoma metastasis progression, forming a circNSUN2/IGF2BP2/HMGA2 RNA-protein ternary complex [24]. In another example, METTL3 and YTHDC1 control the circ-ZNF609 accumulation and back-splicing reaction [25]. Moreover, in colorectal cancer, METTL3 induced the overexpression of circ1662 by binding the flanking sequences through installing its m6A modification [26]. In spite of the evidence showing essential roles of m6A and circRNA on human cancer, the comprehensive analysis of m6A-circRNA at epitranscriptome-wide is still absent. Actually, the direct evidence that in support of the function of m6A-circRNA are worth looking forward.

For the limitations in present, there was one OSCC cell line (SCC25) in the microarray, which is the limitations and insufficient for present research. Besides, despite this sequencing data of m6A-circRNA microarray and MeRIP-Seq, more cellular biochemistry experimental data are also needed.

Conclusion

In conclusion, this research reveals the epitranscriptome-wide mapping of m6A-modified circRNA in OSCC. The findings expand our understanding of circRNAs and enrich the roles of m6A-circRNAs in OSCC, providing potential resolution strategy for OSCC targeted therapy.

Materials and methods

Human m6A-circRNA epitranscriptomic microarray analysis

The m6A-circRNA Epitranscriptomic microarray analysis was performed by Aksomics Inc (KangChen Bio-tech, Shanghai, China). Total RNA from cell samples (SCC25 cells, HNOK cells) was quantified using the NanoDrop ND-1000. The microarray hybridization was performed based on the Arraystar’s standard protocol. In brief, total RNAs were immunoprecipitated with anti-m6A rabbit polyclonal antibody (Synaptic Systems, 202,003). The unmodified RNAs were recovered from the supernatant as ‘Sup’. The m6A-modified RNAs were eluted from the immunoprecipitated magnetic beads as the ‘IP’. Then, ‘Sup’ and ‘IP’ RNAs were administrated with RNase R (Epicentre, Inc.), and subsequently labeled with Cy3 and Cy5 respectively as cRNAs s in separate reactions by Arraystar Super RNA Labeling Kit (Arraystar, AL-SE-005). The cRNAs were combined and hybridized on Arraystar Human circRNA Epitranscriptomic Microarray (8 × 15 K, Arraystar). Slides were washed and the arrays were scanned by an Agilent Scanner G2505C in two-color channels.

Array images were analyzed by Agilent Feature Extraction software (version 11.0.1.1). Raw intensities of Sup (Supernatant, Cy3-labelled) and IP (Immunoprecipitated, Cy5-labelled) were normalized with average of log2-scaled Spike-in RNA intensities. After Spike-in normalization, the probe signals having Present (P) or Marginal (M) QC flags in at least 3 out of 6 samples were retained for further ‘m6A methylation level’ and ‘m6A quantity’ analyses.

‘m6A methylation level’ was calculated for the percentage of modification based on the IP (Cy5-labelled) and Sup (Cy3-labelled) normalized intensities. ‘m6A quantity’ was calculated for the m6A methylation amount based on the IP (Cy5-labelled) normalized intensities. Differentially m6A-methylated circRNAs between two comparison groups were identified by filtering with the fold change and statistical significance (p-value) thresholds. Hierarchical Clustering was performed to show the distinguishable m6A-methylation pattern among samples.

Workflow of m6A-circRNA epitranscriptomic microarray analysis

The workflow of human m6A-circRNA microarray analysis in OSCC cells was presented in the graphical representation (Figure S1), including RNA extraction, quality control (QC), library construction and data analysis. There were 3 pairs of samples for the m6A-circRNA epitranscriptomic microarray analysis, including OSCC cells (3 independent samples, SCC25 cells) and normal cells (3 independent samples, HOK cells). After the m6A immunoprecipitation, the m6A-modified RNAs were eluted from the immunoprecipitated magnetic beads were marked as ‘IP’ group (immunoprecipitated RNAs), and the unmodified RNAs were recovered from the supernatant and marked as ‘Sup’ group (supernatant unmodified RNAs) (Supplementary Info: supplementary Figure S1A). Regarding the raw data, the statistic analysis was performed from two aspects, including m6A methylation level and m6A quantity (Supplementary Info: supplementary Figure S1B). The data from m6A-circRNA epitranscriptomic microarray analysis reveals the landscape of m6A-modified circRNAs in OSCC.

m6A-circRNA data analysis

The ‘m6A methylation level’ for a transcript was calculated as the percentage of modified RNA (modified %) in all RNAs based on the IP (Cy5-labelled) and Sup (Cy3-labelled) normalized intensities:

$$\text{modified}{\%}= \frac{\text{modified RNA}}{\text{Total RNA}} = \frac{\text{IP}}{\text{IP}+\text{Sup}}$$
$$= \frac{{\varvec{I}\varvec{P}}_{Cy5\;normalized\;intensity}}{{\varvec{I}\varvec{P}}_{Cy5\;normalized\;intensity}+{\varvec{S}\varvec{u}\varvec{p}}_{Cy3\;normalized\;intensity}}$$

Postscript: Raw intensities of IP (immunoprecipitated, Cy5-labelled) and Sup (supernatant, Cy3-labelled) were normalized with average of log2-scaled Spike-in RNA intensities.

$$\text{log2}\left({\varvec{I}\varvec{P}}_{Cy5\;normalized\;intensity}\right)=\text{log2}\left({\varvec{I}\varvec{P}}_{Cy5\;raw}\right)-\text{Average}\left[\text{log2}\left({\varvec{I}\varvec{P}}_{spike-in\_Cy5\;raw}\right)\right]$$
$$\text{log2}\left({\varvec{S}\varvec{u}\varvec{p}}_{Cy3\;normalized\;intensity}\right)=\text{log2}\left({\varvec{S}\varvec{u}\varvec{p}}_{Cy3\;raw}\right)-\text{Average}\left[\text{log2}\left({\varvec{S}\varvec{u}\varvec{p}}_{spike-in\_Cy3\;raw}\right)\right]$$

The “m6A quantity” was calculated for the m6A methylation amount of each transcript based on the IP (Cy5-labelled) normalized intensities.

$$\text{Sample m6A quantity}= \text{Sample}\;{\varvec{I}\varvec{P}}_{Cy5\;normalized\;intensity}$$

Postscript: Raw intensities of IP (Cy5-labelled) were normalized by average of log2-scaled Spike-in RNA intensities.

$$\text{log2}\left({\varvec{I}\varvec{P}}_{Cy5\;normalized\;intensity}\right)=\text{log2}\left({\varvec{I}\varvec{P}}_{Cy5\;raw}\right)-\text{Average}\left[\text{log2}\left({\varvec{I}\varvec{P}}_{spike-in\_Cy5\;raw}\right)\right]$$

Cells and culture

OSCC cells (SCC25 cells) were provided by ATCC (American Type Culture Collection, Manassas, VA, USA) and cultured in DMEM Medium supplemented with fetal bovine serum (FBS, 10%), 100 U/ml penicillin, 100 µg/ml streptomycin. Normal cells (Human Oral Keratinocytes, HOK, Catalog No. 2610, ScienCell) were provided by ScienCell (San Diego, California, USA). HOK cells were cultured in Oral Keratinocyte Medium (OKM, Cat. No. 2611, ScienCell) recommended by ScienCell in vitro. Cells were incubated in a 37 °C humidified incubator with 5% CO2.

m6A methylated RNA immunoprecipitation sequencing (MeRIP-seq)

The total RNA was extracted from cell samples, and then divided into two groups, including Input control sample and immunoprecipitation (IP) sample. The RNAs were firstly spliced into ~ 100nt fragments. IP samples provide unbiased measurements of methylated RNA fragments with specific m6A antibodies. Meanwhile, the Input control sample reflected the abundance of basic RNA enrichment. The reference genome was hg38 gencode. Through library construction, high-throughput sequencing and bioinformatics analysis, mapping of whole transcriptome m6A location was generated by Jiayin Biotechnology Ltd. (Shanghai, China).

Statistical analysis

Differentially m6A-methylated RNAs between two comparison groups were identified by filtering with the fold change (> 1.5 fold) and statistical significance (p-value < 0.05) thresholds. Data were shown as means ± standard deviation (SD).

Availability of data and materials

The raw data that support the findings of this study has been deposited into NCBI Gene Expression Omnibus (GEO): accession number GSE198105. Researchers may view the GSE198105 study at: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE198105, and GSE197457 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE197457).

Abbreviations

m6A:

N6-methyladenosine

circRNA:

circular RNA

References

  1. Mendenhall WM, Holtzman AL, Dagan R, Bryant CM, Hitchcock KE, Amdur RJ, Fernandes RP. Current Role of Radiotherapy in the Management of Oral Cavity Squamous Cell Carcinoma. Craniomaxillofacial Trauma Reconstruction. 2021;14(1):79–83.

    Article  Google Scholar 

  2. Qian JM, Schoenfeld JD. Radiotherapy and Immunotherapy for Head and Neck Cancer: Current Evidence and Challenges. Front Oncol. 2020;10:608772.

    Article  Google Scholar 

  3. Wu J, Qi X, Liu L, Hu X, Liu J, Yang J, Yang J, Lu L, Zhang Z, Ma S, et al. Emerging Epigenetic Regulation of Circular RNAs in Human Cancer. Mol Ther Nucleic Acids. 2019;16:589–96.

    Article  CAS  Google Scholar 

  4. Guo J, Su Y, Zhang M. Circ_0000140 restrains the proliferation, metastasis and glycolysis metabolism of oral squamous cell carcinoma through upregulating CDC73 via sponging miR-182-5. Cancer Cell Int. 2020;20:407.

    Article  CAS  Google Scholar 

  5. Chen X, Yu J, Tian H, Shan Z, Liu W, Pan Z, Ren J. Circle RNA hsa_circRNA_100290 serves as a ceRNA for miR-378a to regulate oral squamous cell carcinoma cells growth via Glucose transporter-1 (GLUT1) and glycolysis. J Cell Physiol. 2019;234(11):19130–40.

    Article  CAS  Google Scholar 

  6. Zhang H, Wang Z, Zhang Z. Hsa_circ_0009128 mediates progression of oral squamous cell carcinoma by influencing MMP9. Oral Dis. 2021. https://doi.org/10.1111/odi.14019.

  7. Wang J, Jiang C, Li N, Wang F, Xu Y, Shen Z, Yang L, Li Z, He C. The circEPSTI1/mir-942-5p/LTBP2 axis regulates the progression of OSCC in the background of OSF via EMT and the PI3K/Akt/mTOR pathway. Cell Death Dis. 2020;11(8):682.

    Article  CAS  Google Scholar 

  8. Li B, Wang F, Li X, Sun S, Shen Y, Yang H. Hsa_circ_0008309 May Be a Potential Biomarker for Oral Squamous Cell Carcinoma. Disease Markers. 2018;2018:7496890.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Zhang H, Shen Y, Zhang B, Qian M, Zhang Y, Yang H. Hsa_circ_0003829 serves as a potential diagnostic predictor for oral squamous cell carcinoma. J Int Med Res. 2020;48(9):300060520936880.

    Article  CAS  PubMed  Google Scholar 

  10. Zhao W, Qi X, Liu L, Ma S, Liu J, Wu J. Epigenetic Regulation of m(6)A Modifications in Human Cancer. Mol Ther Nucleic Acids. 2020;19:405–12.

    Article  CAS  Google Scholar 

  11. Liu L, Wu Y, Li Q, Liang J, He Q, Zhao L, Chen J, Cheng M, Huang Z, Ren H, et al. METTL3 Promotes Tumorigenesis and Metastasis through BMI1 m(6)A Methylation in Oral Squamous Cell Carcinoma. Mol Ther. 2020;28(10):2177–90.

    Article  CAS  Google Scholar 

  12. Chen C, Guo Y, Guo Y, Wu X, Si C, Xu Y, Kang Q, Sun Z. m6A Modification in Non-Coding RNA: The Role in Cancer Drug Resistance. Front Oncol. 2021;11:746789.

    Article  Google Scholar 

  13. Zhou C, Molinie B, Daneshvar K, Pondick JV, Wang J, Van Wittenberghe N, Xing Y, Giallourakis CC, Mullen AC. Genome-Wide Maps of m6A circRNAs Identify Widespread and Cell-Type-Specific Methylation Patterns that Are Distinct from mRNAs. Cell Reports. 2017;20(9):2262–76.

    Article  CAS  Google Scholar 

  14. Zhao W, Cui Y, Liu L, Ma X, Qi X, Wang Y, Liu Z, Ma S, Liu J, Wu J. METTL3 Facilitates Oral Squamous Cell Carcinoma Tumorigenesis by Enhancing c-Myc Stability via YTHDF1-Mediated m(6)A Modification. Mol Therapy Nucleic Acids. 2020;20:1–12.

    Article  CAS  Google Scholar 

  15. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.

    Article  CAS  Google Scholar 

  16. Kanehisa M. Toward understanding the origin and evolution of cellular organisms. Protein Sci. 2019;28(11):1947–51.

    Article  CAS  Google Scholar 

  17. Kanehisa M, Furumichi M, Sato Y, Ishiguro-Watanabe M, Tanabe M. KEGG: integrating viruses and cellular organisms. Nucleic Acids Res. 2021;49(D1):D545-d551.

    Article  Google Scholar 

  18. Zhao W, Cui Y, Liu L, Qi X, Liu J, Ma S, Hu X, Zhang Z, Wang Y, Li H, et al. Splicing factor derived circular RNA circUHRF1 accelerates oral squamous cell carcinoma tumorigenesis via feedback loop. Cell Death Different. 2020;27(3):919–33.

    Article  CAS  Google Scholar 

  19. Jiang W, Zhang C, Zhang X, Sun L, Li J, Zuo J. CircRNA HIPK3 promotes the progression of oral squamous cell carcinoma through upregulation of the NUPR1/PI3K/AKT pathway by sponging miR-637. Ann Transl Med. 2021;9(10):860.

    Article  CAS  Google Scholar 

  20. Ai Y, Wu S, Zou C, Wei H: Circular RNA circFOXO3 regulates KDM2A by targeting miR-214 to promote tumor growth and metastasis in oral squamous cell carcinoma. J Cell Mol Med. 2022;26(6):1842–52.

  21. Liu J, Jiang X, Zou A, Mai Z, Huang Z, Sun L, Zhao J. circIGHG-Induced Epithelial-to-Mesenchymal Transition Promotes Oral Squamous Cell Carcinoma Progression via miR-142-5p/IGF2BP3 Signaling. Cancer Res. 2021;81(2):344–55.

    Article  CAS  Google Scholar 

  22. Legnini I, Di Timoteo G, Rossi F, Morlando M, Briganti F, Sthandier O, Fatica A, Santini T, Andronache A, Wade M, et al. Circ-ZNF609 Is a Circular RNA that Can Be Translated and Functions in Myogenesis. Mol Cell. 2017;66(1):22-37.e29.

    Article  CAS  Google Scholar 

  23. Wang X, Ma R, Zhang X, Cui L, Ding Y, Shi W, Guo C, Shi Y. Crosstalk between N6-methyladenosine modification and circular RNAs: current understanding and future directions. Mol Cancer. 2021;20(1):121.

    Article  Google Scholar 

  24. Chen RX, Chen X, Xia LP, Zhang JX, Pan ZZ, Ma XD, Han K, Chen JW, Judde JG, Deas O, et al. N(6)-methyladenosine modification of circNSUN2 facilitates cytoplasmic export and stabilizes HMGA2 to promote colorectal liver metastasis. Nature Commun. 2019;10(1):4695.

    Article  Google Scholar 

  25. Di Timoteo G, Dattilo D, Centrón-Broco A, Colantoni A, Guarnacci M, Rossi F, Incarnato D, Oliviero S, Fatica A, Morlando M, et al. Modulation of circRNA Metabolism by m(6)A Modification. Cell Reports. 2020;31(6):107641.

    Article  Google Scholar 

  26. Chen C, Yuan W, Zhou Q, Shao B, Guo Y, Wang W, Yang S, Guo Y, Zhao L, Dang Q, et al. N6-methyladenosine-induced circ1662 promotes metastasis of colorectal cancer by accelerating YAP1 nuclear localization. Theranostics. 2021;11(9):4298–315.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This work was supported by National Natural Science Foundation of China (No. 82002889, 82104631, 81902740), Tianjin Medical University Stomatological Hospital Foundation (No: 2020YKY01), Hospital Project of Tianjin Medical University Cancer Institute and Hospital (No: B1908), Science & Technology Development Fund of the Tianjin Education Commission for Higher Education (No: 2018KJ079), Tianjin Science and Technology Planning Project (Diversified investment fund projects, general project) (No. 21JCYBJC01150).

Author information

Author notes

  1. Wei Zhao and Jingwen Liu are both first author.

Authors and Affiliations

  1. The School and Hospital of Stomatology, Tianjin Medical University, 300070, Tianjin, China

    Wei Zhao, Jingwen Liu, Jie Wu, Xiaozhou Ma, Xi Wang, Leyu Zhang, Xin Hu & Jiayin Deng

  2. Institute of Orthopaedics, Tianjin Hospital, Tianjin, China

    Zhe Han

  3. Department of Immunology, Key Laboratory of Immune Microenvironment and Disease of the Educational Ministry of China, Tianjin Key Laboratory of Cellular and Molecular Immunology, School of Basic Medical Sciences, Tianjin Medical University, 300060, Tianjin, China

    Jianming Yang

  4. Department of Integrated Traditional & Western Medicine, Tianjin Medical University Cancer Institute and Hospital, National Clinical Research Center for Cancer, 300060, Tianjin, China

    Yameng Cui

Authors
  1. Wei Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Jingwen Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Jie Wu
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Xiaozhou Ma
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Xi Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Leyu Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Zhe Han
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Jianming Yang
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Yameng Cui
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Xin Hu
    View author publications

    You can also search for this author in PubMed Google Scholar

  11. Jiayin Deng
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Wei Zhao, Jingwen Liu performed the experiments and write the paper. Jie Wu, Xi Wang, Leyu Zhang, Xiaozhou Ma, Zhe Han, Jianming Yang acted as the assists. Jiayin Deng, Xin Hu and Cui Yameng are responsible for the designing and funding. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Yameng Cui, Xin Hu or Jiayin Deng.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

NA.

Competing interests

Authors declare no conflict of interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhao, W., Liu, J., Wu, J. et al. High-throughput microarray reveals the epitranscriptome-wide landscape of m6A-modified circRNA in oral squamous cell carcinoma. BMC Genomics 23, 611 (2022). https://doi.org/10.1186/s12864-022-08806-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-022-08806-z

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us