This article is mentioned in:

Introduction

Osteosarcoma (OS) is the most common malignancy of bone in early adolescence (1). Conventional OS, also termed classical OS, is a common type of OS and is universally life-threatening due to its rapid growth, high local aggression and metastatic potential (2). During previous years, considerable progress has been made in identifying the key components in conventional OS, including genes, pathways and microRNAs (miRNAs). For example, during osteoblast differentiation, miRNA (miR)-34 is significantly induced by bone morphogenetic protein 2, and regulates multiple components of the Notch signalling pathway, including Notch1, Notch2 and jagged 1, which affects osteoclast differentiation. This regulatory association may be closely associated with the pathogenesis of OS (3). In addition, phosphatase and tensin homolog (PTEN) has been found to be a potent regulator of the phosphatidylinositol 3-kinase (PI3K) /serine-threonine kinase (Akt) pathway (4), and the loss of PTEN is a common occurrence in conventional OS (5). A previous study has showed that the expression of PTEN can be inhibited by miR-221, which potentiates the PI3K/Akt pathway in the conventional pathogenesis of OS (6). PTEN is also a target of miR-92a, and of members of the miR-17 and miR-130/301 families in OS (7).

In 2010, using genome-wide microarrays, Fritsche-Guenther et al (8) found that the aberrant expression of ephrin receptor A2 (EphA2) and its ligand, EFNA1 in OS can modulate the activation of the mitogen-activated protein kinase (MAPK) pathway. In addition, it was found that the expression of CD52 was higher in OS metastases compared with conventional OS metastases, and CAMPATH-1H, an antibody directed against CD52, reduced the number of viable OS cells (9). In 2013, Luo et al (10) found numerous differentially expressed genes (DEGs) and regulatory associations between transcription factors and DEGs in OS using the microarray data deposited by Fritsche-Guenther et al For example, interleukin 6 can be regulated by tumour protein p53 (TP53), nuclear factor I/C (CCAAT-binding transcription factor), retinoic acid receptor α, and CCAAT/enhancer binding protein β. In 2014, Yang et al (11) also identified a number of DEGs, Gene Ontology (GO) terms, including protein binding, and significant pathways, including focal adhesion, in OS based on a meta-analysis of eight expression profiles, including the one deposited by Fritsche-Guenther (8). However, in these previous studies, the potential miRNAs and regulatory associations between miRNAs and DEGs in OS were not examined.

In the present study, to screen and identify additional DEGs and miRNAs in conventional OS, the microarray data deposited by Fritsche-Guenther (8) were downloaded. Following GO and pathway enrichment analyses, and construction of a protein-protein interaction (PPI) network for the DEGs, the potential miRNAs in the most significant pathway for the upregulated DEGs were identified, and a regulatory network for the miRNAs-DEGs was constructed. The results were expected to assist in elucidating the aetiology of conventional OS, and provide more information to assist in the clinical diagnosis and treatment of this disease.

Materials and methods

Affymetrix microarray data

The GSE14359 (8) gene expression profile data were acquired from the Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/), which was based on the platform of the GPL96 [HG-U133A] Affymetrix Human Genome U133A Array. This dataset contains 10 conventional OS samples from the femur or tibia (two replicates each) from five consenting patients with grade 2–3 conventional OS between 7 and 74 years of age; eight OS lung metastasis tumour samples (two replicates each) from four consenting patients with a grade 1–3 OS lung metastatic tumour; and two non-neoplastic primary osteoblast cell samples with limited life span in vitro from one patient (two replicates). These 10 conventional OS samples and two non-neoplastic primary osteoblast samples were selected for further analysis.

The CEL files and probe annotation files were downloaded, and the gene expression data of all samples were preprocessed via background correction, quantile normalization and probe summarization using the Gene Chip Robust Multi Array algorithm (12) in the Affy software package (version 1.32.0; http://www.bioconductor.org/packages/release/bioc/html/affy.html) (13).

DEG screening

The Linear Models for Microarray Data package (version 3.10.3; http://www.bioconductor.org/packages/2.9/bioc/html/limma.html) (14) of R was used to identify genes, which were significantly differentially expressed in the conventional OS samples. The raw P-value was adjusted by the Benjamin and Hochberg method (15), and only genes meeting the cut-off criteria of |log2 fold-change|>1 and adjusted P<0.01 were selected as DEGs.

Hierarchical clustering analysis of the DEGs

Hierarchical clustering is a common method used to determine clusters of similar data points in multidimensional spaces (16). The pheatmap package (version 1.08; https://cran.r-project.org/web/packages/pheatmap/) (17) was used to perform hierarchical clustering via joint between-within distances for the DEGs in the conventional OS and non-neoplastic primary osteoblasts samples.

GO and pathway enrichment analyses

The Database for Annotation, Visualization and Integrated Discovery (DAVID) provides a set of comprehensive functional annotation tools, which can be used to identify the biological meanings of abundant genes (18). P<0.01 was used as the cut-off criterion for GO and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis using DAVID (version 6.7; https://david-d.ncifcrf.gov/), based on the hypergeometric distribution algorithm.

PPI network construction

The Search Tool for the Retrieval of Interacting Genes database (version 10.0; http://string-db.org/), which provides experimental and predicted interaction information (19), was used to analyse the PPIs for DEGs by calculating the combined score, and a score >0.4 was selected as the cut-off criterion. Subsequently, the PPI network of the upregulated and downregulated DEGs was visualized using Cytoscape (version 3.2.0; http://cytoscape.org/) (20).

Screening and analysis of network modules

The network modules were obtained based on Molecular Complex Detection (MCODE) analysis (21) of the original PPI networks. The default parameters (degree cut-off, 2; node score cut-off, 0.2; K-core, 2) were used as the cut-off criteria for the network module screening.

In order to obtain further information on the gene functions and identify pathways closely associated with the DEGs, functional annotation analysis and subsequent pathway enrichment analysis of the network module with the highest MCODE score were performed using DAVID, with a P<0.01 cut-off.

Integrated miRNA-DEG regulatory network construction

The potential miRNAs for upregulated DEGs in the most significant pathway were predicted using the miRDB database (version 1.24.0; http://www.bioconductor.org/packages/2.8/bioc/html/maDB.html) (22), with a cut-off for the target score of ≥60. The binding sites of miRNAs in the human mRNAs > 800 were abandoned. The integrated miRNA-DEG regulatory network was then visualised with Cytoscape.

Results

Identification of DEGs

Following the data preprocessing, 11,107 probes were obtained. Based on the cut-off criteria, a total of 987 DEGs were screened from the conventional OS samples, including 317 upregulated genes and 670 downregulated genes. The hierarchical cluster analysis of the data revealed that it was possible to use the DEGs to accurately distinguish the conventional OS samples from the non-neoplastic primary osteoblast cell samples (Fig. 1).

Figure 1. Hierarchical clustering analysis of the differentially expressed genes in the conventional OS samples and non-neoplastic primary osteoblasts samples. Each row represents a single gene; each column represents a tissue sample. The gradual color change from green to red represents the change from upregulation to downregulation. R1 and R2 indicate two replicates of one sample. OS, osteosarcoma; Nor, non-neoplastic primary osteoblast sample; OS1-OS5, conventional OS samples.

Enrichment analysis of upregulated and downregulated DEGs

According to the GO functional annotation, the upregulated DEGs were predominantly enriched in GO terms associated with DNA replication, including MCM3, replication factor C (RFC)5, replication protein A3 (RPA3) and flap endonuclease 1 (FEN1), and cell cycle, including cyclin-dependent kinase 1 (CDK1), NDC80 kinetochore complex (NDC80), BUB1 mitotic checkpoint serine/threonine-protein kinase (BUB1) and mitotic arrest deficient 2 like 1 (MAD2L1). A number of downregulated DEGs, including caveolin 1 (CAV1), cadherin 13 (CDH13), vascular endothelial growth factor C (VEGFC) and transforming growth factor β receptor 3 (TGFBR3), were relevant to blood vessel development, whereas epidermal growth factor receptor (EGFR), TP53, VEGFB and MAPK1 were associated with the regulation of cell proliferation (Table IA).

Table I. GO terms and pathways enriched for the upregulated and downregulated DEGs.

Table I.

GO terms and pathways enriched for the upregulated and downregulated DEGs.

A, Top 10 GO terms for the upregulated and downregulated DEGs
Category	Term	Description	P-value	n	Genes
Up	GO:0006259	DNA metabolic process	8.09E-13	38	RPA3, FANCL, FEN1, DTL, CENPF, MCM3, RFC5, RFC3, RFC4, RFC2
	GO:0006260	DNA replication	3.90E-10	21	DTL, CENPF, MCM3, RPA3, RFC5, RFC3, RFC4, RFC2, RRM2, FEN1
	GO:0022403	Cell cycle phase	8.15E-10	30	AURKA, NCAPG, BUB1, CDK1, KIF11, DLGAP5, CENPF, BIRC5, NDC80, MAD2L1
	GO:0000278	Mitotic cell cycle	1.34E-09	28	NCAPG, BUB1, CDK1, KIF11, CENPF, BIRC5, NDC80, CDKN3, MAD2L1, ZWINT
	GO:0051301	Cell division	1.36E-09	25	NCAPG, BUB1, ASPM, CDK1, KIF11, CENPF, BIRC5, NDC80, MCM5, MAD2L1
	GO:0000280	Nuclear division	5.16E-09	21	CDK1, KIF11, CENPF, CDC23, NDC80, BIRC5, SMC4, MAD2L1, NCAPG, BUB1
	GO:0007067	Mitosis	5.16E-09	21	CDK1, KIF11, CENPF, CDC23, NDC80, BIRC5, SMC4, MAD2L1, NCAPG, BUB1
	GO:0000087	M phase of mitotic cell cycle	7.03E-09	21	CDK1, KIF11, CENPF, CDC23, NDC80, BIRC5, SMC4, MAD2L1, NCAPG, BUB1
	GO:0048285	Organelle fission	1.03E-08	21	CDK1, KIF11, CENPF, CDC23, NDC80, BIRC5, SMC4, MAD2L1, NCAPG, BUB1
	GO:0000279	M phase	1.18E-08	25	NCAPG, BUB1, CDK1, KIF11, DLGAP5, CDC23, CENPF, BIRC5, NDC80, MAD2L1
Down	GO:0001568	Blood vessel development	1.57E-06	28	CAV1, THBS1, MMP14, PNPLA6, CDH13, VEGFC, NTRK2, TGFBR3, ENG, TNFAIP2
	GO:0042127	Regulation of cell proliferation	1.63E-06	60	EGFR, CTBP1, TP53, MFGE8, HOXC10, VEGFB, MAPK1, VEGFC, SMAD3, SMAD2
	GO:0007242	Intracellular signaling cascade	2.08E-06	84	RRAS, TP53, CAV1, MAPKAPK3, FHL2, TGFBR3, GRK5, ABL1, CRK, IGFBP5
	GO:0010033	Response to organic substance	2.28E-06	56	TIMP3, STAT6, SRR, PPP3CB, COL6A2, SMAD2, CDH13, ADCY9, SMPD1, TGFBR3
	GO:0051270	Regulation of cell motion	2.47E-06	24	SMAD3, ACTN1, MAPK1, VEGFC, SEMA3F, TGFBR3, RRAS, THBS1, IGFBP3, IGFBP5
	GO:0001944	Vasculature development	2.52E-06	28	CAV1, MYH9, MMP14, PNPLA6, CDH13, VEGFC, NTRK2, TGFBR3, ENG, TNFAIP2
	GO:0008285	Negative regulation of cell proliferation	7.25E-06	34	CAV1, TP53, SMAD3, SMAD2, TGFBR3, ADAMTS1, IGFBP3, ENG, IGFBP5, TOB1
	GO:0040007	Growth	1.22E-05	22	SMAD2, LAMB2, NUPR1, DHCR7, SERPINE1, TGFBR3, BIN3, ADD1, IGFBP5, ERCC2
	GO:0030334	Regulation of cell migration	1.24E-05	21	EGFR, SMAD3, MAPK1, VEGFC, PTP4A1, RRAS, TGFBR3, THBS1, IGFBP3, IGFBP5
	GO:0048514	Blood vessel morphogenesis	3.44E-05	23	CAV1, CDH13, VEGFC, SEMA3C, PLCD1, NR2F2, THBS1, TNFAIP2, ENG, CYR61
B, Pathways enriched for the upregulated and downregulated DEGs
Category	Term	Description	P-value	n	Genes
Up	hsa03030	DNA replication	1.50E-06	9	RFC5, RFC3, RFC4, POLE3, RFC2, MCM3, FEN1, MCM5, RPA3
	hsa04110	Cell cycle	1.77E-04	12	CDC7, CDK1, CDKN1B, MAD2L1, CREBBP, BUB1, PRKDC, CDC23, MCM3, SMC1A
	hsa03040	Spliceosome	1.90E-04	12	RBM22, HNRNPA3, SF3B1, SNRPA1, MAGOH, TRA2B, SNRNP200, LSM5, LSM3, SNRPE
	hsa03430	Mismatch repair	1.96E-03	5	RFC5, RFC3, RFC4, RFC2, RPA3
	hsa03420	Nucleotide excision repair	3.76E-03	6	RFC5, RFC3, RFC4, POLE3, RFC2, RPA3
	hsa04114	Oocyte meiosis	4.63E-03	9	CDK1, MAD2L1, SLK, BUB1, CDC23, AURKA, PPP1CC, SMC1A, PPP1CB
	hsa03410	Base excision repair	9.30E-03	5	HMGB1, POLE3, TDG, POLB, FEN1
Down	hsa04510	Focal adhesion	4.08E-07	29	CAV1, COL6A1, THBS1, FN1, EGFR, FLNC, VEGFB, LAMA2, MAPK1, VEGFC
	hsa04142	Lysosome	1.74E-04	17	SGSH, CLN3, NAGLU, PLA2G15, CLTB, PSAP, CTSA, GLB1, DNASE2, LAMP1
	hsa04520	Adherens junction	3.35E-04	13	EGFR, WASF3, SMAD3, SMAD2, CTNNA1, TCF7L1, CSNK2A2, MAPK1, FYN, MAPK3
	hsa04115	p53 signaling pathway	4.26E-04	12	CCND1, TNFRSF10B, ZMAT3, SERPINE1, DDB2, TP53, FAS, PERP, THBS1, IGFBP3
	hsa05219	Bladder cancer	8.61E-04	9	EGFR, VEGFB, MAPK1, VEGFC, CCND1, MAP2K2, MAPK3, TP53, THBS1
	hsa04540	Gap junction	1.28E-03	13	ADCY3, EGFR, GNAI2, MAP2K2, GNA11, LPAR1, ITPR3, MAPK1, ADCY9, CSNK1D
	hsa00980	Metabolism of xenobiotics by cytochrome P450	2.44E-03	10	GSTM1, AKR1C3, GSTM2, AKR1C2, CYP1B1, ADH5, GSTT1, EPHX1, AKR1C1, ALDH3B1
	hsa05216	Thyroid cancer	2.48E-03	7	MAPK1, CCND1, MAP2K2, RXRA, MAPK3, TP53, TCF7L1
	hsa05212	Pancreatic cancer	2.57E-03	11	EGFR, VEGFB, MAPK1, VEGFC, CCND1, RELA, MAPK3, TP53, SMAD3, SMAD2
	hsa05220	Chronic myeloid leukemia	3.49E-03	11	MAPK1, CTBP1, CCND1, MAP2K2, RELA, MAPK3, TP53, SMAD3, BCL2L1, ABL1

[i] DEGs, differentially expressed genes; GO, Gene Ontology; Up, upregulated; Down, downregulated.

According to the results of the pathway enrichment analysis, the upregulated DEGs were predominantly enriched in seven pathways. In accordance with the GO term analysis, the DNA replication pathway, including RFC2, RFC3, RFC4 and RFC5, and cell cycle pathway, including CDK1, minichromosome maintenance complex component 3 (MCM3) and BUB1, were also enriched in the upregulated genes. The downregulated DEGs were predominantly enriched in the focal adhesion, including CAV1, collagen type VI α1 (COL6A1), thrombospondin 1 (THBS1) and EGFR, and p53 signalling pathways, including TP53, Fas cell surface death receptor (FAS) and TP53 apoptosis effector (PERP), as shown in Table IB.

Construction and analysis of the PPI network

The PPI network for the upregulated and downregulated DEGs consisted of 442 pairs of PPIs. The degrees of DEGs, including CDK1, MAD2L1, NDC80, non-SMC condensin I complex subunit G (NCAPG), BUB1, centromere protein F (CENPF) and kinesin family member 11 (KIF11), were >17 (Table II), indicating that they were important genes in OS.

Table II. Differentially expressed genes with a connectivity degree of ≥10 in the protein-protein interaction network.

Table II.

Differentially expressed genes with a connectivity degree of ≥10 in the protein-protein interaction network.

ID	Degree
CDK1	29
MAD2L1	23
NDC80	20
NCAPG	20
BUB1	19
CENPF	19
KIF11	18
DLGAP5	17
CREBBP	17
BIRC5	17
RFC4	16
RRM2	16
TP53	16
AURKA	16
SF3A2	14
ASPM	14
SNRPG	13
MAPK1	12
HMMR	11
NUP107	11
CDKN3	11
PPP1CC	11
SRSF1	11
RACGAP1	10
NUP160	10
ZWINT	10
SRSF3	10

Analysis of network modules

A total of 10 network modules were obtained from the PPI network using the default criteria, and the module with the highest score contained 16 nodes and 102 edges. In this module, CDK1 interacted with other DEGs, including MAD2L1, BUB1, NCAPG, NDC80 and CENPF (Fig. 2).

Figure 2. Modules with the highest significance in the protein-protein interaction network. The lines between any nodes represent the interrelations of those proteins.

The functional enrichment analysis of the module with the highest score showed that the majority of the DEGs in this module were predominantly associated with the cell cycle. Certain DEGs, including CDK1, MAD2L1, BUB1 and NDC80, were implicated in mitosis and the M phase of mitotic cell cycle; other DEGs, including Rac GTPase-activating protein 1 (RACGAP1) and MAD2L1, were correlated with cell cycle process (Table IIIA). CDK1, MAD2L1, BUB1 and aurora kinase A (AURKA) were significantly enriched in the oocyte meiosis pathway (Table IIIB).

Table III. GO terms and pathways enriched for the DEGs in the most significant module.

Table III.

GO terms and pathways enriched for the DEGs in the most significant module.

A, Top 10 GO terms enriched for the DEGs in the most significant module
Term	Description	P-value	n	Genes
GO:0000280	Nuclear division	9.92E-16	11	CDK1, MAD2L1, KIF11, NCAPG, DLGAP5, BUB1, CENPF, BIRC5, NDC80, AURKA
GO:0007067	Mitosis	9.92E-16	11	CDK1, MAD2L1, KIF11, NCAPG, DLGAP5, BUB1, CENPF, BIRC5, NDC80, AURKA
GO:0000087	M phase of mitotic cell cycle	1.19E-15	11	CDK1, MAD2L1, KIF11, NCAPG, DLGAP5, BUB1, CENPF, BIRC5, NDC80, AURKA
GO:0048285	Organelle fission	1.49E-15	11	CDK1, MAD2L1, KIF11, NCAPG, DLGAP5, BUB1, CENPF, BIRC5, NDC80, AURKA
GO:0000278	Mitotic cell cycle	1.86E-15	12	CDK1, MAD2L1, KIF11, NCAPG, DLGAP5, BUB1, CENPF, BIRC5, NDC80, AURKA
GO:0022402	Cell cycle process	2.11E-15	13	CDK1, KIF11, DLGAP5, CENPF, AURKA, NDC80, BIRC5, MAD2L1, NCAPG, BUB1
GO:0022403	Cell cycle phase	6.44E-15	12	CDK1, MAD2L1, KIF11, NCAPG, DLGAP5, BUB1, BIRC5, NDC80, AURKA, ASPM
GO:0000279	M phase	5.81E-14	11	CDK1, MAD2L1, KIF11, NCAPG, DLGAP5, BUB1, CENPF, BIRC5, NDC80, AURKA
GO:0007049	Cell cycle	9.57E-14	13	CDK1, KIF11, DLGAP5, CENPF, AURKA, NDC80, BIRC5, MAD2L1, NCAPG, BUB1
GO:0051301	Cell division	1.80E-12	10	CDK1, MAD2L1, KIF11, NCAPG, BUB1, CENPF, BIRC5, NDC80, RACGAP1, ASPM
B, Pathways enriched for the DEGs in the most significant module
Term	Description	P-value	n	Genes
hsa04114	Oocyte meiosis	1.88E-04	4	CDK1, MAD2L1, BUB1, AURKA
hsa04914	Progesterone-mediated oocyte maturation	4.06E-03	3	CDK1, MAD2L1, BUB1
hsa04110	Cell cycle	8.43E-03	3	CDK1, MAD2L1, BUB1

[i] DEGs, differentially expressed genes; GO, Gene Ontology.

Analysis of the miRNA-DEG regulatory network

The miRNA-DEG regulatory network contained 63 miRNAs, nine of their corresponding DEGs and 16 DEGs, which interacted with these nine DEGs. DNA polymerase ζ subunit 3 (POLE3) was regulated by 18 miRNAs, including miR-4310, miR-4680-3p, miR-583 and miR-4269; RFC3 was regulated by 16 miRNAs, including miR-802 and miR-649; RFC3 and RFC5 were modulated by miR-522-3p and miR-224-3p. In addition, RFC2, RFC3, RFC4 and RFC5 interacted with each other (Fig. 3).

Figure 3. Regulatory network containing microRNAs and their corresponding differentially expressed genes enriched in the DNA replication pathway. Large circular nodes represent upregulated genes enriched in the DNA replication pathway; small circular nodes represent genes interacting with the genes enriched in the DNA replication pathway; rectangular nodes represent microRNAs. miR, microRNA.

Discussion

In the present study, 317 DEGs were found to be significantly upregulated and 670 were significantly downregulated in the conventional OS samples. A majority of the DEGs were associated with cell cycle. According to the miRNA-DEG regulatory network for the DEGs enriched in DNA replication, RFC2, RFC3, RFC4 and RFC5 were found to interact with each other.

RFC2, RFC3, RFC4 and RFC5 encode members of RFC family, also termed activator 1. These DEGs were enriched in DNA replication, which agreed with the results of previous studies (23,24). DNA replication is an essential event in tumour growth (25). The deregulation of protein complexes involved in the initiation of DNA replication can lead to cancer (26). Several DEGs in the network module, including CDK1, MAD2L1, NDC80 and BUB1, which had higher degrees in the PPI network, were found to interact with RFC4. These four DEGs were predominantly enriched in cell mitosis and cell cycle. Alterations in cell cycle regulation occur in several types of cancer, including OS (27). Cyclin-dependent kinase 1 (CDK1) is an important G2/M checkpoint protein (28), and its inhibitor, SCH 727965 (dinacliclib) can trigger the apoptosis of U-2 OS cells (29). MAD2L1, BUB1 and NDC80 are involved in the spindle checkpoint pathway (30,31). MAD2 has been reported to be commonly overexpressed in human conventional OS (32), and BUB1 has been found to be ectopically expressed in SAOS and U-2 OS cell lines (33). In addition, CDK1, MAD2L1 and BUB1 have been found to be significantly enriched in the pathway of oocyte meiosis, which was found to be markedly altered in high-grade OS cell lines when compared with osteoblasts (34). RFC3 was also modulated by a cluster of miRNAs, including miR-802. The expression of miR-802 has been reported to be upregulated in OS tissues, and to promote cell proliferation by targeting p27 in U-2 OS and MG-63 cells (35). RFC3 and RFC5 are also modulated by miR-224-3p and miR-522-3p. There is no previous evidence indicating that miR-224-3p and miR-522-3p are associated with conventional OS. Therefore, miR-224-3p and miR-522-3p are predicted to be novel biomarkers in conventional OS. Therefore, RFC2-5, together with certain DEGs, including CDK1, MAD2L1, NDC80 and BUB1, and a series of miRNAs, including miR-802, miR-224-3p and miR-522-3p, may be responsible for the initiation and development of conventional OS.

In conclusion, the present study found that the majority of DEGs, including CDK1, MAD2L1, NDC80 and BUB1, were associated with the cell cycle. Other DEGs, including RFC2, RFC3, RFC4 and RFC5, were associated with DNA replication. These, in addition to a number of miRNAs, including miR-802, miR-224-3p and miR-522-3p, may be essential in the pathogenesis of conventional OS, providing novel information to assist in the clinical diagnosis of this disease. However due to limitations in the present study, additional experiments are required to shed light on the molecular mechanisms involved in this life-threatening disease.

References