Transcriptome analysis of osteosarcoma identifies suppression of wnt pathway and up-regulation of adiponectin as potential biomarker

Osteosarcoma (OS) is primary malignant bone tumour with complicated early diagnosis. There are no specific markers currently available for predicting the prognosis and chemosensitivity of OS. In present study we performed transcriptome profiling of single patient tumour tissue with RNA-seq technology. We analysed surgically removed sarcoma sample from single 16 years old male patient. Transciptome analysis was done with RNA-seq technology, bioinformatics with Lifescope and R Bioconductor. Validation experiments were done with quantitative real-time PCR (QRTPCR). After quality and coverage filtering, RNA-seq experiment resulted 29,311,899 mapped reads for sarcoma and 22,099,159 mapped reads for normal bone tissue. 65 genes were differentially expressed with FDR corrected statistical significance below 0.05. Seven genes were down-regulated and 58 genes were up-regulated in sarcoma. The most highly up-regulated gene in sarcoma was adiponectin, ADIPOQ (with adjusted p-value 5.5E-07, log2 fold change was 7.9). Many of the genes we found are related to the adipose tissue metabolism (ADIPOQ, PLIN1, FABP4) and to the Wnt signalling suppression (WIF1, SOST). We also found novel fusion transcript between the genes LMTK2 and ZSWIM5. LMTK2 is lemur tyrosine kinase 2, and it has been shown to be involved in NGF-TrkA signalling. Interestingly, studies support the involvement of LMTK2 in development of prostate cancer. ZSWIM5 is zinc finger SWIM domain protein 5 and its function is not known. Immunohistochemical analysis confirmed positive staining for adiponectin in osteosarcoma. This paper is a good illustration how transcriptome analysis can find new biomarkers and targets for complex diseases.


Introduction
Osteosarcoma (OS) is the most frequent primary bone tumour in children and adolescents [1]. It is malignant mesenchymal sarcoma and if left untreated it progresses to local and systemic disease. OS is devastating disease with poor early diagnosis and a low long-term survival rate despite significant improvements in therapy. Although OS is rare disease, its high mortality makes its overall contribution to cancer mortality high [2]. Occasionally familial involvement in OS has been observed. Therefore, genetic susceptibility plays some role in the disease development and genetic regulation seems to be important in disease progression. P53 has frequently found to be involved in OS development and it is suggested as a target for gene therapy of sarcoma [3,4]. Wnt signalling pathway is involved in the bone formation and also in sarcoma development [5][6][7]. According to current understanding activation of Wnt pathway leads to the development of OS [8]. On the other hand, in another recent study the inactivation of Wnt signalling in high-grade OS biopsies was shown [9]. Activation of Wnt signalling was shown to inhibit proliferation of cells and promote osteogenic differentiation in osteosarcoma cell lines [9]. Therefore, the results on the Wnt signalling involvement in OS development are somewhat controversial. Most likely, this is caused by the differences in models and methods applied in different studies.
In many cases, new fusion transcripts in malignant samples have been described and this reflects genomic complexity of the OS [10]. OS is genomically unstable tumour what makes its chemical and biological therapy complicated [11]. This cancer is very aggressive and develops distant metastases in 45% even with the intensive neoadjuvant treatment with multiple chemotherapeutic drugs [12]. Survival of patients has increased only slightly and no other significant improvement has been made since decades [13][14][15]. Therefore, new drugs are seriously needed for OS. Better understanding of the OS biology and genomics can give us new information and potential targets for drug development.
In present study we performed transcriptome profiling of single patient tumour tissue with RNA-seq technology. We compared tumorous tissue with patient's own normal bone tissue in order to find specific profile in transcriptome characteristic for the sarcoma. RNA-seq is most comprehensive analytical tool for transcriptome profiling and in addition to conventional gene expression profiling it allows finding of transcriptional aberrations, fusions, what can be potentially important for the tumorigenesis. Although this study is exploratory and its results need to be confirmed by further studies, we were able to describe the genetic network behind OS and new potential biomarker.

Subject
Patient, 16 years old male, became ill with complaints of pain in left knee area. History of trauma was missing and GP administered painkillers and vitamins. After 6 months he returned to GP with complaints of pain, swelling and dysfunction in left distal femur and knee area. The swelling line observed in left femoral distal region, the area was thicker and painful during the touch. No changes in skin color were detected. The X-ray investigation showed additional shading and structural change in distal part of left femur. For detailed investigation the MRI was performed and as a result malignant process was suspected. Patient was hospitalized and bone biopsy was taken for histological investigation. The diagnosis of osteosarcoma was confirmed.
Chemotherapy for osteosarcoma started by SSG XIV treatment protocol. After 3 month of chemotherapy surgical removal of tumor (distal part of femoral bone with knee joint) and replacement of knee and lower part of femur with megaprosthesis was performed. Materials for this study were collected from surgically removed tissue.

RNA extraction and RNA-seq
RNA was extracted by using the Trizol combined with columns (RNeasy Mini kit, Qiagen Sciences, Maryland, USA). Up to 100mg of tissue was used for extractions. In brief, mortar was washed with DEPC water (about 20 min) and bone sample was placed into mortar. Liquid nitrogen was added three times. Each time the nitrogen was allowed to boil off and finally, bone was grinded into a powder. Trizol was added to the homogenized bone tissue and the mix was transferred into a polypropylene tube, followed by centrifugation. Supernatant was transferred into a new tube and RNA extraction with RNeasy Mini kit was followed.
The quality of total RNA was evaluated with Agilent 2100 Bioanalyzer and RNA 6000 Nano kit (Agilent Technologies Inc., California, USA). The average RIN number for the samples was at least 7.50 ng of total RNA was amplified by applying Ovation RNA-Seq System V2 (NuGen, Emeryville, CA, USA) after which the resulting cDNAs were pooled in equal amount and the pool was used to prepare the DNA fragment library with SOLiD System chemistry (Life Technologies Corp., Carlsbad, CA, USA). Sequencing was performed using SOLiD 5500xl platform and paired-end DNA sequencing chemistry (Life Technologies Corp., Carlsbad, CA, USA). 75 bp from forward and 35 bp from reverse direction was sequenced and altogether at least 40 million reads per sample was received, which is sufficient for the evaluation of expression pattern of the transcriptome.

Bioinformatics analysis
For bioinformatics analysis Lifescope software was used. This is proprietary software developed by Life Technologies and allows complete analysis of sequencing data from mapping to the annotation and structural analysis. Briefly, raw data were mapped to reference genome and list of gene counts was generated. In addition, exon junctions and gene fusions were detected to find alternative transcriptional modifications.
For statistical analysis DeSeq package for R was used [16]. DeSeq is a package what is specifically developed for RNA-seq or other counts data. The package DeSeq provides methods to test for differential expression by use the negative binomial distribution and a shrinkage estimator for the distribution's variance [16]. DeSeq also allows performing analysis of the data without replicates. In such case, the dispersion is estimated from comparing their counts across conditions as replacement for a proper estimate of the variance across the replicates. Package performs sample comparison and also adjusts P-value to overcome multiple testing problem. DeSeq package uses Benjamini-Hochberg procedure, which controls for false discovery rate (FDR) [17].

Quantitative real-time PCR
For gene-expression analysis first strand cDNA was synthesized by using random primers and SuperScript TM III Reverse Transcriptase (Invitrogen, USA). Gene expression levels of the selected 11 genes were analyzed by applying TaqMan Assay-On-Demand (FAM-labeled MGB-probe) gene expression assay mix (20x, Applied Biosystems, Foster City, CA, USA) and TaqMan Universal PCR Master Mix (Applied Biosystems, Foster City, CA, USA). For quantitative real-time PCR (qRT-PCR) analysis the ABI PRISM 7900HT Fast Real-Time PCR System equipment (PE Applied Biosystems, USA) and the ABI PRISM 7900 SDS 2.2.2 Software were used. Every reaction was made in four replicates to minimize possible errors. All reactions were performed in a final volume of 10 μl, using 50-100 ng of cDNA. Data is presented as 2 -ΔCT calculated in relation to the HPRT1. For HPRT1 analysis, primer for exon 6 (5'-GACTTTGCTTTCCTTGGTCAGG-3') and primer for exon 7 (5'-AGTCTGGCTTATATCCAACACTTCG-3') and VIC-labeled MGB probe was used. All the assays and respective genes used in RT-PCR are listed in (Table 1).

Identification of differentially expressed genes between normal bone and osteosarcoma
After quality and coverage filtering, RNA-seq experiment yielded 29,311,899 mapped reads for sarcoma sample and 22,099,159 mapped reads for normal bone tissue. We had annotated 24,922 genes and 21,825 genes were expressed in one or another sample. Statistical analysis with DeSeq identified that 65 genes were differentially expressed with multiple testing corrected statistical significance below 0.05 ( Table 2). 7 genes were up-regulated in normal tissue and 58 genes were up-regulated in sarcoma. The most highly upregulated gene by statistical significance was adiponectin, ADIPOQ (with adjusted p-value 5.5E -07 ). Adiponectin was almost the most activated gene by fold change, log2 fold change was 7.9, what is about 250 times over-expressed compared to control tissue. By fold change the most up-regulated gene was triadin, TRDN. TRDN has log2 fold change 8.2 ( Table 2).
We made two observations by looking the list of the differentially expressed genes. First, very few genes were overexpressed in the healthy bone compared with sarcoma. This means, sarcoma produces very actively (in high levels) large number of genes. And second, many of the genes we found are related to the adipose tissue metabolism (adiponectin, perilipin 1, fatty acid binding protein 4). We also found typical genes for bone tissue (e.g. Wnt inhibitory factor 1, WIF1 and sclerostin, SOST), but quite many of the genes are characteristic for adipose tissue. There are also genes participating in immune response. WIF1 is a gene what is involved in different malignancies and it regulates cell cycle [18,19]. In recent study its involvement in osteosarcoma was shown [20].
In addition to the counts data, RNAseq enables detection of potential fusion transcripts. In sarcoma sample, we found novel fusion transcript between the genes LMTK2 and ZSWIM5 (Table 3, Figure 1). LMTK2 is lemur tyrosine kinase 2, and it has been shown to be involved in NGF-TrkA signalling [21]. Interestingly, several studies support the involvement of LMTK2 in development of prostate cancer [22][23][24]. ZSWIM5 is zinc finger SWIM domain protein 5 and its function is not known very well.

Confirmation of differentially expressed gene with quantitative RT-PCR
We tested 11 genes with most significant differences from RNAseq data. Quantitative real-time PCR confirmed differences with almost all tested 11 genes (Figure 2). Thus, in case of PANX3 we found significant down-regulation in sarcoma samples and all other genes were up-regulated in sarcoma. Only in case of TRDN we didn't find difference between normal bone and sarcoma sample. This is peculiar, as the difference based in RNAseq analysis was very large. As RNAseq gives information also based on exons, we checked the situation with TRDN. Indeed, we found that in sarcoma sample for exon 1 in TRDN genes there were 361 reads (RPKM 16.96). There were reads also matching to other exons (RPKMs around 3…10). In normal bone sample, no reads were mapped to any exons of TRDN gene. As the real-time PCR assay for TRDN (Hs00952568_m1) detects only junction between 2 and 3 exons, this may explain the discrepancy between RNAseq and RT-PCR data. Therefore, in case of TRDN we have alterative splicing between test samples and this generates some confusion. TRDN is known to have several different transcripts and different alternative forms exist [25]. However, as the function of TRDN is not known, this finding should be taken cautiously.
All genes with statistically significant gene expression levels were used for clustered analysis. This analysis revealed clear differences in gene expression profiles in analysed two samples (sarcoma and normal bone, heatmap in (Figure 3). This indicates correct sample preparation and grouping.

Functional annotation of the RNA-seq data
Functional annotation of RNAseq data was performed with the IPA system. We filtered data and for log2FC values cutoff 4 was used. Most significantly enriched network was network called "Hematological Disease, Metabolic Disease, Endocrine System Development and Function" (Figure 4). This means that the genes with the largest differences are related to the oncogenesis and endocrine regulation.

Confirmation of adiponectin expression with immunohistochemistry
We analysed histological sample of osteosarcoma with adiponectin antibody and found strong staining in the osteoid, matrix and endothelial cells (Figure 5). Therefore, adiponectin expression in osteosarcoma was also confirmed in protein level.

Discussion
Main finding of our study is that compared to normal bone    Table 3. Description of the fusion transcript we found in osteosarcoma sample during RNAseq analysis.
tissue, osteosarcoma expresses quite many specific transcripts. These transcripts can be taken as potential biomarkers or new drug targets for the osteosarcoma. The genes found in our study, seem to have involvement in development of OS or in other cancers. However, more information about their roles in OS development is needed. As this study is also exploratory additional studies should be performed in order to test their validity. We try to analyse the list of genes in the light of more recent findings.
Adiponectin (ADIPOQ) seems to be good biomarker for osteosarcoma studies. While ADIPOQ has not been related to OS in earlier studies, there are some studies where the role of ADIPOQ in cancer development has been addressed. Namely, one very recent meta-study revealed the association of ADIPOQ with risk of cancer development [26]. In addition, adiponectin is able to induce motility of prostate cancer and chondrosarcoma cells [27,28]. Moreover, adiponectin is closely associated with bone remodelling and bone   pathologies [29,30]. Serum adiponectin levels have shown to be inversely associated with bone mineral density [31]. Our immunohistological analysis confirmed high expression of adiponectin in osteoid. Therefore, adiponectin could be useful marker for OS patients as well.
WIF1 is a gene with very clear involvement in oncogenesis. Generally, WIF1 is silenced by epigenetic inactivation during cancer development, but this is not limited to osteosarcoma [19,20]. WIF1 is inhibitor of Wnt signalling and therefore it is involved in development and remodelling of different tissues [32]. What is somewhat controversial is our finding of increased expression of WIF1 in OS patients because it supposed to be silenced or down-regulated in cancer samples [20]. In our study we clearly detected up-regulation of WIF1 in sarcoma patient and this may indicate suppression of Wnt signalling in our sample. This hypothesis needs to get more confirmation with additional studies. On the other hand, increased sclerostin (SOST) expression supports the idea that Wnt signalling is suppressed in the patients in our study. Namely, increased SOST is related to the inhibition of Wnt signalling [33]. Therefore, increased SOST level fits very well with increased WIF1 what decreases activity of Wnt signalling.
SOST seems to be antagonist for Wnt signalling and therefore, up-regulation of WIF1 and SOST seems to fit with each other [33][34][35]. This means that our finding generates functional network describing the restructuring activity in the bone during tumorigenesis. Down-regulation of Wnt signalling pathway in osteosarcoma is also described in some earlier studies. Earlier studies have similarly found clear downregulation of Wnt signalling in 25 high-grade prechemotherapy biopsies [36]. In many studies the comparisons are not made between healthy bone and osteosarcoma, therefore not all published gene expression-profiling studies are comparable with ours [37,38]. Therefore, some controversies exist about the involvement of Wnt signalling pathway in OS development.
Interestingly, in OS sample we found quite many genes related to the adipocyte functioning. Perilipin 1 or PLIN1 is one example. This gene produces protein what coats the lipid storage droplets in adipocytes and protects them until they can be broken down by lipase [39,40]. Similarly, FABP4 is fatty acid binding protein 4 and is described to be specific for adipocytes [41]. Two explanations are possible. This finding can reflect that during tumorigenesis osteoblasts differentiate to adipocytes. On the other hand, this can be also normal phenotype of osteosarcoma cells. According to literature, both variants have been described. Ultrastructural analysis of osteosarcoma cell described lipid droplets in all types of osteosarcoma cells [42]. Osteoblasts can trans-differentiate to  adipocytes by inhibition of gap-junctional communication or by hormonal stimulation [43,44]. Taken together, adipocyte specific transcripts in osteosarcoma sample are quite common and this finding is in accordance of the literature.
Our study has some very clear limitations. For instance, sample size was very small in our study. Analysis based on single sample is very exploratory and does not provide conclusive evidence. We are fully aware, that the results presented here, should be verified in independent, large sample.

Conclusions
In conclusion, our whole transcriptome profiling of osteosarcoma identified down regulation of Wnt signalling pathway and adiponectin as potential biomarker for osteosarcoma relevant for the tumour progression. These changes can potentially explain molecular nature of the osteosarcoma.

Competing interests
The authors declare that they have no competing interests.