miRNA profiling identifies deregulated miRNAs associated with osteosarcoma development and time to metastasis in two large cohorts

Osteosarcoma (OS) is an aggressive bone tumor primarily affecting children and adolescents. The etiology of OS is not fully understood. Thus, there is a great need to obtain a better understanding of OS development and progression. Alterations in miRNA expression contribute to the required molecular alterations for neoplastic initiation and progression. This study is the first to investigate miRNA expression in OS in a large discovery and validation cohort comprising a total of 101 OS samples. We established the signature of altered miRNA expression in OS by profiling the expression level of 752 miRNAs in 23 OS samples using sensitive LNA‐enhanced qPCR assays. The identified miRNA expression changes were correlated with gene expression in the same samples. Furthermore, miRNA expression changes were validated in a second independent cohort consisting of 78 OS samples. Analysis of 752 miRNAs in the discovery cohort led to the identification of 33 deregulated miRNAs in OS. Twenty‐nine miRNAs were validated with statistical significance in the second cohort comprising 78 OS samples. miRNA/mRNA targets were determined, and 361 genes with an inverse expression of the target miRNA were identified. Both the miRNAs and the identified target genes were associated with multiple pathways related to cancer as well as bone cell biology, thereby correlating the deregulated miRNAs with OS tumorigenesis. An analysis of the prognostic value of the 29 miRNAs identified miR‐221/miR‐222 to be significantly associated with time to metastasis in both cohorts. This study contributes to a more profound understanding of OS tumorigenesis, by substantiating the importance of miRNA deregulation. We have identified and validated 29 deregulated miRNAs in the – to our knowledge – largest discovery and validation cohorts used so far for miRNA analyses in OS. Two of the miRNAs showed a promising potential as prognostic biomarkers for the aggressiveness of OS.

Osteosarcoma (OS) is an aggressive bone tumor primarily affecting children and adolescents. The etiology of OS is not fully understood. Thus, there is a great need to obtain a better understanding of OS development and progression. Alterations in miRNA expression contribute to the required molecular alterations for neoplastic initiation and progression. This study is the first to investigate miRNA expression in OS in a large discovery and validation cohort comprising a total of 101 OS samples. We established the signature of altered miRNA expression in OS by profiling the expression level of 752 miRNAs in 23 OS samples using sensitive LNA-enhanced qPCR assays. The identified miRNA expression changes were correlated with gene expression in the same samples. Furthermore, miRNA expression changes were validated in a second independent cohort consisting of 78 OS samples. Analysis of 752 miRNAs in the discovery cohort led to the identification of 33 deregulated miRNAs in OS. Twentynine miRNAs were validated with statistical significance in the second cohort comprising 78 OS samples. miRNA/mRNA targets were determined, and 361 genes with an inverse expression of the target miRNA were identified. Both the miRNAs and the identified target genes were associated with multiple pathways related to cancer as well as bone cell biology, thereby correlating the deregulated miRNAs with OS tumorigenesis. An analysis of the prognostic value of the 29 miRNAs identified miR-221/ miR-222 to be significantly associated with time to metastasis in both cohorts. This study contributes to a more profound understanding of OS tumorigenesis, by substantiating the importance of miRNA deregulation. We have identified and validated 29 deregulated miRNAs in theto our knowledgelargest discovery and validation cohorts used so far for miRNA analyses in OS. Two of the miRNAs showed a promising potential as prognostic biomarkers for the aggressiveness of OS.

Introduction
Osteosarcoma (OS) is an aggressive tumor of the bone most frequently affecting children and adolescents. OS is often diagnosed at a very late stage, as the first indication of the malignancy is pain, which is often mistaken for more common conditions such as 'growing pains' (McCarville, 2009). The late identification of OS leads to 15-25% of the patients exhibiting distant metastases at the time of diagnosis (Bielack et al., 2002). The five-year survival rate of OS patients with no detectable metastases at the time of diagnosis is 65-70% (Bielack et al., 2002), whereas patients with relapsing disease or detectable metastases at the time of diagnosis have only a five-year survival rate of 10-40% (Kager et al., 2003;Meyers et al., 1993).
OS is characterized by diverse cytogenetic alterations and exhibits changes in multiple different pathways with substantial cell-to-cell variation (Bridge et al., 1997). The etiology of OS is not well understood and the only prognostic markers are the absence of metastases at the time of diagnosis and the response to neoadjuvant chemotherapy measured as the extent of tumor necrosis (Wang, 2005). However, the characteristic onset of high-grade neoplasms in young people indicates a consistent, but yet undetected, alteration or group of alterations defining OS.
Cancer initiation and progression is controlled by both genetic and epigenetic events (Hanahan and Weinberg, 2011). Epigenetic changes are particularly important for childhood cancers, which are often characterized by a low mutational load. Molecular subtypes, clinical heterogeneity, and in some cases the prediction of the disease course have been associated with epigenetic changes. This has been demonstrated for ependymomas (Mack et al., 2014) and more recently Ewing sarcoma (Sheffield et al., 2017), which are characterized by no or the same recurrent genetic aberration, respectively. The genetic components of OS have been intensively studied (Bishop et al., 2016;Kansara et al., 2014); however, the elucidation of the epigenetic changes is still in its early stage. Gene expression is regulated by numerous mechanisms and molecules, including microRNAs (miRNAs). These are small noncoding RNA molecules (22-25 nucleotides), regulating gene expression post-transcriptionally by blocking mRNA translation and/or altering the stability of the mRNA by binding perfectly or with mismatches to the 3 0 -untranslated region (UTR) of the mRNA (Stark et al., 2003). Each miRNA can potentially directly or indirectly regulate more than 100 target genes and thereby affect the expression of both tumor suppressors and oncogenes (Brennecke et al., 2005). Furthermore, miRNAs function as classical tumor suppressors and oncogenes/oncomiRs as specific miRNA alterations are associated with carcinogenesis. Deregulation of miRNA expression has been identified in numerous cancers emphasizing the importance of miRNAs in carcinogenesis and their potential as diagnostic, prognostic, and therapeutic means (Calin and Croce, 2006;Negrini et al., 2009). miRNAs are also important regulators of osteogenic signaling pathways and are implicated in normal osteoblast (OB) growth and differentiation (Inose et al., 2009;Lian et al., 2012). Hence, changes in this balance will affect proliferation of bone cells and may contribute to OS development. The impact of alterations in miRNA expression patterns in OS development has been addressed in several studies (Duan et al., 2011;Jones et al., 2012;Lulla et al., 2011;Maire et al., 2011;Namlos et al., 2012;Thayanithy et al., 2012;Won et al., 2013;Zhang et al., 2015). However, none of the results are consistent across the studies, which may reflect different methodologies and normalization strategies, as well as low sample numbers. Most studies have evaluated a single miRNA, but a few studies have established a miRNA expression profile of OS by profiling several hundred miR-NAs (Duan et al., 2011;Jones et al., 2012;Lulla et al., 2011;Maire et al., 2011;Namlos et al., 2012;Won et al., 2013;Zhang et al., 2015). However, all studies have analyzed a limited number of OS samples ranging from 2 to 18.
In this study, we established the signature of miRNA expression alterations in OS by profiling the expression of 752 miRNAs in 23 OS samples using a highly sensitive qPCR technology. Significantly deregulated miRNAs were validated in an independent cohort comprising 78 OS samples. Both the miRNAs and identified target mRNAs were evaluated for their relation to cancer development and progression, as well as bone cell biology.

OS samples
Fresh-frozen diagnostic biopsies of primary tumors from 23 patients diagnosed with OS were collected from the archives of the Institute of Pathology, Aarhus University Hospital (Table 1). The biopsies were immediately snap-frozen and stored at À80°C. 25-50 tissue sections of 10 lm were used for RNA extraction for each sample. Tissue sections were hematoxylin/ eosin (HE)-stained for every fifth 10-lm tissue section cut to confirm the presence of tumor cells.
Formalin-fixed, paraffin-embedded (FFPE) samples from 78 primary tumors from patients diagnosed with OS were collected from the archives of the Institute of Pathology, Aarhus University Hospital (Table 1). A complete outline of the different procedures and analyses performed with each cohort is shown in Fig. 1.

Cell lines
Five OS and three OB cell lines were used in this study. The OS cell lines CRL-1543, CRL-1427, CRL-2098 and the OB cell line CRL-11372 were purchased from American Type Culture Collection (ATCC, Manassas, VA, USA). The OB cell line C-12720 was purchased from PromoCell (Heidelberg, Germany), and the OB cell line HO-f-4610 was purchased from SanBio (Uden, The Netherlands). HTB-96 and HTB-85 were cultured in McCoy's 5a medium modified. CRL-1427 and CRL-1543 were cultured in Eagle's minimum essential medium. CRL-2098 was cultured in RPMI-1640. All OS cell lines were supplemented with 40 mgÁmL À1 streptomycin, 240 mgÁmL À1 penicillin (400 000 UÁmL À1 ), and 10% fetal bovine serum, except for HTB-85, which was supplemented with 15% fetal bovine serum. CRL-11372 was cultured in 1 : 1 mixture of Ham's F12 medium and Dulbecco's modified Eagle's medium without phenol red, supplemented with 0.3 mgÁmL À1 G-418 and 10% fetal bovine serum. C-12720 was cultured in OB growth medium + supplement mix (PromoCell). HO-f was cultured in OB medium (ScienCell Research Laboratories, Carlsbad, CA, USA), and the flasks were coated with poly-L-lysine solution bioreagent 0.01% (Sigma Aldrich, St Louis, MO, USA).
All eight cell lines were cultured in a humidified atmosphere of 5% CO 2 at 37°C, except for CRL-11372, which was cultured at 34°C.

Total RNA extraction from fresh-frozen OS samples and cell lines
Total RNA including miRNAs was extracted from the 23 fresh-frozen OS samples and the eight cell lines. The tissue was homogenized using a TissueLyser bead mill and one tungsten carbide bead (3 mm) per sample. RNA extraction was performed using the miR-Neasy Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions until each sample was separated into three phases. The upper aqueous phase was transferred to a new tube and total RNA  Fig. 1. A detailed outline of the different processes and analyses performed for the OS discovery and validation cohort. QC, quality control; diff., differential; btw., between.
was purified using a QIAcube automated device (Qiagen) according to the manufacturer's instructions including an on-column DNase digest. RNA was quantified using a NanoDrop ND-2000c spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA), and the quality was assessed using the RNA 6000 Nano Kit (Agilent Technologies Inc., Santa Clara, CA, USA) on an Agilent 2100 BioAnalyzer.

RNA extraction from FFPE samples
Four tissue sections of 10 lm from each FFPE sample were used for RNA extraction. RNA was extracted using the miRNeasy FFPE Kit (Qiagen) according to the manufacturer's instructions (Fig. 1). One microlitre Spike-in (UniSp2, UniSp4, and UniSp5), from the RNA Spike-in kit, UniRT (Exiqon, Vedbaek, Denmark), was added to the PKD buffer during the purification.

cDNA synthesis and miRNA expression analysis of fresh-frozen OS samples
For the initial miRNA profiling analysis, 50 ng of total RNA from 23 fresh-frozen OS samples, two OS, and three OB cell lines was reverse-transcribed using the miRCURY LNA TM Universal RT microRNA PCR kit (Exiqon) in a Mastercycler Ò pro S (Eppendorf, Hamburg, Germany) in a 50 lL reaction volume, comprising the Universal cDNA synthesis kit II, RNA Spike-in kit (Exiqon) (Fig. 1). Quality control (QC) was performed according to the manufacturer's instructions using the ExiLENT SYBR Ò Green master mix kit (Exiqon), using 1 : 80 diluted cDNA. Quality of the cDNA synthesis and PCR amplification were evaluated through the detection of miR-103, miR-191-5p, miR-423, and UniSp6 at stable C q levels (C q < 32) (Fig. 1).
The miRNA profile was established using 45 ng cDNA from each sample by qPCR using the micro-RNA Ready-to-Use PCR, Human panel I and panel II v2 analyzing 752 miRNAs (Exiqon) on a Lightcycler Ò 480 V2 Real-Time PCR System (Roche Applied Science, Penzberg, Germany) (Fig. 1). The plates were set up using a Bravo Automated Liquid Handling Platform (Agilent Technologies Inc.).

Validation of miRNA expression in OS FFPE samples
A total of 65 ng RNA (5 ngÁlL À1 ) from the OS FFPE samples and the cell lines was reverse-transcribed in 20 lL reactions using the Universal cDNA Synthesis kit II (Exiqon) and a 2720 Thermal Cycler (Thermo Fisher Scientific) (Fig. 1). Reverse transcription was performed in duplicate for each sample. All the subsequent amplifications were performed using a Lightcycler Ò 480 Real-Time PCR System (Roche Applied Science). QC of the cDNA synthesis and PCR amplification were performed according to the manufacturer's instructions (Exiqon), with slight modifications; the cDNA was diluted 1 : 20 to ensure the detection of lowly expressed miRNAs. Spike-ins were included in the RNA purification and RT step using the RNA Spike-in kit, UniRT (Exiqon). hsa-miR-191-5p, hsa-miR-423-5p, UniSp4, UniSp5, UniSp6, and Cel-miR-39-3p were used for QC of the RNA purification and cDNA synthesis (Fig. 1). UniSp4 and UniSp5 RNA spike-in controls were used for QC of RNA purifications. All primary OS FFPE samples, except one, had C q values < 37 for UniSp5. UniSp6 and Cel-miR-39-3p spike-in control primers were used for QC of the cDNA synthesis and to rule out PCR inhibition, as a higher concentration of cDNA was used in the analysis, than described in the manufacturer's instructions (65 ng RNA instead of 20 ng for the cDNA synthesis and 1 : 20 dilution of the cDNA instead of 1 : 100). These conditions were optimized for these amounts to ensure adequate amplification from the FFPE samples. All primary OS FFPE samples had C q values between 16.5 and 18 for UniSp6. Negative controls included in the cDNA synthesis were a no RNA template and a no RT enzyme. miRNA expression analysis was performed for 78 FFPE primary OS samples, five OS, and the three OB cell lines using 5 lL 1 : 20 diluted cDNA from each sample using a custom Pick&Mix microRNA PCR Panel, 384 well Ready-to-Use (Exiqon) (Fig. 1). Five reference miRNAs, hsa-miR-148b-3p, hsa-miR-185-5p, hsa-miR-191-5p, hsa-miR-423-5p, and hsa-miR-425-3p, were selected based on the results from the discovery cohort, using both geNorm (Vandesompele et al., 2002) and NormFinder (Andersen et al., 2004). All 78 OS FFPE samples had C q values < 30 for the reference miRNAs.

Gene expression analysis
Genomic RNA expression analysis was performed using the SurePrint G3 Human 8X60K One color Microarrays (Agilent), which analyzes 27958 Entrez Gene RNAs and 7419 lincRNAs. Hundred nanogram RNA from each sample of the 23 fresh-frozen OS primary tumors, five OS cell lines, and three OB cell lines was analyzed as per the manufacturer's instructions. The Agilent G3 Microarrays were scanned using the NimbleGen MS200 scanner with Multi-TIFF as File Settings Control, a single scan area was selected, and the entire feature area was captured. Two-micrometre resolution was selected. The red channel was omitted from the analysis. Grid alignment was determined manually. Data were extracted using the Agilent Feature Extraction 11.5.1.1. The raw data were preprocessed using ArrayStar v12.1.0 build 134 (DNASTAR, Madison, WI, USA) and normalized using the quantile normalization method and log 2 -transformed. Genes were considered differentially expressed when the difference between OS samples and OB cell lines was > 2fold and with a Bonferroni-corrected P-value < 0.05. Only genes with an average normalized probe intensity of log 2 > 6 in OB cell lines were considered expressed.

Data analysis
miRNA expression analysis was performed using the GENEX software version 5.4.2.128 (MultiD, G€ oteborg, Sweden). Raw data (C q values) for the 752 miRNAs analyzed using the discovery cohort were imported and merged (Table S1). Each sample has been run with two 384-well plates containing different miRNAs. Hence, the raw data contain numerous NA values for the miR-NAs not present in the respective plates. The mean global expression was used for normalization of the Human Panel I+II v2. The limit for a detectable miRNA was set to C q < 37. Log 2 -normalized data are available in Table S2. Pick&Mix microRNA PCR Panels were normalized using five selected stably expressed reference genes (hsa-miR-148b-3p, hsa-miR-185-5p, hsa-miR-191-5p, hsa-miR-423-5p, and hsa-miR-425-3p) (Fig. 1). Raw data (C q values) for the 33 miRNAs analyzed in the validation cohort are available in Table S3. The raw data have been imported and merged using the GENEX software version 5.4.2.128 (MultiD). Log 2 -normalized data for the validation cohort are available in Table S4.
The statistical software R v.3.0.0 (R Core Team 2013) was used for unsupervised hierarchical clustering using the 'Ward' method as distance measure with an in-house-developed script. Pathway and network analyses were performed using Ingenuity Ò Pathway Analysis (IPA Ò ; Qiagen), the Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa and Goto, 2000;Kanehisa et al., 2014), and DIANA miRPath.v3 (Vlachos et al., 2015). For the IPA core analysis, the KEGG pathway analysis, and the DIANA pathway analysis, human was selected as the species analyzed.
Validated miRNA gene targets were established using miRWalk2.0 (Dweep and Gretz, 2015) with the identification of miRNA binding sites in the 3 0 UTR. The generated gene lists of validated miRNA gene targets were compared with a list of genes comprising Boxplots, Mann-Whitney U-tests, and Kaplan-Meier analysis were made with the GRAPHPAD PRISM 7 Software (GraphPad Software Inc., La Jolla, CA, USA). For Kaplan-Meier analysis, the log-rank (Mantel-Cox) test was used for the determination of P-values. For the prognostic analyses (survival and time to metastasis), median expression was used to subdivide the samples into high and low expression groups.

Data accessibility
All miRNA data generated and analyzed during this study are included in this published article (Tables S1-S4 including both raw and log 2 -normalized data). The normalized log 2 expression data for the 361 genes analyzed in this study are provided in Table S5. The raw gene expression dataset used in the current study is part of another ongoing study and will be published and deposited in a public repository separately. Meanwhile, these data are available from the corresponding author upon request.

Results
In the present study ( Fig. 1), we investigated changes in miRNA expression of 752 miRNAs in primary OS and OS cell lines compared to the miRNA expression levels of OB cell lines. OB cell lines were used as a reference as OBs/osteoprogenitor cells are considered to represent the cell of origin for OS (Abarrategi et al., 2016). The 33 miRNAs identified as differentially expressed in the OS discovery cohort were validated in an independent validation cohort comprising 78 primary OS samples.

Overall miRNA expression
The relative expression levels of 752 miRNAs were evaluated in 23 primary OS samples, two OS, and three OB cell lines. Following normalization, using the global mean, a total of 339 and 202 miRNAs were reliably detected in more than 40% and 90% of the samples, respectively. The frequency of missing data was calculated based on all samples (both OS samples and OB cell lines). Further analyses were performed using the dataset with the expression of 339 miRNAs to ensure analysis of all miRNAs important for OS etiology. If using a 90% threshold, miRNAs not expressed in OB cell lines or miRNAs downregulated/ depleted in OS samples could be removed, thereby excluding biologically relevant miRNAs.
An unsupervised hierarchical clustering was performed using the log 2 -transformed expression values for the 339 miRNAs (Fig. 2). The samples separated into two main clusters, of which one contained 17 OS samples and the two OS cell lines. Even though the OS cell lines were in a cluster with OS samples, the cell lines kept the closest mutual resemblance, indicating a cell line similarity stronger than the cell origin. The other cluster contained six OS samples and the three OB cell lines. There was no difference between tumor content or tumor/patient characteristics in these six OS samples compared with the rest of the OS cohort.
The two observed clusters could not be explained with the available clinical and histopathological data as no separation in terms of the primary tumor site, histology, age, gender, or metastatic status was observed, even though a smaller group containing five OS samples separated from the other OS samples (Fig. S1). These samples had a slightly higher miRNA detection level, which might explain the separation. Of note, these samples were different from the six OS samples grouping apart in the cluster.

Identification of differentially expressed miRNAs in OS
When comparing expression profiles between OS samples and OB cell lines, 76 miRNAs were identified as differentially expressed (P-value < 0.05) ( Table S6). As OB cell lines were used as the normal control, miRNAs with no expression change between OS and OB cell lines were excluded to avoid expression changes caused by cell culture and not as a result of cell origin. This provided a list of 33 miRNAs with deregulated expression levels, of which 26 were significantly downregulated and seven were upregulated in the OS samples compared to OB cell lines (Table 2, Fig. 3 and Fig. S2A,B). Four of these miRNAs (miR-181a-5p, miR-181b-5p, miR-30e-3p, and miR-92b-3p) did not exhibit an expression difference between OS and OB cell lines, but were, nonetheless, selected for further validation as miR-181a-5p, miR-181b-5p, and miR-92b-3p have previously been associated with OS development and miR-30e-3p has been shown to be deregulated in several other cancers (Duan et al., 2011;Jones et al., 2012;Ning et al., 2016;Zhu et al., 2014). Of the 33 differentially expressed miRNAs, 12 miRNAs have not previously been described as differentially expressed in OS (Table S7). Of the previously described miRNAs, 16 were identified with the same direction of the expression change as in this study, while five miRNAs showed an opposite expression change (Table S7). On the other hand, previous studies using microarrays or sequencing to investigate miRNA expression changes in OS identified 14 miR-NAs, which were not part of the 33 validation miRNAs found in this study (each miRNA was identified in three or more studies) (Table S8) investigation of these 14 miRNAs in the OS discovery cohort analyzed in this study showed that six of the miRNAs did exhibit a differential expression in OS samples when compared to OB cell lines, but no difference was observed between OS and OB cell lines; four of the miRNAs did not exhibit a significant difference between OS samples and OB cell lines, and four miR-NAs were not present in our normalized data set.

Validation of differentially expressed miRNAs
The 33 miRNAs identified in the discovery cohort were validated in a larger, independent cohort consisting of 78 primary OS FFPE samples using a custom-designed Pick&Mix microRNA PCR Panel.
Of the 26 miRNAs with decreased expression, 24 were confirmed in the validation cohort, and of the seven miRNAs with increased expression, five were confirmed, thereby corroborating 88% of our initial findings (Table 2, Fig. 3 and Fig. S2A,B).

Pathways associated with the deregulated miRNAs and their target genes
To determine the possible impact of the 29 differentially expressed miRNAs on OS development, enriched biological functions were explored using Ingenuity Pathway Analysis (IPA) and DIANA Tools mir-Path.v3 exploring KEGG pathways. Several biological functions related to cancer development and progression were enriched for the 29 deregulated miRNAs (Table 3). Alterations in miRNA expression may have an extensive effect on gene expression, as one miRNA can directly or indirectly regulate more than 100 target genes by binding to their 3 0 UTR. Genes targeted by at least one of the 29 differentially expressed miRNAs were therefore identified using miRWalk2.0 analyzing only experimentally validated target genes. Genomewide RNA expression changes in OS were established by analyzing the 23 OS samples from the discovery cohort with RNA expression microarrays containing 27958 Entrez Gene RNAs (Andersen et al., in preparation). The identified target genes were correlated with genes identified as differentially expressed in the OS samples compared to OB cell lines (> 2-fold change and P-value < 0.05). This identified 301 genes with upregulated expression in OS samples, targeted by one or more of the identified downregulated miRNAs, as well as 60 genes with decreased expression in the OS samples targeted by one or more of the miRNAs with increased expression in the OS samples (Table S5). The possible influence of the deregulation of miRNA expression on OS development, associated with the affected target genes, was investigated using IPA and KEGG pathway analysis. The most significant pathways, as well as biological functions associated with disease and cellular functions for validated target genes, are reported in Table 4. For IPA, six of the top ten biological functions and three of the top five pathways were associated with cancer development and progression. Furthermore, one of the biological functions, for which deregulated genes were enriched, was 'Skeletal Disorders' and one of the top canonical pathways was 'Role of Osteoblasts, Osteoclasts, and Chondrocytes in Rheumatoid Arthritis'. For the KEGG analysis, nine of the ten highest-ranking pathways were related to cancer development and progression. The 27th pathway for the KEGG analysis was 'Osteoclast Differentiation', thereby linking bone cell biology and cancer development and progression for the target genes with altered gene expression.
An analysis of the biological relations between the target genes identified a network of 21 upregulated Diseases and disorders Cancer c 4.14 9 10 À13 18 Connective tissue disorders 4.14 9 10 À13 12 Organismal injury and abnormalities 4.14 9 10 À13 20 Reproductive system disorders 8.23 9 10 À9 15 Developmental disorder 4.49 9 10 À8 6 Molecular and cellular functions Cell death and survival c 2.84 9 10 À4 5 Cell-to-cell signaling and interaction c 0.00342 1 Cellular development c 0.00342 6 Cellular growth and proliferation c 0.00570 6 Cell cycle c 0.0125 1 (B) Pathway P-value

IPA biological functions P-value
No. of genes associated Diseases and disorders Connective tissue disorders 2.86 9 10 À13 77 Inflammatory disease 2.86 9 10 À13 91 Organismal injury and abnormalities 2.86 9 10 À13 348 Skeletal and muscular disorders b 2.86 9 10 À13 104 Cancer a 4.66 9 10 À11 343 Molecular and cellular functions Cellular growth and proliferation a 1.28 9 10 À11 121 Cellular movement a 1.70 9 10 À10 80 Cellular Development a 1.05 9 10 À9 86 Cell-to-cell signaling and interaction a 5.07 9 10 À8 60 Cell death and survival a 6.10 9 10 À7 102 KEGG pathways and eight downregulated genes significantly associated with 'Cell Death and Survival, Cancer, Hematological System Development and Function' (P-value = 10 À42 ). The intimate relationship of the 29 genes and the seven downregulated and four upregulated miRNAs validated to target these genes is illustrated in Fig. 4. The identified network comprised five genes (CCNG1, MDM2, TNFSF10, TP53, and UCHL1) directly associated with the TP53 pathway. Three additional target genes (ABL1, FASLG, and UBD) not part of the identified network were associated with the TP53 pathway. Four target genes (ABL1, MDM2, STMN1, and TP53) were associated with the RB pathway.

miRNAs associated with the metastatic status and response to chemotherapy
To identify miRNAs with an overall difference in expression based on the metastatic status, the discovery cohort was separated into metastatic and nonmetastatic primary OS (8 and 15 OS samples, respectively). This analysis revealed three miRNAs (miR-29b-3p, miR-29c-3p, and miR-374a-5p) significantly differentially expressed between metastatic and nonmetastatic OS samples (Table S9 and Fig. S3). The response to neoadjuvant chemotherapy was analyzed by subdividing the discovery cohort into poor responders (< 90% necrosis after chemotherapy) and good responders (> 90% necrosis after chemotherapy) (8 and 9 OS samples, respectively). This identified miR-663b as significantly differentially expressed. Also, miR-664a-3p showed a strong trend toward differential expression (Table S9 and Fig. S3). However, when these miRNAs were analyzed in the samples from the validation cohort, none of the three miRNAs associated with the metastatic status nor the two miRNAs associated with the chemotherapeutic response displayed a significant difference (Table S9 and Fig. S3).

Prognostic value of the 29 differentially expressed miRNAs
The only prognostic markers for OS are currently absence of metastases at the time of diagnosis and response to chemotherapy. As 25-30% of OS patients develop distant metastases and the five-year survival rate for patients with relapse is only 10-40%, it is critical to identify new prognostic biomarkers for OS. We therefore investigated the prognostic value in function of survival and time to development of metastatic disease of the 29 differentially expressed miR-NAs. Kaplan-Meier analysis was performed on the discovery cohort divided into two groups based on the median expression (high and low expression groups for each miRNA). For the overall survival analysis, two miRNAs (miR-128-3p and miR-34a-3p) were borderline significant in the discovery cohort (P = 0.0880 and 0.0760, respectively). No significant difference was, however, observed in the larger OS cohort (P = 0.5645 and 0.7507, respectively).
Analysis of the miRNA expression in function of time to metastasis, describing the aggressiveness of the OS tumor, identified miR-222-3p and miR-221-3p as significantly (P = 0.0359) and borderline significantly (P = 0.0863) associated with time to metastasis in the discovery cohort, respectively (Fig. S4A,C). High expression of miR-221-3p and low expression of miR-222-3p correlated with a poor prognostic value and earlier development of distant metastases. The association between expression level and poor prognosis in terms of earlier development of distant metastases was highly significant in the validation cohort (P = 0.0065 and 0.0074, respectively) (Fig. S4B,D) and when combining results from both the discovery and validation cohorts (P = 0.0007 and 0.0015, respectively) (Fig. 5A,  B). When combining the expression level of the two miRNAs into one prognostic signature, there was a clear separation between OS patients with low-miR-221-3p/high-miR-222-3p expression and OS patients with high-miR-221-3p/low-miR-222-3p for both the discovery and validation cohorts and in combination (P = 2.4 9 10 À8 ) ( Fig. 5C and Fig. S4E,F).

Discussion
OS is a complex disease exhibiting alterations in numerous pathways and with major cell-to-cell variation. The understanding of the etiology of OS is limited and no common molecular alteration has been identified. This designates the need for a more fundamental understanding of the underlying mechanisms leading to OS development and progression. Alterations in miRNA expression are one of the mechanisms associated with tumor development. As one miRNA potentially regulates the expression of more than 100 target genes, alterations in miRNA expression can have a vast impact on neoplastic initiation and progression. Hence, the focus of this study was the analysis of miRNA expression changes in OS. Profiling of miRNA expression in OS has been addressed in previous studies (Duan et al., 2011;Jones et al., 2012;Lulla et al., 2011;Maire et al., 2011;Namlos et al., 2012;Thayanithy et al., 2012;Won et al., 2013;. However, the limited number of samples of 2-18 OS samples in the discovery cohorts and 2-12 samples in the respective validation cohorts makes it challenging to draw reliable conclusions from the published data.
We have in this study investigated the expression of 752 miRNAs in 23 OS samples and validated significant findings in 78 additional OS samples. These are, to our knowledge, both the largest discovery and validation cohorts so far used for miRNA expression analysis in OS. The identified miRNA expression signature showed strong statistical significance in both the discovery and the validation cohorts comprising together a total of 101 OS samples. This provides robust evidence of the 29 miRNAs with an altered expression level between OS samples and OB cell lines to be truly deregulated in OS. A difference in miRNA expression between histological subtypes was not identified designating a universal role of the deregulated miRNAs in OS tumorigenesis.
A significant role of miRNAs in OS tumorigenesis was supported by the investigation of biological functions associated with the 29 deregulated miRNAs, revealing a strong correlation with cancer development and progression. This was further substantiated by the investigation of the validated miRNA/mRNA targets, identifying 361 genes with a significantly changed expression level in the OS samples. These genes were associated with multiple pathways and biological functions related to cancer development and progression. Furthermore, both a canonical pathway and biological functions were associated with bone cell pathogenesis linking the cancer-related pathways with OS tumorigenesis.
A relation between the identified pathways associated with the target genes and the TP53 pathway, which is often altered in OS, also supported a correlation between the deregulated miRNAs and OS tumorigenesis. The network analysis and target genes comprised a total of eight genes directly associated with the TP53 pathway (ABL1, CCNG1, FASLG, MDM2, TNFSF10, TP53, UBD, and UCHL1) and four genes (ABL1, MDM2, STMN1, and TP53) associated with the RB pathway. Hence, the altered miR-NAs identified in this study fit well with the current view of OS pathogenesis that the TP53 and RB pathways are the two major pathways involved in OS development and progression. Many potential functional miRNA/mRNA relationships were identified in this study. These need to be further investigated to determine their function and involvement in OS tumorigenesis.
The prognostic value of the 29 deregulated miRNAs was also investigated in function of the overall survival and time to metastasis by subdividing the OS samples into high/low expression groups. This identified two miRNAs, miR-221-3p and miR-222-3p, in which the expression level showed a significant association with time to metastasis. When combining the high/low expression scores for the two miRNAs for each OS patient, the separation became even more distinct and significant, signifying the robustness of the miR-221/ 222 signature as a prognostic biomarker. Increased expression of miR-221 and miR-222 has also been shown to be significantly associated with lymphatic metastasis in gastric, colorectal, and breast cancer (Falkenberg et al., 2013;Fu et al., 2014;Liu et al., 2012b;Sun et al., 2011). Furthermore, these two miR-NAs have been shown to be a promising biomarker in breast cancer in function of time to developing distant metastasis, although a high expression of both miR-NAs was associated with poor prognosis in breast cancer (Falkenberg et al., 2013). The discrepancy of miR-222 expression between these studies and the OS samples in this study may be explained by the ability of a miRNA to target several hundred genes, thereby promoting tumor growth in some malignancies and repressing it in others. This is further substantiated by the fact that miR-222 was shown to target MMP1, which is associated with increased aggressiveness and metastatic risk in oral tongue squamous cell carcinoma (Liu et al., 2009) (MMP1 showed a 10-fold increase in expression in the OS samples in this study). High expression of MMP1 has also been associated with the development of metastasis in colorectal and breast cancer (Liu et al., 2012a;Sunami et al., 2000). miR-221 targets among other genes BCL2L11, which exhibits pro-apoptotic properties, thereby suppressing tumor growth (Aichberger et al., 2009), as well as the gatekeeper PTEN . In different cancers, high expression of miR-221/222 has also been shown to facilitate epithelial/mesenchymal transition, a prerequisite for metastasis formation (Li et al., 2015;Stinson et al., 2011). Furthermore, modulation of miR-221 levels in OS cell lines demonstrated its role in carcinogenesis with high levels leading to increasing invasiveness and migration of the cells (Zhu et al., 2015). Hence, both the deregulation itself and the target genes of miR-221 and miR-222 indicate a prominent role of these two miRNAs in the ability of a neoplasm to spread to other organs and in the aggressiveness of specific cancer types.
Of the 29 deregulated miRNAs identified in this study, expression changes in nine miRNAs have not previously been associated with OS in miRNA profiling studies, while 20 miRNAs previously associated with OS development and progression were uncovered (Duan et al., 2011;Jones et al., 2012;Lulla et al., 2011;Maire et al., 2011;Namlos et al., 2012;Thayanithy et al., 2012;Won et al., 2013;Zhang et al., 2015). However, as OS is characterized by many different alterations, confirmation of previous findings is important for this particular disease. Furthermore, many of the formerly reported miRNA changes in OS have been identified in cohorts comprising less than 10 OS samples. Previous studies profiling miRNA expression changes in OS revealed 14 miRNAs significantly changed in these studies, which were not among the 29 miRNAs identified in our study. This may be due to differences in the cohorts, but also due to low sample number in the published studies or the use of normal bone as the control tissue, comprising a lower number than the OS cohorts (comprising 2-12 normal bone samples) (Jones et al., 2012;Namlos et al., 2012;Thayanithy et al., 2012;Won et al., 2013;Zhang et al., 2015). The cell of origin for OS has not yet been completely clarified. Both mesenchymal stem cells (MSCs) and derived osteoprogenitors such as OBs have been analyzed for their ability to form OS. Several studies have shown that OB precursors presented a higher incidence of OS compared to early mesenchymal progenitors (Abarrategi et al., 2016), supporting the model of OS originating from cells with OB commitment rather than immature MSCs. OB cell lines have therefore been used as the normal controls in this study, as well as in several previous studies (comprising 1-3 OB cell lines) (Duan et al., 2011;Lulla et al., 2011;Maire et al., 2011). Only one study has so far been using both normal bone tissue and OB cell lines as controls, but analyzed only 15 miRNAs previously identified to be differentially expressed between normal bone and OS cell lines. Nine of these were validated in primary OS samples when compared to normal bone. However, for some of these miRNAs, the authors observed opposite expression changes in OS when compared to either OB cell lines or normal bone (Namlos et al., 2012). This may be due to OS being at an intermediate state between the undifferentiated OB cells and the fully differentiated bones comprising several different cell types. Furthermore, normal bone was obtained from cancer patients suggesting a different age distribution between OS samples and bone samples, which could lead to differential miRNA expression profiles, and raise the question of the possible presence of cancer-related epigenetic alterations in these samples. These points, as well as the small number of miRNAs analyzed in primary OS, bone and OB cell lines, make it difficult to draw general conclusions from this study. Even though OB cells can be considered to be the cell of origin of OS, the use of cell lines as normal control holds the limitation of identifying miRNAs with altered expression due to differences between cell lines and human tissue and not due to the cell origin. Unfortunately, no normal bone tissue was available for the current study. Therefore, future studies should comprise pediatric OB cells as normal control as well as normal bone to obtain a more comprehensive description of the role of the identified miRNAs in OS tumorigenesis.
Several of the 24 significantly downregulated miR-NAs identified in this study have previously been identified as tumor suppressor miRNAs in other malignancies, among others miR-100-5p, miR-125b-5p, miR-127-3p, miR-370-3p, miR-335-5p, and miR-411-5p (Gao et al., 2016;He et al., 2013;Jiang et al., 2014;Li et al., 2014;Luan et al., 2016;Zhang et al., 2016). Of these, miR-100-5p, miR-127-3p, and miR-335-5p have previously been identified as downregulated in OS studies analyzing 2-18 samples (Jones et al., 2012;Maire et al., 2011;Thayanithy et al., 2012;Zhang et al., 2015). However, miR-127-3p and miR-370-3p have also been detected upregulated in two OS samples in one study . miR-125b-5p and miR-411-5p have not previously been identified having an altered expression in any OS miRNA profiling studies. miR-125b has been shown to be an important regulator of both proliferation and differentiation in different cell types (Scott et al., 2007;Sempere et al., 2004). Moreover, miR-125b has been shown to inhibit normal OB proliferation in mouse cells (Mizuno et al., 2008), indicating a role of this miRNA in bone development and thereby in OS tumorigenesis. miR-335-5p has also been identified to be important in the differentiation process of bone cells, as a positive regulator of Wnt signaling through the inhibition of Wnt antagonists (Zhang et al., 2011). Wnt signaling is elevated when osteoprecursor cells differentiate into OBs. However, Wnt signaling is inhibited when the OB cells differentiate further toward the terminal stage of bone-forming cells, which may be a functional explanation for the downregulation of miR-335-5p in OS.
Future functional studies for each identified miRNA and their target genes in OS development and progression have a strong potential for discovering novel diagnostic and prognostic markers as well as developing new treatment strategies, thereby improving the longterm survival rate for these young patients. Circulating miRNAs detected in serum and plasma of patients with cancer are being widely investigated for their potential as both diagnostic and prognostic markers of cancer, including OS (Ram Kumar et al., 2016). This might also provide a tool of diagnosing OS at an earlier time point when individuals complain repeatedly about 'growing' pains. Furthermore, miRNAs hold the promise of novel therapeutic targets, either by blocking the expression of an oncomiR or by substituting the loss of a tumor suppressor miRNA.

Conclusion
In conclusion, this study substantiates the importance of miRNA deregulation in OS and its association with OS tumorigenesis. We have identified and validated 29 deregulated miRNAs in, to our knowledge, the largest discovery and validation cohorts published, comprising a total of 101 OS samples. Both the miRNAs and their identified target genes are associated with multiple pathways and biological functions, related to cancer development and progression, as well as bone cell biology, thereby associating the deregulated miRNAs with OS tumorigenesis. Interestingly, a promising potential as prognostic biomarkers for the aggressiveness of OS was seen for two of these miRNAs (miR-221-3p and miR-222-3p).

Supporting information
Additional Supporting Information may be found online in the supporting information tab for this article:   Boxplots of the expression level of three miR-NAs associated with the metastatic potential of OS; two miRNAs associated with the chemotherapeutic response of OS. Fig. S4. Prognostic value of the expression of two miRNAs in function of time to metastasis. Table S1. Raw expression data (C q values) for 752 miRNAs analyzed with the microRNA Ready-to-Use PCR, Human panel I and panel II v2 (Exiqon, Denmark) using the OS discovery cohort. Table S2. Log 2 normalized data for the 752 miRNAs analyzed for the OS discovery cohort. Table S3. Raw expression data (C q values) for 33 miR-NAs to be validated and the five reference miRNAs analyzed with the Pick&Mix microRNA PCR Panel, 384 well Ready-to-Use (Exiqon, Denmark) using the OS validation cohort. Table S4. Log 2 normalized data for the 33 miRNAs analyzed in the OS validation cohort. Table S5. Expression of 301 upregulated and 60 downregulated genes identified to be targeted by the 29 deregulated miRNAs. Table S6. Expression of 76 miRNAs identified as differentially expressed in OS samples when compared to OB cell lines. Table S7. Previous miRNA profiling studies describing the 33 miRNAs identified in the discovery cohort. Table S8. Comparison of profiling studies investigating miRNA expression changes in OS and their expression in the OS investigation cohort in this study. Table S9. Expression of three miRNAs associated with the metastatic potential of OS and two miRNAs associated with the chemotherapeutic response of OS.