Identification of the miRNA-mRNA regulatory network of small cell osteosarcoma based on RNA-seq

Small cell osteosarcoma (SCO) is a rare subtype of osteosarcoma characterized by highly aggressive progression and a poor prognosis. The miRNA and mRNA expression profiles of peripheral blood mononuclear cells (PBMCs) were obtained in 3 patients with SCO and 10 healthy individuals using high-throughput RNA-sequencing. We identified 37 dysregulated miRNAs and 1636 dysregulated mRNAs in patients with SCO compared to the healthy controls. Specifically, the 37 dysregulated miRNAs consisted of 27 up-regulated miRNAs and 10 down-regulated miRNAs; the 1636 dysregulated mRNAs consisted of 555 up-regulated mRNAs and 1081 down-regulated mRNAs. The target-genes of miRNAs were predicted, and 1334 negative correlations between miRNAs and mRNAs were used to construct an miRNA-mRNA regulatory network. Dysregulated genes were significantly enriched in pathways related to cancer, mTOR signaling and cell cycle signaling. Specifically, hsa-miR-26b-5p, hsa-miR-221-3p and hsa-miR-125b-2-3p were significantly dysregulated miRNAs and exhibited a high degree of connectivity with target genes. Overall, the expression of dysregulated genes in tumor tissues and peripheral blood samples of patients with SCO measured by quantitative real-time polymerase chain reaction corroborated with our bioinformatics analyses based on the expression profiles of PBMCs from patients with SCO. Thus, hsa-miR-26b-5p, hsa-miR-221-3p and hsa-miR-125b-2-3p may be involved in SCO tumorigenesis.


INTRODUCTION
Osteosarcoma (OS) is a common primary bone tumor in young adults and adolescents [1], which is characterized by a high metastatic potential: approximately 20-25% of newly diagnosed patients harbor detectable lung-related metastases [2,3].Moreover, the distal femur and proximal tibia are frequent metastatic locations for OS.Despite recent improvements in therapeutic strategies, including surgical, radiotherapy and neoadjuvant chemotherapy, the 5-year survival of patients with OS remains below 70% [4].
Small-cell osteosarcoma (SCO) is a rare subtype of OS and accounts for less than 1% of OS cases [5].SCO is highly aggressive and frequently found in the metaphyses of long bones, which results in a poor prognosis.Histologically, SCO presents as small, round, primitive, and undifferentiated cells with osteoid production [6].Moreover, SCO does not respond to radiotherapy and is more aggressive than other subtypes of OS [7].Thus, novel therapeutic strategies for patients with SCO are urgently needed.
The etiologic factors and mechanism underlying the pathogenesis of OS are currently unclear.An increasing number of studies show that microRNAs (miRNAs) play important roles in the tumorigenesis and development of OS. miRNAs are small non-coding RNA molecules (18-25 nt) that increase or repress the expression of target genes [8].Specifically, the down-regulation of miR-133a www.impactjournals.com/oncotargetand miR-539 are associated with an unfavorable prognosis for patients suffering from OS [9].Moreover, miR-664 acts as an oncogene and promotes the proliferation of OS cells by down-regulating the expression of FOXO4.In addition, miR-664 promotes OS cell invasion and migration by suppressing the expression of SOX7 [10,11].Osteopontin is a phosphorylated glycoprotein involved in the invasion of OS cells, and the suppression of miR-4262 in OS cells promotes osteopontin-mediated cancer invasion [12].Furthermore, the down-regulation of miR-26a and up-regulation of miR-27a reportedly contribute to the aggressiveness of OS, and low miR-26a expression and high miR-27a expression are associated with advanced TNM stage, tumor grade, and distant metastasis in patients with OS [13].
SCO is an exceedingly rare subtype of OS associated with a poor prognosis.However, the mechanism underlying the pathogenesis of SCO is unclear, and previous studies of SCO are case reports due to the rarity of this disease.In this study, high-throughput RNA sequencing was performed to obtain the miRNA and mRNA expression profiles of peripheral blood mononuclear cells (PBMCs) from patients with SCO to identify genes that are differentially expressed in these patients compared with healthy individuals and explore potential diagnostic biomarkers for SCO.

Genes dysregulated in patients with SCO compared with healthy individuals
Blood samples from 3 patients with SCO and 10 healthy control individuals were subjected to highthroughput RNA sequencing; mRNA reads were used to align to the UCSC human reference genome (hg.19), and miRNA reads were matched to miRBase database.Differences in the expression levels of miRNA and mRNA between samples from patients with SCO and healthy individuals were then analyzed.
Compared with healthy individuals, 37 dysregulated miRNAs, including 27 up-regulated and 10 downregulated miRNAs, were identified in the PBMCs of patients with SCO patients, as shown in Table 1.Moreover, 1636 dysregulated mRNAs, including 555 up-regulated mRNAs and 1081 down-regulated mRNAs, were also identified in SCO.The top 20 significantly upregulated and top 20 down-regulated mRNAs are listed in Table 2.

KEGG pathway annotation of dysregulated genes in SCO
KEGG is a knowledge database for the systematic functional analysis of genes to connect genomic information with higher functional information by computerizing current knowledge on cellular processes and standardizing gene annotations.Thus, KEGG analyses were performed to assign pathways and functionally classify genes that are dysregulated in SCO.An FDR < 0.05 was defined as significant KEGG pathway enrichment.As shown in Table 3, pathways in cancer, mTOR signaling, oxidative phosphorylation, cell cycle, apoptosis, and DNA replication were significantly enriched KEGG pathways.

Verification of the expression level of dysregulated genes in tumor tissues and peripheral blood samples from patients with SCO
To verify the RNA-seq data, the expression levels of differentially expressed miRNAs and mRNAs in tumor tissues and peripheral blood samples from patients with SCO were quantified by qRT-PCR.As shown in Figure 3A-3C, hsa-miR-221-5p, hsa-miR-26b-5p and hsa-miR-21-5p were significantly up-regulated in both tumor tissues and peripheral blood samples from patients with SCO.The expression of miR-5706 was significantly up-regulated in peripheral blood samples from patients with SCO patients with healthy individuals, but its expression did not significantly differ between SCO tumors and adjacent nontumor tissues (Figure 3D).As shown in Figure 3E and 3F, miR-656-3p and RIF1 were significantly down-regulated in both tumor tissues and peripheral blood samples from patients with SCO, whereas FAM89A was significantly up-regulated in peripheral blood samples from patients with SCO compared with healthy individuals, but its  expression did not significantly differ between SCO tumors and adjacent non-tumor tissues (Figure 3G).In general, qRT-PCR data were consistent with our bioinformatics analyses.

DISCUSSION
miRNAs are endogenous, small non-coding RNAs that participate in the regulation of diverse cellular processes, including cell growth, cell differentiation, apoptosis and the immune response, by binding to the 3' untranslated region of certain subsets of messenger RNAs.A growing body of evidence demonstrates that the abnormal expression of miRNAs might serve as a potential diagnostic and prognostic biomarker for cancers [14].
In our study, the mRNA and miRNA expression profiles of peripheral blood samples from patients with SCO were obtained via RNA sequencing.Dysregulated mRNAs and miRNAs were identified in silico, and the expression levels of candidate mRNAs and miRNAs in both tumor tissues and peripheral blood samples from patients with SCO were verified by qRT-PCR.In general, the expression status of candidate mRNAs and miRNAs in peripheral blood samples of SCO patients were consistent with our bioinformatics analysis based on the expression profiling of PBMCs of SCO patients.Moreover, the expression levels of candidate mRNAs and miRNAs were also detected in SCO tumor tissues and adjacent non-tumor tissues.qRT-PCR results of candidate mRNAs and miRNAs had the similarity expression status in tumor tissues of OS patients based on bioinformatics analyses.
hsa-miR-26b-5p was significantly up-regulated in the PBMCs from patients with SCO compared to healthy control patients.In the miRNA-target gene interaction network, miR-26b-5p exhibited the highest connectivity with miRNAs and regulated the expression of 197 targetgenes.Specifically, miR-26b-5p is reportedly downregulated in OS tissue and cells [15,16] and inhibits OS cell proliferation, migration, invasion and apoptosis by downregulating the expression of PFKFB3 [15].Moreover, miR-26b inhibits OS cell metastasis by suppressing the expression of CTGF and Smad1 [17].Huaier (Trametes robiniophila murris) is a fungus used in traditional Chinese medicine and suppresses cell proliferation and induces apoptosis in human lung cancer cells by upregulating miR-26b-5p [18].miR-26b plays important roles in the proliferation and metastasis of various cancer types, such as epithelial ovarian carcinoma, hepatocellular carcinoma and prostate cancer [19][20][21].FASLG is one of top 20 down-regulated genes in SCO and targeted by miR-26b.Moreover, FASLG encodes Fas ligand, which is a member of tumor necrosis factor superfamily.FAS then triggers cell apoptosis by binding to FASLG.The FAS/ FASLG signaling pathway is essential for immune system regulation, including activation-induced T cell death and cytotoxic T lymphocyte-induced cell death.In our study,  SCO: small-cell osteosarcoma; FDR: false discovery rate; No.: number of genes in pathway.miR-26b-5p was up-regulated in the PBMCs of patients with SCO and targeted FASLG, which may play a key role in the immune response to SCO.Furthermore, hsa-miR-221-3p was significantly down-regulated in the PBMCs from patients with SCO.In the miRNA-mRNA regulatory network, miRNA-221 exhibited high connectivity with several mRNAs and targeted 18 mRNAs, and the expression levels of miR-221 in tumor tissue, cell lines and the sera from patients with OS are reportedly dramatically up-regulated.Thus, serum miR-221 may serve as a diagnostic and prognostic marker for OS [22,23].Moreover, miR-221 enhances the proliferation, invasion and migration of OS cells by suppressing PTEN expression [22].TIMP3 is one of top 20 up-regulated genes in SCO and targeted by miR-221.TIMP3 encodes TIMP metallopeptidase inhibitor 3, an inhibitor of matrix metalloproteinases, which are involved in the degradation of the extracellular matrix.TIMP3 is hypermethylated in OS tumor tissues compared with normal tissues [24] and has been identified as a tumor suppressor that plays essential roles in the inhibition of tumor angiogenesis.The down-regulation of miR-221/222 significantly increases the expression of TIMP3 and enhances the sensitivity of breast cancer cells to tamoxifen [25].TIMP3 expression also predicts favorable survival in HCC [26].In our study, miR-221-3p was down-regulated and targeted TIMP3 in the PBMCs of patients with SCO, which suggested that miR-221-3p and TIMP3 play essential roles in the immune response to SCO.
hsa-miR-125b-2-3p is one of the top 2 significantly down-regulated miRNAs and targeted 13 mRNAs, including PTGFRN (one of the top 20 up-regulated mRNAs) in the miRNA-mRNA regulatory network of SCO.PTGFRN encodes prostaglandin F2 receptor inhibitor, and the expression of PTGFRN correlates with the metastatic status of human OS [27].However, the biological functions of miR-125b-2-3p and PTGFRN are unclear and need to be elucidated.
Pathways in cancer, mTOR signaling and oxidative phosphorylation were significantly enriched in the KEGG pathway analysis of dysregulated mRNAs identified in the PBMCs of patients with SCO compared to healthy individuals.Dysregulated pathways in cancer contribute to aberrant cell growth, cell death and cell motility in various cancer types, including colorectal cancer, pancreatic cancer, glioma, basal cell carcinoma, renal cell carcinoma, prostate cancer and melanoma [28,29].Moreover, a series of articles report that abnormal mTOR signaling plays key roles in OS tumorigenesis and therapy [30], and arsenic sulfide promotes cell apoptosis and autophagy by suppressing Akt/ mTOR signaling in OS [31].In addition, cucurbitacin E also inhibits OS cell growth and invasion by suppressing PI3K/ AKT/mTOR signaling [32], whereas a lack of oxidative phosphorylation prevents cell apoptosis in colorectal cancer and OS.Accordingly, PSB-603 (A2b adenosine receptor antagonist) increases colorectal cancer cell death by promoting oxidative phosphorylation and ROS production [33].A lack of oxidative phosphorylation also decreases susceptibility to apoptosis in OS [34].In addition, cell cycle, apoptosis, and DNA replication pathways were significantly enriched in our study.Based on the aforementioned studies, the identified abnormally expressed mRNAs might play vital roles in OS tumorigenesis by regulating these enriched KEGG pathways.
The expression levels of candidate miRNAs and mRNAs in peripheral blood samples from patients with SCO were verified by qRT-PCR.In general, the qRT-PCR results were consistent with the RNA sequencing and bioinformatics analyses.However, our study was also subject to limitations.First, the SCO tumor sample size for the RNA sequencing and qRT-PCR validation analyses were small due to the rarity of SCO.Second, the expression levels of candidate miRNAs contradicted previously published results.Thus, the biological roles of these miRNAs need to be explored in vivo and in vitro to elucidate their roles in the tumorigenesis of SCO.
In conclusion, we constructed an miRNA-mRNA regulatory network by identifying miRNAs and target genes that are differentially expressed in the PBMCs of patients with SCO compared to healthy individuals.In this regulatory network, hsa-miR-26b-5p, hsa-miR-221-3p and hsa-miR-125b-2-3p were significantly dysregulated and exhibited high connectivity with target genes.Thus, these miRNAs may play essential roles in the initiation and development of SCO, but their biological roles in SCO tumorigenesis require further exploration.Moreover, a clinical study based on a large cohort of patients with SCO is essential to validate the diagnostic value of these miRNAs.

Sample isolation and characterization
Thirteen patients treated at the Third Affiliated Hospital of Kunming Medical University hospital, including 3 patients with SCO and 10 healthy control individuals, were enrolled in our study.Ten milliliters of peripheral blood was obtained from each of patients, and PMBCs were isolated.The total RNA was then extracted from the PMBCs using TRIzol reagent (Invitrogen, Carlsbad, CA, USA).Our study was approved by the Third Affiliated Hospital of Kunming Medical University hospital, and informed written consent was obtained from all patients.Our study was conducted in accordance with the Declaration of Helsinki.

Transcriptome library preparation and sequencing
The Illumina TruSeq RNA sample Prep Kit (Illumina, Inc., San Diego, CA, USA) was used to prepare the mRNA library according to the manufacturer's protocol.Oligo(dT) www.impactjournals.com/oncotargetbeads were used to isolate polyA mRNA from total RNA, and the mRNA was cleaved into short fragments for random hexamer priming to synthesize first-strand cDNA, followed by the second-strand cDNA synthesis and end repair.Next, the short fragments were connected with sequencing adapters.After PCR amplification, the enriched cDNA libraries were sequenced using an Illumina HiSeq 2500 (Illumina, Inc., San Diego, CA, USA).For miRNA sequencing, the miRNAs were isolated, and adapters were added.The library was constructed by using TruSeq Small RNA Library Preparation Kits (Illumina, Inc., San Diego, CA, USA).

Data analyses
The raw RNA-sequencing reads were stored in FASTQ format.The raw reads containing adapter sequences and low-quality sequences (reads with ambiguous bases 'N') were removed using FASTx-tool SeqPrep (https://github.com/jstjohn/SeqPrep)and Sickle (https://github.com/najoshi/sickle). Clean reads were aligned with the UCSC human reference genome (build hg19) using TopHat v1.3.1 [35].Sequences were matched to miRNAs using the miRDeep2 tool and miRBase (release 21) (http://www.mirbase.org/).Aligned read files were then processed by Cufflinks v1.2.1 [36], which measured the relative expression of the genes with the normalized RNA-Seq fragment counts.Fragments per kilobase of exon model per million mapped reads (FPKM) were used to present the expression level of gene.The transcripts differentially expressed in SCO were identified, and P-values < 0.05 and < 0.01 indicated differentially expressed mRNA and miRNA, respectively.

Functional annotation of differentially expressed genes
The biological roles of differentially expressed genes were predicted using the GeneCoDis3 (http:// genecodis.cnb.csic.es/analysis)database to describe the gene ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) pathways [37,38].An FDR < 0.05 was set as the cut-off for selecting significantly enriched GO terms and KEGG pathway.

Construction of miRNA-target gene network
The target-genes of differentially expressed miRNAs were predicted using the miRWalk database (http://www.umm.uni-heidelberg.de/apps/zmf/mirwalk/), in which the correlations between target-genes and miRNAs have been experimentally confirmed in vivo and in vitro [39].We used 6 algorithms, RNA22, miRanda, miRDB, miRWalk, PICTAR2 and TargetScan, to predict the target-genes of miRNAs; if more than 4 of 6 algorithms predicted the same gene for an miRNA, the gene was identified as a target-gene of that miRNA [39].The interaction network of differentially expressed miRNAs and target-genes were constructed using the Cytoscape software (http:// cytoscape.org)[40].

Quantitative real-time polymerase chain reaction (qRT-PCR)
To validate the expression levels of dysregulated genes in SCO, tumor tissues and peripheral blood samples were obtained from other 3 patients with SCO and 3 healthy individuals.Total RNA was extracted from tissues samples and peripheral blood samples using TRIzol reagent (Invitrogen, CA, USA) according to the manual instructions.The ReverTra Ace qPCR RT Master Mix Kit (TOYOBO, Shanghai, China) was used to synthesize cDNA, and qRT-PCR reactions were performed using SYBR ® FAST qPCR Kits (KAPA bio, Boston, USA) on a LightCycler 480 (Roche Indianapolis, IN, USA).U6 and GAPDH were used as internal controls for miRNA and mRNA detection, respectively.The relative expression levels of target genes were calculated using ΔCT method, and the mean ± standard deviation and independent samples t-test were used in the statistical analysis.P < 0.05 was considered significant.The PCR primers used are shown in Supplementary Table 1.

Figure 1 :
Figure 1: The Gene Ontology term enrichment of dysregulated mRNAs in SCO.FDR: false discovery rate.(A): Biological processes of Gene Ontology terms; (B): molecular functions of Gene Ontology terms; (C): cellular components of Gene Ontology terms.

Figure 2 :
Figure 2: miRNA-mRNA interaction network in SCO.(A): Interaction network between up-regulated miRNAs and downregulated mRNAs; (B): interaction network between down-regulated miRNAs and up-regulated mRNAs.Circular nodes represented mRNAs, and diamond nodes represented miRNAs.Green and red colors represented down-regulation and up-regulation, respectively.Solid lines indicated interaction associations between miRNAs and mRNAs.

Table 1 : Differentially expressed miRNAs in SCO compared with normal controls miRNA log 2 FC P-value Up-regulated miRNAs
SCO: small-cell osteosarcoma; FC: fold change.

Table 2 : Differentially expressed mRNAs in SCO compared with normal controls
SCO: small-cell osteosarcoma.