Background/Aims: Osteosarcoma (OS) is the most common primary malignant bone tumor in children and adolescents. However, the molecular mechanisms regulating osteosarcoma tumorigenesis and progression are still poorly understood. Circular RNAs (circRNAs) have been identified as microRNA sponges and are involved in many important biological processes. This study aims to investigate the global changes in the expression pattern of circRNAs in osteosarcoma and provide a comprehensive understanding of differentially expressed circRNAs. Methods: Microarray based circRNA expression was determined in osteosarcoma cell lines and compared with hFOB1.19, which was used as the normal control. We confirmed the microarray data by real time-qPCR in both osteosarcoma cell lines and tissues. The circRNA/microRNA/mRNA interaction network was predicted using bioinformatics. Gene Ontology analysis and 4 annotation tools for pathway analysis (KEGG, Biocarta, PANTHER and Reactome) were used to predict the functions of differentially expressed circRNAs. Results: We revealed a number of differentially expressed circRNAs and 12 of them were confirmed, which suggests a potential role of circRNAs in OS. Among these differentially expressed circRNAs, hsa_circRNA_103801 was up-regulated in both osteosarcoma cell lines and tissues, while hsa_circRNA_104980 was down-regulated. The most likely potential target miRNAs for hsa_circRNA_103801 include hsa-miR-370-3p, hsa-miR-338-3p and hsa-miR-877-3p, while the most potential target miRNAs of hsa_circRNA_104980 consist of hsa-miR-1298-3p and hsa-miR-660-3p. Functional analysis found that hsa_circRNA_103801 was involved in pathways in cancer, such as the HIF-1, VEGF and angiogenesis pathway, the Rap1 signaling pathway and the PI3K-Akt signaling pathway, while hsa_circRNA_104980 was related to some pathways such as the tight junction pathway. Conclusions: This study has identified the comprehensive expression profile of circRNAs in osteosarcoma for the first time. And the ceRNA network prediction and bioinformatics functional analysis could provide a comprehensive understanding of hsa_circRNA_103801 and hsa_circRNA_104980, which may be involved in the initiation and progression of osteosarcoma. The present study indicates that circRNAs may play important roles in osteosarcoma and thus serve as biomarkers of osteosarcoma diagnosis and treatment.

Osteosarcoma(OS) is the most common primary malignant bone tumor, which mainly affects children and adolescents [1]. The prognosis of patients with osteosarcoma can be improved with the combination of surgery and chemotherapy, and the 5-year survival rate has reached 60% to 70% in recent years [2]. However, a considerable number of patients respond poorly to chemotherapy and have a high risk of local relapse or lung metastasis even after normative chemotherapy and curative resection of the primary lesion [3]. And the precise molecular mechanisms involved in the tumorigenesis, chemotherapy resistance and lung metastasis of osteosarcoma are still poorly understood. Moreover, there are almost no effective diagnostic markers or therapeutic targets for osteosarcoma patients, which strongly limits improvements in prognosis. Therefore, one strategy to avoid chemo-resistance, reduce lung metastasis and improve clinical outcomes is to identify effective diagnostic biomarkers or therapeutic targets that play important roles in osteosarcoma.

Circular RNAs (circRNAs), are a class of non-coding RNAs (ncRNAs) that regulate transcriptional and posttranscriptional gene expression [4]. Unlike traditional linear RNAs, circRNAs are circular in shape, where the 3' and 5' ends are covalently linked and form a closed continuous loop containing multiple microRNA (miRNA) binding sites capable of sequestering miRNAs with their own miRNA response elements (MREs) [5, 6]. Intracellular circRNAs with competing endogenous RNAs (ceRNAs) activity may act as miRNA sponges that sequester miRNAs by binding miRNAs with MREs, which strongly suppress miRNA activity and result in increased levels of miRNA target genes [7, 8]. Therefore, circRNAs have been considered important biological regulators for understanding the molecular mechanisms of disease and identifying effective diagnostic biomarkers or therapeutic targets.

Recent studies have shown and emphasized the importance of circRNAs in regulating cancer-related signaling pathways [9, 10]. Additionally, circRNAs may be associated with tumor types and serve as risk factors for some types of cancer [11, 12]. To our knowledge, however, little is known about the expression and function of circRNAs in osteosarcoma. In this study, we report, for the first time, the comprehensive expression profile of circRNAs in osteosarcoma. We also performed a systemic bioinformatics analysis to identify circRNAs that are essential for the biological processes of osteosarcoma, which may provide potential targets for the development of novel diagnostic and therapeutic strategies against osteosarcoma.

Cell culture

The human osteoblast hFOB1.19 and the human osteosarcoma cell lines U2OS, MG63, HOS and 143B were obtained from the American Type Culture Collection (ATCC). hFOB1.19 was cultured by following the culture method described by ATCC. The ZOS and ZOS-M cell lines have been described previously [13]. U2OS/MTX300 cells, a methotrexate-resistant derivative of the U2OS human osteosarcoma cell line, were provided by Dr. M. Serra (Istituti Ortopedici Rizzoli, Bologna, Italy) and were continuously cultured in the presence of 300µg/L MTX [14]. All the other cells were cultured in DMEM (Gibco, Grand Island, NY, USA) and supplemented with 10% fetal bovine serum (Gibco, Grand Island, NY, USA), penicillin (10,000 U/L), and streptomycin (100 mg/L) at 37°C in a 5% CO2 humidified incubator.

RNA extraction and RNA Sample QC

Total RNA was extracted using the TRIZOL reagent (Invitrogen, NY, USA) in accordance with the manufacturer's protocol. And then total RNA from each sample was quantified using the NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA). RNA integrity and gDNA contamination of each sample was assessed by denaturing agarose gel electrophoresis.

Labeling and hybridization

The sample preparation and microarray hybridization were performed based on Arraystar’s standard protocols. Briefly, total RNA was digested with Rnase R (Epicentre, Inc.) to remove linear RNAs and enrich circular RNAs. Then, the enriched circular RNAs were amplified and transcribed into fluorescent cRNA utilizing a random priming method (Arraystar Super RNA Labeling Kit; Arraystar). The labeled cRNAs were hybridized onto the Arraystar Human circRNA Array (8xl5K, Arraystar). After having washed the slides, the arrays were scanned by the Agilent Scanner G2505C. The workflow of microarray expression profiles of circRNAs is shown in Fig. 1.

Fig. 1.

Experimental workflow of microarray expression profiles of circular RNAs in osteosarcoma.

Fig. 1.

Experimental workflow of microarray expression profiles of circular RNAs in osteosarcoma.

Close modal

Microarray data collection and analysis

Agilent Feature Extraction software (version 11.0.1.1; Agilent Technologies, Santa Clara, USA) was used to analyze the acquired array images. Quantile normalization and subsequent data processing were performed using the R software limma package. Differentially expressed circRNAs between two samples were identified through fold change filtering. Hierarchical clustering was performed to show the distinguishable circRNAs expression pattern among samples. Specifically, the R software limma package was used for normalization. And data transformation was performed by using the basic functions of R software. The package gplots and function heatmap2 in R software were used for mapping. And the distance metric was euclidean distance.

The full set of raw data from this study was deposited in NCBI’s Gene Expression Omnibus (GEO) and is accessible through the GEO Series accession number GSE96964.

Real-time quantitative PCR validation of microarray data

Real-time quantitative PCR was used to confirm the circRNA expression profiles obtained from the microarray data. Total RNA was extracted using the TRIZOL reagent (Invitrogen, NY, USA) in accordance with the manufacturer's protocol. And then total RNA from each sample was quantified using the NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA). Then cDNA was prepared by reverse transcription from the total RNA of the human osteoblast hFOB1.19 and five human osteosarcoma cell lines (U2OS, U2OS/MTX300, 143B, ZOS and ZOSM). The method for preparation cDNA for qPCR is described below.

Step one is to prepare the annealing mixture, combine the following:

incubate the tube for 5min under 65°C water bath temperature, then put the tube on the ice for 2min.

Step two is to prepare the RT reaction mixture, combine the following:

incubate the tube for 1min at 37°C followed by incubating at 50°C for 1h, then terminate at 70°C for 15min to inactivate the enzymes (Gene Amp PCR System 9700, Applied Biosystems).

In total, 6 up-regulated circRNAs and 6 down-regulated circRNAs were analyzed using the ABI PRISM® 7500 Sequence Detection System (Applied Biosystems). Specific primer sequences were presented in Table 1. The relative expression of circRNAs was determined by the 2-△△Ct method, and β-actin expression was used to normalize the data.

Table 1.

Primers used in this study

Primers used in this study
Primers used in this study

Detection of circRNA expression profiles in tissue samples

Fifteen pairs of fresh frozen tissue samples (15 tumor tissues and 15 paired normal tissues) obtained from patients with OS were used in this study. Total RNA was extracted by using the TRIZOL reagent (Invitrogen, NY, USA). And cDNA was prepared by reverse transcription of the total RNA from tissue samples. The relative expression of the circRNAs that were confirmed by cell line samples was determined using the ABI PRISM® 7500 Sequence Detection System (Applied Biosystems).

Annotation for circRNA/miRNA interaction

The circRNA/miRNA interaction was predicted with Arraystars homemade miRNA target prediction software based on TargetScan [15] and miRanda [16], and the differentially expressed circRNAs within all the comparisons were annotated in detail with the circRNA/miRNA interaction information. The competing endogenous RNA (ceRNA) network of circRNAs was characterized by the use of an approach that was termed mutually targeted MRE enrichment (MuTaME) [17]. And the analysis of the circRNA/microRNA/mRNA regulatory networks was conducted according to the target genes of circRNA-targeting miRNAs.

Bioinformatics analysis

Gene Ontology (GO) analysis was performed to explore the functional roles of circRNA-targeting genes in terms of biological processes, cellular components and molecular functions. Biological pathways defined by the Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/), Biocarta (https://cgap.nci.nih.gov/Pathways/BioCarta_Pathways), PANTHER (http://www.pantherdb.org) and Reactome (http://www reactome.org) were used to explore the pathways related to circRNA-targeting genes.

Statistics

Statistical analyses were performed using SPSS software, version 17.0 (SPSS Inc.). CircRNA expression profiles in 15 pairs of fresh frozen tissue samples were analyzed using Paired t test. Correlations between the relative expression levels of circRNAs and their ceRNAs were analyzed using Pearson’s correlation method. A P value of less than 0.05 was considered statistically significant.

Study approval

This study was approved by the medical ethical committee of the First Affiliated Hospital of Sun Yat-Sen University, and written informed consent was obtained from the patients or their guardians before sample collection.

Analysis of differentially expressed circRNAs

In total, the expression of 4660 human circRNAs was detected by the Arraystar Human circRNA Array. Hierarchical clustering was performed to show all expressed circRNAs (Fig. 2). Box plots show that the distributions of circRNAs for the eight compared samples were nearly the same after normalization (Fig. 3A). Based on the statistical analysis of circRNA expression data, differentially expressed circRNAs were discriminated between the human osteoblast hF0B1.19 (control) and the human osteosarcoma cell lines group. Differentially expressed circONAs were sorted by their fold change (FC ≥ 2.0), which are located above the top green line and below the bottom green line in the scatter plots shown in Fig. 3B. A total of 252 circRNAs were differentially expressed, consisting of 71 up-regulated circRNAs and 181 down-regulated circRNAs. And a heat map was drawn to show the top 20 up-regulated circRNAs and top 20 down-regulated circRNAs from all differentially expressed circRNAs (Fig. 3C). According to the annotation of human circRNAs using the RefSeq database and a catalogue of non-coding RNAs, genomic origin of human circRNAs includes exons, introns, UTRs, antisense to known transcripts, intergenic regions and unannotated regions of the genome [18]. We summarized the classification of differentially expressed circRNAs (Fig. 3D). Among the up-regulated circRNAs, there were 62 exonic, 8 intronic and 1 origin from UTRs. Among the down-regulated circRNAs, there were 172 exonic, 8 intronic and 1 antisense.

Fig. 2.

A hierarchical clustering analysis shows the distinguishable circRNAs expression pattern among samples. The R software limma package was used for normalization. And data transformation was performed by using the basic functions of R software. The package gplots and function heatmap2 in R software were used for mapping. And the distance metric was euclidean distance.

Fig. 2.

A hierarchical clustering analysis shows the distinguishable circRNAs expression pattern among samples. The R software limma package was used for normalization. And data transformation was performed by using the basic functions of R software. The package gplots and function heatmap2 in R software were used for mapping. And the distance metric was euclidean distance.

Close modal
Fig. 3.

Analysis of differentially expressed circRNAs in osteosarcoma. (A) Box plots show that the distributions of circRNAs for the eight compared samples were nearly the same after normalization. (B) CircRNAs in the Scatter-Plot located above the top green line and below the bottom green line indicated more than a 2.0-fold change of circRNAs. (C) The heat map shows the top 20 up-regulated circRNAs and top 20 down-regulated circRNAs. Red represents high relative expression, and green represents low relative expression. (D) The classification of differentially expressed circRNAs based on the genomic origin are shown.

Fig. 3.

Analysis of differentially expressed circRNAs in osteosarcoma. (A) Box plots show that the distributions of circRNAs for the eight compared samples were nearly the same after normalization. (B) CircRNAs in the Scatter-Plot located above the top green line and below the bottom green line indicated more than a 2.0-fold change of circRNAs. (C) The heat map shows the top 20 up-regulated circRNAs and top 20 down-regulated circRNAs. Red represents high relative expression, and green represents low relative expression. (D) The classification of differentially expressed circRNAs based on the genomic origin are shown.

Close modal

Real-time qPCR validation of 12 differentially expressed circRNAs

We selected 12 differentially expressed circRNAs, including 6 up-regulated circRNAs (hsa_circRNA_102049, hsa_circRNA_103309, hsa_circRNA_103801, hsa_circRNA_101113, hsa_circRNA_102485 and hsa_circRNA_100241) and 6 down-regulated circRNAs (hsa_circRNA_104692, hsa_circRNA_104892, hsa_circRNA_102678, hsa_circRNA_103110, hsa_circRNA_100258 and hsa_circRNA_104980), for validation in human osteoblast hFOB1.19 and five human osteosarcoma cell lines (U2OS, U2OS/MTX300, 143B, ZOS and ZOSM) (Fig. 4A-L). There was a great consistency between the real-time qPCR results and microarray analysis data. All the selected circRNAs were confirmed, which demonstrated the high reliability of the microarray expression data.

Fig. 4.

Real-time qPCR validation of 12 differentially expressed circRNAs and detection of circRNA expression profiles in 15 pairs of fresh frozen tissue samples. Fold changes were determined by the 2-△△Ct method. (A-F) The validation results of 6 up-regulated circRNAs are shown. (G-L) The validation results of 6 down-regulated circRNAs are shown. (M-N) The expression of hsa_circRNA_103801 was up-regulated in osteosarcoma tissue samples. (O-P) The expression of hsa_circRNA_104980 was down-regulated in osteosarcoma tissue samples. Data represent the mean ± SD of 3 independent experiments. *p<0.05, **p<0.01 by Paired t test, SPSS 17.0 (N and P).

Fig. 4.

Real-time qPCR validation of 12 differentially expressed circRNAs and detection of circRNA expression profiles in 15 pairs of fresh frozen tissue samples. Fold changes were determined by the 2-△△Ct method. (A-F) The validation results of 6 up-regulated circRNAs are shown. (G-L) The validation results of 6 down-regulated circRNAs are shown. (M-N) The expression of hsa_circRNA_103801 was up-regulated in osteosarcoma tissue samples. (O-P) The expression of hsa_circRNA_104980 was down-regulated in osteosarcoma tissue samples. Data represent the mean ± SD of 3 independent experiments. *p<0.05, **p<0.01 by Paired t test, SPSS 17.0 (N and P).

Close modal

Detection of circRNA expression profiles in tissue samples

To determine which circRNAs played the most important roles in osteosarcoma, realtime qPCR was used to detect the expression of the selected 12 differentially expressed circRNAs noted above in 15 pairs of fresh frozen tissue samples (15 tumor tissues and 15 paired normal tissues) obtained from patients with OS (Fig. 4, M and O). We noted that the expression profile of two circRNAs was consistent between cell lines and tissue samples, while the expression profile of other circRNAs was not quite consistent (data not shown). The expression level of hsa_circRNA_103801 was up-regulated in both osteosarcoma cell lines and osteosarcoma tissue samples (Fig. 4C and 4N). And hsa_circRNA_104980 was down-regulated in both osteosarcoma cell lines and osteosarcoma tissue samples (Fig. 4L and 4P).

Prediction and validation of the circRNA/microRNA/mRNA interaction network

Recent studies on circRNA/microRNA interaction have indicated that circRNAs play a key role in the regulation of gene expression by interacting with their target miRNAs [7, 18]. To find the potential target miRNAs that may be associated with osteosarcoma, two experimentally confirmed circRNAs (hsa_circRNA_103801 and hsa_circRNA_104980) were selected for further analysis and the circRNA/microRNA interaction was predicted with Arraystars homemade miRNA target prediction software based on TargetScan and miRanda.

The most likely potential target miRNAs for hsa_circRNA_103801 include hsa-miR-370-3p, hsa-miR-338-3p and hsa-miR-877-3p (Fig. 5A). The results of sequence analysis of microRNA response elements (MREs) are shown in Fig. 5B. The 2D structure shows the MRE sequence, the target miRNA seed type (8mer) and the 3 pairing sequence (nucleotides 13-16). The precise base positions are shown in the alignments in the upper left and right corners. The seed type 8mer represents an exact match to positions 2-8 of the target miRNA (the seed + position 8) followed by a base ‘A’. The AU content 30nt upstream and downstream of the seed sequence is shown in the "Local AU" column. The red bars stand for A/U and high accessibility, while the black bars stand for G/C and low accessibility of the seed. And the extent of the accessibility is shown by the height of the bars. The position column shows the most likely relative MRE position on the linear presentation of hsa_circRNA_103801.

Fig. 5.

Prediction of the circRNA/microRNA/mRNA interaction network of hsa_circRNA_103801. (A) The most likely potential target miRNAs for hsa_circRNA_103801 include hsa-miR-370-3p, hsa-miR-338-3p and hsa-miR-877-3p. (B) The results of the sequence analysis of microRNA response elements (MREs) are shown. (C) The top 10 potential hsa_circRNA_103801 ceRNAs are shown. (D) A network diagram shows the circRNA/microRNA/mRNA interaction network based on the predicted target genes of hsa_circRNA_103801-targeting miRNAs. The blue triangle in the center represents the target circRNA (hsa_circRNA_103801). The purple squares represent potential target miRNAs for hsa_circRNA_103801. The green triangles (circRNAs) and red rounds (mRNAs) represent the potential hsa_circRNA_103801 ceRNAs. (E) hsa_circRNA_103801 expression was correlated with the expression of TANC1, one of the putative targets of the putative target miRNAs in 10 human osteosarcoma tissues (Pearson’s correlation analysis, P = 0.0002).

Fig. 5.

Prediction of the circRNA/microRNA/mRNA interaction network of hsa_circRNA_103801. (A) The most likely potential target miRNAs for hsa_circRNA_103801 include hsa-miR-370-3p, hsa-miR-338-3p and hsa-miR-877-3p. (B) The results of the sequence analysis of microRNA response elements (MREs) are shown. (C) The top 10 potential hsa_circRNA_103801 ceRNAs are shown. (D) A network diagram shows the circRNA/microRNA/mRNA interaction network based on the predicted target genes of hsa_circRNA_103801-targeting miRNAs. The blue triangle in the center represents the target circRNA (hsa_circRNA_103801). The purple squares represent potential target miRNAs for hsa_circRNA_103801. The green triangles (circRNAs) and red rounds (mRNAs) represent the potential hsa_circRNA_103801 ceRNAs. (E) hsa_circRNA_103801 expression was correlated with the expression of TANC1, one of the putative targets of the putative target miRNAs in 10 human osteosarcoma tissues (Pearson’s correlation analysis, P = 0.0002).

Close modal

To identify and characterize the competing endogenous RNA (ceRNA) network of hsa_ circRNA_103801, we used an approach that was termed mutually targeted MRE enrichment (MuTaME)[17]. By using MuTaME, we sought to identify mRNAs that are targeted by hsa_ circRNA_103801-targeting microRNAs. The rna22 microRNA target prediction algorithm was used to generate MuTaME scores for the entire human protein-coding transcriptome [19]. And some authoritative studies have shown a low rate of false prediction for rna22 [20, 21]. Using MuTaME, we identified 3115 candidate ceRNAs for hsa_circRNA_103801. The top 10 potential hsa_circRNA_103801 ceRNAs are shown in Fig. 5C. A network diagram is drawn to show the circRNA/microRNA/mRNA interaction network based on the predicted target genes of hsa_circRNA_103801-targeting miRNAs (Fig. 5D). And qPCR validation of the predicted ceRNAs for hsa_circRNA_103801 in human osteosarcoma tissues was performed. The result indicated that hsa_circRNA_103801 expression was correlated with expression of some predicted ceRNAs, one of which is TANC1 (Fig. 5E). However, the expression of some predicted ceRNAs showed little correlation with hsa_circRNA_103801(data not shown), which may require further validation.

For hsa_circRNA_104980, the most likely potential target miRNAs are hsa-miR-1298-3p and hsa-miR-660-3p (Fig. 6A). The results of the sequence analysis of MREs are shown (Fig. 6B). The 2D structure shows the MRE sequence, the target miRNA seed type (7mer-m8) and the 3 pairing sequence (nucleotides 13-16). The precise base positions are shown in the alignments in the upper left and right corners. The seed type 7mer-m8 represents an exact match to positions 2-8 of the target miRNA (the seed + position 8). As mentioned above, the AU content 3 Ont upstream and downstream of the seed sequence is shown in the "Local AU" column. The position column shows the most likely relative MRE position on the linear presentation of hsa_circRNA_104980.

Fig. 6.

Prediction of the circRNA/microRNA/mRNA interaction network of hsa_circRNA_104980. (A) The most likely potential target miRNAs for hsa_circRNA_104980 include hsa-miR-1298-3p and hsa-miR-660-3p. (B) The results of the sequence analysis of microRNA response elements (MREs) are shown. (C) The top 10 potential hsa_circRNA_104980 ceRNAs are shown. (D) A network diagram shows the circRNA/microRNA/mRNA interaction network based on the predicted target genes of hsa_circRNA_104980-targeting miRNAs. The blue triangle in the center represents the target circRNA (hsa_circRNA_104980). The purple squares represent potential target miRNAs for hsa_circRNA_104980. The green triangles (circRNAs) and red rounds (mRNAs) represent the potential hsa_circRNA_104980 ceRNAs. (E) hsa_circRNA_104980 expression was correlated with the expression of NFATC2IP, one of the putative targets of the putative target miRNAs in 10 human osteosarcoma tissues (Pearson’s correlation analysis, P = 0.0787).

Fig. 6.

Prediction of the circRNA/microRNA/mRNA interaction network of hsa_circRNA_104980. (A) The most likely potential target miRNAs for hsa_circRNA_104980 include hsa-miR-1298-3p and hsa-miR-660-3p. (B) The results of the sequence analysis of microRNA response elements (MREs) are shown. (C) The top 10 potential hsa_circRNA_104980 ceRNAs are shown. (D) A network diagram shows the circRNA/microRNA/mRNA interaction network based on the predicted target genes of hsa_circRNA_104980-targeting miRNAs. The blue triangle in the center represents the target circRNA (hsa_circRNA_104980). The purple squares represent potential target miRNAs for hsa_circRNA_104980. The green triangles (circRNAs) and red rounds (mRNAs) represent the potential hsa_circRNA_104980 ceRNAs. (E) hsa_circRNA_104980 expression was correlated with the expression of NFATC2IP, one of the putative targets of the putative target miRNAs in 10 human osteosarcoma tissues (Pearson’s correlation analysis, P = 0.0787).

Close modal

We also identified mRNAs that are targeted by hsa_circRNA_104980-targeting microRNAs by using MuTaME. In total, we identified 1819 candidate ceRNAs of hsa_ circRNA_104980. And the top 10 potential hsa_circRNA_104980 ceRNAs are shown in Fig. 6C. A network diagram was drawn to show the circRNA/microRNA/mRNA interaction network based on the predicted target genes of hsa_circRNA_104980-targeting miRNAs (Fig. 6D). Further qPCR validation of the predicted ceRNAs in human osteosarcoma tissues showed that hsa_circRNA_104980 expression was correlated with the expression of some predicted ceRNAs, one of which is NFATC2IP. Although the P value was higher than 0.05, the result has shown the positive correlation between hsa_circRNA_104980 and NFATC2IP (Fig. 6E). It may require further validation in more tissue samples.

Bioinformatics analysis

In the present study, GO analysis and 4 annotation tools for pathway analysis (KEGG, Biocarta, PANTHER and Reactome) were used to investigate the biological functions of two experimentally confirmed circRNAs (hsa_circRNA_103801 andhsa_circRNA_104980), which were based on mRNAs targeted by circRNAs-targeting microRNAs.

As for hsa_circRNA_103801, the biological process analysis (Fig. 7A) revealed that its target genes were mainly involved in biological regulation, lipid phosphorylation and regulation of molecular function, among other biological processes. The cellular component analysis (Fig. 7B) indicated that the target genes were mainly involved in the plasma membrane and Golgi apparatus, among other cellular components. While the molecular function analysis (Fig. 7C) showed that hsa_circRNA_103801-related genes took part in phosphatidylinositol kinase activity, guanyl-nucleotide exchange factor activity and Ras guanyl-nucleotide exchange factor activity, among other activities.

Fig. 7.

GO analysis and KEGG pathway annotation of hsa_circRNA_103801. (A) Biological process analysis of hsa_circRNA_103801-targeting genes. (B) Cellular component analysis of hsa_circRNA_103801-targeting genes is shown. (C) Molecular function analysis of hsa_circRNA_103801-targeting genes is shown. (D) KEGG pathway analysis shows the top 10 pathways related to hsa_circRNA_103801. (E) The related pathways were ranked according to GeneRatio (Selection Counts/Selection Size). (F) Pathway analysis using Biocarta database shows the top 10 pathways related to hsa_circRNA_103801. (G) Pathway analysis using PANTHER database shows the top 10 pathways related to hsa_circRNA_103801. (H) Pathway analysis using Reactome database shows the top 10 pathways related to hsa_circRNA_103801.

Fig. 7.

GO analysis and KEGG pathway annotation of hsa_circRNA_103801. (A) Biological process analysis of hsa_circRNA_103801-targeting genes. (B) Cellular component analysis of hsa_circRNA_103801-targeting genes is shown. (C) Molecular function analysis of hsa_circRNA_103801-targeting genes is shown. (D) KEGG pathway analysis shows the top 10 pathways related to hsa_circRNA_103801. (E) The related pathways were ranked according to GeneRatio (Selection Counts/Selection Size). (F) Pathway analysis using Biocarta database shows the top 10 pathways related to hsa_circRNA_103801. (G) Pathway analysis using PANTHER database shows the top 10 pathways related to hsa_circRNA_103801. (H) Pathway analysis using Reactome database shows the top 10 pathways related to hsa_circRNA_103801.

Close modal

KEGG analysis (Fig. 7D) showed that there were 19 pathways related to hsa_circRNA_103801, including breast cancer, proteoglycans in cancer, pathways in cancer, and transcriptional misregulation in cancer. When the related pathways were ranked according to GeneRatio (Selection Counts/Selection Size), the top three related pathways were pathways in cancer, the PI3K-Akt signaling pathway and proteoglycans in cancer, respectively (Fig. 7E). Biocarta analysis (Fig. 7F) indicated that the most relevant pathway was VEGF, hypoxia and angiogenesis. And PANTHER analysis (Fig. 7G) showed that the top three related pathways were interleukin signaling pathway, FGF signaling pathway and hypoxia response via HIF activation. The Reactome analysis (Fig. 7H) also suggested that hsa_circRNA_103801 was involved in FGFR signaling pathway and PI3K-Akt signaling pathway.

For hsa_circRNA_104980, the biological process analysis (Fig. 8A) found that its target genes were mainly related to cellular protein modification processes, protein modification processes and macromolecule modification, among other similar processes. And the cellular component analysis (Fig. 8B) showed that the target genes were mainly involved in organelles, intracellular organelles and coated vesicles, among other similar cellular components. While the molecular function analysis (Fig. 8C) revealed that hsa_circRNA_104980-related genes took part in binding phosphoric ester hydrolase activity, and carbohydrate derivative binding, among other similar processes. KEGG analysis (Fig. 8D) found that there were 22 pathways related to hsa_circRNA_104980, and the top 10 pathways are shown. And the related pathways were ranked in order of GeneRatio (Fig. 8E). Biocarta analysis (Fig. 8F) indicated that the most relevant pathway was the information-processing pathway at the IFN-beta enhancer. And PANTHER analysis (Fig. 8G) showed that the related pathways of hsa_circRNA_104980 included Nicotinic acetylcholine receptor signaling pathway and p38 MAPK pathway. The Reactome analysis (Fig. 8H) also suggested that hsa_circRNA_104980 was related to asparagine N-linked glycosylation and cytokine signaling in immune system.

Fig. 8.

GO analysis and KEGG pathway annotation of hsa_circRNA_104980. (A) Biological process analysis of hsa_circRNA_104980-targeting genes. (B) Cellular component analysis of hsa_circRNA_104980-targeting genes is shown. (C) Molecular function analysis of hsa_circRNA_104980-targeting genes is shown. (D) KEGG pathway analysis shows the top 10 pathways related to hsa_circRNA_104980. (E) The related pathways were ranked according to GeneRatio (Selection Counts/Selection Size). (F) Pathway analysis using Biocarta database shows the top 10 pathways related to hsa_circRNA_104980. (G) Pathway analysis using PANTHER database shows the top 10 pathways related to hsa_circRNA_104980. (H) Pathway analysis using Reactome database shows the top 10 pathways related to hsa_circRNA_104980.

Fig. 8.

GO analysis and KEGG pathway annotation of hsa_circRNA_104980. (A) Biological process analysis of hsa_circRNA_104980-targeting genes. (B) Cellular component analysis of hsa_circRNA_104980-targeting genes is shown. (C) Molecular function analysis of hsa_circRNA_104980-targeting genes is shown. (D) KEGG pathway analysis shows the top 10 pathways related to hsa_circRNA_104980. (E) The related pathways were ranked according to GeneRatio (Selection Counts/Selection Size). (F) Pathway analysis using Biocarta database shows the top 10 pathways related to hsa_circRNA_104980. (G) Pathway analysis using PANTHER database shows the top 10 pathways related to hsa_circRNA_104980. (H) Pathway analysis using Reactome database shows the top 10 pathways related to hsa_circRNA_104980.

Close modal

Recent studies have demonstrated that circRNAs are an abundant, stable and conserved class of RNA molecules that can function as miRNA sponges in gene regulation to affect disease initiation and progression [8, 22-24]. Other studies have also revealed important functions for circRNAs, including protein sequestration [25], transcriptional regulation [23], and potential functions in cancer [10]. Therefore, circRNAs have been considered important biological regulators for understanding the molecular mechanisms of disease and identifying effective diagnostic biomarkers or therapeutic targets.

In this study, we performed a comprehensive analysis of the global changes in the expression pattern of circRNAs in osteosarcoma. This study reports for the first time the comprehensive expression profile of circRNAs in osteosarcoma. And systemic bioinformatics analysis was conducted to identify circRNAs essential for the biological processes of osteosarcoma, which may provide potential targets for the development of novel diagnostic and therapeutic strategies against osteosarcoma.

A total of 4660 human circRNAs were detected. There were 252 differentially expressed circRNAs, consisting of 71 up-regulated and 181 down-regulated circRNAs. Among them, 12 differentially expressed circRNAs were selected for validation. The great consistency between the real-time qPCR results and microarray data confirmed the high reliability of the microarray expression data, which supports the further functional analysis of the differentially expressed circRNAs.

Following circRNA detection within the tissue samples, we found that hsa_ circRNA_103801 was up-regulated in both osteosarcoma cell lines and osteosarcoma tissue samples, while hsa_circRNA_104980 was down-regulated. This interesting result drove us to focus on these two circRNAs, but not on others. The result indicates that hsa_circRNA_103801 and hsa_circRNA_104980 may play an important role in the biological processes of osteosarcoma, such as cell proliferation, differentiation, and invasion, among other similar processes. Therefore, both prediction of the circRNA/microRNA/mRNA interaction network and bioinformatics analysis were conducted for a more comprehensive understanding of these two circRNAs.

In recent years, there is increasing amount of experimental evidence supporting the hypothesis that circRNAs may possess ceRNAs activity, which means circRNAs may have an irreplaceable role in regulatory RNA networks. And circRNAs may take an important role in regulating the biological behaviors of cancer cells through competitive binding with their target miRNAs.

The most likely potential target miRNAs for hsa_circRNA_103801 include hsa-miR-370-3p, hsa-miR-338-3p and hsa-miR-877-3p. Among them, hsa-miR-338-3p was found to be down-regulated in esophageal squamous cell carcinoma (ESCC) tissues compared with adjacent non-tumor tissues [26]. And the aberrant expression of hsa-miR-338-3p increased the risk of esophageal cancer. This result suggests that hsa-miR-338-3p may be a tumor suppressor gene, which is associated with the incidence and development of cancer. In our study, hsa_circRNA_103801 was up-regulated in both osteosarcoma cell lines and osteosarcoma tissue samples, which indicates that hsa_circRNA_103801 may be an oncogene in osteosarcoma and may promote the incidence and development of osteosarcoma by suppressing hsa-miR-338-3p activity.

Intracellular circRNAs with competing endogenous RNAs (ceRNAs) activity may compete with linear RNAs by binding miRNAs with MREs, which strongly reduce the ability of miRNAs to bind their target mRNAs and result in increased expression levels of these genes [7, 8]. In this study, we have found that hsa_circRNA_103801 expression was significantly correlated with the expression of TANC1, one of the putative targets of its putative target miRNAs (Fig. 5E). This result suggests that hsa_circRNA_103801 may increase the expression of TANC1 through competitive binding with hsa-miR-877-3p, whose target gene is TANC1. Avirneni-Vadlamudi U indentified misregulated myoblast fusion caused by ectopic TANC1 expression as a Rhabdomyosarcoma (RMS) neoplasia mechanism [27]. And hsa_circRNA_103801 may act as an oncogene and promote the incidence and development of osteosarcoma by hsa_circRNA_103801/hsa-miR-877-3p/TANCl pathway, which highly reflects the effect of ceRNA regulatory network.

At the same time, KEGG analysis showed that there were 19 pathways related to hsa_circRNA_103801, including breast cancer, proteoglycans in cancer, pathways in cancer, and transcriptional misregulation in cancer. When the related pathways were ranked according to GeneRatio, the most related pathway was pathways in cancer, which strongly suggests that hsa_circRNA_103801 may be of importance in the incidence and development of osteosarcoma. And we found that hsa_circRNA_103801 was also involved in the HIF-1 signaling pathway, the Rapl signaling pathway and the PI3K-Akt signaling pathway. It is a meaningful finding that all the four annotation tools for pathway analysis (KEGG, Biocarta, PANTHER and Reactome) suggested the important role of hsa_circRNA_103801 in HIF-1, VEGF and angiogenesis pathways. Recently, El-Naggar AM et al. [28]. demonstrated that translational activation of HIF1α may promote sarcoma metastasis, which indicates that hsa_circRNA_103801 may promote osteosarcoma metastasis through the HIF-1, VEGF and angiogenesis pathways. And the Rap1 signaling pathway has been reported to regulate tumorigenesis and tumor progression [29, 30]. There is plenty of evidence showing the important role of the PI3K-Akt signaling pathway in osteosarcoma [31-33]. All these results strongly indicate that hsa_circRNA_103801 may be closely related to the initiation and progression of osteosarcoma.

The predicted ceRNA network can help to clarify the molecular mechanisms of hsa_ circRNA_103801 in osteosarcoma. We would focus on the biological functions and molecular mechanisms of hsa_circRNA_103801 in follow-up studies.

For hsa_circRNA_104980, the most likely potential target miRNAs consist of hsa-miR-1298-3p and hsa-miR-660-3p. But the biological functions of hsa-miR-1298-3p and hsa-miR-660-3p in cancer have not been reported. In the present study, hsa_circRNA_104980 was down-regulated in both osteosarcoma cell lines and osteosarcoma tissue samples, indicating hsa_circRNA_104980 may function as a tumor suppressor gene in osteosarcoma. GO analysis revealed that hsa_circRNA_104980-related genes took part in binding phosphoric ester hydrolase activity, and carbohydrate derivative binding, among other similar activities. KEGG analysis found that there were 22 pathways related to hsa_circRNA_104980, which were ranked in order of GeneRatio. The other annotation tools (Biocarta, PANTHER and Reactome) also suggested the related pathways of hsa_circRNA_104980, which may provide some useful information for understanding its biological functions. These bioinformatics analysis data could help to explore the role of hsa_circRNA_104980. Since the function of hsa_circRNA_104980-targeting miRNAs remains unknown, it is worth further study to determine the role of hsa_circRNA_104980 in osteosarcoma.

In summary, the present study identified the comprehensive expression profile of circRNAsinosteosarcomaforthefirsttime.TheceRNAsnetworkprediction and bioinformatics analysis could provide a comprehensive understanding of hsa_circRNA_103801 and hsa_ circRNA_104980, which may be involved in the initiation and progression of osteosarcoma. The comprehensive expression profile of circRNAs in osteosarcoma may contribute to the further study of potential diagnostic biomarkers or therapeutic targets.

The authors declare no conflict of interest.

This research was supported by National Natural Science Foundation of China (81472506 to Jingnan Shen) and Natural Science Foundation of Guangdong Province (2016A030313227 and S2013010016847 to Junqiang Yin) and Young teachers cultivating project of Sun Yat-Sen University (14ykpy16 to Junqiang Yin).

1.
Biermann JS, Adkins DR, Agulnik M, Benjamin RS, Brigman B, Butrynski JE, Cheong D, Chow W, Curry WT, Frassica DA, Frassica FJ, Hande KR, Hornicek FJ, Jones RL, Mayerson J, McGarry SV, McGrath B, Morris CD, ODonnell RJ, Randall RL, Santana VM, Satcher RL, Siegel HJ, von Mehren M, Bergman MA, Sundar H: National complrehensive cancer n: Bone cancer. J Natl Compr Cane Netw 2013; 11: 688-723.
2.
Heare T, Hensley MA, Dell’Orfano S: Bone tumors: osteosarcoma and Ewings sarcoma. Curr Opin Pediatr 2009; 21: 365-372.
3.
Bielack SS, Kempf-Bielack B, Delling G, Exner GU, Flege S, Helmke K, Kotz R, Salzer-Kuntschik M, Werner M, Winkelmann W, Zoubek A, Jurgens H, Winkler K: Prognostic factors in high-grade osteosarcoma of the extremities or trunk: an analysis of 1, 702 patients treated on neoadjuvant cooperative osteosarcoma study group protocols. J Clin Oncol 2002; 20: 776-790.
4.
Kramer MC, Liang D, Tatomer DC, Gold B, March ZM, Cherry S, Wilusz JE: Combinatorial control of Drosophila circular RNA expression by intronic repeats, hnRNPs, and SR proteins. Genes Dev 2015; 29: 2168-2182.
5.
Barrett SP, Salzman J: Circular RNAs: analysis, expression and potential functions. Development 2016; 143: 1838-1847.
6.
Wu HJ, Zhang CY, Zhang S, Chang M, Wang HY: Microarray Expression Profile of Circular RNAs in Heart Tissue of Mice with Myocardial Infarction-Induced Heart Failure. Cell Physiol Biochem 2016; 39: 205-216.
7.
Hansen TB, Jensen TI, Clausen BH, Bramsen JB, Finsen B, Damgaard CK, Kjems J: Natural RNA circles function as efficient micro RNA sponges. Nature 2013; 495: 384-388.
8.
Qu S, Yang X, Li X, Wang J, Gao Y, Shang R, Sun W, Dou K, Li H: Circular RNA: A new star of noncoding RNAs. Cancer Lett 2015; 365: 141-148.
9.
Beermann J, Piccoli MT, Viereck J, Thum T: Non-coding RNAs in Development and Disease: Background, Mechanisms, and Therapeutic Approaches. Physiol Rev 2016; 96: 1297-1325.
10.
Guarnerio J, Bezzi M, Jeong JC, Paffenholz SV, Berry K, Naldini MM, Lo-Coco F, Tay Y, Beck AH, Pandolfl PP: Oncogenic Role of Fusion-circRNAs Derived from Cancer-Associated Chromosomal Translocations. Cell 2016; 165: 289-302.
11.
Nair AA, Niu N, Tang X, Thompson KJ, Wang L, Kocher JP, Subramanian S, Kalari KR: Circular RNAs and their associations with breast cancer subtypes. Oncotarget 2016; 7: 80967-80979.
12.
Xu L, Zhang M, Zheng X, Yi P, Lan C, Xu M: The circular RNA ciRS-7 (Cdrlas) acts as a risk factor of hepatic microvascular invasion in hepatocellular carcinoma. J Cancer Res Clin Oncol 2017; 143: 17-27.
13.
Zou CY, Wang J, Shen JN, Huang G, Jin S, Yin JQ, Guo QC, Li HM, Luo L, Zhang M, Zhang LJ: Establishment and characteristics of two syngeneic human osteosarcoma cell lines from primary tumor and skip metastases. Acta Pharmacol Sin 2008; 29: 325-332.
14.
Serra M, Scotlandi K, Manara MC, Maurici D, Lollini PL, De Giovanni C, Toffoli G, Baldini N: Establishment and characterization of multidrug-resistant human osteosarcoma cell lines. Anticancer Res 1993; 13: 323-329.
15.
Enright AJ, John B, Gaul U, Tuschl T, Sander C, Marks DS: MicroRNA targets in Drosophila. Genome Biol 2003; 5:R1.
16.
Pasquinelli AE: MicroRNAs and their targets: recognition, regulation and an emerging reciprocal relationship. Nat Rev Genet 2012; 13: 271-282.
17.
Tay Y, Kats L, Salmena L, Weiss D, Tan SM, Ala U, Karreth F, Poliseno L, Provero P, Di Cunto F, Lieberman J, Rigoutsos I, Pandolfl PP: Coding-independent regulation of the tumor suppressor PTEN by competing endogenous mRNAs. Cell 2011; 147: 344-357.
18.
Memczak S, Jens M, Elefsinioti A, Torti F, Krueger J, Rybak A, Maier L, Mackowiak SD, Gregersen LH, Munschauer M, Loewer A, Ziebold U, Landthaler M, Kocks C, le Noble F, Rajewsky N: Circular RNAs are a large class of animal RNAs with regulatory potency. Nature 2013; 495: 333-338.
19.
Miranda KC, Huynh T, Tay Y, Ang YS, Tam WL, Thomson AM, Lim B, Rigoutsos I: A pattern-based method for the identification of MicroRNA binding sites and their corresponding heteroduplexes. Cell 2006; 126: 1203-1217.
20.
Hammell M, Long D, Zhang L, Lee A, Carmack CS, Han M, Ding Y, Ambros V: mirWIP: microRNA target prediction based on microRNA-containing ribonucleoprotein-enriched transcripts. Nat Methods 2008; 5: 813-819.
21.
Ritchie W, Flamant S, Rasko JE: Predicting microRNA targets and functions: traps for the unwary. Nat Methods 2009; 6: 397-398.
22.
Hansen TB, Kjems J, Damgaard CK: Circular RNA and miR-7 in cancer. Cancer Res 2013; 73: 5609-5612.
23.
Li Z, Huang C, Bao C, Chen L, Lin M, Wang X, Zhong G, Yu B, Hu W, Dai L, Zhu P, Chang Z, Wu Q, Zhao Y, Jia Y, Xu P, Liu H, Shan G: Exon-intron circular RNAs regulate transcription in the nucleus. Nat Struct Mol Biol 2015; 22: 256-264.
24.
Cortes-Lopez M, Miura P: Emerging Functions of Circular RNAs. Yale J Biol Med 2016; 89: 527-537.
25.
Du WW, Yang W, Liu E, Yang Z, Dhaliwal P, Yang BB: Foxo3 circular RNA retards cell cycle progression via forming ternary complexes withp21 and CDK2. Nucleic Acids Res 2016; 44: 2846-2858.
26.
Yang M, Liu R, Sheng J, Liao J, Wang Y, Pan E, Guo W, Pu Y, Yin L: Differential expression profiles of microRNAs as potential biomarkers for the early diagnosis of esophageal squamous cell carcinoma. Oncol Rep 2013; 29: 169-176.
27.
Avirneni-Vadlamudi U, Galindo KA, Endicott TR, Paulson V, Cameron S, Galindo RL: Drosophila and mammalian models uncover a role for the myoblast fusion gene TANC1 in rhabdomyosarcoma. J Clin Invest 2012; 122: 403-407.
28.
El-Naggar AM, Veinotte CJ, Cheng H, Grunewald TG, Negri GL, Somasekharan SP, Corkery DP, Tirode F, Mathers J, Khan D, Kyle AH, Baker JH, LePard NE, McKinney S, Hajee S, Bosiljcic M, Leprivier G, Tognon CE, Minchinton AI, Bennewith KL, Delattre O, Wang Y, Dellaire G, Berman JN, Sorensen PH: Translational Activation of HIFlalpha by YB-1 Promotes Sarcoma Metastasis. Cancer Cell 2015; 27: 682-697.
29.
Ishida D, Kometani K, Yang H, Kakugawa K, Masuda K, Iwai K, Suzuki M, Itohara S, Nakahata T, Hiai H, Kawamoto H, Hattori M, Minato N: Myeloproliferative stem cell disorders by deregulated Rapl activation in SPA-1-deficient mice. Cancer Cell 2003; 4: 55-65.
30.
Okada T, Sinha S, Esposito I, Schiavon G, Lopez-Lago MA, Su W, Pratilas CA, Abele C, Hernandez JM, Ohara M, Okada M, Viale A, HeguyA, Socci ND, Sapino A, Seshan VE, Long S, Inghirami G, Rosen N, Giancotti FG: The Rho GTPase Rndl suppresses mammary tumorigenesis and EMT by restraining Ras-MAPK signalling. Nat Cell Biol 2015; 17: 81-94.
31.
Moriarity BS, Otto GM, Rahrmann EP, Rathe SK, Wolf NK, Weg MT, Manlove LA, LaRue RS, Temiz NA, Molyneux SD, Choi K, Holly KJ, Sarver AL, Scott MC, Forster CL, Modiano JF, Khanna C, Hewitt SM, Khokha R Yang Y, Gorlick R, Dyer MA, Largaespada DA: A Sleeping Beauty forward genetic screen identifies new genes and pathways driving osteosarcoma development and metastasis. Nat Genet 2015; 47: 615-624.
32.
Li YJ, Dong BK, Fan M, Jiang WX: BTG2 inhibits the proliferation and metastasis of osteosarcoma cells by suppressing the P13K/AKT pathway. Int J ClinExp Pathol 2015; 8: 12410-12418.
33.
Zhao J, Zhang ZR, Zhao N, Ma BA, Fan QY: VEGF Silencing Inhibits Human Osteosarcoma Angiogenesis and Promotes Cell Apoptosis via PI3K/AKT Signaling Pathway. Cell Biochem Biophysiol 2015; 73: 519-525.

W. Liu, J. Zhang and C. Zou contributed equally to this work.

Open Access License / Drug Dosage / Disclaimer
This article is licensed under the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License (CC BY-NC-ND). Usage and distribution for commercial purposes as well as any distribution of modified material requires written permission. Drug Dosage: The authors and the publisher have exerted every effort to ensure that drug selection and dosage set forth in this text are in accord with current recommendations and practice at the time of publication. However, in view of ongoing research, changes in government regulations, and the constant flow of information relating to drug therapy and drug reactions, the reader is urged to check the package insert for each drug for any changes in indications and dosage and for added warnings and precautions. This is particularly important when the recommended agent is a new and/or infrequently employed drug. Disclaimer: The statements, opinions and data contained in this publication are solely those of the individual authors and contributors and not of the publishers and the editor(s). The appearance of advertisements or/and product references in the publication is not a warranty, endorsement, or approval of the products or services advertised or of their effectiveness, quality or safety. The publisher and the editor(s) disclaim responsibility for any injury to persons or property resulting from any ideas, methods, instructions or products referred to in the content or advertisements.