Profiling and Bioinformatics Analyses of Differential Circular RNA Expression in Glioblastoma Multiforme Cells Under Hypoxia

The hypoxia microenvironment is highly associated with GBM’s malignant phenotypes. CircRNAs were reported involved in GBM’s biological characteristics and regulated by HIF-1α. However, the differential expression profile and role of circRNAs in GBM cells under hypoxia are still unclear. The expression profiles of circRNAs in LN229 and T98G under hypoxia were explored via circRNA sequencing analysis. Those circRNAs significantly dysregulated both in LN229 and T98G and could be found in circBase were selected and validated by qRT-PCR, RNase R digestion reaction, and Sanger sequencing. Normal cell line and fresh GBM tissues were also used for qRT-PCR validation. The roles of differentially expressed circRNAs were evaluated by bioinformatics analyses. There were 672 dysregulated circRNAs in LN229 and 698 dysregulated circRNAs in T98G. GO analysis indicated that the alteration of circRNA expression related to GBM cell’s biogenesis and metabolism. KEGG analysis demonstrated that TGF-β signaling pathway, HIF-1 signaling pathway, and metabolism-related signaling pathway were closely associated with differentially expressed circRNAs under hypoxia. These results were confirmed by GSEA analysis. The 6 selected and dysregulated circRNAs both in LN229 and T98G including hsa_circ_0000745, hsa_circ_0020093, hsa_circ_0020094, hsa_circ_0000943, hsa_circ_0004874, and hsa_circ_0002359 were validated by qRT-PCR. Inhibition of hsa_circ_0000745 inhibited GBM cell’s proliferation, migration, and invasion. HIF-1α centered circRNA-miRNA-mRNA networks analysis showed that the 6 validated circRNAs could cross-talk with 11 related miRNAs. The circRNA expressions are dysregulated in GBM cell under hypoxia. The 6 validated circRNAs could participate in GBM’s development and progression when hypoxia occurs. They might be the candidates for prognostic markers and adjuvant therapeutics of GBM in the future.


Introduction
Glioblastoma multiforme (GBM) is the most malignant brain tumor in the central nervous system (CNS) (Furnari et al. 2007;Ostrom et al. 2013). Currently, the most effective treatment of GBM is gross total surgical resection and standard postoperative concurrent chemoradiotherapy. However, this treatment still delivers an insufficient therapeutic effect for GBM (Bartek et al. 2012;Dunn et al. 2012). Hypoxia is a common microenvironment for solid tumors including GBM that can occur under both physiological and pathological conditions. The hypoxia microenvironment has been shown to be highly associated with tumorigenesis, angiogenesis, invasion, and resistance to chemoradiotherapy which are the main causes of GBM patients' death. Some studies reported that hypoxia is the driving force of GBM and could be used as a novel treatment target (Brahimi-Horn et al. 2007;Kalkan 2015;Yang et al. Zheng Chen, Shaohua Su and Min Yang contributed equally to this work. 1 3 2012; Mudassar et al. 2020). In addition, a new type of endogenous noncoding RNA, circular RNA (circRNA), is becoming an area of great interest in the field of RNA research recently. CircRNAs are covalently contiguous closed loops without 5′ and 3′ ends and are structurally more stable than linear RNAs and less susceptible to degradation by nucleic acid exonucleases. CircRNAs are a group of non-coding RNAs (ncRNAs) which can play crucial roles in regulating gene expression to affect the progression of various tumors. Accumulating studies have demonstrated that circRNAs can serve as a competing endogenous RNA (ceRNA) to target microRNAs (miRNAs), thus regulating the biological activities of tumor cells, and they can be regulated by hypoxia-inducible factor-1α (HIF-1α) in several kinds of tumor cells including glioma (Boeckel et al. 2015;Su et al. 2019;Ding et al. 2020;Jin et al. 2018;Wang et al. 2018). Therefore, circRNAs have great potential to be used as therapeutic targets of GBM under hypoxia.
To clarify the expression profile and role of circRNAs in GBM under hypoxia, we first used circRNA sequencing analysis to evaluate the circRNA expression profiles in two GBM cell lines under hypoxia and normoxia. A step-wise bioinformatics analysis was performed to predict functions and pathways in differentially expressed circRNAs. The same dysregulated circRNAs in two GBM cell lines were then validated and confirmed by quantitative real-time polymerase chain reaction (qRT-PCR).

Cell Culture
GBM cell lines LN229, T98G, and HEB (as normal control group) were purchased and authenticated by the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). The culture medium for LN229 and T98G was composed of DMEM (Life Technologies/Gibco, USA) and 10% FBS (Life Technologies/Gibco, USA), 100 U/mL penicillin, and 100 μg/mL streptomycin (Gibco, USA). For normoxic culture, the cells were incubated at 37 °C with 95% air, 5% CO 2 , and 100% humidity. For hypoxic culture, the cells were incubated at 37 °C with 1% O 2 , 5% CO 2 and 100% humidity for 72 h. Three pairs of LN229 cells under hypoxia (LNT) for 72 h and normoxia (LNC), 3 pairs of T98G cells under hypoxia (T98GT) for 72 h and normoxia (T98GC), and 3 pairs of HEB cells under hypoxia for 72 h and normoxia were prepared for further experiments.

circRNA Sequencing Analysis
The total RNA was isolated from both normoxic and hypoxic cultured GBM cell lines using TRIzol reagent (Invitrogen, USA) according to the manufacturer's instructions. The RNA purity and concentration were determined with Nan-oDrop ND-1000 (Thermo, USA). The RNA integrity of the samples was assessed by denaturing agarose gel electrophoresis. The Ribo-Zero rRNA Removal Kit (Illumina, USA) was used to remove the rRNA. The high-throughput whole transcriptome sequencing and subsequent bioinformatics analysis were performed by Cloud-Seq Biotech (Shanghai, China) as previously reported. Paired-end reads were acquired from the Illumina HiSeq 4000 sequencer. Edger software (v3.16.5) was used to identify the differentially expressed circRNAs. A differential expression of the cir-cRNA between 2 sets of samples is calculated by using a standardized number of reads. Any circRNA exhibiting the fold change (FC) ≥ 2.0, p value < 0.05, can be considered as significantly differential expression.

Functional Analysis of Different circRNAs
Gene Ontology (GO) functional analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were used to analyze differentially expressed circRNA-associated genes. In order to understand the potential functional mechanism of differentially expressed circRNA-associated genes better and in more details, the genes of two different cell lines were separately divided into upregulated group and downregulated group for GO analysis. GO functions were classified into 3 subcategories: Biological Process (BP), Cellular Component (CC) and Molecular Function (MF). The 10 most enriched GO terms were ranked by Enrichment Score. The KEGG pathway analysis was performed to explore the biological pathways which differentially expressed circRNAs might be involved in. Similarly, the associated genes of two different cell lines were also divided into upregulated group and downregulated group for KEGG analysis. Besides that, gene set enrichment analysis (GESA) was also used to confirm the GO terms and KEGG pathways in the GBM cells between under hypoxia and normoxia.

Analysis of circRNA-miRNA-mRNA Network and Related Prediction
Differentially expressed circRNAs were used to predict the potential miRNA response elements and the binding sites of miRNAs with CloudSeq's in-house software based on miRanda and TargetScan (CloudSeq Inc., China). The network between circRNAs and miRNAs was also constructed by Cytoscape based on the binding sites of the differentially expressed circRNAs and miRNAs. Different nodes represent circRNAs and miRNAs, respectively. Solid lines between two nodes represent potential binding. Because this study was related to the hypoxic microenvironment of GBM cells, the network was designed and presented HIF-1α centered.

qRT-PCR of circRNAs
Among all of the circRNAs identified, the same upregulated and downregulated circRNAs both in LN229 and T98G under hypoxia were firstly selected, and then the circRNAs which could not be found in circBase (http:// www. circb ase. org/) were excluded. Next, the remaining circRNAs which were verified significantly differentially expressed in circBase were finally selected for qRT-PCR validation. Glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was chosen to be a reference for normalization. The primers of all selected circRNAs are shown in Supplementary Table 1. Total RNA was reverse-transcribed into complementary DNA using a PrimeScript RT Reagent Kit (TaKaRa, Japan) and subjected to qRT-PCR on an Applied Biosystems 7500 Fast Real-Time PCR System. Three independent assays were performed for all samples. The expression was determined by using the threshold cycle, and relative expression levels were calculated via the 2 −△△Ct method. Likewise, circRNAs' expression in fresh GBM tissue (the central part of tumor) and fresh para-tumor tissue (peritumoral edema tissue containing GBM cells) were also verified by qRT-PCR (this study was approved by the Ethics Committee of Xinhua Hospital, School of Medicine, Shanghai Jiaotong University. All patients (or their guardians) who enrolled in the study had signed informed consents.). The expression of circRNAs in HEB cell line under hypoxia and normoxia were examined by qRT-PCR as well.

Ribonuclease R Digestion
RNase R digestion reaction was conducted as previously described. Total RNA (2 μg) was incubated at 37 °C for 30 min with 4 U/μg of RNase R (Epicentre Inc., USA) to enrich circRNAs, and after digestion, the total RNA was incubated at − 80 °C overnight to stop the reaction, and then, the RNA was subjected to cDNA synthesis and PCR. The amplification products were assessed by agarose gel electrophoresis.

Sanger Sequencing
To confirm the back-splice junction site of circRNAs, we designed divergent primers for performing PCR. The full length of PCR amplification products was determined by inserting them into a T-vector for the Sanger sequencing. The primer synthesis and Sanger sequencing were performed by the Biosune Company.

Silencing of hsa_circ_0000745
Three small interfering RNAs (siRNAs) and a negative siRNA control (si-NC) were synthesized (GenePharma, China

Cell Viability Assay
Cell viability and proliferation was, respectively, evaluated using the cell counting kit-8 solution (CCK-8) assay (Dojindo, Japan). Cells were collected and seeded at a density of 5 × 10 3 cells/well in 96-well plates for CCK-8 assay. After various treatments, cells were incubated with 10 μL CCK-8 solution for 1 h following the manufacturer's specifications. The data was assessed by measuring the optical density (OD) at 450 nm using a microplate reader (BioTek, USA).

Wound Healing Assay
Approximately 1 × 10 6 cells with different treatment were seeded into 6-well plates and incubated at 37 °C until cells reached a confluence of at least 90%. Wounds were created by scratching cell monolayers with a 200 μL plastic pipette tip and then incubated in fresh medium containing 1% FBS and different treatment for 24 h. Photographs were taken to estimate the mean number of migrating cells per field.

Transwell Invasion Assay
Cell invasion assays were performed using 24-well transwell plates (Corning, USA). Approximately 1 × 10 5 cells with different treatment were seeded in the upper chamber with serumfree medium in triplicate. Medium containing 10% FBS (300 μL) was added to the lower chamber as a chemo-attractant. After incubation for 24 h, the cells above the upper chambers were removed by cotton swab, and the cells below the membrane were fixed by methanol, stained with 0.1% crystal violet for 10 min, and counted from 5 randomly chosen fields for each well.

Statistical Analysis
Data are presented as means ± standard error (SE). SPSS 17.0 (SPSS Inc., USA) was used for the data analysis. The artworks were created by GraphPad Prism 5.0 (GraphPad Software, USA). The two-tailed t test was used for comparing the circRNA expression levels in GBM cell lines under hypoxia versus under normoxia. p value < 0.05 was considered statistically significant.

RNA-seq Profiling of circRNAs in GBM Cell Lines Under Hypoxia and Normoxia
A total of 12,539 circRNAs were detected in GBM cell lines. The proportion of upregulated and downregulated circRNAs in LN229 and T98G were showed in Fig. 1. The proportion of newly discovered circRNAs was also shown. There were a variety of catalogs of circRNAs, and most of them were exonic circRNAs. All of these circRNAs were widely distributed across all human chromosomes ( Fig. 1 and Supplementary Table 2). Differential expression of circRNA in LN229 or T98G under hypoxia could be visualized by a scatter plot. Furthermore, a volcano plot showed the differentially expressed circRNAs by the filtering criteria (FC ≥ 2.0, p value < 0.05) in LN229 or T98G under hypoxia. Among the 672 dysregulated circRNAs identified in LN229 under hypoxia, 349 circRNAs were upregulated and 323 circRNAs were downregulated (Supplementary Table 3). Among the 698 dysregulated cir-cRNAs identified in T98G under hypoxia, 338 circRNAs were upregulated, and 360 circRNAs were downregulated (Supplementary Table 4). Cluster heatmap clearly showed differentially dysregulated circRNAs in GBM cell lines between under hypoxia and normoxia (Fig. 2).

Prediction of circRNA Functions in GBM Cell Lines Under Hypoxia
The top 10 (ranked by Enrichment Score) BP, CC, and MF terms are shown in Fig. 3 and Supplementary Table 5. In LN229, the top 3 identified BP terms corresponding to the upregulated circRNAs were organelle organization, cellular macromolecule metabolic process, and cellular component organization. The top 3 identified CC terms were related to intracellular, nucleoplasm, and cytosol. The top 3 identified MF terms were catalytic, acting on a protein and transferase activity. KEGG analysis showed that there were 10 pathways related to the upregulated cir-cRNAs, including transforming growth factor-β (TGF-β) signaling pathway, adherens junction, cell cycle, ubiquitinmediated proteolysis, viral carcinogenesis, lysine degradation, HIF-1 signaling pathway, RNA degradation, thyroid hormone signaling pathway, and p53 signaling pathway (Fig. 4, Supplementary Fig. 1 and Supplementary Table 6). In T98G, the top 3 identified BP terms corresponding to the upregulated circRNAs were regulation of biological quality, cell cycle, and response to cholesterol. The top 3 identified CC terms were related to intracellular, cytoplasm, intracellular organelle, and organelle. The top 3 identified MF terms were protein binding, enzyme binding, and purine nucleotide binding ( Fig. 3 and Supplementary  Table 7). KEGG analysis showed that the top 3 pathways related to the upregulated circRNAs were lysine degradation, TGF-β signaling pathway, and morphine addiction ( Fig. 4 and Supplementary Table 8). The downregulated circRNAs were analyzed in the same way. In LN229, the top 3 identified BP terms corresponding to the downregulated circRNAs were organelle organization, cellular component organization, or biogenesis and cellular component organization. The top 3 identified CC terms were related to intracellular, organelle, and intracellular organelle. The top 3 identified MF terms were protein binding, binding, cytoskeletal protein binding, and regulatory RNA binding ( Fig. 3 and Supplementary Table 5). The top 3 identified KEGG pathways were ubiquitin-mediated proteolysis, glycosaminoglycan biosynthesis-heparan sulfate/heparin, and lysine degradation ( Fig. 4 and Supplementary Table 6). In T98G, the top 3 identified BP terms corresponding to the downregulated circRNAs were cellular component organization, cellular component organization, or biogenesis and regulation of cellular metabolic process. The top 3 identified CC terms were related to intracellular, intracellular organelle, and organelle. The top 3 identified MF terms were binding, GTPase regulator activity, and protein binding ( Fig. 3

Validation of the Differentially Expressed circRNA in GBM Cell Lines, Tissues, and Normal Cell Line Under Hypoxia
We selected the same dysregulated circRNAs both in LN229 and T98G under hypoxia based on circRNA sequencing analysis results (Supplementary Table 9) and circBase data screening for validation. Divergent primers for those circR-NAs were designed (Supplementary Table 1), and the appearance of one single peak in the melting curve confirmed the specificity of the amplified products for each circRNA. The expression levels of 6 selected circRNAs were determined by using qRT-PCR. The expression levels of hsa_circ_0000745 were significantly higher both in LN229 and T98G under hypoxia than under normoxia. Meanwhile, hsa_circ_0020093, hsa_circ_0020094, hsa_circ_0000943, hsa_circ_0004874, and hsa_circ_0002359 expression levels were significantly lower both in LN229 and T98G under hypoxia compared with normoxic control (Fig. 5). We then conducted Sanger Fig. 1 Expression profiles of circRNAs in GBM cell lines under hypoxia. A Proportion of dysregulated circRNAs and newly discovered circRNAs in LN229 under hypoxia. The proportion of upregulated circRNAs (348, 52%) was slightly higher than downregulated circRNAs (324, 48%). Most of the dysregulated circRNAs (516, 77%) can be found in circBase database. Other 156 dysregulated cir-cRNAs (23%) were newly discovered. B Proportion of dysregulated circRNAs and newly discovered circRNAs in T98G under hypoxia. The proportion of downregulated circRNAs (359, 52%) was slightly higher than upregulated circRNAs (334, 48%). Most of the dysregulated circRNAs (534, 77%) can be found in circBase database. Other 159 dysregulated circRNAs (23%) were newly discovered. C The distribution of circRNAs from different classifications based on the genomic origin in LN229 under hypoxia. Most of the dysregulated circRNAs were exonic circRNAs. In downregulated circRNAs, the proportion of exonic circRNAs (77%) was lower than in upregulated circRNAs (84%). The proportions of intronic and sense overlapping circRNAs in downregulated circRNAs (9% and 14%) were lower than in upregulated circRNAs (6% and 10%). D The distribution of cir-cRNAs from different classifications based on the genomic origin in T98G under hypoxia. Most of the dysregulated circRNAs were exonic circRNAs. The proportions of exonic (82% vs 81%), intronic (6% vs 6%), sense overlapping (12% vs 13%), intergenic (0% vs 0%), and antisense (0% vs 0%) circRNAs between in upregulated and downregulated circRNAs were almost the same. E The distribution of dysregulated circRNAs in different chromosomes in LN229 under hypoxia. The differences of circRNAs' distribution chromosomes including chr8 (12 vs 4), chr17 (26 vs 8), chr18 (5 vs 12), chrX (2 vs 13) and chrM (0 vs 8) between in upregulated and downregulated groups were more than doubled. The upregulated circRNAs were most distributed on chr3 (30). The downregulated circRNAs were most distributed on chr1 (34). F The distribution of dysregulated circRNAs in different chromosomes in T98G under hypoxia. The differences of circR-NAs' distribution chromosomes including chr2 (14 vs 29), chr13 (2 vs 15), chr16 (13 vs 6), chr18 (7 vs 3), and chrM (12 vs 0) between in upregulated and downregulated groups were more than doubled. Both of the upregulated and downregulated circRNAs were most distributed on chr1 (33 and 46). In sum, there were significant differences in the chromosome distributions of circRNAs between the upregulated group and the downregulated group both in LN229 and T98G, and there was even a wide gap in some distribution positions, indicating that these differences may be specific and worthy of further study Compared with the upregulated group, CC in downregulated group had more cytosol and cytoplasm components, and MF in downregulated group had more RNA binding-related function. C GO analysis results of upregulated circRNAs' host genes in T98G under hypoxia. D GO analysis results of downregulated circRNAs' host genes in T98G under hypoxia. Compared with the upregulated group, BP in downregulated group had more cellular metabolic process and less sterolrelated biological process. The MF in downregulated group had more GTPase-related activity, ubiquitin-related activity, and less nucleotide binding-related function than in upregulated group sequencing to confirm the head-to-tail splicing junction in circRNAs PCR product amplified with divergent primers. In addition, we designed divergent primers for the 6 selected circRNAs and convergent primers for their corresponding mRNA to amplify them. The PCR amplification products of circRNAs were resistant to RNase R digestion, whereas the divergent primers amplified no products in genomic DNA (gDNA) (Fig. 5). In contrast, the PCR products of corresponding mRNA amplified with convergent primers disappeared after digestion by RNase R. In addition, results showed that when under hypoxia, the expression of hsa_circ_0000745 was significantly higher in 2 GBM cell lines than in HEB cell line. The expression of hsa_circ_0020093, hsa_circ_0020094, hsa_circ_0000943, hsa_circ_0004874, and hsa_circ_0002359 were significantly lower in GBM cell lines than in HEB cell line. The expression of hsa_circ_0000745 was significantly higher in GBM tissue than in para-tumor tissue. The expression of hsa_circ_0020093, hsa_circ_0020094, hsa_ circ_0000943, hsa_circ_0004874, and hsa_circ_0002359 were significantly lower in GBM tissue than in para-tumor tissue. There were not any significant expression changes between HEB cell line under hypoxia and normoxia (Fig. 6).

Construction of circRNA-miRNA-mRNA Interaction Networks
Since this study focused on hypoxic condition, we revolved around HIF-1α and utilized the 6 qPCR validated circRNAs to construct a predicted circRNA-miRNA-mRNA interaction network (Fig. 6). The target genes and potential miRNA targets of these 6 circRNAs were also listed in this study (Supplementary Table 10).

Inhibition of hsa_circ_0000745 Inhibits GBM Cell's Proliferation, Migration, and Invasion
In GBM cell lines, specifically designed siRNAs were used to knockdown hsa_circ_0000745 expression under hypoxia. Hsa_circ_0000745 expression was significantly reduced in cells transfected with siRNA (si-circ_0000745 #2 ). This construct was used in all subsequent analyses (Fig. 7A). Compared with the negative-control treatment, hsa_circ_0000745 knockdown obviously suppressed the proliferation ability of LN229 and T98G under hypoxia, as examined using CCK-8 and clone formation assays (Fig. 7B, C). GBM cells showed Compared with downregulated group, the upregulated group was more closely related to TGF-β signaling pathway, viral carcinogenesis, HIF-1 signaling pathway, RNA degradation, and p53 signaling pathway. It was worth noting that HIF-1 signaling pathway is directly related to hypoxia. Our result suggested that these upregulated circRNAs in LN229 may play an important role in GBM's development under hypoxia. C Top 10 pathways shown by KEGG analysis in upregulated circRNAs' host genes in T98G under hypoxia. D Top 10 pathways shown by KEGG analysis in downregulated circRNAs' host genes in T98G under hypoxia. Compared with downregulated group, the upregulated group was more closely related to TGF-β signaling pathway. This result was the same as in LN229 which indicated that TGF-β signaling pathway was not only activated in GBM cells under hypoxia but also participated in the upregulation of circRNAs decreased invasion through transwell assay after downregulation of hsa_circ_0000745 under hypoxia (p value < 0.01, Fig. 7D). There was a similar trend detected in the migration experiments (p value < 0.01, Fig. 7E). Based on these results, knockdown of hsa_circ_0000745 noticeably inhibited the proliferative, invasive and migratory ability of GBM cells under hypoxia.

Discussion
Glioma is the most common primary intracranial malignant tumor, among which GBM with the highest malignancy is the main pathological type (Louis et al. 2016;Nabors et al. 2017;Reni et al. 2017). Currently, total surgical resection combined with radiotherapy and chemotherapy is the most standard treatment method for GBM. However, the median survival period of GBM is only 14 to 16 months. For this reason, it is urgent for us to find more effective adjuvant treatment strategies for GBM (Bartek et al. 2012;Dunn et al. 2012;Hu et al. 2018;Jiang et al. 2016). The hypoxia microenvironment of GBM has been shown to be closely related to malignant progression, therapy resistance, and prognosis. HIF-1α was a master transcriptional regulator of the adaptive response to hypoxia. It was reported that HIF-1α could drive glioma from low-grade astrocytoma to high-grade GBM, and higher HIF-1α expression in GBM has been correlated with worse prognosis (Chen et al. 2019a;Huang et al. 2016;Jensen 2006;Korkolopoulou, et al. 2004;Lo Dico et al. 2018;Zagzag et al. 2000). Recently, circRNAs have become more attention as new prognostic or therapeutic targets for tumors. CircRNAs were found to be differentially expressed in GBM (Cen et al. 2021;Rezaei et al. 2020;Shahzad, et al. 2021;Zhang et al. 2020a;Balandeh et al. 2021;Guo and Piao 2021;Lu et al. 2021;Luo et al. 2021). In addition, previous research showed that circRNA was regulated by HIF-1α under hypoxia condition. In glioma, circDENND2A/ miR-625-5p axis was reported being associated with HIF-1α (Su et al. 2019). Therefore, as hypoxic microenvironmentrelated circRNAs participate in diverse biological characteristics of tumor cells, circRNAs have great potential to be used as therapeutic targets for GBM under hypoxia.
GO analysis results indicated that most target genes of dysregulated circRNAs founded in this study functioned in GBM cells under hypoxia. BP and MF results of GO analysis showed that the alteration of circRNAs expression related to cell biogenesis and metabolism. This meant that change of circRNAs expression may influence GBM cells' metabolism, further affecting the biological characteristics of tumor cells under hypoxia. The hypoxia microenvironment of GBM was highly associated with tumor invasion, chemoradiotherapy, and prognosis. Some studies indicated that circRNAs played important roles in GBM cell's malignant phenotypes. However, few studies reported the changes of circRNAs expression in GBM cells under hypoxia. Our results demonstrated that under hypoxia condition, the cir-cRNA expression was changed, and this change may affect GBM's metabolism through regulating target genes in tumor cells. GSEA analysis showed that the GO terms including metabolism and biogenesis were positively correlated with GBM cells under hypoxic condition. This was basically consistent with the results of general GO analysis we mentioned above. Further, KEGG analysis results showed the most relevant signaling pathways of differentially expressed circR-NAs in GBM cells under hypoxia. Activation of the HIF-1 pathway is a common feature of glioma and may explain the intense vascular hyperplasia often seen in GBM (Krock et al. 2011). In LN229, HIF-1 signaling pathway was found closely related to upregulated circRNAs. This result proved again that hypoxic condition could regulate circRNAs' regulation. In addition, TGF-β signaling is highly active in high grade gliomas, and elevated TGF-β activity has been associated with poor clinical outcomes in this deadly disease. It was reported that activated TGF-β1/Smad2 signaling in response to hypoxia in GBM contributed to the upregulation of left-right determination factor (LEFTY) whose overexpression with Nodal could increase GBM cells' proliferation rates and the expression of epithelial-mesenchymal transition (EMT)/glioma stem cell (GSC)-related markers (Matsumoto et al. 2020). Researches showed that integrin inhibition may downregulate TGF-β pathway activation in GBM. Anti-TGF-β strategies may indirectly inhibit vascular endothelial growth factor (VEGF) pathway activation (Seystahl et al. 2015). It was worth noting that irradiation and also hypoxia-induced attraction of haematopoietic stem and progenitor cells (HPC) by glioma cells involved a TGFβ-dependent, HIF-1α-mediated release of CXC chemokine ligand 12 (CXCL12). Induction of transcriptional activity of HIF-1α by hypoxia or irradiation requires an intact TGF-β signaling cascade in glioma. This delineated a novel stress signaling cascade in glioma cells involving TGF-β, HIF-1α and CXCL12 (Tabatabai et al. 2006). In our study, TGF-β signaling pathway was found associated with upregulated circRNAs under hypoxia both in LN229 and T98G. This result indicated that TGF-β pathway was not only activated ◂ in GBM cells under hypoxia but also probably participated in the dysregulation of circRNAs. During hypoxic conditions, tumor cells upregulated crucial metabolic enzymes that help them adapt to the demand for nutrients and changes in their redox status. Adaption to the hypoxic environment was partly mediated by the activation and stabilization of HIF-1α (Eales et al. 2016). Coincidentally, GSEA analysis also showed that TGF-β signaling pathway was positively correlated with GBM cells under hypoxia. This result proved once again that TGF-β signaling pathway played an important role in the development of GBM cells under hypoxia. Other signaling pathways such as ubiquitin-mediated proteolysis, purine metabolism, pyrimidine metabolism, glycosaminoglycan biosynthesis, and lysine degradation were found correlated with dysregulated circRNAs in this study which meant that dysregulation of circRNAs under hypoxia could further influence GBM cells' metabolism and biochemical reaction.
Validation experiments confirmed the reliability of our circRNA sequencing analysis data. We finally selected and identified 6 circRNAs that were dysregulated both in LN229 and T98G under hypoxia. Their expression difference whatever in GBM cell lines, tissues, and normal cell line all indicated that they might have regulatory roles such A Results showed that when under hypoxia, the expression of hsa_ circ_0000745 was significantly higher in GBM cell lines than in HEB cell line. The expressions of hsa_circ_0020093, hsa_circ_0020094, hsa_circ_0000943, hsa_circ_0004874, and hsa_circ_0002359 were significantly lower in GBM cell lines than in HEB cell line. B Results showed that the expression of hsa_circ_0000745 was significantly higher in GBM tissue than in para-tumor tissue. The expres-sion of hsa_circ_0020093, hsa_circ_0020094, hsa_circ_0000943, hsa_circ_0004874, and hsa_circ_0002359 were significantly lower in GBM tissue than in para-tumor tissue. C Results showed that there were not any significant expression changes between HEB cell line under hypoxia and normoxia. D The predicted HIF-1α centered cir-cRNA-miRNA-mRNA interaction networks. The green diamonds are circRNAs and the red triangles are miRNA. (*p value < 0.05, **p value < 0.01, ***p value < 0.001, HEB-Treated: HEB under hypoxia, HEB-Control: HEB under normoxia, LNT: LN229 under hypoxia, T98GT: T98G under hypoxia.) as tumor promoting circRNAs and hypoxia stress response circRNAs in biological processes of GBM under hypoxia, which could be verified by further experiments. Recently, it was demonstrated that circRNAs could serve as new prognostic markers and therapeutic targets. Sponging miRNAs in circRNA-miRNA-mRNA regulatory axes is one of the well-recognized functions of circRNAs in glioma. Because the circRNA sequencing analysis in this study was under hypoxia, bioinformatics analysis here only provided several HIF-1α centered possible pathways for each network. Previous literature indicated that circARHGAP35 was associated with poor survival in cancer patients. CircARHGAP35 depletion remarkably reduced cell migration and invasion abilities, while an increase in cell motility was observed in the linear ARHGAP35 knockdown tumor cells (Li et al. 2021). It was found in this study that hsa_circ_0000943 whose target gene was ARHGAP35 significantly downregulated both in LN229 and T98G under hypoxia. This result indicated that hsa_circ_0000943 may played a role in GBM cell's malignant phenotype under hypoxic condition. Hsa_circ_0000745 (circ-SPECC1), a novel discovered circRNA, was reported to promote tumorigenesis of gastric Fig. 7 Silencing of hsa_circ_0000745 and its related functional experimental results under hypoxia. A Three small interfering RNA (siRNAs) and a negative siRNA control (si-NC) were transfected into GBM cells under hypoxia. The interference efficiency was measured using qRT-PCR after 48 h. The most effective siRNA (si-circ_0000745. #2 ) was selected for all subsequent experiments. B The results of CCK8 assay showed that silencing of hsa_circ_0000745 decreased the proliferation of LN229 and T98G under hypoxia.
C The results of clone formation assay showed that silencing of hsa_circ_0000745 decreased the proliferation of LN229 and T98G under hypoxia. D The results of transwell assay showed that silencing of hsa_circ_0000745 decreased the invasion of LN229 and T98G under hypoxia. E The results of wound healing assay showed that silencing of hsa_circ_0000745 decreased the migration of LN229 and T98G under hypoxia. (*p value < 0.05, **p value < 0.01, ***p value < 0.001, ****p value < 0.0001) cancer cells via miR-526b/KDM4A/YAP1 axis. It modulated TGF-β2 and autophagy under oxidative stress by sponging miR-33a to promote hepatocellular carcinoma (Chen et al. 2019b;Zhang et al. 2020b). It was found in our study that the expression of circ-SPECC1 was upregulated significantly in two kinds of GBM cell lines under hypoxia. Knockdown of circ-SPECC1 noticeably inhibited the proliferative, invasive, and migratory ability of GBM cells under hypoxia. However, its definite mechanism needed to be further researched. Interestingly, the expression difference of circRNA under hypoxia in GBM cells has been studied in the past. Although the results may be different from our study due to the different cell lines (U87 and A172), this study also proved that the differentially expressed circRNA under hypoxia in GBM cells was closely related to HIF-1α. In addition, the study verified that dysregulated circDENND2A can affect the migration and invasion of tumor cells through functional experiments (Su et al. 2019).
This is a preliminary study, and there are still a few limitations. First, different GBM cell lines may respond differently to hypoxia condition. In this research, to avoid errors as far as possible, only dysregulated circRNAs expressed both in 2 cell lines were selected for study. Though the parallelism of the same set of repeat samples was already high for circRNA research. However, sample size here is still small, and there are many other GBM cell lines such as U373 and U87 that are not studied whose response to hypoxia condition remains unknown. Therefore, we still consider planning to increase sample size and other cell lines for difference analysis of circRNA in further study. More verification at basic research level and subsequent mechanism including ceRNA will also be studied in future. Second, 11 miRNA were showed directly related to HIF-1α in our network illustration. Their specific functions and mechanisms were still unclear in GBM's biological characteristics under hypoxia. They also deserve to be further studied.
In conclusion, our study reveals that the circRNA expressions are dysregulated in GBM cells under hypoxia. It means that circRNAs, especially the 6 validated circRNAs, could participate in GBM's development and progression when hypoxia occurs. They might be the candidates for prognostic markers and adjuvant therapeutics of GBM in the future. However, based on these findings, much work is still needed to uncover the underlying molecular mechanisms further.