Integrated Analysis of Competing Endogenous RNAs Network Reveals Potential Signatures in Osteosarcoma Development.

The purpose of this work was to extract key players such as mRNAs and long non-coding RNA (lncRNAs) in the etiopathogenesis of osteosarcoma (OS). The sequencing analyses (mRNAs and lncRNAs) of OS were conducted followed by differentially expressed mRNAs and lncRNAs (DEmRNAs and DElncRNAs) identification between U-2OS cells with has-miR-590-5p overexpression and negative control cells. Following this, the co-expression and functional enrichment analyses of DEmRNAs and DElncRNAs were carried out. Also, the miRNAs-DElncRNAs-DEmRNAs regulatory network was constructed with DElncRNAs-miRNAs and DElncRNAs-DEmRNAs pairs after the target gene analysis of miRNA. In addition, the ceRNA-has-miR-590-5p was further extracted based on the has-miR-590-5p-DElncRNAs and DElncRNAs-DEmRNAs interactions. Finally, the results of the bioinformatics analysis was verified by reverse-transcription polymerase chain reaction (RT-PCR). Totally, 980 DEmRNAs (539 up-regulated DEmRNAs and 441 down-regulated DEmRNAs) and 682 DElncRNAs (352 up-regulated DElncRNAs and 330 down-regulated DElncRNAs) were extracted between cells with hsa-miR-590-5p overexpression and normal cells. The functional analyses suggested that up-regulated genes were significantly enriched in several GO terms such as signal transduction and cytokine-cytokine receptor interaction pathway while down-regulated genes (SCUBE3, HIST1H4E and EDIL3) were associated with calcium ion binding, cell surface function and nucleosome assembly. Additionally, the miRNAs-DEmRNAs-DEmRNAs network represented 220 pairs among 41 miRNAs, 38 DElncRNAs and 61 DEmRNAs. Furthermore, the ceRNA-hsa-miR-590-5p network consisted of 70 interaction pairs including hsa-miR-590-5p-SCUBE3-CTB-113D17.1, hsa-miR-590-5p-EDIL3-CTB-113D17.1 and hsa-miR-590-5p-HIST1H4E-CTB-113D17.1) among hsa-miR-590-5p, 30 DEmRNAs and 4 down-regulated DElncRNAs. Meanwhile, the RT-PCR results incidated that compared with the blank (KB) and negative control (NC) group, the mRNA expression of SCUBE3, HIST1H4E, and EDIL3 were significantly descreased in mimics group (P value <0.05). The lncRNA CTB-113D17.1 might implicate with OS development probably via serving as a hsa-miR-590-5p sponge to regulate gene targets (SCUBE3, EDIL3 and HIST1H4E), which will facilitate the deep understandings of OS progression.

Introduction

Osteosarcoma (OS) is as one of aggressive bone malignant tumor and primarily influences children and adolescents.[1] Existing investigations have reported that OS had an increasing incidence and a higher recurrence as well as metastasis.[2] More specifically, nearly 15% to 20% of patients with OS have recognized metastases in clinical practice.[3] In recent years, numerous studies have concentrated on exploring the feasible approaches such as surgery for OS management.[4] However, this did not produce favorable clinical outcomes and accompanied a higher recurrence. Encouragingly, increasing researchers have investigated possible molecular mechanisms of OS initiation and development to seek novel therapeutic targets to OS.

Notably, noncoding RNAs transcripts such as microRNAs (miRNAs) and long noncoding RNAs (lncRNAs) has been reported to play significant roles on molecular mechanism of cancers. The competing endogenous RNAs (ceRNAs) hypothesis suggests different RNA transcripts participate in pathological processes mainly via completely targeting the binding sites of shared miRNAs. Andersen et al examined the expression levels of miRNAs by analyzing the miRNA profiling data of patients suffering from OS.[5] Consequently, they found and verified 29 de-regulated miRNAs in OS tumorigenesis. Liao et al argued that lncRNA SNHG16 was closely associated with the development of OS probably through sponging miR-98-5p. Interestingly, many researchers have conducted the ceRNA network analysis using a bioinformatics method to decipher the underlying mechanisms of OS. For instance, Zhang et al suggested that the lower expression of 2 lncRNAs LINC00028 and LINC00323 exhibited a better survival outcome and higher levels of RAP1B, ATF2 and PPM1B were predominately linked with OS recurrence according to ceRNA regulatory analyses.[6] In addition, Wang et al recently also stated that HOTAIR as a ceRNA via sponging miR-217 modulated OS development by up-regulating ZEB1 level.[7] Although existing evidence has implied that several noncoding RNA transcripts possibly contributed to OS progression, the detailed molecular mechanisms of OS has not been fully understood.

Multiple research groups have reported that has-miR-590 was a novel signature in diagnosis and treatment of OS target. For instance, Cai et al indicated that miR 590 5p suppresses OS cell proliferation and invasion via targeting KLF5.[8] However, the precise regulatory mechanisms of miR-590 in the OS carcinogenic process have not been investigated. Therefore, we performed a bioinformatics analysis based on mRNAs and lncRNAs profile data from OS cell lines with has-miR-590-5p overexpression (OSM) and has-miR-590-5p negative control (NC) cell. The differentially expressed mRNAs (DEmRNAs) and differentially expressed lncRNAs (DElncRNAs) between OSM and NC cells were identified followed by functional enrichment analysis. Afterward, miR-590-5p-miRNA-mRNA regulatory network was constructed to screen the pivotal ceRNAs involved in OS progression, which will provide better understandings and promising therapeutic targets for OS.

Materials and Methods

Cell Culture and Transfection

U-2OS cells were cultured in RPMI-1640 medium (including 1% penicillin/streptomycin) with 10% fetal calf serum (GIBCO, Newyork, USA) at 37°C and kept in a humidified atmosphere of 5% CO2. Subsequently, the logarithmic growth phase cells were seeded in 6 well plate (1.5 × 105 / well). They were divided into 3 groups, such as miRNA mimics group, miRNA inhibitor group, and miRNA negative control. When cells grew into 60%, miR-590-5p vectors were transfected into U-2OS cells through Lipofactamine 2000 according to a standard procedure at room temperature.

Sequencing

The total RNA was extracted from cells using Trizol reagent following the manufacturer’s instructions (Invitrogen, USA). After the quality control of RNA, mRNA was enriched by Dynabeads® oligo dT magnetic beads (Thermo Fisher Scientific, USA), and then broke into shot fragments by fragmentation buffer (Agilent Technologies, California, USA). Afterward, the RNA fragments were reverse transcribed into the first strand cDNA with random hexamers. The second strand cDNA was compounded by adding into buffer, dNTPs, RNase H and DNA polymerase I. The final cDNA library was constructed after double strands cDNA were purified and repaired. The concentration of cDNAs in the library was attenuated into 1 ng/µL with a Qubit 2.0 fluorometer, and then cDNAs were detected using the Agilent Bioanalyzer 2100 (Agilent Technologies, California, USA). The libraries were pooled according to the data size and effective cDNA concentration. Finally, the cDNA libraries were sequenced on an Illumina HiSeq™ 3500.

Data Acquisition and Pre-Processing

To decipher the molecular mechanisms of miR-590-5p in the initiation and development of osteosarcoma, the mRNA and lncRNA sequencing were conducted respectively and corresponding data was also obtained, which was derived from 3 control osteosarcoma cells and 3 miR-590-5p-overexpressed cell lines. Then raw reads were screened by filtering unreliable reads to generate clean reads as follows: 1) adaptors; 2) reads without inserted fragment; 3) low quality bases (QC < 20) and the whole sequence involving bases with QC < 10; 4) reads with N ratio > 10%; 5) sequence with the length < 20 bp. After that, the quality control of clean reads was performed by fastqc[9] software (v 0.11.5, http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and then aligned to human reference genome (Hg38) using TopHat[10] (v 2.1.0, http://ccb.jhu.edu/software/tophat/) according to default parameters. The featureCounts[11] (v 1.6.0, http://subread.sourceforge.net/) was employed to annotate genome sequence using human genome annotation file (gencode) as a reference, and subsequently, lncRNA and mRNA data were determined. Finally, the Pearson’s correlation coefficient (P) between 2 samples was measured with the cor function (https://stat.ethz.ch/R-manual/R-devel/library/stats/html/cor.html) in R 3.5.2 and the P value approximates to 1, the expression pattern between samples is more similar.

Identification of Differentially Expressed mRNAs and lncRNAs

The normalization of 6 samples read count data (mRNA and lncRNA) was processed with trimmed mean of M values (TMM) method provided by edgeR[12,13] software in R. Afterward, the DEmRNAs and DElncRNAs between control (NC) and miR-590-5p overexpression (X590) group were identified by quasi-likelihood (QL) F-test in edgeR with the cutoffs of | log2 fold change (FC) | > 1 and P value < 0.05. The heatmaps of DEmRNAs and DElncRNAs were constructed respectively using pheatmap (version: 1.0.10, https://cran.r-project.org/web/packages/pheatmap/index.html) software in R.

Co-Expression Analysis of DEmRNAs and DElncRNAs

The co-expression pairs of DEmRNA and the DElncRNA were established on the basis of the pairwise Pearson correlation coefficient method[14] with the thresholds of | r | > 0.9 and P value < 0.05. Then the P value in the lncRNA and mRNA co-expression analysis was corrected by Benjamini & Hochberg (1995) (“BH” or its alias “FDR”) with the threshold of FDR <0.05 and | r | >0.99. Of these, the positive correlation between mRNA and lncRNA was determined with r > 0.9 while r < −0.9 was considered as positive correlation lncRNA-mRNA pairs.

Gene Ontology (GO) Functional and Pathway Enrichment Analyses

The GO and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted by DAVID[15] (v 6.8, https://david.ncifcrf.gov/home.jsp) to further investigate the functional characteristics of DEmRNAs. The significance cutoff of hyper-geometric test (P values < 0.05) and count ≥2 were regarded as the thresholds for significant enrichment. In addition, KEGG pathway analysis of positively correlated DElncRNAs in co-expression network was also carried out by clusterprofiler[16] using the same criteria mentioned above.

Prediction of miRNAs Targets

Admittedly, miRNAs exert crucial roles in multiple biological processes primarily via binding to the 3’ un-translated region (UTRs) of mRNAs or lncRNAs. Here, we utilized miRanda[17] to predict the binding sites between hsa-miR-590-5p and 3’ UTRs of DEmRNAs and DElncRNAs based on -sc140 and -strict parameters. Notably, the cell lines with hsa-miR-590-5p overexpression were used in this investigation, and therefore, we concentrated on down-regulated DEmRNAs and DElncRNA in co-expression network.

CeRNA Network Construction

The interactive correlation among mRNAs, lncRNAs and their miRNA targets can be determined and evaluated by ceRNA network. Herein, DElncRNAs-miRNA and DElncRNAs-DEmRNA pairs were screened and used for ceRNA network construction. Then the KEGG pathway analysis was performed on the lncRNAs obtained in ceRNA network by the clusterprofiler package of R language, and the P value was corrected by BH with the threshold of P value < 0.05. Moreover, DElncRNAs-miRNA interactions related to hsa-miR-590-5p were further extracted. Afterward, DElncRNAs-hsa-miR-590-5p-DEmRNA network construction was carried out by combining with DElncRNAs- hsa-miR-590-5p and positive co-expression DElncRNAs-DEmRNA pairs. Finally, the ceRNA network was visualized with cytoscape[18] software. Here, the Counts per Million (CPM) was represented the absolute expressed values of the gene, and the logFC was represented the differently expression multiple.

Verification

The reverse-transcription polymerase chain reaction (RT-PCR) was used to verify the expression levels of key mRNAs. Total RNA of 5-10 × 106 cell samples (U-2OS cells with has-miR-590-5p overexpression and negative control cells) were extracted using TRIzol (9109, TaKaRa) and the extracted RNA quality was determined with a microplate reader (Infinite M100 PRO, TECAN). cDNA was synthesized through reverse transcription using 5x primeScript RT Master MIX (perfect Real Time) (RR036A, TaKaRa). The groups were as follows: blank (KB) group, negative control (NC) group and mimics group. The RT-PCR amplification was executed using SYBR Green PCR master mix (A25742, Thermo) and fluorescence quantitative PCR instrument (7900HTFAST, ABI). All samples were normalized to the corresponding expression of internal control GAPDH. and the relative expression level was calculated through the 2-ΔΔCt analysis method.

Statistical Analysis

All data were statistically analyzed using the Graphpad Prism 5 software (Graphpad Software, San Diego, CA) and presented as the means ± standard deviation. Statistical significance of the differences among the various experimental groups were determined by 1-way ANOVA and Newman-Keuls Multiple Comparison Test. P value < 0.05 was considered statistically significant.

Results

Correlation Analysis Between Samples and Identification of DEmRNAs and DElncRNAs

After data mining and annotation, the pairwise correlation analysis between samples was conducted. Our results showed that the Pearson’s correlation coefficient was more than 0.98, suggesting the rationality of sample selection and reliability of this investigation (Figure 1). Furthermore, a total of 980 DEmRNAs and 682 DElncRNAs were extracted with | log2 FC| >1 and P value < 0.05 according to the method mentioned above, which was comprised of 539 up-regulated genes, 411 down-regulated genes, 352 up-regulated lncRNAs and 330 down-regulated lncRNAs. The heatmaps of DEmRNAs and DElncRNAs were displayed in Figure 2 A and B.

Figure 1.

The heatmap of correlation analysis between samples based on mRNA expression level. NC_1: sample 1 in control group; NC_2: sample 2 in control group; NC_3: sample 3 in control group; X-590_1: sample 1 in tumor group with hsa-miR-590-5p overexpression; X-590_2: sample 2 in tumor group with hsa-miR-590-5p overexpression; X-590_3: sample 3 in tumor group with hsa-miR-590-5p overexpression.

Figure 2.

(DPI: 300) Identification and functional analyses of differentially expressed mRNAs and lncRNAs (DEmRNAs and DElncRNAs). (A) the heatmap of DEmRNAs; (B) the heatmap of DElncRNAs; (C) the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of up-regulated DEmRNAs; (D) GO and KEGG enrichment analyses of down-regulated DEmRNAs. The horizontal axis represents the number of DEmRNAs and the vertical axis showed names of GO terms and KEGG enrichment pathways.

Co-Expression Analysis

After the P value was corrected by BH with the threshold of FDR <0.05 and | r | >0.99, the co-expression correlation between DElncRNAs and DEmRNAs was investigated and the findings revealed that there were 4691 positively co-expressed DElncRNA-DEmRNA pairs, containing 337 DElncRNAs and 634 DEmRNAs (Supplementary file 1).

Functional Analyses of DEmRNAs and DElncRNAs

The GO enrichment analysis of up-regulated DEmRNAs revealed that the significantly enriched GO-biological process (BP) terms included signal transduction and cell adhesion while the extracellular region and proteinaceous extracellular matrix were 2 dramatically enriched GO-cellular components (CC) terms. Furthermore, we also found that these genes played essential roles in actin binding and actin filament binding according to the GO-molecular function (MF) analysis. Meanwhile, the KEGG pathway analysis suggested that up-regulated genes were responsible for cytokine-cytokine receptor interaction and platelet activation pathway (Figure 2C and Table 1). For down-regulated genes, we noted that these DEmRNAs such as SCUBE3 (Signal Peptide, CUB Domain and EGF Like Domain Containing 3) and EDIL3 (EGF Like Repeats And Discoidin Domains 3) were primarily associated with calcium ion binding in GO-MF term and cell surface function in GO-CC term. Simultaneously, several genes including HIST1H4E (Histone Cluster 1 H4 Family Member E) and AKAP5 (A-Kinase Anchoring Protein 5) were enriched in multiple GO-BP terms such as nucleosome assembly and chemical synaptic transmission. Additionally, the most enriched KEGG pathways were pathways in cancer and alcoholism (Figure 2 D and Table 1).

Table 1.

Top 5 Gene Ontology (GO) Terms and Pathways of Differentially Expressed mRNAs.

Category	Term	Count	Gene
GOTERM_BP	GO:0007268∼chemical synaptic transmission	17	KCNMB4, SYT10, OPRL1, KCNA1, GRIN2A, VIPR1, GPR1, GABRR2, SLC1A3, SYPL1, SLC1A6, CACNA1G, AKAP5, SLC18A3, NPFFR1, KCNQ2, PTGDR2
GOTERM_BP	GO:0017158∼regulation of calcium ion-dependent exocytosis	6	SYT10, SYT2, SYT11, SYT8, SYTL5, TRPV6
GOTERM_BP	GO:0042475∼odontogenesis of dentin-containing tooth	7	BMP4, SOSTDC1, JAG2, LEF1, CA2, PITX2, FGF4
GOTERM_BP	GO:0006334∼nucleosome assembly	10	HIST1H1E, HIST1H2BC, HMGB2, H2BFS, HIST1H2BK, ANP32E, HIST1H3B, HIST1H4E, HIST1H4I, HIST1H3H
GOTERM_BP	GO:0045668∼negative regulation of osteoblast differentiation	6	HOXA2, ID2, SFRP1, GDF10, PTCH1, CITED1
GOTERM_CC	GO:0098793∼presynapse	8	SYNDIG1, SYT10, SYT11, SYNJ1, SYT8, SYTL5, WNT7A, CHAT
GOTERM_CC	GO:0032809∼neuronal cell body membrane	5	UNC5A, KCNA2, CX3CR1, SLC4A8, KCNE3
GOTERM_CC	GO:0060077∼inhibitory synapse	4	SLC32A1, SYT11, NPTN, IQSEC3
GOTERM_CC	GO:0000786∼nucleosome	8	HIST1H1E, HIST1H2BC, HIST1H2BK, HIST1H2AG, HIST1H3B, HIST1H4E, HIST1H4I, HIST1H3H
GOTERM_CC	GO:0009986∼cell surface	23	HAVCR2, PPFIA2, ARSB, LPL, SCUBE3, RTN4RL1, ADAMTS15, KCNA1, GRIN2A, IQGAP2, ACKR3, ALPP, ADGRG2, SDC2, FZD6, AMBP, SLC1A3, SFRP1, ITGA7, NPTN, FGFBP1, WNT7A, GHR
GOTERM_MF	GO:0005509∼calcium ion binding	33	SYT2, NELL1, FAM20C, NELL2, JAG2, SYT8, CABP7, PCDHGC4, MYL10, EDIL3, DCHS1, CDH6, PCDHAC1, EFHD1, ANXA8, VWA2, IHH, MATN2, TCHH, MATN3, SCUBE3, SYT10, SYT11, CAPSL, NID2, CLGN, ANXA10, BNIP2, ANXA8L1, SYTL5, DSC2, NCAN, ADGRL3
GOTERM_MF	GO:0005544∼calcium-dependent phospholipid binding	8	ANXA8, ANXA10, SYT10, ANXA8L1, SYT2, SYT11, SYT8, SYTL5
GOTERM_MF	GO:0008201∼heparin binding	11	BMP4, LPL, ADAMTS8, SFRP1, FGF9, RSPO3, ADAMTS15, PTCH1, FGFBP1, PCOLCE2, FGF4
GOTERM_MF	GO:0005201∼extracellular matrix structural constituent	6	MATN3, HAPLN1, LUM, COL3A1, COL1A1, NCAN
GOTERM_MF	GO:0030276∼clathrin binding	5	SYT10, SYT2, SYT11, SYT8, SYTL5
KEGG_PATHWAY	hsa05034: Alcoholism	12	DRD1, HIST1H2BC, HIST1H2BK, HIST1H2AG, HIST1H3B, GNG13, HIST1H4E, GRIN2A, GNG3, GNG4, HIST1H4I, HIST1H3H
KEGG_PATHWAY	hsa05217: Basal cell carcinoma	6	BMP4, LEF1, PTCH1, WNT7A, FZD6, GLI1
KEGG_PATHWAY	hsa04725: Cholinergic synapse	8	GNG13, SLC18A3, GNG3, GNG4, KCNQ2, PIK3R3, CHRNA3, CHAT
KEGG_PATHWAY	hsa03320: PPAR signaling pathway	6	LPL, ACOX1, OLR1, ACSBG2, ACAA1, ACSBG1
KEGG_PATHWAY	hsa05200: Pathways in cancer	17	BMP4, E2F2, FGF9, SPI1, BIRC7, GNG13, LEF1, RAD51, FZD6, GLI1, CCNE2, PTCH1, GNG3, GNG4, PIK3R3, WNT7A, FGF4

KEGG: Kyoto Encyclopedia of Genes and Genomes; BP: biological process; CC: cellular components; MF: molecular function.

The KEGG enrichment analysis of DElncRNAs in co-expression network was also conducted and the results revealed that these DElncRNAs were mainly correlated with osteoclast differentiation, rheumatoid arthritis, and proteoglycans in cancer (Figure 3). DElncRNAs including AC003104, AC147651 and CTB-193M12 were critical players implicated with these biological processes.

Figure 3.

Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of DElncRNAs in co-expression network. The horizontal axis indicates the name of DElncRNAs and the number in brackets represents the number of annotated genes in KEGG database. The vertical axis was names of KEGG pathways. GeneRatio represents the ratio of the number of co-expressed genes and the number of annotated genes in each KEGG pathway. The different color showed the different p value and the color ranging from blue to red indicates the increasing significance.

Construction of ceRNA Network and Sub-Network Extraction

According to the prediction analyses of miRNA targets, miRNAs-DElncRNAs pairs were identified and the integrative ceRNA network, containing 220 DElncRNA-miRNA-DEmRNA interactions among 41 miRNAs, 38 DElncRNAs and 61 DEmRNAs, was constructed as showed in Figure 4A. The top 10 nodes with the higher degree were respectively exhibited in Table 2. In addition, there were 10 lncRNAs involved in the KEGG pathway (Table 3). Besides, the sub-ceRNA network was further identified and visualized based on the down-regulated DElncRNAs-hsa-miR-590-5p pairs and the positive co-expression analysis of these DElncRNAs, which represented 70 interactions such as hsa-miR-590-5p-SCUBE3, hsa-miR-590-5p-EDIL3-CTB-113D17.1 and hsa-miR-590-5p-HIST1H4E pairs among hsa-miR-590-5p, 30 DEmRNAs and 4 down-regulated DElncRNAs (Figure 4B; Table 4). In addition, the results shown that SCUBE3 (CPM: the average of X590 group was 4.513, and the average of NC group was 9.881, logFC= -1.336); EDIL3 (CPM: the average of X590 group was 53.520, and the average of NC group was 183, logFC= -1.793); HIST1H4E (CPM: the average of X590 group was 1.260, and the average of NC group was 1.604, logFC= -1.20679) (Supplementary file 2).

Figure 4.

The ceRNA regulatory network and ceRNA sub-network. (A) miRNAs-lncRNAs-mRNAs (ceRNA) network: the orange nodes show up-regulated genes, green nodes represent down-regulated genes, rhomboid nodes indicate lncRNAs, and circle nodes represent mRNAs. The dotted line represents lncRNAs-mRNAs interactions and solid arrow shows the miRNA-mRNA pairs and miRNA-lncRNA interactions. (B) the hsa-miR-590-5p-lncRNAs-mRNAs network: the green nodes show down-regulated genes, arrow nodes indicate hsa-miR-590, and circle nodes represent mRNAs. The dotted line represents lncRNAs-mRNAs interactions and solid arrow shows the hsa-miR-590-5p -mRNA regulatory pairs and hsa-miR-590-5p–lncRNA interactions.

Table 2.

Top 10 miRNAs, mRNAs and lncRNAs in ceRNA Regulatory Network.

Description	Name	Degree	Description	Name	Degree	Description	Name	Degree
mRNA_down	AURKA	6	lncRNA_down	NORAD	13	miRNA	hsa-miR-619-5p	17
mRNA_up	SCN4A	6	lncRNA_up	RP11-284F21.10	13	miRNA	hsa-miR-8485	10
mRNA_up	KSR2	5	lncRNA_up	RP11-301G19.1	9	miRNA	hsa-miR-6089	9
mRNA_down	TCF4	4	lncRNA_up	CTD-3088G3.6	8	miRNA	hsa-miR-1285-3p	5
mRNA_up	LPCAT1	4	lncRNA_down	RP11-404E16.1	7	miRNA	hsa-miR-3620-5p	5
mRNA_down	LEPROT	4	lncRNA_down	USP2-AS1	6	miRNA	hsa-miR-3960	5
mRNA_up	PAPLN	4	lncRNA_up	RP11-488P3.1	5	miRNA	hsa-miR-5189-5p	5
mRNA_up	TTLL11	4	lncRNA_down	AC008982.2	4	miRNA	hsa-miR-10400-5p	5
mRNA_up	MYO16	4	lncRNA_up	RP11-474B16.1	4	miRNA	hsa-miR-877-3p	4
mRNA_down	ADAMTS15	4	lncRNA_down	RP11-28G8.1	4	miRNA	hsa-miR-4739	4

Table 3.

The KEGG Pathway Analysis of lncRNAs in the ceRNA Network.

Cluster	Description	pvalue	p.adjust	Count	Gene
AC147651.4	Rheumatoid arthritis	2.86E-05	0.001314	4	MMP1/ATP6V0A1/MMP3/CCL3
RP11-488P3.1	Rheumatoid arthritis	0.000131	0.009047	4	MMP1/CCL2/MMP3/ATP6V0A1
RP11-488P3.1	IL-17 signaling pathway	0.000136	0.009047	4	MMP1/MAPK13/CCL2/MMP3
RP11-403A3.1	Rheumatoid arthritis	4.81E-05	0.004001	5	CCL2/MMP3/CCL3/MMP1/ATP6V0A1
RP11-403A3.1	IL-17 signaling pathway	5.06E-05	0.004001	5	CCL2/MMP3/MMP1/MAPK13/FOSB
RP11-403A3.1	Human cytomegalovirus infection	0.00039	0.020531	6	PLCB2/CCL2/CCL3/IL6R/MAPK13/NFATC2
RP11-403A3.1	Chagas disease (American trypanosomiasis)	0.000981	0.038757	4	PLCB2/CCL2/CCL3/MAPK13
RP11-403A3.1	TNF signaling pathway	0.00139	0.043927	4	CCL2/CASP7/MMP3/MAPK13
RP11-403A3.1	Yersinia infection	0.001793	0.047225	4	NLRP3/CCL2/MAPK13/NFATC2
THAP7-AS1	Huntington disease	0.025003	0.025003	1	DNAH2
RP11-214K3.24	Glucagon signaling pathway	0.013318	0.026637	1	GCGR
RP11-214K3.24	Neuroactive ligand-receptor interaction	0.042719	0.042719	1	GCGR
RP11-150O12.6	Galactose metabolism	0.00079	0.028794	2	PGM1/GAA
RP11-150O12.6	Starch and sucrose metabolism	0.001066	0.028794	2	PGM1/GAA
RP3-412A9.16	Primary immunodeficiency	0.009527	0.04457	1	BLNK
RP3-412A9.16	B cell receptor signaling pathway	0.020501	0.04457	1	BLNK
RP3-412A9.16	TGF-beta signaling pathway	0.023483	0.04457	1	AMHR2
RP3-412A9.16	NF-kappa B signaling pathway	0.025469	0.04457	1	BLNK
RP3-412A9.16	Osteoclast differentiation	0.031908	0.044671	1	BLNK
RP11-531A24.5	Hippo signaling pathway	0.019349	0.037542	1	PARD6B
RP11-531A24.5	Tight junction	0.021234	0.037542	1	PARD6B
RP11-531A24.5	Axon guidance	0.022742	0.037542	1	PARD6B
RP11-531A24.5	Rap1 signaling pathway	0.026385	0.037542	1	PARD6B
RP11-531A24.5	Endocytosis	0.031285	0.037542	1	PARD6B
RP11-531A24.5	Human papillomavirus infection	0.041462	0.041462	1	PARD6B
RP11-95M15.3	Glycosphingolipid biosynthesis—lacto and neolacto series	0.003392	0.003392	1	GCNT2
CH17-353B19.1	Choline metabolism in cancer	0.012313	0.012313	1	SLC44A5

Table 4.

The Interaction Pairs in ceRNA-hsa-miR-590-5p Regulatory Network.

MiRNAs/LncRNA	LncRNAs/mRNAs	mRNAs	LncRNAs
hsa-miR-590-5p	LncRNA LVCAT1	NCEH1	LVCAT1
hsa-miR-590-5p	LncRNA AP002954.3	VGLL2	LVCAT1
hsa-miR-590-5p	LncRNA CTB-113D17.1	BARX2	LVCAT1
hsa-miR-590-5p	LncRNA RP11-817I4.2	NPFFR1	LVCAT1
hsa-miR-590-5p	C8orf33	NCEH1	AP002954.3
hsa-miR-590-5p	NCEH1	VGLL2	AP002954.3
hsa-miR-590-5p	VGLL2	EDIL3	CTB-113D17.1
hsa-miR-590-5p	BARX2	C8orf33	CTB-113D17.1
hsa-miR-590-5p	NPFFR1	DRD1	CTB-113D17.1
hsa-miR-590-5p	EDIL3	NEK7	CTB-113D17.1
hsa-miR-590-5p	DRD1	DGKK	CTB-113D17.1
hsa-miR-590-5p	NEK7	VPS36	CTB-113D17.1
hsa-miR-590-5p	DGKK	AKAP5	CTB-113D17.1
hsa-miR-590-5p	VPS36	RP2	CTB-113D17.1
hsa-miR-590-5p	AKAP5	HIPK3	CTB-113D17.1
hsa-miR-590-5p	RP2	CACNA1G	CTB-113D17.1
hsa-miR-590-5p	HIPK3	SHISA9	CTB-113D17.1
hsa-miR-590-5p	CACNA1G	KCNA1	CTB-113D17.1
hsa-miR-590-5p	SHISA9	RMND5A	CTB-113D17.1
hsa-miR-590-5p	KCNA1	PWP1	CTB-113D17.1
hsa-miR-590-5p	RMND5A	PCGF5	CTB-113D17.1
hsa-miR-590-5p	PWP1	SCUBE3	CTB-113D17.1
hsa-miR-590-5p	PCGF5	MBTPS2	CTB-113D17.1
hsa-miR-590-5p	SCUBE3	NCEH1	CTB-113D17.1
hsa-miR-590-5p	MBTPS2	GPR180	CTB-113D17.1
hsa-miR-590-5p	GPR180	DNAJB9	CTB-113D17.1
hsa-miR-590-5p	DNAJB9	DSC2	CTB-113D17.1
hsa-miR-590-5p	DSC2	VGLL2	CTB-113D17.1
hsa-miR-590-5p	ARRDC3	ARRDC3	CTB-113D17.1
hsa-miR-590-5p	MAP3K1	MAP3K1	CTB-113D17.1
hsa-miR-590-5p	PPP1R3D	PPP1R3D	CTB-113D17.1
hsa-miR-590-5p	RSAD2	NPFFR1	CTB-113D17.1
hsa-miR-590-5p	SRD5A2	RSAD2	CTB-113D17.1
hsa-miR-590-5p	HIST1H4E	SRD5A2	CTB-113D17.1
LncRNA LVCAT1	C8orf33	HIST1H4E	RP11-817I4.2

The primer sequence of RT-PCR was shown in Table 5. As illustrated in Figure 5, compared with the KB and NC group, the mRNA expression of SCUBE3, HIST1H4E, and EDIL3 were significantly descreased in mimics group (P value < 0.05), which was consistent with the result of the above bioinformatics analysis.

Table 5.

Primers Used for RT-PCR.

Gene	Primer-sequences
SCUBE3	Forward, 5’-ATGCCGTCCTGGCTTTGAG-3’
	Reserve, 5’-CCGTCGGTATGGAGCACAAA-3’
HIST1H4E	Forward, 5’-GCAAAGGCGGAAAGGGACT-3’
	Reserve, 5’-CGGATGGCAGGCTTGGTAAT-3’
EDIL3	Forward, 5’-GATGGAAAGACTTGGGCAATG-3’
	Reserve, 5’-GGCTCAGAACAACCCGACAG-3’
GAPDH	Forward, 5’-TGACAACTTTGGTATCGTGGAAGG-3’
	Reserve, 5’-AGGCAGGGATGATGTTCTGGAGAG-3’

Figure 5.

The expression of SCUBE3, HIST1H4E and EDIL3 by RT-PCR. *represents the mimics group compare with KB group, # represents the mimics group compare with NC group. P value < 0.05 was defined as a significant difference. **P value < 0.01, *P value < 0.05, #P < 0.05. KB represents blank group, NC represents negative control group.

Discussion

OS is reported to be a malignant bone tumor of mesenchymal origin and characterized by a high mortality and poor prognosis.[19] Although increasing researchers have focused on dissecting and elaborating the underlying molecular etiology of OS over few years, precise mechanisms of OS initiation and development are still poorly understood. Notably, ceRNA hypothesis, which claimed that lncRNAs can exert regulator roles via competing binding sites of shared miRNAs with mRNAs, has been frequently applied in the illuminating the pathogenic mechanisms of OS.[20] In the current research, a comprehensive bioinformatics analysis was performed on the basis of the RNA sequencing data to screen pivotal modulators involved in the progression of OS. Our findings suggested that 980 DEmRNAs (539 up-regulated DEmRNAs and 441 down-regulated DEmRNAs) and 682 DElncRNAs (352 up-regulated DElncRNAs and 330 down-regulated DElncRNAs) were identified between cell lines with hsa-miR-590-5p overexpression and NC cell groups. Furthermore, the co-expression analysis revealed that there were 4691 positively expressed DEmRNAs-DElncRNAs pairs. The functional analyses indicated that up-regulated genes were significantly enriched in several GO terms (signal transduction and actin binding) and pathways such as cytokine-cytokine receptor interaction while down-regulated genes such as SCUBE3, HIST1H4E and EDIL3 were predominated with calcium ion binding, cell surface function and nucleosome assembly. Additionally, the ceRNA network construction was implemented via integrating screened DElncRNAs-DEmRNAs pairs with DElncRNAs-miRNAs interactions, which represented 220 pairs including 41 miRNAs, 38 DElncRNAs and 61 DEmRNAs. Furthermore, the ceRNA- hsa-miR-590-5p network was also constructed and contained 70 interaction pairs (such as hsa-miR-590-5p-SCUBE3-CTB-113D17.1, hsa-miR-590-5p-EDIL3-CTB-113D17.1 and hsa-miR-590-5p-HIST1H4E-CTB-113D17.1) among hsa-miR-590-5p, 30 DEmRNAs and 4 down-regulated DElncRNAs.

SCUBE3, a member of vertebrate SCUBE family, was initially identified from human vascular endothelial cells and predominately enriched in primary osteoblasts.[21] Existing evidence has implied that this gene can encode secreted cell surface-correlated membrane glycoproteins and composes 5 motifs including 9 epidermal growth factor-like (EFG) repeat domains, a N-terminal signal peptide sequence, a spacer region with several glycosylation sites, 3 duplicated cysteine domains and a C-terminal CUB (complement components of C1r/C1 s, EGF-associated sea urchin protein and bone morphogenetic protein 1) domain.[22] Recently, multiple investigators emphasized that SCUBE played essential roles in a wide variety of cancers. Wu et al demonstrated that SCUBE was highly expressed in lung cancer cells and acted as an endogenous transforming growth factor (TGF)-beta receptor ligand to regulate lung cancer metastasis.[23] Moreover, a previous study conducted the survival analysis for those patients with lung cancer and found that the over-expression of this gene exhibited a worse prognosis.[24] More notably, Liang et al found that SCUBE3 had a higher expression level in OS cells compared with those normal osteoblasts and SCUBE3 suppression could significantly retarded the proliferation of OS cells,[25] suggesting that SCUBE3 was also probably associated with OS occurrence and progression. Herein, our functional analysis showed that SCUBE3 participated in biological processes in cell surface, especially in calcium ion pathway. Numerous studies have reported that calcium pathway severed as a key regulatory mechanism in human carcinogenesis,[26] which further implied that the molecular interactions between SCUBE3 and calcium ions might be a novel therapeutic target for OS treatment. Here, lncRNA CTB-113D17.1 might function as a ceRNA in regulating SCUBE3 expression of OS through competitively binding to hsa-miR-590-5p. Overwhelming evidence has unraveled that lncRNA severed as functional modulators in tumor development.[27] Taken together, we could speculate that lncRNA CTB-113D17.1 contributed to the OS progression possibly via targeting the SCUBE3-hsa-miR-590-5p axis, which still needed to be confirmed in further bioinformatics analysis and experimental evidence.

HIST1H4E, another prominent regulator in ceRNA-hsa-miR-590-5p network, was found to be closely linked with the ribosome assembly process and also exhibited a significant correlation with CTB-113D17.1. Reports pointed out that the ribosome biosynthesis was a key prerequisite for gene expression and numerous controls in cellular processes mediated ribosome production and cell proliferation.[28] Sulima et al highlighted that there was strong correlations between ribosome defects and the risk of cancers,[29] which provided direct support for our findings. In addition, down-regulated EDIL3 also strongly interacted with hsa-miR-590-5p and CTB-113D17.1. Moreover, the functional enrichment analysis indicated this gene was significantly enriched in calcium ion binding pathway. Accumulating evidence has hinted that EDIL3 were responsible for a wide variety of cancers. Xia et al found that EDIL3 could induce the activation of TGF-beta signaling to trigger angiogenesis and recurrence of liver cancer.[30] However, few investigations have explored the influence of EDIL3 on OS development and progression, therefore, the detailed and integrated study required to be carried out to verify our results.

Although we explored the potential molecular mechanisms of HCC occurrence and development based on a ceRNA network analysis using a bioinformatics method, there were still several limitations in this research. For instance, an integrated bioinformatics analysis with a larger sample size still was needed to further illuminate the molecular interactions among different RNA transcripts.

In conclusion, 980 DEmRNAs and 682 DElncRNAs were screened between hsa-miR-590-5p-overexpression OS cells and normal cells. Moreover, lncRNA CTB-113D17.1 was possibly involved in OS pathogenesis by sponging hsa-miR-590-5p-SCUBE3/HIST1H4E/EDIL3 axis. However, our findings still required to be validated in future studies.

Click here for additional data file.

Supplementary_file_1 for Integrated Analysis of Competing Endogenous RNAs Network Reveals Potential Signatures in Osteosarcoma Development by Lisong Heng, Zhen Jia, Jian Sun, Yitong Zhao, Kun Zhang, Yangjun Zhu and Shemin Lu in Technology in Cancer Research & Treatment

Click here for additional data file.

Supplementary_file_2 for Integrated Analysis of Competing Endogenous RNAs Network Reveals Potential Signatures in Osteosarcoma Development by Lisong Heng, Zhen Jia, Jian Sun, Yitong Zhao, Kun Zhang, Yangjun Zhu and Shemin Lu in Technology in Cancer Research & Treatment