Identification and characterization of N6‐methyladenosine modification of circRNAs in glioblastoma

Abstract This research systematically profiled the global N6‐methyladenosine modification pattern of circular RNAs (circRNAs) in glioblastoma (GBM). Based on RNA methylation sequencing (MeRIP sequencing or N6‐methyladenosine sequencing) and RNA sequencing, we described the N6‐methyladenosine modification status and gene expression of circRNAs in GBM and normal brain tissues. N6‐methyladenosine–related circRNAs were immunoprecipitated and validated by real‐time quantitative PCR. Bioinformatics analysis and related screening were carried out. Compared with those of the NC group, the circRNAs from GBM exhibited 1370 new N6‐methyladenosine peaks and 1322 missing N6‐methyladenosine peaks. Among the loci associated with altered N6‐methyladenosine peaks, 1298 were up‐regulated and 1905 were down‐regulated. The N6‐methyladenosine level tended to be positively correlated with circRNA expression. Bioinformatics analysis was used to predict the biological function of N6‐methyladenosine–modified circRNAs and the corresponding signalling pathways. In addition, through PCR validation combined with clinical data mining, we identified five molecules of interest (BUB1, C1S, DTHD1, F13A1 and NDC80) that could be initial candidates for further study of the function and mechanism of N6‐methyladenosine–mediated GBM development. In conclusion, our findings demonstrated the N6‐methyladenosine modification pattern of circRNAs in human GBM, revealing the possible roles of N6‐methyladenosine–mediated novel noncoding RNAs in the origin and progression of GBM.


| INTRODUC TI ON
Glioma is a common brain tumour that accounts for nearly 80% of all primary brain neoplasms. Among them, glioblastoma (GBM) is a life-threatening tumour with a worse survival outcome. Despite the use of multiple aggressive treatments, such as surgery and/or chemoradiotherapy, the survival rate of GBM patients is relatively low. [1][2][3] Therefore, to develop a better treatment strategy, it is important to understand the molecular features of GBM occurrence. Recent epigenetic studies have found that RNA posttranscriptional modification plays an essential role in regulating cell growth and metabolism as well as the biological behaviour of tumours. Over 60% of RNA modifications belong to N6-methyladenosine, which is also a prevalent epigenetic modification in eukaryotic mRNA. The modification process is reversible and is completed by 'writers', 'erasers' and 'readers'. [4][5][6] Interestingly, in addition to mRNAs, the well-known noncoding RNAs-circRNAs (a class of covalently linked single-stranded closed circRNAs without a 3′-end poly(A) tail or 5′-end cap)-also possess extensive N6-methyladenosine modification sites. 7 CircRNAs are closely related to the occurrence and development of glioma. 8,9 Moreover, a previous study suggested that different sets of principles can affect N6-methyladenosine biogenesis in circRNAs and mRNA because a number of N6-methyladenosine-circRNAs are produced from exons whose corresponding mRNAs do not contain N6-methyladenosine peaks. Despite this, there is a lack of studies on N6-methyladenosine modification of circRNAs, and the role of N6-methyladenosine modification of circRNAs in GBM pathogenesis remains unclarified. [10][11][12] In this study, we reported for the first time the circRNA-based analysis of N6-methyladenosine modification in GBM tissue and normal brain tissue, which proved that there was a high degree of difference and diversity in the N6-methyladenosine modification patterns between GBM and control groups. Meanwhile, abnormal N6-methyladenosine modification of circRNA in GBM has been shown to be involved in transcriptional regulation and cancer-associated pathways. We hope that this study will help further investigate the potential effects of N6methyladenosine modification on GBM pathogenesis.

| Patients & samples
Glioma tissues (confirmed as GBM by postoperative pathological diagnosis) and normal cerebral cortex tissues (traumatic brain contusion and laceration combined with cerebral hernia, requiring internal decompression) were collected intraoperatively. After the samples were isolated, they were rinsed with normal saline and immediately transferred to a 1.8-ml RNA-free cryopreserved tube, which was stored in a refrigerator at −80°C for RNA isolation. Five clinical samples from the GBM group and normal control group (NC group) were selected for N6-methyladenosine sequencing and RNA sequencing, and the remaining samples were stored for verification (Table S1). Ethical approval for this study was obtained from the Ethics Committee of the Affiliated Hospital of Hebei University, and written informed consent was received from all participants.

| Bioinformatics analysis
The paired-end sequences were acquired from an Illumina NovaSeq 6000 sequencer and subjected to quality control (Q30).
Subsequently, 3′ adapter trimming and low-quality sequence removal were performed by cutadapt 1.9.3 software. 13 The alignment of clean sequences was performed against a reference genome (hg19, UCSC) using STAR software. 14 The identification of circRNAs was carried out by DCC software based on the aligned sequences. 15 Data normalization and differential circRNA expression analysis were performed using EdgeR software (v3. 16.5), 16 and the high-quality sequences from all libraries were aligned against the reference genome using HiSat2 software (v2.0.4). 17 MACS software was used to screen potential methylated sites on RNAs (peaks), 18 whereas diffReps was of Hebei University of Chinese Medicine, Grant/Award Number: KTZ2019019 revealing the possible roles of N6-methyladenosine-mediated novel noncoding RNAs in the origin and progression of GBM.

K E Y W O R D S
CircRNA, decoration pattern, GBM, glioblastoma, N6-methyladenosine modification, ncRNA, noncoding RNA used to identify differentially methylated sites. 19 The peaks found by these two programmes as overlapping exons of mRNA and cir-cRNA were identified and extracted by homemade scripts. Gene ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were utilized to determine the source genes of differentially methylated circRNAs and differentially expressed circRNAs.

| Gene-specific N6-methyladenosine qPCR validation
Five genes with differentially methylated sites according to N6methyladenosine sequencing and RNA sequencing were subjected to reverse transcription (RT)-qPCR. A portion of the fragmented RNA was saved for use as the input control. The remaining RNA was bound to anti-N6-methyladenosine antibody-coupled beads, and the N6-methyladenosine-containing RNA was then immunoprecipitated and eluted from the beads. Both the N6-methyladenosine-IP sample and input control were subjected to RT-qPCR with genespecific primers. The primer sequences are listed in Table S2.

| Statistical analysis
The mean ±standard deviation (SD) was calculated from the data of 3 independent experiments. Statistical tests were conducted using SPSS 25.0 and GraphPad Prism 7.0 software. Paired Student's t tests were performed between GBM group and NC group samples.
One-way ANOVA was utilized to compare the differences among three or more groups. Data were considered statistically significant as follows: *P-value <.05, **P-value <.01, ***P-value <.001 and ****P-value <.0001. All experiments were repeated three times independently.

| Overall N6-methyladenosine modification pattern in mRNAs from both groups
Human GBM tissue and normal brain tissue (n = 5) were selected for RNA sequencing and transcriptome-wide N6-methyladenosine sequencing assays. A total of 29 161 N6-methyladenosine peaks were identified by MACS in the GBM group, representing transcripts of 11 637 genes. In the NC group, 28 733 N6-methyladenosine peaks were identified, which corresponded to 11 096 gene transcripts ( Figure 1A,B). Among them, 20 335 individual N6-methyladenosine peaks in 9732 N6-methyladenosine-modified genes were detected in the two groups. Notably, the GBM group had 8826 new peaks and 8398 missing peaks compared to the NC group, revealing that the global N6-methyladenosine modification patterns were markedly different between the GBM and NC groups ( Figure 1A-C).
The N6-methyladenosine methylomes were further mapped by HOMER software. The top consensus motif in the 37,559 identified N6-methyladenosine peaks was GGACU ( Figure 1D).
By analysing the N6-methyladenosine-modified peaks of each gene, we found that 57% of all modified genes (4384/7757) had a unique N6-methyladenosine modification peak. The majority of genes (6877/7757) had one to three N6-methyladenosinemodified sites ( Figure 1E). In particular, genes with GBM-unique N6-methyladenosine tended to have more N6-methyladenosinemodified sites than genes with NC-unique N6-methyladenosine (genes with two N6-methyladenosine-modified sites: 14.00% vs 13.00%; genes with three or more N6-methyladenosine-modified sites: 6.50% vs 6.00%; Figure 1F 260 from the 3′UTR. In general, N6-methyladenosine peaks tend to occur in CDS regions, which means that N6-methyladenosine modification is likely to play a crucial role in regulating encoded proteins, but the mechanism needs further study.

| Introduction to basic retouching patterns
Genome-wide maps of N6-methyladenosine-modified circRNAs in the GBM and NC groups (n = 5 per group) were constructed. A total of 2,997 N6-methyladenosine-circRNA peaks overlapped in the control and GBM groups, whereas 1,322 N6-methyladenosine-circRNA peaks were only found in the control group but not in the GBM group, and 1,370 N6-methyladenosine-circRNA peaks were only detected in the GBM group but not in the control group ( Figure 2A).
Based on the motif analysis of 2000 circRNA peaks with the top enrichment score (−log10, P), a consensus sequence (GGACU) was identified in both the control and GBM groups ( Figure 2B), suggesting that the data are reproducible. As demonstrated in Figure 2C,D, N6-methyladenosine circRNA expression was relatively lower in the GBM group than in the control group, which was also true for non-N6-methyladenosine circRNAs, indicating that the downregulation of circRNAs in the GBM group seems to be unrelated to the assembly of N6-methyladenosine. Most N6-methyladenosine-circRNAs and non-N6-methyladenosine-circRNAs were frequently composed of a single one or more exons ( Figure 2E).  Table 1. In addition, the N6-methyladenosine-circRNA sites were significantly differentially expressed between the GBM and control groups (fold-change≥2.0, P≤0.00001; Figure 3A). The comparative analysis of differentially expressed N6-methyladenosine-circRNAs revealed that the most significant N6-methyladenosine peaks were frequently composed of exonic sequences ( Figure 3B). It has been reported that the majority of circRNAs are derived from protein-coding genes (PCGs) that span 2-3 exons. In this work, we found that most differentially methylated circRNAs originated from PCGs spanning a single exon, and the length of one exon in N6-methyladenosine-circRNAs was greater than that in two or more exons of N6-methyladenosine-circRNAs ( Figure 3C). In addition, the distribution patterns of modified N6-methyladenosine peaks in GBMs revealed that the abnormal N6-methyladenosine peaks were ascribed to all chromosomes, but chromosomes 1, 2 and 6 were more prominently represented ( Figure 3D). Among these chromosomes, the 3 chromosomes harbouring the highest number of differentially methylated N6-methyladenosine peaks were chromosomes 1 (251), 2 (246) and 6 (235).

| Bioinformatic analyses of the altered N6methyladenosine circRNAs
To determine the pathophysiological role of N6-methyladenosine modification in GBM, GO enrichment and KEGG pathway analyses were conducted on the modified N6-methyladenosine peaks. As demonstrated in Figure 4A, the peaks that were up-regulated in GBM were markedly associated with cellular component organiza- In contrast, the down-regulated peaks were obviously related to glutamatergic synapses and morphine addiction ( Figure 4D).  Figure S1D). In addition, the distribution patterns of circRNAs in GBMs demonstrated that the altered circRNAs were attributed to all chromosomes, but chromosomes 1, 2 and 3 were overrepresented ( Figure S1E).
In addition, the top 20 modified circRNAs are summarized in Table 2. GO and KEGG pathway analyses also highlighted that the  We wondered whether the number of N6-methyladenosine peaks in each gene was related to gene expression levels. As shown in Figure 1E, different genes had different numbers of N6-methyladenosine-modified sites. By analysing the relative expression levels of these genes, it was observed that a smaller number of N6-methyladenosine-modified regions in each gene was related to elevated gene expression. Genes with two and four N6methyladenosine-modified sites tended to have less circRNA abundance than those with unique N6-methyladenosine modification peaks ( Figure 5D).  (Tables S3-S6). It also indicates that for a circRNA that only has one N6-methyladenosine modification site, which is likely to be a potentially important molecule with significant changes in both methylation and expression levels. As gene expression is regulated by various factors, the impact of differential N6-methyladenosine modifications on gene expression is worth further investigation.

| Effects of N6-methyladenosine modification on circRNA expression in GBM
To determine whether N6-methyladenosine methylation could affect circRNA expression levels, we examined the differential However, we found that among the circRNAs with up-regulated N6-methyladenosine levels, more circRNAs exhibited increased expression levels than exhibited decreased expression levels (11.00% vs 6.00%) ( Figure 6A,B). Among the genes with increased circRNA expression in the GBM group, more circRNAs showed increased N6-methyladenosine levels than showed decreased N6-methyladenosine levels (11.00% vs 7.00%). More N6-methyladenosine circRNAs (18%) were detected among the up-regulated circRNAs than among the down-regulated circRNAs (17%) (Figure 6C,D). This result indicated that N6methyladenosine modification tends to exhibit a significant correlation with circRNA expression in GBM.

| Validation of molecules of interest (sites) based on MeRIP-PCR
In combination with the RNA sequencing and N6-methyladenosine sequencing data, we developed screening criteria: all circRNAs (including methylation sites) conforming to the criteria (fold-change >2, P < .001) were ranked from high to low in terms of the size of foldchange value, and then, the intersection of the two data was taken. In addition, the expression levels of the above circRNAs were detected in GBM and NC samples, and the results showed that the change trend of circRNA expression levels was similar to that of N6-methyladenosine methylation levels ( Figure 7A,B). In conclusion, GBM samples had unique N6-methyladenosine modification patterns that are distinct from those of normal tissues at both the transcriptome-wide and gene-specific scales. In addition, although there is limited evidence at the gene database level to find direct associations between circRNAs and GBM prognosis, we identified several miRNAs that competitively bind cir-cRNAs based on the molecular sponge function (ceRNA) of circRNAs and GBM disease association. And the data mining of GBM correlation was carried out. The specific screening process was as follows:

| Correlation analysis of the clinical prognosis of key molecules
There were about 2500 miRNAs in the miRbase, and about 600 miRNAs related to GBM in the miRCancer database. We selected 4 circRNAs with consistent sequencing and validation results. These 600 were used for circRNA-miRNA analysis, and the top5 of each circRNA was selected to display for the results.
We performed patient-based data analysis for each circRNAtargeted binding miRNA and presented survival analysis for the four miRNAs most closely associated with prognosis ( Figure S4A-E). The results showed that these miRNAs were significantly positively correlated with the prognosis of glioma patients (P < .05), and played a

| D ISCUSS I ON
Many studies have suggested that epigenetic modification plays a crucial role in the pathogenesis of GBM. 20,21 Recently, N6methyladenosine modification has attracted extensive attention, but the specific mechanism of this novel RNA modification in the occurrence and development of GBM has not been fully studied, especially in circRNA. 11 We used RNA N6-methyladenosine sequencing to explore the state of N6-methyladenosine-circRNA mod-  Of the five circRNAs of interest we selected, BUB1 has been considered a novel therapeutic target for glioblastoma and plays a key role in promoting tumour proliferation and radiation resistance in GBM in a forkhead box protein M1 (FOXM1)-dependent manner. 36 In addition, ZWINT is significantly positively correlated with the mitochondrial protein NDC80, the serine/threonine protein kinase PLK1 (PLK1), and complex spindle and mitochondrial-related subunit 1 (SKA1) and is associated with the regulation of mitosis and the cell cycle in GBM. 37 In primary glioblastoma, the low tumour copy number of the F13A1 gene fragment is associated with poor survival, 38 which is inconsistent with our sequencing results; this complex mechanism needs further study. These key molecules can be used as the preferred genes for future research on the function and mechanism of N6-methyladenosine-mediated GBM development.
This study also has some limitations.

| CON CLUS ION
Our findings provide the first map of human circRNA N6methyladenosine modification in GBM and identify differentially expressed circRNA transcripts associated with hypermethylated and hypomethylated modifications, revealing the potential association between abnormal N6-methyladenosine modifications and cancer-related gene expression and function, which is helpful to further study the mechanism of N6-methyladenosine-mediated gene expression regulation. We hope that it can provide a road map for discovering the mechanism of action of N6-methyladenosine modification in noncoding RNAs and the development of new therapeutic strategies for GBM or provide new ideas for further research by regulating N6-methyladenosine modification transcripts or N6methyladenosine-related genes.

ACK N OWLED G EM ENTS
We would like to give special thanks to NewCore BiodataStudio in Projects of Hebei University (hbu2020bs003).

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflicts of interest. Writing-review & editing (equal).

DATA AVA I L A B I L I T Y S TAT E M E N T
The data sets used and analysed during the present study are available from the corresponding author on reasonable request.