Screening and functional prediction of differentially expressed circular RNAs in human glioma of different grades

Circular RNAs (circRNAs) have a critical regulatory function in human glioma. However, novel circRNAs related to different pathological grades of glioma and their crucial potential function are worth screening and prediction. CircRNA expression profiling was performed for 6 paired high- and low-grade glioma tissues and 5 adjacent normal brain tissues through next-generation sequencing. Quantitative real-time PCR (qRT-PCR) was conducted to validate circRNA expression. Bioinformatics analysis was performed, and circRNA-miRNA-mRNA networks were constructed. The expression and survival data of miRNAs and target genes were examined by GEPIA, Chinese Glioma Genome Atlas (CGGA), ONCOMINE, and cBioPortal databases. The RNA binding proteins (RBPs), open reading frames (ORFs) and N6-methyladenosine (m6A) modifications of the identified circRNAs were also predicted. Through multilevel research screening, 4 circRNAs (hsa_circ_0000915, hsa_circ_0127664, hsa_circ_0008362, and hsa_circ_0001467) were associated with glioma of different pathological grades and could be preferred candidates for subsequent functional analysis. Therefore, circRNAs are associated with the different pathological grades of glioma and reveal their potential critical regulatory function. CircRNAs might provide vital molecular biomarkers and potential therapeutic targets for glioma.

AGING because of its highly infiltrative growth characteristics. Low-grade gliomas (LGGs), such as diffuse astrocytic and oligodendroglial tumors, have low-grade malignancy and slow growth rates, resulting in a long survival time for patients. There are also several reports of progression of some LGG into HGG in the clinic [2]. However, which factors influence glioma grade is unclear, and the mechanisms underlying the progression, etiology, pathogenesis and molecular characteristics of the different grades of glioma require further study.
As a novel kind of endogenous noncoding RNA (ncRNA), circular RNAs (circRNAs) possess a covalently closed structure with the absence of a 5′ and 3′ end [3]. CircRNAs are specifically expressed in cells or tissues, evolutionarily conserved across species, and relatively stable as a result of resistance to RNase R [4]. CircRNAs have been associated with how malignant tumors are regulated, including gliomas, as indicated by the latest studies. CircRNAs have crucial functions in regulating the cell cycle, apoptosis, angiogenesis, proliferation, invasion, migration, metastasis and radioand chemoresistance in glioma cells [5][6][7]. Recent reports have claimed that certain circRNAs are closely related to the pathological grade of gliomas. For example, circFOXO3 [8], circNT5E [9], and hsa_circ_0001730 [10] are differentially expressed in HGG and LGG. CircRNA expression is disordered in HGG compared with matching normal tissues or cells [11]. Currently, as a detection method, next-generation sequencing is a crucial method to explore novel circRNAs. However, the differentially expressed circRNAs of different grades of human glioma have not been fully investigated. Novel circRNAs related to different pathological grades of glioma and their specific crucial functions are worth screening and prediction.
The biological function of circRNAs is mainly reflected in the regulation of gene expression at the transcriptional or posttranscriptional level. Regarding their regulatory mechanisms, circRNAs have been confirmed to act as miRNA sponges to regulate the expression of miRNA target genes by competitively binding miRNA response elements (MREs). Acting as a miRNA sponge or as a competing endogenous RNA (ceRNA) is still one of the most recognized mechanisms. CircRNAs could compete with mRNAs, lncRNAs, and pseudogenes that also contain MREs to adsorb miRNAs, relieve the negative regulation of miRNAs on target gene expression, and improve the function and expression levels of target genes. For instance, circCCDC9 could act as a miR-6792-3p sponge to relieve the negative regulation of miR-6792-3p on its target gene CAV1 and then inhibit the progression of gastric cancer [12]. CircMAPK4 could function as an oncogene and suppress apoptosis in glioma by sponging miR-125a-3p [13]. Therefore, to explore more circRNAs closely related to multiple cancers, including glioma, the mechanism of miRNA sponges is still a direction worthy of attention and research.
Apart from their function as miRNA sponges, circRNAs can also function as RBP sponges or protein scaffolds and can be translated into either peptides or proteins [14]. According to one report, an internal ribosome entry site (IRES)-driven open reading frame (ORF) is present in circ-SHPRH, which has been reported to be translated into SHPRH-146aa [15]. This finding suggests that the potential translation function may be an important research direction for circRNAs. N6methyladenine (m 6 A) modification is the most abundant form of posttranscriptional RNA modification in eukaryotes. In recent years, the role of m 6 A in cancers including gliomas has been gradually recognized [16,17]. It has been reported that N6-methyladenosine (m 6 A) modification could also initiate the translation of circRNAs [18]. In linear mRNAs, m 6 A is usually located around the start and stop codons, and it also appears to be enriched in the junction sequences of circRNAs [19]. Recently, it has been newly reported that m 6 A can also mediate the biogenesis of coding circRNAs [20], and circRNA can not only be modified by m 6 A but also regulate the process of m 6 A modification by binding to the enzyme modified by m 6 A [21]. However, these potential vital functions of circRNAs have not been fully investigated and clarified. Hence, the prediction and analysis of RBP sponges and m 6 A modification of circRNAs requires deep future research.
In this research, we innovatively performed RNA-seq analysis of gliomas of different grades and analyzed and predicted the multiple crucial functions of the identified circRNAs. Our results demonstrated that 4 differentially expressed circRNAs are associated with the pathological characteristics of glioma, and their multilevel functions merit systematic analysis. The potential utility of circRNAs as vital molecular biomarkers or potential therapeutic targets in glioma has been indicated.

CircRNA expression profiles
CircRNA expression profiles were performed on 6 HGG and 6 LGG patient tissues as well as 5 adjacent normal tissues. After the raw data were quantile normalized, the circRNA expression profiles were AGING observed. A fold change (FC) > 2 and p < 0.05 were the cutoff criteria to identify differentially expressed circRNAs between each pair of groups. Eventually, a total of 153799 circRNAs were discovered in the highlow grade group. Among them, 296 exhibited upregulation and 332 showed downregulation. With regard to the high-normal grade group, 118649 circRNAs were found. Among them, 718 showed upregulation, and 269 exhibited downregulation. For the low-normal grade group, 128334 circRNAs were discovered, among which 945 showed upregulation and 252 exhibited downregulation ( Figure 1A). As demonstrated by the hierarchical clustering heatmap ( Figure 1B), a distinctive mode of circRNA expression was observed among the three groups. Based on how the different circRNA types were distributed, a large proportion of the circRNAs originate from exons. In comparison, other circRNAs stem from intronic and/or intergenic genomic regions ( Figure 1C). How the expression of circRNAs was differentiated to a significant extent between each group is presented in the form of volcano plots ( Figure 2A). The Circos plots show all the differentially expressed circRNAs compared with HGG, LGG tissue and adjacent normal brain tissue ( Figure 2B), which led to the finding that every single chromosome showed differences in the expression of circRNAs.

Validation of the circRNA-seq data via qRT-PCR
To verify the circRNA-seq analysis data, we focused on the downregulated circRNAs with certain circBase IDs. Five circRNAs were first selected as examples in further analysis from the most downregulated circRNAs. A design was developed for the divergent primers spanning the back-splice junction sites. The qRT-PCRamplified products of each candidate circRNA were identified by Sanger sequencing to validate the sequence of head-to-tail splice junctions. The primer sequences are presented in Supplementary Table 1.
Samples were obtained from 15 HGG tissues, 15 LGG tissues, and 15 adjacent normal tissues for qRT-PCR analysis. The analysis revealed that 4 circRNAs were significantly downregulated, showing the same trend as the RNA-seq results, except 1 circRNA (hsa_circ_0000199) was not significantly differentially expressed. In conclusion, the qRT-PCR results suggested that the RNA-seq was reliable and that most of the screened circRNAs merit further investigation. As shown in Figure 3A, the expression of hsa_circ_0000915 and hsa_circ_0127664 was downregulated in the HGG group compared with the LGG group, and the expression of hsa_circ_0001467 and hsa_circ_0008362 was downregulated in the HGG group compared with the normal group ( Figure 3B). Sanger sequencing was used to validate the head-to-tail splice junctions of the above 4 circRNAs ( Figure 3C). Through characterization of the identified circRNAs, we found that hsa_circ_0000915 (position: chr19: 18650181-18650530) is encoded from exons 3 and 4 of FKBP8. Its splicing length is 259 bp. Moreover, the characterization information for hsa_circ_0000915, hsa_circ_0127664, hsa_circ_0008362, and hsa_circ_ 0001467 were drawn via circPrimer software [22] ( Figure 3D).

Diagnostic accuracy analysis of candidate circRNAs
To explore the potential diagnostic value of the identified circRNAs, receiver operating characteristic (ROC) curves were analyzed. Relative to low grade, HGG patients showed a downregulation of hsa_circ_0000915 and hsa_circ_0127664 expression. To obtain the diagnostic value of these 2 circRNAs in HGG and LGG patients, a ROC curve analysis was conducted. The AUC was 0.8267 (95% CI: 0.6819-0.9714, p < 0.01) for hsa_circ_0000915 and 0.8044 (95% CI: 0.6457-0.9632, p < 0.001) for hsa_circ_0127664. Moreover, the diagnostic value of hsa_circ_0008362 and hsa_circ_0001467 was also evaluated by ROC curve analysis in HGG and normal tissue. The AUCs for hsa_circ_0008362 and hsa_circ_0001467 were 0.9467 (95% CI: 0.8728-1.000, p < 0.0001) and 0.9644 (95% CI: 0.9088-1.000, p < 0.0001), respectively. The potential diagnostic value of these circRNAs in glioma patients is displayed in Figure 3E.

GO enrichment and Reactome pathway analysis
By conducting Gene Ontology (GO) enrichment analysis in the biological process (BP), cellular component (CC), and molecular function (MF) categories, the host genes (linear counterparts) of the differentially regulated circRNAs in high-low grade were investigated, based on which their potential roles were deduced. The 10 most enriched BP, CC and MF GO terms of dysregulated circRNAs from the high-low grade glioma group are shown in Figure 4A. The regulation of neuron projection development, GTPase binding, and postsynaptic density were the most enriched GO entries in the BP, MF, and CC categories, respectively. Then, we performed Reactome pathway analysis, and 20 ranking predominant pathways are displayed in Figure 4B. According to the Reactome analysis, the top three enriched and meaningful pathways were signaling by Rho GTPases (Reactome, R-HSA-194315; p= 0.000312945), signaling by receptor tyrosine kinases (Reactome, R-HSA-9006934; p= 0.001761105) and transcriptional regulation by TP53 (Reactome, R-HSA-3700989; p= 0.001084954).

AGING
Additionally, GO enrichment and Reactome pathway analysis of the downstream target genes in high-low grade was conducted. Figure 5A presents the 10 most enriched BP, CC and MF GO terms for the target genes in high-low grade. PremiRNA processing, ion channel activity, and RISC complex were the most enriched GO entries in the BP, MF, and CC categories, respectively. In addition, the top 20 predominant pathways were (E) ROC curve analysis of four candidate circRNAs. hsa_circ_0000915 and hsa_circ_0127664 in HGG patients were analyzed against low-grade glioma patients via ROC curves, and hsa_circ_0008362 and hsa_circ_0001467 in HGG were analyzed against normal tissue using ROC curves. AGING displayed by Reactome pathway analysis ( Figure 5B). According to the Reactome analysis, the top three pathways were cellular responses to stress (Reactome, R-HSA-2262752; p= 0.015340627), transcriptional regulation by RUNX1 (Reactome, R-HSA-8878171; p=0.021079352) and transcriptional regulation by MECP2 (Reactome, R-HSA-8986944; p=0.000234161).

CircRNA-miRNA-mRNA network analysis
Emerging studies have revealed that circRNAs could be miRNA sponges, which are effective in regulating how mRNA is expressed [23,24]. In this paper, miRNAs that could potentially bind circRNAs were predicted using miRanda software (total score ≥ 150, total energy ≤ −20). A total of 129 miRNAs were found to be potential targets of the 4 identified circRNAs (Supplementary Table 2). Based on the bioinformatics analysis and combined with the research reports about miRNAs related to glioma in PubMed, we found that the most likely potential target miRNAs for hsa_circ_0000915 included hsa-miR-6765-3p and hsa-miR-330-5p. The most likely potential target miRNAs for hsa_circ_0127664 included hsa-miR-3945 and hsa-miR-99b-3p. hsa-miR-4656 and hsa-miR-217-5p were also predicted to be potential targets of hsa_circ_0008362. hsa-miR-6827-5p and hsa-miR-302d-5p were predicted to be potential targets of hsa_circ_0001467. The partially screened detailed sequence analysis of MREs for hsa_circ_0000915, hsa_circ_0127664, hsa_circ_0008362 and hsa_circ_ 0001467 is presented in Figure 6A.
The top target genes of miRNAs were predicted by TargetScan and miRDB (Criteria: Cumulative weighted context++ score ≤ −0.24 for TargetScan, Target Score ≥ 90 for miRDB). By analysis, we found that FAIM2, DLGAP2, ATP1B1, and RALYL could be the targets of hsa-miR-330-5p, hsa-miR-99b-3p, hsa-miR-217-5p and hsa-miR-302d-5p, respectively. We compared the expression levels of FAIM2, DLGAP2, ATP1B1 and RALYL using ONCOMINE databases and found that these 4 genes were significantly downregulated in brain and CNS cancer compared with the normal group ( Figure 7A). The expression data for these target genes in patients with GBM and LGG were examined by GEPIA. Compared with the normal group, FAIM2, DLGAP2, ATP1B1 and RALYL were downregulated in GBM and LGG patients ( Figure 7B, 7C).
The survival analysis was performed with the CGGA database. The group with low FAIM2, DLGAP2, ATP1B1 and RALYL expression showed a reduced length of survival (Supplementary Figure 1A). We analyzed the FAIM2, DLGAP2, ATP1B1 and RALYL alterations using the cBioPortal online tool for Glioblastoma Multiforme (TCGA, Firehose Legacy). Four genes were altered in 28 (21%) queried patients/samples (136 samples) (Supplementary Figure  1B). These results indicated that these genes were likely to be glioma-related tumor suppressor genes. The circRNA-miRNA-mRNA network was diagrammed with the assistance of Cytoscape software (Supplementary Figure 2) (Supplementary Table 3).

RBP binding and m 6 A modification function prediction
According to the data from the Cancer-Specific CircRNA Database (CSCD), RBPs binding with 4 identified circRNAs were drawn via Cytoscape software ( Figure 8A). Twenty-two vital RBPs, IGF2BP1, IGF2BP3, EIF4A3, hnRNPC, and AGO2, can bind to circRNAs and may play important biological functions in glioma. According to the data from CSCD, the structure of 4 circRNAs was determined. Diagrams containing information about RBPs and ORFs are presented in Figure 8B. In addition, an interaction map of the difference in gene expression of RBPs between GBM and LGG and cancer pathway activity (activation and inhibition) was analyzed by GSCALite (Supplementary Figure 3). For example, we found that IGF2BP2, which is closely related to GBM and LGG, can be closely related to the mesenchymal transition (EMT) pathway of cancer progression.
According to data from the circBank database, the m 6 A levels for the 4 circRNAs were obtained and are shown in Table 1. We then analyzed their m 6 A modification sites using the SRAMP predictor tool. Only high confidence m 6 A sites are listed ( Table 2). It was found that hsa_circ_0127664 and hsa_circ_0008362 had a relatively high m 6    . The prediction of coding proteins of identified circRNAs was analyzed by circRNADb, and hsa_circ_0000915 and hsa_circ_0001467 had IRES as well as ORFs, revealing their potential functions in coding proteins (Table 3).

DISCUSSION
Circular RNAs (circRNAs) have a critical regulatory function in human glioma [25]. As a detection method, next-generation sequencing is a crucial method to explore novel circRNAs. However, regarding circRNAseq data for glioma, most previous studies have focused on HGG, particularly glioblastoma (GBM, WHO IV), and thus, dysregulated circRNAs have only been analyzed between glioma and normal tissues [9]. The differentially expressed circRNAs of different grades of human glioma have not been fully investigated. Novel circRNAs related to different pathological grades of glioma and their crucial function are also worth   Subsequently, to verify the RNA-seq data, the circRNA expression level in tissues was evaluated by qRT-PCR.
To determine the potential clinical diagnostic value of circRNAs with differential expression, ROC curves were analyzed. In addition, to preliminarily explore the multilevel functions as well as potential mechanisms associated with malignant development of glioma, GO enrichment, Reactome pathway, and circRNA-miRNA-mRNA network analyses were performed. The expression and overall survival data for the target miRNAs and mRNAs were analyzed and examined using the GEPIA, CGGA, ONCOMINE, and cBioPortal databases. We screened 4 important circRNAs and predicted their potential functions. As a novel direction worth exploring for circRNAs, m 6 A modification levels, RBPs, and ORFs were also predicted for further analysis.
We found that hundreds of circRNAs, with upregulation and downregulation, were significantly differentially expressed (FC ≥ 2 or ≤ 0.5 and p < 0.05) in different grade glioma patients. As revealed by the RNA-seq analysis carried out in the high-normal and low-normal groups, the number of upregulated circRNAs exceeded that of circRNAs exhibiting downregulation. The same phenomenon was observed in a study published by Xia et al. [26]. However, in the high-low group, we found more downregulated circRNAs, which aroused our interest. In addition, a large proportion of the circRNAs identified in our study originated from exons. In comparison, other circRNAs stemmed from intronic and/or intergenic genomic regions, which was also observed in a study published by Zhu et al. [27].
AGING Despite many studies focusing on upregulated circRNAs in multiple cancer types, other studies have suggested that a majority of circRNAs with active expression play a significant role in normal physiological functions [28]. Relative to colorectal and ovarian cancer tissues, a high proportion of circRNAs are abundant in normal human tissues [29]. Many studies have revealed that certain circRNAs are downregulated in tumors compared with normal tissues and cells [30][31][32][33], which indicates that these circRNAs may function as suppressor genes and have a protective function in inhibiting cancer. Hence, we chose specific downregulated circRNAs for further research. Five downregulated circRNAs were first selected from the most downregulated circRNAs, and the RNA-seq reliability was evaluated with qRT-PCR. In 15 paired HGG and LGG tissues and 15 adjacent normal tissues, the expression of these 5 circRNAs was validated via RT-qPCR. Among the 5 circRNAs, four exhibited a clear conformance with RNA-seq, indicating that the RNA-seq data were reliable. Furthermore, we found that the selected representative circRNA, hsa_circ_0001467, was distinctly downregulated in glioma tissues, which is related to its diagnostic value for glioma patients and shows its potentially important clinical significance. Thus, downregulated circRNAs should not be neglected because they may serve as a novel tumor suppressive factor and potential target for new therapies in human cancer as well as other diseases [34][35][36].
The differential expression of circRNAs is reasonably related to their parental genes because they are mainly encoded from exons and/or introns of host genes [3,37]. To explore the roles of circRNAs through this potential mechanism, the biological functions and identified predominant pathways of their host gene as well as target genes were conducted through GO enrichment and Reactome pathway analyses. According to GO enrichment analysis, the deregulation of circRNAs was related to the regulation of neuron projection development, GTPase binding, and postsynaptic density. According to the Reactome analysis of host genes, the top three pathways were signaling by Rho GTPases, signaling by receptor tyrosine kinases and transcriptional regulation by TP53. However, the top three pathways were cellular responses to stress, transcriptional regulation by RUNX1, and transcriptional regulation by MECP2 in the Reactome pathway analyses of target genes. We noted that the Rho GTPase pathway has been shown to be associated with infiltrative growth and invasion of glioblastoma [38,39]. In addition, RUNX1 is considered to a key factor participating in the malignant biological behavior of glioma cells [40]. It is highly expressed and could contributes to the mesenchymal subtype of GBM in a TGFβ pathway-dependent manner [41]. Although many studies have analyzed the host genes of dysregulated circRNAs based on GO enrichment and pathway analyses [27,42], we support that GO and pathway analyses aimed at target genes can also better predict circRNA function.
Increasing evidence suggests that circRNAs are related to various diseases that humans can develop, particularly carcinomas, such as glioma [10], gastric cancer [43,44], hepatocellular carcinoma [45][46][47], pancreatic cancer [48], thyroid carcinoma [49] and breast cancer [50]. Thus, it is speculated that circRNAs might assist with the diagnosis and treatment of these human diseases. Our results also revealed that these circRNAs were suitable and sensitive diagnostic markers in identifying glioma. Additionally, many studies have indicated that circRNA is capable of acting as a miRNA sponge, which leads to interference with miRNAs or other varieties of molecules associated with the target genes to be inhibited [23]. Considering the significance of miRNAs in the pathogenesis of glioma, it was suspected that the identified circRNAs might cause malignant development of glioma through interactions with miRNAs. Therefore, miRNAs targeted by these differentially expressed circRNAs were predicted. For instance, hsa_circ_0000915 has the potential to bind to hsa-miR-330-5p. As demonstrated in a prior study, hsa-miR-330-5p is related to the stemness and development of glioma [11]. In addition, according to our analysis based on the CGGA database, hsa-miR-330-5p also shows important clinical utility. A study on hsa-miR-302d-5p and glioblastoma showed that hsa-miR-302d-5p could function as a tumor inhibitor inhibiting glioblastoma via targeting NF-κB as well as FGF2 [51]. In addition, miR-217 could target YWHAG and promote the viability, proliferation, migration, invasion and mitosis of glioblastoma cells both in vitro and in vivo [52]. Accordingly, the prediction results of circRNAs and miRNAs in this paper deserve further study.
RBPs can act as activators or inhibitors of circRNA formation and regulate the expression level of circRNAs [53]. Recently, binding to RBPs and the potential coding of circRNAs have attracted the attention of researchers; for example, circ-Amotl1 could bind to PDK1 as well as AKT1, facilitating the cardioprotective nuclear translocation of pAKT [54]. In our research, the potential RBPs of these 4 circRNAs were analyzed. Twenty-two vital RBPs were found, including IGF2BP1, IGF2BP3, EIF4A3, hnRNPC, and AGO2. For instance, it was recently found that E2F1 and EIF4A3 might promote the biogenesis of circSEPT9 and the occurrence and development of triple-negative breast cancer by acting on the miR-637/LIF axis [55]. The interaction of the difference in gene expression of RBPs between GBM and LGG and cancer pathway activity was also analyzed. We found that some RBP genes were closely related to GBM and LGG, and they were predicted to have a close correlation with multiple cancer pathways. In our research, IGF2BP2 was predicted to be closely related to the EMT pathway of cancer progression. A recent study also confirmed this conclusion. It has been reported that miR-138 may function as a tumor inhibitor by directly inhibiting IGF2BP2 and suppressing EMT during the progression of LGG [56]. Hence, RBPs of circRNAs in this paper merit further study.
A large proportion of circRNAs, which are derived from coding genes and have ORFs, show potential functions of coding proteins [20]. In one study, circRNAs showed many m 6 A consensus motifs, and a m 6 A site was enough to initiate circRNA translation in a cap-independent fashion involving the m 6 A reader YTHDF3 and the translation initiation factors eIF4G2 and eIF3A. circRNA translation efficiency can be affected by the level of m 6 A modification [18,57]. Interestingly, after searching the literature and GEPIA databases, it was found that eIF4G2 and YTHDF3 were upregulated in glioma compared with normal brain tissues [58] (GEPIA2 screening condition: log2foldchange cutoff=1, p-value cutoff=0.01). In another study, m 6 A modification was also discovered to be common in circRNAs [59]. Recently, m 6 A was reported to mediate the biogenesis of coding circRNAs [20] and could modulate cytoplasmic export and promote CRC liver metastasis [60]. The functional crosstalk and their function in tumor regulation between non-coding RNAs (miRNA, lncRNA, and circRNA) and m 6 A have been explored [16,61]. Hence, the analysis and validation of m 6 A in circRNAs requires further research. We analyzed the m 6 A modification levels of candidate circRNAs. For example, the m 6 A modification level of hsa_circ_0127664 is 0.6784, which was provided by circBank. The coding proteins prediction of hsa_circ_0001467 suggests that it might function as a coding protein. Thus, this modification in circRNAs also merits further analysis.

CONCLUSIONS
Through next-generation sequencing and bioinformatics analysis, the important dysregulated circRNAs were determined. Our findings indicated that 4 circRNAs (hsa_circ_0000915, hsa_circ_0127664, hsa_circ_ 0008362, and hsa_circ_0001467) are associated with the pathological characteristics of glioma and revealed critical regulatory functions at multiple levels. The results presented here might provide vital molecular biomarkers and potential therapeutic targets for glioma.

Clinical samples
HGG and LGG patient tissues and adjacent normal tissues were collected from the Second Affiliated Hospital of Hebei Medical University and Affiliated Hospital of Hebei University. None of the recruited patients received preoperative radiotherapy or chemotherapy before biopsy. After resection from glioma patients, the tissues were instantly immersed in liquid nitrogen and then stored at -80°C for later use. The Ethics Committee of the Second Affiliated Hospital of Hebei Medical University and Affiliated Hospital of Hebei University authorized this study. Written informed consent was obtained under a grant from every glioma patient. Histopathological examination was independently confirmed by two pathologists.

Total RNA isolation and quality control
In line with the instructions from the manufacturer, a HiPure Total RNA Mini Kit (Magen) was employed to extract total RNA from tissues. In addition, a Eukaryote Total RNA Nano 6000 assay was used to assess RNA quality and integrity on an Agilent Technologies 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). RNA integrity numbers (RINs) were considered integrity values for RNA measurements. The results that showed RINs≥7 were kept for further analysis. After isolation, RNA was stored at -80°C for subsequent experiments.

RNA-seq analysis
To remove rRNA, the total RNA samples were treated with a RiboErase kit (Human/Mouse/Rat) (Kapa Biosystems, Inc., Woburn, MA). Before the RNA-seq libraries were constructed, linear RNAs were eliminated for subsequent treatment using RNase R (Epicenter Technologies, Madison, WI, USA). Briefly, as specified by the manufacturer, an RNA Hyper Prep Kit (Kapa Biosystems, Inc., Woburn, MA) was employed to construct strand-specific RNA-seq libraries. Following fragmentation, the digested RNA samples were applied to synthesize first-and second-strand cDNAs with the assistance of random hexamer primers; dTTP was substituted by dUTP in the second-strand synthesis AGING reaction. Subsequent to cDNA synthesis, end repair was carried out for the double-stranded products, a single 'A' base was added, and the cDNA products were ligated to adaptors. With the aid of magnetic particles, cDNA fragments were obtained for subsequent treatment using uracil DNA glycosylase to achieve removal of the second-strand cDNA. The first-strand cDNA was purified and then PCR-amplified, while a 2100 Agilent Bioanalyzer (Agilent, Santa Clara, CA, USA) was applied for quality control of the libraries. RNA-seq was performed with an Illumina HiSeq X Ten (X10) platform (Illumina, San Diego, CA, USA) via standard protocols on a 150-bp paired-end run. The information of relevant patients participating in the RNA-seq is presented in Supplementary Table 4.

Quantitative real-time PCR (qRT-PCR)
In line with the manufacturer's instructions, TRIzol reagent (Life Technologies, CA, USA) was employed to separate total RNA from tissues. Electrophoresis on a denaturing agarose gel was employed to assess RNA integrity. In addition, cDNA was acquired with a High-Capacity cDNA Reverse Transcription kit (Geneseed® II First Strand cDNA Synthesis Kit, Guangzhou, China). For qRT-PCR on an ABI 7500 real-time PCR system (Applied Biosystems, Foster City, CA, USA), a SYBR Green kit (Geneseed® qPCR SYBR® Green Master Mix, Guangzhou, China) was employed. For collection of fluorescent signals, the PCR program was operated for 5 minutes at 95°C prior to 40 thermal cycles, each with 10 sec at 95°C and 32 sec at 60°C. Glyceraldehyde phosphate dehydrogenase (GAPDH) (Sangon Biotech, Shanghai, China) was the endogenous control used to normalize circRNA expression. After PCR, the primer specificity was assessed by melt curve analysis. Using the average of the GAPDH-normalized 2 -ΔΔCt values, the circRNA expression levels were comparatively quantified through calculation [62]. The circRNA sequences were obtained from circBase (http://www.circbase.org/) [63] using divergent primers spanning the back-splice junction sites rather than convergent primers. A 1% (w/v) agarose gel was applied to separate the amplified PCR products, which were subjected to confirmation by Sanger sequencing.

Statistical analysis
GraphPad Prism 8 (GraphPad Software, Inc., La Jolla, CA, USA) was used for statistical analysis. Comparisons of two variables in the RNA-seq data were performed using the Student's t-test. Fold change (FC) ≥ 2 or ≤ 0.5 and p < 0.05 were considered statistically significant in the circRNA-seq analysis. To adjust the pvalue in RNA-seq analysis, the false discovery rate was determined by calculations. For each circRNA, the expression level is indicated as the fold change determined by applying the 2 -ΔΔCt method for qRT-PCR analysis. All data are denoted as the mean ± S.D. of a minimum of three separate experiments. The Mann-Whitney test was used for data comparison. p-values < 0.05 indicate statistical significance.

Data Availability Statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
XC G, YH Z, and Q L were the major contributors to the writing and revision of the manuscript. W X, XM L collected and prepared the samples. XC G, YH Z, and Q L contributed to performing the experiments and data interpretation. L S and WT Y collected the related references and participated in discussions. SG S and H W made substantial contributions to the conception or design of the work. SG S and H W approved the final version of the manuscript. All authors read and approved the final manuscript.