Overexpression of miR-206 in osteosarcoma and its associated molecular mechanisms as assessed through TCGA and GEO databases.

Osteosarcoma (OS) is a primary malignant bone tumor that predominantly occurs in adolescents. Different types of OS tumor are highly malignant, associated with a poor prognosis and are invasive with blood-vessel dissemination characteristics, thus affected patients are prone to early lung metastasis. MicroRNAs (miRNAs/miR) are small non-coding RNA molecules that act as oncogenes or tumor suppressors during tumor development. The present study investigated the role of miR-206 in OS development. Bioinformatics analysis demonstrated that miR-206 was upregulated in OS and thus may serve as a risk factor for cancer prognosis. Subsequently, in response to miR-206 overexpression, differentially expressed genes were screened and analyzed using the Database for Annotation, Visualization and Integrated Discovery, Gene Ontology enrichment analysis, the Kyoto Encyclopedia of Genes and Genomes pathways and protein-protein interaction network construction, in order to identify key miR-206 targets. The results demonstrated that high miR-206 expression inhibited OS cell proliferation, which was associated with a good patient prognosis. Thus, miR-206 may serve as a potential target for OS treatment, in order to improve early disease diagnosis. Copyright: © Xu et al.

Introduction

In the United States, osteosarcoma (OS) is a highly prevalent primary bone tumor, which accounted for 0.2% of all human solid tumor malignancies from 1973 to 2004 (1). OS includes metaphysis of the proximal tibia, distal femur and other long tubular bones (2). Despite recent development of novel adjuvant chemotherapy techniques and surgical methods, which have increased the 5-year survival rate for OS to ~70% (3), OS mortality and metastasis rates remain high (4). Regardless of the identification of several anticancer drugs and tumor suppressors, the underlying molecular mechanisms in OS tumorigenesis remain unclear (5,6).

MicroRNAs (miRNAs/miRs) are a class of endogenous small non-coding RNAs that directly bind to the 3′-untranslated region (UTR) of target mRNAs and regulate gene expression (7). An increasing number of studies have reported that miRNA expression profiles are altered in various cancer cells and tissues, suggesting their value as biomarkers for cancer and targets for therapy (8,9).

miR-206 has been demonstrated to play vital roles in different types of cancer. For example, miR-206 is a well-known tumor suppressor in human breast cancer, which regulates estrogen receptor-α expression during normal breast epithelial cell development (10). Furthermore, miR-206 has been reported in ovarian (11), gastric (12,13), colorectal (14), laryngeal (15), cervical (16), lung (17) and liver (18) cancer. The present study investigated miR-206 expression levels in both primary and metastatic OS tissues compared with normal tissues, in order to determine whether miR-206 has the potential to serve as a biomarker for the diagnosis and treatment of patients with OS.

Materials and methods

Data sources

Datasets were retrieved from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). A total of 224 datasets were downloaded and screened, and common gene expression profiles were selected from the GSE65071 (19). The GSE65071 dataset, published in January 2015 and based on the GPL 19631 platform (G-U133A) (http://www.exiqon.com/mirna-pcr-panels) Exiqon human V3 miRNA PCR panel I+II, includes data from 15 normal controls and 20 patients with OS (10 primary OS and 10 metastatic OS).

A total of 64 datasets were retrieved from the GEO database in order to identify target OS genes. The data regarding the differentially expressed genes (DEGs) were obtained from the GSE89074 dataset (Han et al unpublished data). Gene expression profiles of two OS cell lines with miR-206 overexpression, two empty vector controls and two normal control cell lines were selected from the GSE89074 dataset. The GSE89074 dataset was published in October 2016 and based on the GPL570 platform [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array chip data (Affymetrix; Thermo Fisher Scientific, Inc.).

Screening the DEGs

GEO2R (20) was used to analyze the data from the miRNA expression profiles of selected OS cell lines and the DEGs derived from miR-206 overexpression. GEO2R is a web-based tool based on the limma R package (version 3.10; http://www.bioconductor.org). DEGs in OS were screened, with P<0.05 and |log fold-change (FC)|>1 set as the threshold values. DEGs from the primary OS and normal control cell lines were screened, as well as the metastatic OS and normal control cells lines, and DEGs sets of the two groups were obtained. A Venn diagram was used to cross DEGs between the two groups, and the intersecting genes were considered to be OS-associated.

Functional annotation of DEGs

The Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.ncifcrf.gov) was used to perform functional and pathway enrichment analyses. DAVID is a systematic and integrative functional annotation tool that allows researchers to unravel the biological meaning behind large lists of genes (21). Gene Ontology analysis, including the cellular component (CC), molecular function (MF) and biological process (BP) (22), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis (23) were performed for the upregulated and downregulated genes, respectively. P<0.05 was considered to indicate a statistically significant difference.

Protein-protein interaction (PPI)

The Search Tool for the Retrieval of Interacting Genes/Proteins (STRING; http://string-db.org) is a biological database and web resource (24), which was used to construct a PPI network of the DEGs. Based on the STRING database, PPIs of DEGs were selected with scores ≥0.9 (highest confidence), and the PPI networks were visualized using Cytoscape software (version 3.6.1; http://cytoscape.org).

Screening the hub genes

A plugin cyto-Hubba (version 3.6.1) (25) analysis was performed within Cytoscape to detect hub genes with the strongest interactions between other genes (26). A total of 10 genes with high degree scores were identified and selected in the PPI network.

Survival analysis in Gene Expression Profiling Interactive Analysis (GEPIA)

Following collection of the research subjects from The Cancer Genome Atlas (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga), the online database, GEPIA (http://gepia.cancer-pku.cn) was used to assess the association between gene expression and prognosis (27). The effect of the genes on the prognosis of patients with OS was evaluated and key genes that influence OS prognosis were screened. P<0.05 was considered to indicate a statistically significant difference.

Dataset validation

The association between the expression levels of the five key genes and pulmonary metastasis of osteosarcoma was validated using the GSE14359 dataset within the GEO database.

Results

Screening for differentially expressed miRNAs

miRNAs differentially expressed in the primary OS tissues were screened within the GSE65071 dataset. A total of 277 differentially expressed miRNAs were identified, of which 66 were downregulated and 211 were upregulated. miRNAs differentially expressed in the OS tissues with lung metastasis were subsequently screened, which identified 265 differentially expressed miRNAs (58 downregulated and 207 upregulated). A total of 253 differentially expressed miRNAs were identified (Fig. 1). Previous studies have highlighted miR-206 as an important cancer-associated miRNA (28,29); however, to the best of our knowledge, its role in OS development following upregulation in OS tissues has not yet been investigated. The results of the present study demonstrated that miR-206 expression was increased in both the primary OS tissues and metastatic OS tissues compared with normal tissues, respectively. The difference was statistically significant in both cases (P<0.001; Fig. 2).

Figure 1.

Differentially expressed genes between primary OS and metastatic tumors compared with normal tissues in the GSE65071 dataset. OS, osteosarcoma.

Figure 2.

miR-206 expression in each sample of the GSE65071 dataset. ***P<0.001 vs. normal. miR-206, microRNA-206; Normal, normal tissue; Primary, primary OS; Metastatic, metastatic OS; OS, osteosarcoma.

Screening for DEGs in cells overexpressing miR-206 within the GSE89074 dataset

A total of 2,057 DEGs were obtained from cells overexpressing miR-206, including 1,540 upregulated and 517 downregulated genes. All DEGs are presented in the volcanic map (Fig. 3).

Figure 3.

Volcano plot presenting the differentially expressed genes in osteosarcoma cells overexpressing microRNA-206 in the GSE89074 dataset. Red, upregulated genes; green, downregulated genes; black, no significant differences.

Enrichment analysis

Functional pathway enrichment analysis was performed for the DEGs that were upregulated and downregulated in response to miR-206 overexpression. Functional enrichment analysis demonstrated that the upregulated genes were significantly enriched in 127 BPs, 65 CCs and 34 MFs. Subsequently, the present study identified the 10 most significantly enriched BPs, CCs and MFs in the upregulated genes (Fig. 4). Conversely, functional enrichment analysis demonstrated that the downregulated genes were significantly enriched in 152 BPs, 16 CCs and 44 MFs. Fig. 5 indicates the 10 most significantly enriched BPs, CCs and MFs in the downregulated genes.

Figure 4.

Enrichment analysis of the upregulated genes in the top 10 BPs, CCs and MFs. BP, biological process; CC, cell composition; MF, molecular function. The length of the histogram represents a significant degree of enrichment.

Figure 5.

Enrichment analysis of the downregulated genes in the top 10 BPs, CCs and MFs. BP, biological process; CC, cell composition; MF, molecular function.

Analysis of KEGG pathways

The present study analyzed signal pathway enrichment of the DEGs in cells overexpressing miR-206. The upregulated genes were notably enriched in 32 signaling pathways, of which 10 were associated with OS (Fig. 6; Table I). The downregulated genes were notably enriched in 35 signaling pathways, of which 10 were associated with OS (Fig. 7; Table II).

Figure 6.

Top 10 most significantly enriched signaling pathways for the miR-206 upregulated genes.

Table I.

A total of 10 enriched signaling pathways for microRNA-206 upregulated genes.

Term	Signaling pathway	Count	P-value
hsa05146	Amoebiasis	20	9.64×10−6
hsa04510	Focal adhesion	28	5.83×10−5
hsa05205	Proteoglycans in cancer	26	5.72×10−4
hsa05200	Pathways in cancer	38	2.44×10−3
hsa04110	Cell cycle	16	5.55×10−3
hsa05222	Small cell lung cancer	12	1.01×10−2
hsa04151	PI3K-Akt	32	1.04×10−2
hsa04919	Thyroid hormone	14	1.62×10−2
hsa04068	FoxO	15	2.42×10−2
hsa04310	Wnt	15	3.02×10−2

Figure 7.

Top 10 most significantly enriched signaling pathways for the miR-206 downregulated genes.

Table II.

A total of 10 enriched signaling pathways for microRNA-206 downregulated genes.

Term	Signaling pathway	Count	P-value
hsa04010	MAPK	27	1.05×10−7
hsa04630	JAK-STAT	15	4.50×10−5
hsa04620	Toll-like receptor	12	6.78×10−4
hsa04668	TNF	10	2.79×10−3
hsa04064	NF-κB	8	1.02×10−2
hsa04151	PI3K-Akt	18	1.50×10−2
hsa04014	Ras	13	2.31×10−2
hsa05206	MicroRNAs in cancer	15	2.76×10−2
hsa04068	FoxO	9	3.26×10−2
hsa04115	p53	6	3.78×10−2

Network maps and hub gene screening

A total of 1,540 upregulated and 517 downregulated DEGs were uploaded onto STRING to obtain the PPI data. Samples with PPI scores ≥0.9 were selected to construct the PPI network. The PPI network of the upregulated genes consisted of 1,129 nodes and 1,862 edges (Fig. 8). The central nodes of this network were as follows: IL6, FOS, JUN, IRF7, EGR1, OAS1, OAS2, MX1, XAF1 and IFIT3 (Fig. 9). The PPI network of the downregulated genes consisted of 144 nodes and 545 edges (Fig. 10). The central nodes of this network were as follows: PDGFB, NEK2, CHEK1, CCNB1, RBBP4, ANAPC4. HACE1, FBXL5, HERC2 and VPRBP (Fig. 11).

Figure 8.

Protein-protein interaction network of the upregulated genes.

Figure 9.

Upregulated hub genes in the protein-protein interaction network. The different colors represent the difference in significance of the hub genes, in terms of their degree of connectivity. An increased shade gradient indicates a more notable association with other genes.

Figure 10.

Protein-protein interaction network of the downregulated genes.

Figure 11.

Downregulated hub genes in the protein-protein interaction network. The different colors represent the difference in significance of the hub genes, in terms of their degree of connectivity. An increased shade gradient indicates a more notable association with other genes.

Survival analysis

The association between the upregulated and downregulated genes, and the survival of patients with OS was analyzed. A total of five genes (CCNB1, CHEK1, IL6, NEK2 and RBBP4) demonstrated a significant prognostic value. Progression-free survival (PFS) time of patients with high CCNB1 expression was lower than those with low CCNB1 expression (Fig. 12A). Furthermore, PFS (Fig. 12B) and overall survival time (Fig. 12D) were lower in patients with high CHEK1 expression than those with low CHEK1 expression. Similarly, PFS (Fig. 12C) and overall survival time (Fig. 12F) were lower in patients with high NEK2 expression than those with low NEK2 expression. Overall survival time of patients with high RBBP4 expression was lower than those with low RBBP4 expression (Fig. 12G); however, overall survival of patients with high IL6 expression was higher than those with low IL6 expression (Fig. 12E).

Figure 12.

Kaplan-Meier survival curves. (A-C) Progression-free survival analysis with high and low expression levels of CCNB1, CHEK1 and NEK2. (D-G) Overall survival analysis with high and low expression levels of CHEK1, IL6, NEK2 and RBBP4. HR, hazard ratio.

The present study validated the association between gene expression and OS type, in primary OS or pulmonary metastasis of OS, using the GSE14359 dataset (30) within the GEO database. The dataset contained 10 OS samples with pulmonary metastasis and 8 primary OS samples. The mRNA expression levels of the five genes were upregulated in OS lung metastasis compared with primary OS (Fig. 13), with the largest upregulation difference in NEK2 (Log2 FC, 2.28639273; P=2.12×10−3), and the smallest upregulation difference in RBBP4 (Log2 FC, 0.69641798; P=1.81×10−2). Similarly, increased expression levels of CCNB1 (Log2 FC, 1.85398694; P=1.99×10−3), CHEK1 (Log2 FC, 1.13599575; P=7.20×10−3) and IL6 (Log2 FC, 1.23455773; P=4.49×10−3) were observed in OS lung metastasis compared with primary OS. The results suggest that the five genes identified for their prognostic value are closely associated with the metastasis and prognosis of OS.

Figure 13.

Expression levels of the five key genes in primary OS and OS lung metastasis in the GSE14359 dataset. **P<0.05. OS, osteosarcoma; metastasis: Metastatic OS of the lung; primary, primary OS.

Discussion

Increasing evidence demonstrates the role of miRNAs in OS tumorigenesis and tumor development (31,32). miRNAs and their target genes represent potential novel therapeutic biomarkers for OS (33,34). Previous studies have reported downregulated miR-206 expression in OS cells (35–38). However, miR-206 expression was demonstrated to be upregulated in human OS tissues in the present study. A possible reason for the discrepancies observed may be due to the opposing roles miR-206 plays at different stages of OS occurrence and development. For example, miR-206 expression was downregulated in the plasma of patients with early OS, while expression was upregulated in advanced OS. The databases screened in the present study contained data from plasma samples of patients with advanced OS.

miR-206 is transcribed by RNA polymerase II to produce pri-miRNA transcripts (pri-miR-206). Pre-miR-206 precursors with stem-ring structures are produced by processing pri-miR-206 in the nucleus with the RNA endonuclease III, Drosha. Subsequently, the pri-miR-206 is transported to the cytoplasm by the Exportin-5 protein, and further processed by the secondary RNA endonuclease III, Dicer, in order to produce mature double-stranded RNA molecules. One of the mature strands is inserted into the RNA-induced silencing complex, which binds to the 3′-UTRs of target genes and cleaves target RNAs (39). miR-206 has been reported to inhibit the expression of multiple target genes, which are also regulated by multiple miRNAs (40).

The identification of target genes is critical to understanding the role of miRNAs during tumorigenesis (41). In the present study, overexpression of miR-206 resulted in downregulation of the CCNB1 and NEK2 genes, and upregulation of the IL6 gene. CCNB1 is associated with mitosis (42), whereby its aberrant cell cycle regulation is a major cause of excessive cell proliferation and tumorigenesis (43). CCNB1 is closely associated with tumor progression, where its overexpression in tumor cells and tissues leads to uncontrolled phosphorylation and dysregulation of the maturation promotion factor (MPF). The MPF is activated following DNA damage and the affected cells progress through mitosis, proliferating to form different types of tumor (44). Thus, the CCNB1 gene is considered an oncogene and tumor antigen (45). The present study constructed a PPI network map and analyzed key genes, which demonstrated that CCNB1 was downregulated in OS tissue. It is believed that CCNB1 may play a role in the occurrence and development of OS as an anti-oncogene, which can be targeted for the treatment of OS.

Previous studies have demonstrated that NEK2 expression is upregulated in several types of human cancer, including non-small cell lung carcinoma (NSCLC) (46,47), myeloma (48), ovarian cancer (49), breast cancer (50,51), prostate cancer (52), colorectal cancer (53), malignant peripheral neurilemmoma (54), renal cell carcinoma (50) and pancreatic ductal adenocarcinoma (55), compared with the corresponding normal tissues. NEK2 mediates the separation of chromosomes into two daughter cells by regulating centrosome separation and spindle formation. Aberrant NEK2 expression is associated with unregulated cell division through the premature separation of immature centrosomes, abnormal spindle formation, excessive centrosome duplication and abnormal chromosome segregation, these abnormalities promote chromosome aneuploidy and instability, aberrant NEK2 expression is believed to be a key driving force for cellular deterioration in cancer (48). The present study demonstrated that NEK2 expression was downregulated in OS tissues, thus it may serve as an anti-oncogene that can be targeted for OS prevention and treatment.

The Janus kinase 2/signal transducer and activator of transcription 3 (JAK2/STAT3) signaling pathway is one of the major signaling pathways by which IL6 exerts its biological effects (56). The results of the present study demonstrated that the JAK-STAT signaling is significantly enriched in downregulated genes. The occurrence and development of several types of human tumor are closely associated with abnormal IL6 expression (57). Furthermore, a previous study confirmed the role of miR-206 and IL6 in NSCLC via STAT3 signaling (58).

IL6 promotes tumor cell proliferation and angiogenesis by inducing epithelial-to-mesenchymal transition, promoting the expansion and recruitment of myeloid inhibitory cells, altering the inherent biological characteristics of tumor cells and optimizing the external growth environment of several types of tumor (59). IL6 activates STAT3 in OS cells, promoting proliferation and migration. STAT3 activation stimulates the expression of genes associated with cell proliferation, anti-apoptosis, hypoxia, metastasis and angiogenesis (60). These genes include CCND1, cell division cycle protein 2, B-cell lymphoma 2, hypoxia-inducible factor 1-α, heat shock protein 90 and vascular endothelial growth factor (VEGF) (61). STAT3 signaling plays an important role in OS progression. LLL12, a STAT3 inhibitor, notably inhibits the expression of VEGF, matrix metallopeptidase 9 and fibroblast growth factor-1 in OS cells, effectively hindering angiogenesis in vivo and in vitro (62). Furthermore, LLL12 promotes OS cell apoptosis and simultaneously impairs cell adhesion and migration. In addition, OS growth in nude mice is markedly inhibited (63). As IL6 promotes OS development via STAT3 signaling, it has the potential to be used as a target for the prevention and treatment of OS.

The present study confirmed that miR-206 is highly expressed in OS. miR-206 promotes OS development by regulating target gene networks via specific signaling pathways. The potential target genes and biological function of miR-206 provide novel insight into the DEGs of OS. Overall, the results indicate that miR-206 may be used as a novel target for the early diagnosis and treatment of OS.