Hsa-circ-0064636 regulation of the target gene VDAC1/UBE4A by hsa-miR-326/hsa-miR-503-5 in osteosarcoma

Background: Circular RNAs (circRNAs) are a subclass of noncoding RNAs that play a critical role in the regulation of gene expression in eukaryotic organisms. Recent studies have revealed the critical role of circRNAs in cancer progression. Yet, little is not well understood of hsa-circ-0064636 in osteosarcoma (OS). Methods: The differential expression of hsa-circ-0064636 in OS cell lines was veried by quantitative real-time PCR (qRT-PCR). Differentially expressed mRNAs and miRNAs were screened in OS mRNA and miRNA expression datasets. miRNAs that interacted with hsa-circ-0064636 were predicted by RNAhybrid, TargetScan and miRanda. and were further detected using RNAhybrid, TargetScan and miRanda. miRWalk, miRMap, and miRNAMap were used to perform target gene prediction on the intersected miRNAs to construct a circ-miRNA-mRNA interactor network. The target genes were then subjected to survival analysis using PROGgeneV2, which resulted in a circ-miRNA-mRNA interaction subnetwork with prognostic impact. Results: The qRT-PCR experiments successfully veried that hsa-circ-0064636 was signicantly overexpressed in the OS cell line. Hsa-mir-326(miR-326) and hsa-mir-503-5p(miR-503-5p) are target miRNAs of hsa-circ-0064636 in the target genes obtained from the miR-326 and miR-503-5p screens. UBE4A and VDAC1 had a signicant effect on prognosis. UBE4A is a target gene for miR-326, while VDAC1 is a target gene for miR-503-5p . Conclusion: hsa-circ-0064636 may be involved in OS development through acting as a sponge to inhibit miR-326 and miR-503-5p , thus regulating the expression of VDAC1 and UBE4A.


Introduction
Osteosarcoma (OS) is one of the most common primary malignant bone tumors, most commonly seen in children and adolescents, and is characterized by high aggressiveness, high metastatic potential, rapid progression, and chemoresistance 1 . In the past few decades, with extensive resection surgery and multidrug adjuvant chemotherapy, 5-year survival rates have improved to 70-80% 2 . Regrettably, the 5year survival rate for patients with chemical resistance is signi cantly reduced to less than 20% due to poor sensitivity to chemotherapy. OS patients may be better off with additional therapies, such as small molecule targeted agents, but these strategies are not currently making breakthroughs in clinical trials 3,4 . Therefore, there is an urgent need for further research in the search for new OS therapies, especially complex gene regulatory networks.
In recent decades, studies have successively identi ed new classes of noncoding RNAs (ncRNAs), for instance circular RNAs (circRNAs), microRNAs (miRNAs), and long ncRNAs (lncRNAs). They all have different regulatory functions, can effectively regulate essential protein effectors of cellular function and play an important role in the development of osteosarcoma 5,6 . Among them, circRNA is a highly stable and conserved special RNA with covalent closed-loop properties that is not easily degraded by endonucleases and is widely found in various tissues and organs 7 . A growing body of research shows that circRNAs are widely participation in disease development in many ways and have important roles in cancer 8,9 . However, the potential molecular mechanisms and regulatory relationships of circRNAs in osteosarcoma are still not de ned.
In the present study, we mainly focused on the differential expression of hsa-circ-0064636 in OS and further studied its regulatory mechanism by bioinformatics to clarify its expression and prognostic signi cance and to provide a new experimental basis and research ideas for future OS molecular marker studies. serum with the addition of G418 (0.3 mg/mL) as well as penicillin (100 U/mL) and streptomycin (100 U/mL). OS cell lines were all cultured with DMEM/F12 medium containing 10% fetal bovine serum with the addition of streptomycin (100 U/mL) and penicillin (100 U/mL). All cells were incubated in 5% CO 2 at 37°C in a constant temperature incubator. Fresh complete medium was replaced every 2 or 3 days for further incubation, and when the cell density reached 90% or more, cells were then passaged to the next generation at a ratio of 1:3.

Materials And
2.3. RNA Extraction and qRT-qPCR. Total cellular RNA was extracted by TRIzol, and cDNA was synthesized using a reverse transcription kit. qRT-qPCR was performed using SYBR green PCR mix and divergent primers from hsa-circ-0064636 with the following reaction conditions: predenaturation at 95°C for 10 min, denaturation at 95°C for 15 s, annealing at 60°C for 30 s, and 45 cycles. The glyceraldehyde phosphate dehydrogenase (GAPDH) and U6 expression was used as a control. The relative expression of the nal screened genes was calculated by a comparative method 2 -ΔCt . The primers required in this study are listed in Table 1. 2.4. CircRNA-miRNA-mRNA network construction. The miRNAs interacting with hsa-circ-0064636 were predicted via the circRNA target miRNA prediction tools RNAhybrid 10 (http://bibiserv.techfak.unibielefeld.de/rnahybrid/), TargetScan 11 http://www.targetscan.org/ and miRanda 12 (http://www.microrna.org/microrna/home.do), which intersected to identify miRNAs that were simultaneously targeted by the four prediction tools. For miRNA targeting, the gene prediction tools RNAhybrid, TargetScan, miRanda, miRWalk 13 http://zmf.umm.uni-heidelberg.de/apps/zmf/mirwalk2/index.html , miRMap 14 (http://mirnamap.mbc.nctu.edu.tw/), and miRNAMap 15 http://mirnamap.mbc.nctu.edu.tw/ were used to make target gene predictions on the intersected miRNAs. They were compared with the screened miRNAs to identify seven intersects from differentially expressed miRNAs. The nal selected circRNA-miRNA and miRNA-mRNA were used to construct a visual regulatory network for circRNA-miRNA-mRNA using Cytoscape 3.8.0.
2.5 Wayne diagramming. Wayne diagrams are plotted using the R package 'Venn'.
2.6 Differential expression analysis. The miRNA expression dataset GSE65071 (9 normal samples, 14 OS samples) and the mRNA expression dataset GSE1608816 (9 normal samples, 14 OS samples) were downloaded from the GENE EXPRESSION OMNIBUS (https://www.ncbi.nlm.nih.gov/geo/, GEO ) database and the Bioconductor affy was used. The Limma package was used to identify differentially expressed genes between normal and tumor samples. log2|foldchange(FC) | > 1 and adj.P-values < 0.05 were considered to be statistically signi cantly different, and the resulting differentially expressed genes were used for volcano mapping.
2.7. Survival analysis. PROGgeneV2 16 (http://genomics.jefferson.edu/proggene) is a tool that can be used along with publicly available data to investigate the prognostic signi cance of genes. The previously obtained common target genes were entered into a database, the dataset for osteosarcoma was selected, and the overall survival map was constructed based on the median expression of a given gene that could be classi ed into high and low expression groups (Kaplan-Meier, KM plots). PROGeneV2 was used for hypothesis testing using the SURVIVAL package for R. The data were analyzed using the log-rank test. Statistical analysis was performed using the log-rank test, and the threshold for a meaningful survival prognosis was a p-value < 0.05. Based on the survival analysis, a new circRNA-miRNA-mRNA visual regulatory subnetwork was obtained.
2.8 Statistical analysis. GraphPad Prism 8 was used to statistically analyze the experimental data, and the statistical test data were normally distributed. A two-tailed value of p<0.05 was considered statistically signi cant, and a t-test was used for comparisons between the two groups.

Result
3.1 Analysis of variance expression. In the miRNA expression pro le data, 114 miRNAs were upregulated and 117 miRNAs were downregulated in osteosarcoma, and the differential expression volcanoes of miRNAs were plotted ( Figure 1A). The differential analysis of mRNA expression spectrum data yielded 716 mRNAs for downregulated genes and 1924 mRNAs for downregulated genes, and the differential expression volcanoes of mRNAs were mapped ( Figure 1B).
3.2 Hsa_circ_006463 targeting miRNA prediction. In order explore the regulatory mechanism of hsa_circ_006463 in OS species, hsa_circ_006463 target miRNAs were predicted by four databases. TargetScan, miRanda and RNAhybrid predicted 498, 965 and 226 target miRNAs, respectively, while differential analysis yielded 88 target miRNAs. All four intersections were used to obtain common target miRNAs for miR-326 and miR-503-5p ( Figure 2).
3.3 Prediction and analysis of miRNA targeted genes. The mechanism of the regulatory action of hsa-miR-326 and miR-503-5p in OS was studied and the possible binding of target genes was further predicted. miR-326 target genes were predicted by miRWalk, miRanda, miRMap, miRNAMap, RNAhybrid, and TargetScan databases; respectively, 5071, 3364, 6616, 6716, 16373, and 4830 target genes were obtained. Intersection with the 1924 upregulated differential genes obtained from the analysis yielded 31 target genes, which were plotted on a Wayne diagram and labeled with the names of these 31 target genes ( Figure 3A). miR-503-5p target genes were analyzed by miRWalk, miRanda, miRMap, miRNAMap, RNAhybrid and TargetScan databases, which yielded 2034, 844, 6692, 1496, 12416 and 2196 target genes, respectively. These were intersected with the 1924 upregulated differential genes from the previous analysis to obtain 31 target genes, which were plotted on a Wayne plot and labeled with the six target gene names ( Figure 3B).
3.4 CircRNA-miRNA-mRNA network construction. Six databases were used for hsa-miR-326 and miR-503-5p to predict binding targets and differential gene heptad intersections to obtain 31 and 9 target gene relationships, respectively, to construct relationship networks (Figure 4). hsa-circ-0064636 may be  Table 2. Unfortunately, HNRNPA2B1, FOXO3, MYOF, C9orf3, and CUX1 could not be retrieved in the database. In the survival analysis of 31 target genes, only two differences were signi cant (<0.05), namely, UBE4A and VDAC1, and their survival curves were obtained, as shown in Figure 5. UBE4A is a target gene of miR-326 , while VDAC1 is a target gene of miR-503-5p.
3.6 CircRNA to miRNA to mRNA network construction with prognostic impact signi cance. Based on the original circRNA-miRNA-mRNA ceRNA interaction axes , the mRNA retained only the two target genes(UBE4A and VDAC1) with prognostic impact on survival (circRNA and miRNA were also retained) and removed the other target genes, thus constructing a new subnetwork with prognostic implications, shown in Figure 6

Discussion
OS is one of the common malignant tumor of bone, what's more the treatment of OS is still a major problem for human health. The current treatment for OS is mainly surgery and chemotherapy, but the prognosis of patients is still not satisfactory. Previous studies have shown that many genes are associated with the development and progression of OS, however the emphasis has been on proteincoding genes or miRNAs 17 . The molecular mechanisms of OS are still unknown, as well as in the recent past, some noncoding RNAs have been shown to play an important role in the development of OS 18,19 . For example, it has been found that the lncRNA(TMPO-AS1), by regulating the miR-199a-5p to WNT7B axis, acts as a ceRNA that promotes osteosarcoma tumorigenesis 20 . miRNA-1236-3p is signi cantly overexpressed in the OS, thereby inhibiting cell proliferation and inducing apoptosis by targeting KLF8. 21 .
In contrast to other non-coding RNAs, for instance lncRNAs as well as miRNAs, circRNAs are highly conserved and highly stable in mammalian cells. These properties make circRNAs potentially ideal biomarkers and potential therapeutic targets 22 . In recent years, a growing body of research reveals that the important role of circRNAs in different biological processes, especially in tumorigenesis, development, and metastasis, such as gastric cancer 23 , colorectal cancer 24 , lung cancer 25 , and cervical cancer 26 . A previous study reported that circTADA2A is differentially expressed between OS cell lines and corresponding noncancer cell lines 27 . Nevertheless, There are also many circRNAs whose role in the development of osteosarcoma is still largely unknown.In this study, we focused on the signi cantly upregulated expression of hsa-circ-0064636 in OS cell lines and normal tissue cell lines and further predicted the regulatory network between its targeted regulatory miRNAs and corresponding target genes.
In this study, miR-326 and miR-503-5p were predicted to be the target miRNAs of hsa-circ-0064636 by multiple databases of circRNA-miRNA interactions. miR-326 and hsa-miR-503-5 were also used in the GSE65071 dataset and were underexpressed in the osteosarcoma samples compared to normal samples in the difference-in-differences analysis to verify their expression in osteosarcoma. As mentioned in a previous study, hsa-circ-0064636 is signi cantly up-regulated in OS, regulates down-regulation of miR-326 and promotes cervical cancer progression through up-regulation of ELK1. MiR-326 overexpression transfection experiments veri ed that it can lead to inhibition of proliferation and invasion and induce apoptosis and autophagy in cervical cancer cells 28 . MiR-326 is also signi cantly underexpressed in prostate cancer 29 and correlated with prognosis. However, miR-503-5p has been reported to inhibit tumorigenesis, angiogenesis, and lymphangiogenesis in colon cancer through direct downregulation of VEGF-A 30, and in hepatocellular carcinoma studies, miR-503-5p has been reported to regulate epithelial to mesenchymal transformation, metastasis, and prognosis of hepatocellular carcinoma cells through inhibition of WEE1 31 . It follows that miR-326 and miR-503-5p play a key role in cancer suppression in many cancers, as well as their expression is low in cancer. Therefore, hsa-circ-0064636 may lead to the development of OS by suppressing the expression of miR-326 and miR-503-5p in OS. In the existing reports, no query has been made about miR-326 and miR-503-5p in OS.
This study also further used the miRNA prediction mRNA database with analysis of signi cantly differentially expressed genes in OS to screen the potential target genes of miR-326 and miR-503-5p to obtain a circRNA-miRNA-mRNA ceRNA Interaction Network. On the basis of this, we further screened the mRNAs with prognostic impact to obtain a circRNA-miRNA-mRNA regulatory subnetwork. ube4A was a potential direct target of miR-326, and VDAC1 was a potential target of miR-503-5p. We also found signi cant prognostic differences between UBE4A and VDAC1 in the survival analysis of OS patients. The results of survival analysis plots showed the prognosis of UBE4A and VDAC1 was signi cantly worse in the high expression group compared to the low expression group.
UBE4A is belonging to the U-box ubiquitin ligase class of enzymes. The encoded protein is involved in polyubiquitin chain assembly and plays a key role in chromosome condensation and polyubiquitination segregation via securin. Autoantibodies against the encoded protein Potential to become a marker for scleroderma and Crohn's disease 32,33 . UBE4A is signi cantly upregulated in comparison between ovarian plasmatic cystic carcinoma and adjacent normal tissues. The VDAC1 voltage-dependent anion channel protein is a major component of the outer mitochondrial membrane 34 . VDAC1 is expressed in all compartments including mitochondria, the plasma membrane and other cells. The protein regulates major metabolic and energetic functions of the cell, including Ca 2+ homeostasis, oxidative stress, and mitochondria-mediated apoptosis. Because overexpression of VDAC1 triggers cell death, It may be related to the destruction of nerve cells 35 . VDAC1 was identi ed as a tumor promoter in cervical cancer, as well as knock-down of VDAC1 in cervical cancer Cell lines increase the rate at which apoptosis occurs, which was partially reversed by overexpression of VDAC1 36 .
Importantly, we veri ed by qRT-PCR that hsa-circ-00063636, VDAC1 and UBE4A were highly expressed in osteosarcoma cell lines, and miR-326 and miR-503-5p were lowly expressed in osteosarcoma cell lines.The results of qRT-PCR were consistent with the results of our bioinformatics analysis. The disadvantage is that the combination of the regulatory network nally obtained still needs to be further veri ed in the subsequent experiments.

Conclusion
The target binding miRNAs of hsa-circ-0064636 were miR-326 and miR-503-5p. Further predicting the mRNAs to which the miRNAs might bind and taking the intersection results with differentially expressed mRNAs showed that miR-326 had 31 potential targets. miR-503-5p had six potential binding coding genes, and a circRNA-miRNA-mRNA regulatory network was derived according to this. Survival analysis results showed that UBE4A and VDAC1 had signi cant impact signi cance, resulting in a survivalsigni cant subregulatory network of hsa-circ-0064636-miR-326 /miR-503-5p-UBE4A/VDAC1. It was veri ed by qRT-PCR experiments that hsa-circ-0064636, UBE4A and VDAC1is signi cantly differentially overexpressed in OS cell lines while miR-326 and miR-503-5p were under-expressed in osteosarcoma cell lines. In conclusion, hsa-circ-0064636 may be involved in OS development by acting as a sponge, suppressing miR-326 and miR-503-5p to Facilitate up-regulation of VDAC1 and UBE4A. Tables   Table 1. qRT-qPCR primer sequences.
Gene   Schematic representation of hsa-circ-0064636 targeted miRNAs. TargetScan, miRanda, and RNAhybrid database predicted hsa-circ-0064636 target miRNAs and miRNA Wayne plots were derived from the differential analysis, with the two common miRNAs labeled on the right side of the gure.

Figure 3
MiRNA target gene predictions and differential gene Wayne plots. (A) Wayne plot of the intersection of miR-326 using six databases to predict binding target genes and differentially expressed genes. The 31 miR-326 target mRNA names are shown on the right side of (A). (B) Six databases were used for miR-503-5p to predict the intersection of binding target genes, and differentially expressed genes are shown in a Wayne diagram. The names of the six miR-503-5p target mRNAs are shown on the right side of (B).