Construction of a circRNA-miRNA-mRNA network based on competitive endogenous RNA reveals the function of circRNAs in osteosarcoma

Background Osteosarcoma (OS) is a common primary malignant bone tumour. Growing evidence suggests that circular RNAs (circRNAs) are closely related to the development of tumours. However, the function of circRNAs in OS remains unknown. Here, we aimed to determine the regulatory mechanisms of circRNAs in OS. Methods The expression profiles of OS circRNA (GSE96964), microRNA (GSE65071) and mRNA (GSE33382) were downloaded from the Gene Expression Omnibus (GEO) database to identify differentially expressed circRNAs, miRNAs and mRNAs in OS. A ceRNA network was constructed based on circRNA-miRNA pairs and miRNA-mRNA pairs. MRNAs with significant prognostic differences were identified by the TARGET database in the network. Functional and pathway enrichment analyses were performed, and interactions between proteins were predicted using Cytoscape. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed to elucidate the possible functions of these differentially expressed circRNAs. Results A total of 15 downregulated circRNAs, 136 upregulated miRNAs and 52 downregulated mRNAs were identified in OS. Finally, a circRNA-miRNA-mRNA network was constructed in OS based on 14 circRNAs, 24 miRNAs, and 52 mRNAs. GO and KEGG pathway analyses suggested that the mRNAs in the network may be involved in the pathogenesis and progression of OS. Four mRNAs identified by the TARGET database were significantly associated with OS survival prognosis. A circRNA-miRNA-mRNA subnetwork was constructed based on these four mRNAs. Conclusion Our results provide a deeper understanding of the regulatory mechanisms by which circRNAs compete for endogenous RNAs in OS.


Background
Osteosarcoma (OS) is one of the most common bone malignancies [1]. It mainly occurs during adolescence, and its incidence decreases with age [2,3]. At present, the treatment of OS mainly involves surgery and adjuvant chemotherapy, but the prognosis of patients is still not optimistic [4,5]. Studies have found that the occurrence and development of OS is closely related to gene mutations and gene expression disorders [6][7][8]. Currently, gene-targeted therapy has been widely studied in various cancers and may be a great potential treatment strategy [9][10][11]. Therefore, studying the molecular mechanism of OS regarding its occurrence and development processes

Open Access
Cancer Cell International *Correspondence: qibaochuang@163.com 1 Department of Orthopaedics, First People's Hospital of Yibin, Yibin 644000, Sichuan, China Full list of author information is available at the end of the article Qiu et al. Cancer Cell Int (2020) 20:48 and finding effective molecular targets are crucial for the diagnosis and treatment of OS. Circular RNAs (circRNAs) are a new class of noncoding RNAs with covalently closed loop structures [12,13]. Due to the lack of a 5′ cap structure and a 3′ poly (A) tail, circRNAs are resistant to exonuclease and are more stable than linear RNAs [14]. A growing body of research indicates that circRNAs play an important regulatory role in the development of cancer in vivo. Zhou et al. found that by sponging miR-203, circRNA-0008717 inhibits the proliferation and invasion of OS cells [15]. Song et al. found that hsa_circ_0001564 acts as a miR-29c-3p sponge to cause tumourigenesis and may serve as a potential biomarker for OS [16]. In addition, a related study found that circRNA CDR1as/miR-7 signalling promotes tumour growth in OS, with potential therapeutic and diagnostic value [17]. These studies indicate that cir-cRNAs play an important role in the development and progression of OS.
In this study, we collected the expression profiles of circRNAs, miRNAs and mRNAs in OS from the Gene Expression Omnibus (GEO) database and the TARGET database. Differentially expressed circRNAs, miRNAs and mRNAs were identified. CircRNA-miRNA pairs and miRNA-mRNA pairs were screened, and a circRNA-miRNA-mRNA network was constructed. The functional pathways of the mRNAs in this network were analysed by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment. To better understand the pathogenesis of OS, TARGET data were used to screen mRNAs related to OS prognosis, and a circRNA-miRNA-mRNA subnetwork related to OS prognosis was constructed. Our research provides new insights into the mechanisms of the occurrence and development of OS. The flowchart for this procedure is shown in Fig. 1.

Data acquisition
GEO (http://www.ncbi.nlm.nih.gov/geo) is a public functional genomic database. Three datasets (GSE96964, GSE65071, and GSE33382; analysed on Agilent-069978 Arraystar Human CircRNA microarray V1, Exiqon human V3 microRNA PCR panel I + II, and Illumina human-6 v2.0 expression beadchip, respectively) were downloaded from GEO. The probes were converted to the corresponding gene symbols based on the annotation information on the platform. The GSE96964 dataset contained a human osteoblast cell line (hFOB1.19) and 7 human osteosarcoma cell lines (U2OS, MTX300, HOS, MG63,143B, ZOS, and ZOSM). The GSE65071 dataset contained 15 healthy people and 20 patients with osteosarcoma. The GSE33382 dataset contained 3 human osteoblasts and 83 osteosarcoma samples. The TARGET database (https ://ocg. cance r.gov/progr ams/targe t) is a database of childhood tumours that seeks to determine molecular changes in the development and progression of refractory childhood cancer using a comprehensive genomic approach. This database contains gene expression profiles and clinical information for patients with osteosarcomas. We obtained 101 osteosarcoma samples with complete clinical information from the TARGET database. Approval by the Ethics Committee was not necessary because all data were collected from publicly available databases (GEO and TARGET).

Identification of differential expression of circRNAs, miRNAs and mRNAs
The raw data were processed by background correction and normalization by using the affy package of R/Bioconductor. The Limma and edgeR packages were used to identify differentially expressed circRNAs, miRNAs and mRNAs between normal samples and tumour samples. Log [fold change (FC)] > 2 and adj. P value < 0.05 were considered to indicate a statistically significant difference.

Construction of the ceRNA network
Information about circRNAs is available on the circBase website (http://www.circb ase.org/). The cancer-specific circRNA database (CSCD) (https ://gb.whu.edu. cn/CSCD/) was used to predict target miRNAs. These target miRNAs were further screened by miRNAs differentially expressed in GSE65071. Then, the miRDB, TargetScan and miRTarBase databases were employed to predict specific miRNAs. Only mRNAs recognized by the 3 databases were considered to be candidate targets and were intersected with the identified differentially expressed mRNAs to identify the differentially expressed mRNAs targeted by the differentially expressed miRNAs. Finally, the circRNA-miRNA-mRNA regulatory network was constructed by screening the circRNA-miRNA and miRNA-mRNA groups. The data were visualized using Cytoscape 3.7.1.

GO and KEGG functional enrichment analysis
Database for Annotation, Visualization and Integrated Discovery (DAVID; version 6.7; http://david .ncifc rf.gov) is an online bioinformatics database that integrates biological data and analysis tools and provides gene annotation information and protein data. KEGG is a database resource for understanding advanced functions and biological systems from large-scale molecular data generated by high-throughput experimental techniques. GO is a major bioinformatics tool for annotating genes and analysing their biological processes (BPs), molecular functions (MFs) and cellular components (CCs).
To determine the function of differentially expressed mRNAs in the network, analysis was performed using DAVID. An adj. P-value < 0.05 was considered to indicate a statistically significant difference.

Survival analysis
To analyse the mRNAs selected from the ceRNA network, the OS mRNA expression profile data and clinical sample data in the TARGET database were downloaded. The overall survival analysis was performed using Kaplan-Meier curves. The log-rank test was used for statistical analysis. The threshold for survival prognosis significance was a P-value < 0.05. A subnetwork was constructed based on these mRNAs.

Identification of differentially expressed circRNAs, miRNAs and mRNAs
The basic information regarding the expression of circR-NAs, miRNAs and mRNAs in OS and control tissues in the three microarray datasets (GSE96964, GSE65071 and GSE33382) is listed in Table 1. In the circRNA expression profile data, a total of 15 circRNAs were downregulated in OS (Fig. 2, Additional file 1: Figure S1). The basic information for the 15 circRNAs is listed in Table 2. One of the circRNAs (hsa_circ_0000144) was not found in the CSCD database. From the CSCD database, we predicted that 503 miRNAs might target these 14 circRNAs (Additional file 2: Table S1). The basic structural pattern of the 14 circRNAs is shown in Fig. 3. In the miRNA expression profile data, a total of 136 miRNAs were upregulated in OS (Additional file 3: Table S2). Twenty-four intersecting miRNAs were obtained and are shown in Additional file 4: Table S3 (Fig. 4a). A total of 324 potential target genes (Additional file 5: Table S4) were predicted by the TargetScan, miRTarBase and miRDB databases. A total of 52 mRNAs were downregulated in OS in the mRNA expression profile data (Additional file 6: Table S5). There were 52 intersecting mRNAs (Fig. 4b), as shown in Additional file 7: Table S6.

Survival analysis and subnetwork construction of related mRNAs
We downloaded OS mRNA expression profile data (101 samples) and clinical sample data from the TARGET database. We performed survival analysis for the 52 mRNAs in the constructed circRNA-miRNA-mRNA network. Four mRNAs (ARID5B, ELL2, PHC2, and STAT4) were significantly associated with survival prognosis in patients with OS (Fig. 7). We constructed a circRNA-miRNA-mRNA subnetwork with these four mRNAs (Fig. 8).

Discussion
CircRNAs are a class of circular noncoding RNAs that have been observed in a variety of cancers [18][19][20], but their function is still unclear. Due to their special circular structure, circRNAs have high stability and great potential as new tumour markers [21]. Recent studies have found that circRNAs are involved in the development and progression of various cancers [22][23][24], such as gastric cancer, breast cancer, and colorectal cancer. Current evidence suggests that circRNAs can target miRNAs, often referred to as "miRNA sponges, " to reduce the levels of miRNAs and release their targeted inhibition of mRNAs, thereby regulating the expression of protein-coding genes [25][26][27]. These studies have found that circRNAs may play a key role in the development and progression of cancers.
OS is a common bone malignancy [1]. Current treatments for OS are mainly surgery and chemotherapy, but patient prognosis remains unsatisfactory [4]. Due to their stable structure, circRNAs may be a new strategy in the targeted therapy of OS [21]. A related study found that the circRNA LRP6 negatively regulates the levels of KLF2 and APC to promote the development of OS [28]. Dysregulated circRNA_100876 suppresses the proliferation of OS cancer cells by targeting microRNA-136 [29]. The upregulation of the circular RNA circ_0001721 predicts unfavourable prognosis in OS and facilitates cell cycle progression by sponging miR-569 and miR-599 [30]. Although these studies have found that circRNAs are dysregulated in OS and are involved in the development of cancer, the specific mechanism of circRNAs in OS remains unclear.
To investigate the role of circRNAs in OS, a circRNA-miRNA-mRNA regulatory network was constructed based on bioinformatics predictions and transcriptome data. GO enrichment analysis showed that genes in this network were mainly involved in the regulation of protein serine/threonine kinase activity, the transferase complex and ubiquitin protein ligase binding. Many studies have found that serine/threonine kinases play an important role in the occurrence and development of diseases. Insulin-like growth factor I (IGF-1) induces BNIP3 expression through a known AKT serine/ threonine kinase 1 (AKT1)-mediated inhibitory phosphorylation of glycogen synthase kinase 3β (GSK-3β), which promotes cell growth and tumourigenesis [31].
Cooperative recruitment of Arl4A and Pak1 to the plasma membrane contributes to sustained Pak1 activation for cell migration [32]. Furthermore, the ubiquitination of proteins has been found to be involved in the cancer process. Suppression of the ubiquitin pathway by small molecule binding to ubiquitin enhances the doxorubicin sensitivity of cancer cells [33]. Epigenetic silencing of ubiquitin-specific protease 4 by snail1 contributes to macrophage-dependent inflammation and therapeutic resistance in lung cancer [34]. KEGG pathway analysis showed that the network was enriched in the PI3K-Akt signalling pathway, cell cycle and FoxO signalling, which have been confirmed by a large number of studies to participate in cancer progression [35][36][37][38]. These results suggest that the circRNAs in this network may play a key role in the occurrence and development of osteosarcoma. This result is worthy of further investigation.
To further study which circRNAs affect the prognosis of OS patients, we analysed the correlation between the expression levels of 52 mRNAs in this network and overall survival in OS patients based on TARGET data. Our study found that four mRNAs (ARID5B, ELL2, PHC2, and STAT4) may affect the prognosis of patients with OS. We constructed a circRNA-miRNA-mRNA subnetwork. There were 8 circRNAs (hsa_ circ_0002052, hsa_circ_0088214, hsa_circ_0088227, hsa_circ_0008792, hsa_circ_0061265, hsa_circ_0002017, hsa_circ_00054442 and hsa_circ_0077765) and four miRNAs (miR-33a-3p, miR-205-3141, miR-576-3p) in the network. Studies have found that STAT4 is closely related to the prognosis of a variety of cancers, such as breast, liver, and stomach cancers [39][40][41]. By the targeting of ELL2 by microRNA-299 in glioblastoma multiforme, the upregulation of microRNA-299 is associated with a poor prognosis [42]. A previous study found that Phc2 controls the mobilization of bone marrow haematopoietic stem cells and progenitor cells by inhibiting Vcam1 expression [43]. ARID5B influences antimetabolite drug sensitivity and prognosis of acute lymphoblastic leukaemia [44]. Although we have identified four prognostic genes that have not been reported in OS, our predictions are consistent with studies in other cancers. Whether the dysregulation of these four genes will affect the prognosis of OS patients requires more indepth research in the future.
(See figure on next page.) Fig. 4 The differentially expressed miRNAs and mRNAs in the microarray datasets. a Twenty-four miRNAs were obtained from the intersection of miRNAs predicted from the website and the differentially expressed miRNAs identified in the miRNA dataset (GEO dataset). The expression levels of 24 miRNAs in the miRNA dataset are shown in a. b A total of 52 differentially expressed mRNAs were obtained by analysing the mRNA dataset and predicting the targeted binding of genes by the 24 miRNAs. The expression levels of these 52 mRNAs in the mRNA dataset are shown in b. (|log 2 (fold change (FC))| > 2, and an adjusted P value < 0.05 were the cut-off criteria). *P < 0.05, **P < 0.01, ***P < 0.001 Qiu et al. Cancer Cell Int (2020) 20:48 At present, there are few reports on the mechanism of circRNAs in osteosarcoma. The novelty of this study is not only that the circRNA-miRNA-mRNA network was constructed through the GEO database but is also related to the prognosis of OS through the TARGET database. Previous studies identified OS survival prognostic genes based only on the GEO database, which contains few clinical samples of osteosarcoma. The TARGET database contains more osteosarcoma samples with clinical information, which can provide us with richer data for identifying genes related to survival prognosis. However, this study is also limited because it is mainly based on the analysis of sequencing data and may require further experimental exploration in the future.

Conclusion
We identified differentially expressed circRNAs, miR-NAs and mRNAs from publicly available microarray data to construct a circRNA-related ceRNA network. In addition, four mRNAs related to prognosis were identified based on the TARGET database. With these four mRNAs, we constructed a circRNA-miRNA-mRNA subnetwork that may be associated with prognosis in OS. The findings of this study may provide new insights into the pathogenesis of OS and suggest potential therapeutic targets that warrant further investigation.