- Primary research
- Open Access
Whole exome sequencing of a single osteosarcoma case—integrative analysis with whole transcriptome RNA-seq data
Human Genomicsvolume 8, Article number: 20 (2014)
Osteosarcoma (OS) is a prevalent primary malignant bone tumour with unknown etiology. These highly metastasizing tumours are among the most frequent causes of cancer-related deaths. Thus, there is an urgent need for different markers, and with our study, we were aiming towards finding novel biomarkers for OS.
For that, we analysed the whole exome of the tumorous and non-tumour bone tissue from the same patient with OS applying next-generation sequencing. For data analysis, we used several softwares and combined the exome data with RNA-seq data from our previous study.
In the tumour exome, we found wide genomic rearrangements, which should qualify as chromotripsis—we detected almost 3,000 somatic single nucleotide variants (SNVs) and small indels and more than 2,000 copy number variants (CNVs) in different chromosomes. Furthermore, the somatic changes seem to be associated to bone tumours, whereas germline mutations to cancer in general. We confirmed the previous findings that the most significant pathway involved in OS pathogenesis is probably the WNT/β-catenin signalling pathway. Also, the IGF1/IGF2 and IGF1R homodimer signalling and TP53 (including downstream tumour suppressor gene EI24) pathways may have a role. Additionally, the mucin family genes, especially MUC4 and cell cycle controlling gene CDC27 may be considered as potential biomarkers for OS.
The genes, in which the mutations were detected, may be considered as targets for finding biomarkers for OS. As the study is based on a single case and only DNA and RNA analysis, further confirmative studies are required.
Osteosarcoma (OS) is a most prevalent primary malignant bone tumour and mostly occurs in children and adolescents—75% of patients with OS are 15 to 25 years old. The etiology is unknown; however, a genetic predisposition has been suggested ,. Reviewed in , these tumours have high potential to metastasize and are one of the most frequent causes of cancer-related deaths. The survival rate increased up to 70% after chemotherapy became available . However, no further improvements have been made in the last decades in terms of survival. Thus, the survival plateau forces scientists to look for new biomarkers (diagnostic, disease monitoring, response, resistance markers, drug targets), which could lead to, i.e. applying new therapeutic agents. While OS is rare and very heterogeneous (inter-patient, inter-tumour and intra-tumour heterogeneity), the clinical study progress is slow; thus, the preclinical studies are vital. Furthermore, finding the biomarkers and detecting the potential targets for new drugs are essential to improve the present situation.
There are several next-generation sequencing (NGS) and genome-wide association studies (GWAS) about OS, which associate different genes and pathways with pathogenesis of OS -. With whole exome sequencing (WES) and whole genome sequencing (WGS) studies, TP53, PTEN and PRB2 are found to be mutated in significant frequency . High mutation rate in TP53 has also demonstrated in OS cell lines. Additionally, deletion of CDKN2A/B locus and amplification of MDM2 were detected . With GWAS studies, a single nucleotide variant (SNV) in GRM4 was detected as potential biomarker for OS . Gene expression studies reveal that, i.e. WNT inhibitory factor (WIF1) has a loss of expression in OS cell lines ; however, we found in our previous work that the expression has increased significantly . Thus, as demonstrated, the expression pattern of WNT pathway genes in different OS cases may not be similar. When correlating the expression patterns of miRNA/mRNA pairs, miRNAs regulating TGFBR2, IRS1, PTEN and PI3K have been detected . In addition, several serine/threonine kinases (mechanistic target of rapamycin (mTOR)) or tyrosine kinases (SRC, IGF1R, PDGFR, KIT) are considered as targets in OS treatment ,,.
When observing the related pathways, the WNT/β-catenin pathway is one of the most thoroughly studied among bone malignancies. For example, the tumour growth is regulated through this pathway and the overexpression of BMP9 suppresses its activity ,. Furthermore, PI3K/AKT/mTOR signalling pathway was brought forward as a potential target for therapy, and also, pathways associated to TP53 may be altered . Hypoxia-HIF-1α-CXCR4 pathway plays a crucial role during the migration of human osteosarcoma cells . These are just a few examples—the network of associated genes and pathways is complex.
OS has a very unstable genome—it may contain aberrant number of chromosomes, and in most cases, these chromosomes display major structural abnormalities including amplification, deletions and translocations. For example, several studies have demonstrated the gain of chromosomal arms 6p, 8q and 17p in the case of OS ,. To be more precise, i.e. VEGFA amplification and LSAMP deletion have been detected in OS ,. Thus, it is suggested that genomic instability is linked to the development of this tumour ,-. Furthermore, the genomic aberrations are more frequent in metastases than in primary tumours . The genes responsible for cell cycle regulation are suggested to be associated to DNA breakage and genomic instability, i.e. CDC5L overexpression and mutations in TP53 gene are correlated to the high genomic instability in OS ,.Moreover, the chromothripsis event is characteristic to OS—it generates new fusion products. This may explain the sudden onset of OS and the complexity and heterogeneity of OS genome . All these changes make it difficult to find biomarkers suitable for targeting OS, as there are so many different subtypes.
In the present work, we analysed the whole exome of the tumorous and non-tumour bone tissue from the same patient with osteosarcoma. We used next-generation sequencing to study how the coding region of the tumour genome has altered. Additionally, we analysed together the WES genotyping and RNA expression data (from our previous RNA-seq analysis).
Materials and methods
The protocols and informed consent form used in this study were approved by the Ethical Review Committee on Human Research of the University of Tartu. The patient signed a written informed consent, which also includes the acceptance of the report to be published. A 16-year-old Caucasian male patient with an OS diagnosis was studied. In more detail, the patient became ill with complaints of pain in the left knee area. History of trauma was missing, and GP administered painkillers and vitamins. After 6 months, the patient returned to GP with complaints of pain, swelling and dysfunction in the left distal femur and knee area. The swelling line was observed in the left femoral distal region, and the area was thicker and painful to touch. No changes in skin colour were detected. The X-ray investigation showed additional shading and structural change in the distal part of the 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 Scandinavian Sarcoma Group (SSG) XIV treatment protocol. The patient responded well to the therapy—the histological analysis confirmed the necrotic tissue in tumour. After 3 month of chemotherapy, surgical removal of tumour (distal part of femoral bone with knee joint) and replacement of the knee and the lower part of the femur with megaprosthesis was performed. Pathologist confirmed that resection line was without tumour cells and OS was referred as NAS (Not Further Specified). After the patient had recovered from surgery, the SSG XIV chemotherapy treatment protocol was followed. Materials for this study were collected from the surgically removed tissue.
The genomic DNA (gDNA) was extracted from two bone samples from different locations—one sample from tumour area and another sample from the uninvolved normal bone tissue as a control. For gDNA extraction, the tissue was homogenized applying liquid nitrogen and a mortar, and after that, the PureLink Genomic DNA kit (Life Technologies Corp., Carlsbad, CA, USA) was used according to manufacturer’s protocol. The Target Seq Exome Enrichment System and SOLiD 5500 barcoded adaptors (Life Technologies Corp., Carlsbad, CA, USA) were used to prepare the libraries. The SOLiD 5500xl platform and paired-end DNA sequencing chemistry (75 bp forward and 35 bp reverse direction) were applied to sequence the samples.
The data analysis
Offline cluster was used for data processing and analysis. For bioinformatic analysis, LifeScope version 2.5 was applied. LifeScope performed colour space mapping and pairing. Tertiary analysis consisted of SNV discovery (diBayes algorithm) and detection of small indels. Hg19 (GRCh37.p13) was used as a reference, and before mapping, the multifasta file was verified in order to increase the mapping quality.
The SNVs and small indel .gff3 files were used as input in ANNOVAR software (AS; www.openbioinformatics.org/annovar/)  and Ingenuity Variant Analysis (IVA; http://www.ingenuity.com) QIAGEN, Redwood City, MD, USA) software. Applying refGene hg19, dbSNP135 and dbCOSMIC67 databases, AS annotated and predicted the effects of SNVs and small indels we detected in our study samples. AS also provides other prediction tools in order to get prediction scores (PolyPhen-2, SIFT, ljb2 etc.) -. Comparative distribution of SNVs and small indels between different samples was performed with Galaxy software bundle ,. IVA provided tools to annotate SNVs and small indels, which may be associated to cancer. The tumour and control samples were compared, and the lists for diseases, processes and pathways related to cancer were received as output.
The .bam and .bai files were used as input in CEQer software (CS) (www.ngsbicocca.org/html/ceqer.html), which is a tool for analysing copy number variants (CNVs) and loss of heterozygosity (LOH).
About the RNA-seq data analysis, please see our previous article, where we used the bone samples from the same patient .
For comparing the tumour tissue and non-tumour tissue (control tissue) from the same individual, different approaches were applied. After mapping the data to a reference genome, we used several tools to perform the tertiary analysis.
Sequencing statistics from LifeScope software
In the case of the tumour tissue, over 130 million (58%) mappable reads were in target and the enrichment fold was 48%. Eighty-five percent of the detected targets were covered over 20 times, and the average coverage was 185.5. In the case of the control tissue, over 154 million (61%) mappable reads were in target and the enrichment fold was 51%. Eighty-three percent of the detected targets were covered over 20 times, and the average coverage was 157.
SNVs, small indels and CNVs
1)Results from ANNOVAR software
Using refGene hg19 database, AS was able to annotate 37,990 SNVs and 1,484 small indels. In the case of SNVs, we considered the data reliable, if the coverage was over 20; thus, 25,914 SNVs remained. In the case of SNVs, there were 23,767 germline mutations (9,067 in homozygous form and 14,700 in heterozygous form) and 2,147 somatic mutations (in the tumour tissue—116 in homozygous form and 2,031 in heterozygous form) (Table 1, Additional file 1). Furthermore, there were 896 germline small indels (278 in homozygous form and 618 in heterozygous form) and 588 somatic indels (in the tumour tissue—177 in homozygous form and 411 in heterozygous form).
Applying dbSNP135, we were able to annotate 5,281 SNVs and 239 small indels. With dbCOSMIC67, we annotated 2,569 SNVs and 59 small indels—none of these were noted to be associated to bone cancer. Applying ljb2 database, we found 469 SNVs to potentially cause a disease (average ljb2 score over 0.918), including 31 germline mutations and 4 somatic mutations (ESX1: c.A578G/p.K193R; CDC27: c.A17G/p.E6G; TMEM120B: c.G274A/p.D92N; TMEM131: c.C3947T/p.P1316L) in homozygous form in the tumour tissue.
2) Results from Ingenuity Variant Analysis software
Altogether, 207 cancer driver variants (CD-SNVs) were found in 123 genes according to IVA (Additional file 2). Fourteen CD-SNVs potentially gain and 186 lose the gene function. Only seven SNVs may have no drastic effect on gene function in the tumour tissue. Furthermore, according to IVA, none of these 207 SNVs affect the gene functionality in the control tissue. Thirteen of the CD-SNVs were homozygous in the tumour tissue (Table 2). There were no cancer-associated homozygous mutations present in the control tissue; thus, the homozygous CD-SNVs in the tumour tissue are all somatic.
According to IVA, six cancer-associated small indels were found (Table 2). Four of them are homozygous and two are heterozygous in the tumour tissue—the effect is most probably the loss of gene function. These indels are predicted to have no effect in the control tissue.
In most of the genes brought front by IVA, one CD-SNV was found in coding region in heterozygous form. However, some of the genes have more CD-SNVs in coding regions: MUC4 had even 22, ZNF717 had 8, CTBP2 had 7 and OR4C3 had 5 CD-SNVs, whereas these were not present in the control tissue (data not shown). When observing from a slightly different angle—the gene complexes, we can see that the mucin complex has the highest significance—three genes and 27 CD-SNVs are considered (Table 3). There are also other gene complexes, which are potentially associated to cancer processes, and in different complexes, the CD-SNVs are either somatic or germline (Table 3).
In the case of cancer-associated small indels, the statistically most significant results were with complexes related to RELA gene—NFKB1-RELA and RELA-REL complexes both had p value 7.56E-4.
IVA provided the first 100 cancer-associated processes and diseases related to CD-SNVs and small indels. Seventy-three genes and 135 CD-SNVs were found associated to process named as “disorder of genitourinary system” (Table 4). These findings were present in both the tumour and control tissues. There were also two processes associated to bone “myelopoiesis of bone marrow” (associated genes NPM1, RARA) and “quantity of trabecular bone” (associated genes CREBBP, SMO)—these findings were present only in the tumour tissue. In the case of small indels, all the findings were somatic and ALK and RELA genes were associated to “outgrowth of bone marrow cells” and “inflammatory response of bone marrow-derived macrophages”, respectively.
IVA found 111 genes with 202 germline CD-SNVs associated to cancer (Table 5). Fifteen genes, which had 43 somatic CD-SNVs were associated to “bone marrow cancer and tumours”. In the case of small indel, all six genes, with a finding, are associated to cancer and the found small indels are all somatic. The disease named as “tumourigenesis of bone tumour” was associated to small indel in ALK gene and was present only in the tumour tissue.
With the osteosarcoma patient’s tumour and control tissue, WES data IVA found six pathways associated to CD-SNVs and six to cancer driver small indels (Table 6). All the mutations considered here were somatic. In the case of CD-SNVs, the statistically most significant association was between tumour and WNT/β-catenin signalling pathway. In the case of small indels, associations with different cytokine pathways were found. Also, a pathway directly linked to the bone tissue—“RANK signalling in osteoclasts” was brought front.
3) Results from CEQer software
We applied CS to analyse CNVs in tumour and non-tumour tissue exomes. Compared to the control tissue, in the tumour tissue, the loss of coding sequences was found in 6 chromosomes and 183 genes and gain of coding sequences in 4 chromosomes and 65 genes (Figure 1). The loss or gain of coding sequences was altogether in 8 chromosomes, and the most altered were chromosomes 2 and 19 (193,701 bp and 115,358 bp, respectively; Figure 2). The loss of heterozygosity was detected altogether in 68 regions in 37 genes, located in 15 different chromosomes (Additional file 3).
The integrative analysis narrows down the large list of findings from NGS data. When combining the results from WES data (AS, IVA, CS) and RNA-seq data , we found some interesting and rather logical associations, which we would like to emphasize.
SNVs, small indels and RNA expression
To reduce down the complexity of data we received from AS, we decided to perform as follows. In the case of SNV data, we observed both somatic and germline SNVs, which are homozygous in the tumour tissue and should have an effect on translation (nonsynonymous, stopgain, stoploss findings). Thus, we got 527 homozygous germline SNVs (in 392 genes) and 8 homozygous somatic SNVs (in 7 genes), which are located in genes with altered expression in the tumour tissue compared to that in the control tissue. If also considering the ljb2 database scores, seven homozygous SNVs with high disease-causing probability remained (Table 7).
In the case of small indels detected with AS, we observed the somatic and germline indels, which were homozygous in the tumour tissue. There was 52 germline and 26 somatic indels in introns of the genes, which expression pattern has also changed (data not shown). Furthermore, there was five germline and three somatic indels in exons of the genes with altered expression. Thus, we found altogether three frameshift small indels, which possibly have an effect on translation (frameshift insertions and deletion in exons) (Table 7).
In the case of homozygous cancer driver SNVs and small indels found with IVA (Table 2), only four genes have altered expression pattern in the tumour tissue compared to that in the control tissue. The mRNA expression was increased in the case of PLAT (log fold change (logFC) = 3.65, false discovery rate (FDR; corrected statistical significance) = 8.27E-27), AGTPBP1 (logFC = 0.91, FDR = 0.039) and LRRC37A3 (logFC = 1.14, FDR = 0.0072) and decreased in the case of SFB1 (logFC = −1.33, FDR = 0.0037).
CNVs, LOHs and RNA expression
When analysing the CNV results together with RNA expression results, we found that with gained copy numbers, there were altogether 22 genes, with altered expression profile—20 genes with increased and 2 genes with decreased mRNA expression. In the case of loss copy of number, 74 genes’ expression profile had changed—11 genes with increased and 63 genes with decreased mRNA expression. In Table 8, the genes with the lowest FDR values for gene expression results are presented. Here, we would emphasize that the INSR, which has copy number loss in area covering 174,552 bp has also a remarkable decrease in mRNA expression (3.36 times; FDS = 9.67E-31). However, there are also several genes with CNVs, which could be associated to cancer.
Combining the LOH and mRNA expression data, we found that in the tumour tissue, the expression of four genes with LOH has increased significantly and expression of five genes with LOH has decreased significantly (Table 9). The rest of the genes with LOHs had no significant changes in mRNA expression level, and two genes were not detected with RNA-seq (FLJ20518, MANSC4) .
For additional information, please see the supplementary material as separate files for AS, IVA and CS combined with RNA-seq data.
In this study, the exome profiles of the osteosarcoma patient’s tumour and normal bone tissue were compared. Additionally, the RNA-seq data from our previous work was used . For WES data analysis, several softwares were applied and possibly some of them are better in detecting some mutations and not so effective in detecting others. Still, we think it is more beneficial to use different approaches and we believe it is easier to follow, if we discuss separately the results gained from each software.
The ANNOVAR software annotated a large amount of genes with SNVs and small indels, applying refGene hg19 database. Over 2,700 somatic SNVs and small indels were detected specifically in the tumour tissue, from which almost 300 are homozygous. These findings are located all over the exome. This demonstrates that the changes in OS genome are not concentrated into a single or few areas but are rather distributed.
When using ljb2 database, AS detected four homozygous somatic mutations in the tumour tissue, which could potentially cause a disease. These nonsynonymous mutations were located in ESX1, CDC27, TMEM120B and TMEM131. Additionally, in the case of TMEM120B and TMEM131, the mRNA expression has decreased substantially in the tumour tissue compared to that in the control tissue ; however, further studies are needed to confirm the possible associations between found mutations and gene expression level. Available data about the possible associations between OS and these genes is very limited. In TMEM120B, a gene with an unclear function, the mutation COSM1599921 has been previously detected in glioma . The CDC27 is a gene possibly controlling the timing of mitosis and may have an important role in tumour cell division . In addition to the somatic mutation, the CDC27 had 33 heterozygous germline disease-causing mutations (nonsynonymous) (data not shown). In the case of breast cancer, the CDC27 has been demonstrated to be a promising biomarker in predicting the disease progression and prognostication . Thus, these somatic mutations may have some effect on OS pathogenesis. Especially the abundant changes in CDC27 may be important in terms of regulating OS tumour cell division.
In the tumour tissue, we detected homozygous somatic small indels causing the frameshift in five genes—EI24, ALG1L2, TIGD6, GPATCH4 and SSPO. None of these genes have previously been associated to OS, and according to our RNA-seq data, only EI24 of these five genes has altered mRNA expression—it has decreased in the tumour tissue , which could be due to the insertion in exon 9. The EI24 encodes a tumour suppressor and is an immediate-early induction target of TP53-mediated apoptosis—it binds to antiapoptotic BLC2. Furthermore, the EI24 has found to be highly mutated in the case of aggressive breast cancer and is rather associated to tumour invasiveness than development of the primary tumour -. In the present case, we found no mutations in TP53 nor was the expression altered ; thus, according to this data, we may suggest that the TP53 is functional in the tumour tissue. However, the TP53 pathway may still be suppressed due to mutated and downregulated EI24. Moreover, the aggressive nature of OS is correlated to this finding.
Appling Ingenuity Variant Analysis software, we found over 200 cancer driver variants and 93% of these possibly cause the loss of gene function. Thirteen homozygous somatic CD-SNVs were detected in different genes—RGPD3, PRDM9, FOXK1, CCZ1, PLAT, AGTPBP1, SARDH, FAH, CDC27, SBF1, LRRC37A3, ARL17A and LILRB3. The mRNA expression of PLAT, AGTPBP1 and LRRC37A3 has increased and of SFB1 has decreased significantly . We found no previous data about the associations between OS and these genes, except SBF1. With previous OS studies, another missense mutation (p.E1539K) has detected in SBF1 . SBF1 is a SET (a nuclear oncogene) binding factor 1 and may inhibit the cell division . The decreased expression in the tumour tissue may be responsible for the increased cell proliferation. Some other associations, which might be interesting—PLAT gene is important for cell migration and tissue remodelling and the overexpression might cause hyperfibrinolysis , which has not previously described in the case of OS. Two mutations in ARL17A have detected in chondrosarcoma cells . In the case of CDC27, the same mutation (p.E6G) was also brought front by AS as potentially disease causing, which is discussed above. Thus, it is highly likely that at least some of these genes participate in some level of OS pathogenesis.
Additionally, with IVA four homozygous somatic small indels were detected in the tumour tissue. These were in noncoding regions of genes CTCFL, PRR23C, CDCA7L and ALK; thus, the effect might be post-transcriptional. CTCFL is a genetic paralog of CTCF; latter is an important methylation pattern regulator. In the case of CTCF, it has previously demonstrated that in the OS tissue, the changes in its methylation pattern may also cause loss of imprinting of IGF2 and H19 genes, which further alters their expression pattern . In our OS patient’s tumour tissue, the mRNA expression of both IGF2 and H19 has increased significantly (FDR = 3.46E-15 and FDR = 0.0015, respectively) . Thus, the association may be valid here also. In PRR23C, one missense mutation (p.R190W) has detected previously in the OS tissue . ALK encodes a receptor tyrosine kinase and is rearranged, mutated or amplified in several tumours. However, in the case of OS, there are only few reports about ALK ,. In addition, two heterozygous somatic small indels were detected in DSPP and RELA exons; however, we found no previous data about these findings and associations to OS. The small indels might have an effect on the expression of these genes both pre- and post-transcriptional level; however, these suggestions need to be further studied.
According to IVA, there were several genes with more than one mutation—in MUC4, there were even 22 somatic mutations in exons and 44 in introns, although they all were heterozygous. Thus, we found MUC4 locus to be the most altered in the tumour tissue compared to that in the control tissue. This might explain why its mRNA expression in the tumour tissue has decreased (FDR = 0.012) . Mucin 4 is among major constituents of mucus, and it has demonstrated that primary bone tumours rarely express MUC4 protein , which correlates to our finding. Furthermore, with IVA, we found mucin complex (MUC2, MUC4, MUC6) to have a highest significance in OS among others. However, there are also other mucin genes (MUC16, MUC17, MUC20) with somatic heterozygous CD-SNVs. The expression pattern of all other detected mucins has not changed significantly. Thus, mucins may have a role in OS pathogenesis, but we dear not to make any further conclusions.
With IVA, there was four bone-related processes brought front only in the case of the tumour tissue—“myelopoiesis of bone marrow” (NPM1, RARA), “quantity of trabecular bone” (CREBBP, SMO), “outgrowth of bone marrow cells” (ALK) and “inflammatory response of bone marrow-derived macrophages” (RELA). Furthermore, in disease list, 16 genes with over 40 somatic variations were associated to “bone marrow cancer” and “bone tumour”; however, there were also over 200 germline CD-SNVs associated to cancer. Thus, here, we would like to emphasize that in the case of both cancer-associated processes and diseases, the ones associated with bone are somatic mutations; however, the findings possibly promoting cancer are germline mutations. This is one of the phenomena, which we would like to observe in our future studies.
The most significant pathway found with IVA was “WNT/β-catenin signalling pathway” (altered genes: CREBBP, RARA, SMO, SOX10, SOX7, TCF3). Reviewed in , the pathway is required for bone development and has demonstrated to be altered in pathogenesis of OS—overexpression of numerous WNT pathway components including WNT ligands, FZDs and LRP receptors and epigenetic silencing of the pathway inhibiting genes, i.e. WIF1. However, in our previous study, we found WNT7B and WNT11 to be downregulated and WNT2B and WNT5B upregulated; FZD4 and FZD8 upregulated and LRP8 and LRP12 downregulated and DVL3 downregulated and WIF1 and SOST upregulated. Additionally, genes with CD-SNVs—RARA, SMO and SOX7 were upregulated . Thus, our results are rather controversial to several previous studies demonstrating the WNT/β-catenin pathway to be upregulated -. However, there are also studies correlating to our findings ,. As our study is based on a single case, we dear not to conclude, why the WNT/β-catenin pathway is rather downregulated here, but we suggest the controversial results may occur due to major heterogeneity of OS. Nevertheless, the present study demonstrates that in addition to altered expression patter, the genes involved in WNT/β-catenin signalling pathway carry the CD-SNVs.
In the case of small indels, the IVA brought front the pathways associated to RELA and these are mostly cytokine signalling pathways (Table 6). Previously, it has demonstrated that interaction of IL17A and IL17AR promotes metastasis in OS cells. Furthermore, IL17 stimulates osteoclast resorption . In our previous study, we found IL17AR to be significantly upregulated . Osteoclasts are important in pathogenesis of OS—the more active they are, the more aggressive the tumour is . RELA is demonstrated to enhance the osteoclast differentiation . As IVA predicts the loss of RELA functionality (at least partially, as the small indel is heterozygous), in the present case, the OS might not have been as aggressive as it usually would.
Previously, it has demonstrated that chromothripsis event is common to early stage of OS—hundreds of genomic rearrangements will appear in a single instability event . In the present case, the CEQer software detected nearly 2,400 gain and loss events in 8 chromosomes involved, which should qualify as the chromothripsis. However, the initiating cause of this massive rearrangement is unknown, as there were no traumas or other environmental causes we are aware of.
In general, the gain of copy number should increase the mRNA expression and loss of copy number should decrease the expression . In present work, this pattern was valid in the case of 86.5% of the genes with CNVs and altered expression. One of the strongest findings here was the amount of CNVs in INSR, which expression has decreased remarkably (Table 8). The main physiological role of the insulin receptor appears to be metabolic regulation . However, together with IGF1R it forms a hybrid receptor for IGF1, latter together with IGF2 is thought to have a key role in driving the proliferation and survival of sarcoma cells . Furthermore, the growth hormone and IGF1 axis controls the growth and bone modelling/remodelling . Additionally, the IRS1, which is phosphorylated by the INSR, is important for both metabolic and mitogenic pathways . In the present case, the mRNA expression of both IGF1 and IGF2 has increased (FDR = 4.65E-35 and FDR = 3.46E-15, respectively); however, the expression of IGF1R remained the same in the tumour tissue compared to that in the control tissue . Furthermore, in IGF1R we found a heterozygous germline nonsynonymous mutation (p.G1117R) with AS, which according to ljb2 database is a disease causing (data not shown). Similarly to INSR, the mRNA expression of IRS1 is decreased in the tumour tissue compared to that in the control tissue (FDR = 2.62E-10) . Thus, in the present case it seems, the proliferation of tumour cells might be rather supported by increased effect of IGF1, IGF2 and IGF1R homodimer associations, than IGF1, IGF2 and INSR-IGF1R heterodimer associations or INSR effects on IRS1.
The loss of heterozygosity has been reported to be extensive in OS exomes . In the present case, we did not detect whole chromosome or gene region loss; however, we did detect the loss of heterozygosity in smaller regions. The genes with LOH findings and increased mRNA expression—MS4A14, DSC2, IL7R and PCNXL2 have not associated to OS previously. However, in the case of DSC2, the overexpression has demonstrated to be inversely correlated to bone metastasis-free survival . The mutations in IL7R exon 6 have been demonstrated to be present in leukaemia patients’ bone marrow samples but not associated to other solid tumours . The five genes with LOHs and decreased mRNA expression—RPS4X, RPS23, HILPDA (C7orf68), HRNR and MUC4 also do not have previous information associated to OS. Nonetheless, also the LOH analysis brought forward different genes in mucin family. In addition to MUC4, there were also other genes with LOHs but with insignificant mRNA expression changes in the tumour tissues—MUC2, MUC6 and MUC17. Thus, these results also support the idea that mucins might have a role in pathogenesis of osteosarcoma.
In summary, the present case has several characteristics previously demonstrated in OS. The wide genomic arrangements have appeared—SNVs and small indels all over the genome and CNVs in some chromosomes; and in several cases, these rearrangements may have an effect on gene expression. Furthermore, the germline mutations seem to be associated to cancer in general and somatic mutations to bone tumours. The most significant pathway was the one probably most thoroughly studied in the case of OS—the WNT/β-catenin signalling pathway. We found several genes in this pathway carrying the cancer driver variances. Additionally, the IGF1/IGF2 and IGF1R homodimer signalling might have an essential effect on OS pathogenesis. Which also needs to be emphasized is that according to our data (based on DNA and RNA studies), there is no evidence of a nonfunctional TP53; however, the TP53 pathway might be suppressed in further levels—the downregulation of EI24. In addition, with this study, we found associations between different genes and OS pathogenesis, which have not demonstrated before in earlier studies. We found the MUC4 locus to be the most altered in the tumour tissue compared to that in the control tissue; furthermore, several other mucin genes are also possibly associated to OS. The somatic mutation in CDC27 was brought front by two different data analysis softwares and might have a role in OS pathogenesis.
All genes, in which the mutations were detected, may be considered as potential targets for additional studies (i.e. functional, histopathological, clinical studies) for finding OS biomarkers. The present study brought front the WNT pathway genes, IGF1/IGF2 and IGF1R homodimer signalling pathway genes, TP53 together with EI24, MUC4 together with other mucin genes and CDC27 as potential biomarkers for OS. Finally, as this study is based on a single case and only DNA and RNA analysis, these data may not be taken as conclusive evidence and further studies are needed to confirm the present findings.
Picci P: Osteosarcoma (osteogenic sarcoma). Orphanet J Rare Dis. 2007, 2: 6-10.1186/1750-1172-2-6.
Ottaviani G, Jaffe N: The etiology of osteosarcoma. Cancer Treat Res. 2009, 152: 15-32. 10.1007/978-1-4419-0284-9_2.
Botter SM, Neri D, Fuchs B: Recent advances in osteosarcoma. Curr Opin Pharmacol. 2014, 16C: 15-23. 10.1016/j.coph.2014.02.002.
Allison DC, Carney SC, Ahlmann ER, Hendifar A, Chawla S, Fedenko A, Angeles C, Menendez LR: A meta-analysis of osteosarcoma outcomes in the modern medical era. Sarcoma. 2012, 2012: 704872-10.1155/2012/704872. http://www.ashg.org/2013meeting/programguide/files/assets/basic-html/page258.html
Kiezun A, Janeway K, Tonzi P, Mora J, Aguiar S, Mercado G, Melendez J, Garraway L, Rodriguez-Galindo C, Orkin S, Golub T, Getz G, Yunes JA: Next generation sequencing of osteosarcoma identifies the PI3K/mTOR pathway as a unifying vulnerability to be exploited for targeted therapy. In.: ASHG meeting 2013; 2013.
Kuijjer ML, Hogendoorn PC, Cleton-Jansen AM: Genome-wide analyses on high-grade osteosarcoma: making sense of a genomically most unstable tumor. Int J Cancer. 2013, 133 (11): 2512-2521.
Savage SA, Mirabello L, Wang Z, Gastier-Foster JM, Gorlick R, Khanna C, Flanagan AM, Tirabosco R, Andrulis IL, Wunder JS, Gokgoz N, Patiño-Garcia A, Sierrasesúmaga L, Lecanda F, Kurucu N, Ilhan IE, Sari N, Serra M, Hattinger C, Picci P, Spector LG, Barkauskas DA, Marina N, de Toledo SR, Petrilli AS, Amary MF, Halai D, Thomas DM, Douglass C, Meltzer PS, et al: Genome-wide association study identifies two susceptibility loci for osteosarcoma. Nat Genet. 2013, 45 (7): 799-803. 10.1038/ng.2645.
Ottaviano L, Schaefer KL, Gajewski M, Huckenbeck W, Baldus S, Rogel U, Mackintosh C, de Alava E, Myklebost O, Kresse SH, Meza-Zepeda LA, Serra M, Cleton-Jansen AM, Hogendoorn PC, Buerger H, Aigner T, Gabbert HE, Poremba C: Molecular characterization of commonly used cell lines for bone tumor research: a trans-European EuroBoNet effort. Genes Chromosomes Cancer. 2010, 49 (1): 40-51. 10.1002/gcc.20717.
Kansara M, Tsang M, Kodjabachian L, Sims NA, Trivett MK, Ehrich M, Dobrovic A, Slavin J, Choong PF, Simmons PJ, Dawid IB, Thomas DM: Wnt inhibitory factor 1 is epigenetically silenced in human osteosarcoma, and targeted disruption accelerates osteosarcomagenesis in mice. J Clin Invest. 2009, 119 (4): 837-851. 10.1172/JCI37175.
Märtson A, Kõks S, Reimann E, Prans E, Erm T, Maasalu K: Transcriptome analysis of osteosarcoma identifies suppression of wnt pathway and up-regulation of adiponectin as potential biomarker. Genomic Discovery. 2013, 1: 1-9. 10.7243/2052-7993-1-1.
Namlos HM, Meza-Zepeda LA, Baroy T, Ostensen IH, Kresse SH, Kuijjer ML, Serra M, Burger H, Cleton-Jansen AM, Myklebost O: Modulation of the osteosarcoma expression phenotype by microRNAs. PLoS One. 2012, 7 (10): e48086-10.1371/journal.pone.0048086.
Hingorani P, Zhang W, Gorlick R, Kolb EA: Inhibition of Src phosphorylation alters metastatic potential of osteosarcoma in vitro but not in vivo. Clin Cancer Res. 2009, 15 (10): 3416-3422. 10.1158/1078-0432.CCR-08-1657.
Yap TA, Arkenau HT, Camidge DR, George S, Serkova NJ, Gwyther SJ, Spratlin JL, Lal R, Spicer J, Desouza NM, Leach MO, Chick J, Poondru S, Boinpally R, Gedrich R, Brock K, Stephens A, Eckhardt SG, Kaye SB, Demetri G, Scurr M: First-in-human phase I trial of two schedules of OSI-930, a novel multikinase inhibitor, incorporating translational proof-of-mechanism studies. Clin Cancer Res. 2013, 19 (4): 909-919. 10.1158/1078-0432.CCR-12-2258.
Lv Z, Wang C, Yuan T, Liu Y, Song T, Liu Y, Chen C, Yang M, Tang Z, Shi Q, Weng Y: Bone morphogenetic protein 9 regulates tumor growth of osteosarcoma cells through the Wnt/beta-catenin pathway. Oncol Rep. 2014, 31 (2): 989-994.
Cai Y, Cai T, Chen Y: Wnt pathway in osteosarcoma, from oncogenic to therapeutic. J Cell Biochem. 2014, 115 (4): 625-631. 10.1002/jcb.24708.
Guo M, Cai C, Zhao G, Qiu X, Zhao H, Ma Q, Tian L, Li X, Hu Y, Liao B, Ma B, Fan Q: Hypoxia promotes migration and induces CXCR4 expression via HIF-1alpha activation in human osteosarcoma. PLoS One. 2014, 9 (3): e90518-10.1371/journal.pone.0090518.
Raymond A, Ayala A, Knuutila S: Conventional osteosarcoma, genetics. 2002, IARC Press, Lyon
Lau CC, Harris CP, Lu XY, Perlaky L, Gogineni S, Chintagumpala M, Hicks J, Johnson ME, Davino NA, Huvos AG, Meyers PA, Healy JH, Gorlick R, Rao PH: Frequent amplification and rearrangement of chromosomal bands 6p12-p21 and 17p11.2 in osteosarcoma. Genes Chromosomes Cancer. 2004, 39 (1): 11-21. 10.1002/gcc.10291.
Kresse SH, Ohnstad HO, Paulsen EB, Bjerkehagen B, Szuhai K, Serra M, Schaefer KL, Myklebost O, Meza-Zepeda LA: LSAMP, a novel candidate tumor suppressor gene in human osteosarcomas, identified by array comparative genomic hybridization. Genes Chromosomes Cancer. 2009, 48 (8): 679-693. 10.1002/gcc.20675.
Pasic I, Shlien A, Durbin AD, Stavropoulos DJ, Baskin B, Ray PN, Novokmet A, Malkin D: Recurrent focal copy-number changes and loss of heterozygosity implicate two noncoding RNAs and one tumor suppressor gene at chromosome 3q13.31 in osteosarcoma. Cancer Res. 2010, 70 (1): 160-171. 10.1158/0008-5472.CAN-09-1902.
Tarkkanen M, Karhu R, Kallioniemi A, Elomaa I, Kivioja AH, Nevalainen J, Bohling T, Karaharju E, Hyytinen E, Knuutila S: Gains and losses of DNA sequences in osteosarcomas by comparative genomic hybridization. Cancer Res. 1995, 55 (6): 1334-1338.
Tarkkanen M, Elomaa I, Blomqvist C, Kivioja AH, Kellokumpu-Lehtinen P, Bohling T, Valle J, Knuutila S: DNA sequence copy number increase at 8q: a potential new prognostic marker in high-grade osteosarcoma. Int J Cancer. 1999, 84 (2): 114-121. 10.1002/(SICI)1097-0215(19990420)84:2<114::AID-IJC4>3.0.CO;2-Q.
Overholtzer M, Rao PH, Favis R, Lu XY, Elowitz MB, Barany F, Ladanyi M, Gorlick R, Levine AJ: The presence of p53 mutations in human osteosarcomas correlates with high levels of genomic instability. Proc Natl Acad Sci U S A. 2003, 100 (20): 11547-11552. 10.1073/pnas.1934852100.
Yen CC, Chen WM, Chen TH, Chen WY, Chen PC, Chiou HJ, Hung GY, Wu HT, Wei CJ, Shiau CY, Wu YC, Chao TC, Tzeng CH, Chen PM, Lin CH, Chen YJ, Fletcher JA: Identification of chromosomal aberrations associated with disease progression and a novel 3q13.31 deletion involving LSAMP gene in osteosarcoma. Int J Oncol. 2009, 35 (4): 775-788.
Lu XY, Lu Y, Zhao YJ, Jaeweon K, Kang J, Xiao-Nan L, Ge G, Meyer R, Perlaky L, Hicks J, Chintagumpala M, Cai WW, Ladanyi M, Gorlick R, Lau CC, Pati D, Sheldon M, Rao PH: Cell cycle regulator gene CDC5L, a potential target for 6p12-p21 amplicon in osteosarcoma. MCR. 2008, 6 (6): 937-946. 10.1158/1541-7786.MCR-07-2115.
Stephens PJ, Greenman CD, Fu B, Yang F, Bignell GR, Mudie LJ, Pleasance ED, Lau KW, Beare D, Stebbings LA, McLaren S, Lin ML, McBride DJ, Varela I, Nik-Zainal S, Leroy C, Jia M, Menzies A, Butler AP, Teague JW, Quail MA, Burton J, Swerdlow H, Carter NP, Morsberger LA, Iacobuzio-Donahue C, Follows GA, Green AR, Flanagan AM, Stratton MR, et al: Massive genomic rearrangement acquired in a single catastrophic event during cancer development. Cell. 2011, 144 (1): 27-40. 10.1016/j.cell.2010.11.055.
Wang K, Li M, Hakonarson H: ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010, 38 (16): e164-10.1093/nar/gkq603.
Adzhubei IA, Schmidt S, Peshkin L, Ramensky VE, Gerasimova A, Bork P, Kondrashov AS, Sunyaev SR: A method and server for predicting damaging missense mutations. Nat Methods. 2010, 7 (4): 248-249. 10.1038/nmeth0410-248.
Gonzalez-Perez A, Lopez-Bigas N: Improving the assessment of the outcome of nonsynonymous SNVs with a consensus deleteriousness score. Condel Am J Human Genetics. 2011, 88 (4): 440-449. 10.1016/j.ajhg.2011.03.004.
Kumar P, Henikoff S, Ng PC: Predicting the effects of coding non-synonymous variants on protein function using the SIFT algorithm. Nat Protoc. 2009, 4 (7): 1073-1081. 10.1038/nprot.2009.86.
Blankenberg D, Von Kuster G, Coraor N, Ananda G, Lazarus R, Mangan M, Nekrutenko A, Taylor J: Galaxy: a web-based genome analysis tool for experimentalists. Current protocols in molecular biology/edited by Frederick M Ausubel [et al.] 2010, Chapter 19:Unit 19.10.1:1–21.,
Giardine B, Riemer C, Hardison RC, Burhans R, Elnitski L, Shah P, Zhang Y, Blankenberg D, Albert I, Taylor J, Miller W, Kent WJ, Nekrutenko A: Galaxy: a platform for interactive large-scale genome analysis. Genome Res. 2005, 15 (10): 1451-1455. 10.1101/gr.4086505.
Yost SE, Pastorino S, Rozenzhak S, Smith EN, Chao YS, Jiang P, Kesari S, Frazer KA, Harismendy O: High-resolution mutational profiling suggests the genetic validity of glioblastoma patient-derived pre-clinical models. PLoS One. 2013, 8 (2): e56185-10.1371/journal.pone.0056185.
Topper LM, Campbell MS, Tugendreich S, Daum JR, Burke DJ, Hieter P, Gorbsky GJ: The dephosphorylated form of the anaphase-promoting complex protein Cdc27/Apc3 concentrates on kinetochores and chromosome arms in mitosis. Cell Cycle. 2002, 1 (4): 282-292. 10.4161/cc.1.4.139.
Talvinen K, Karra H, Pitkanen R, Ahonen I, Nykanen M, Lintunen M, Soderstrom M, Kuopio T, Kronqvist P: Low cdc27 and high securin expression predict short survival for breast cancer patients. APMIS: acta pathologica, microbiologica, et immunologica Scandinavica. 2013, 121 (10): 945-953. 10.1111/apm.12110.
Zhao YG, Zhao H, Miao L, Wang L, Sun F, Zhang H: The p53-induced gene Ei24 is an essential component of the basal autophagy pathway. J Biol Chem. 2012, 287 (50): 42053-42063. 10.1074/jbc.M112.415968.
Gu Z, Flemington C, Chittenden T, Zambetti GP: ei24, a p53 response gene involved in growth suppression and apoptosis. Mol Cell Biol. 2000, 20 (1): 233-241. 10.1128/MCB.20.1.233-241.2000.
Zhao X, Ayer RE, Davis SL, Ames SJ, Florence B, Torchinsky C, Liou JS, Shen L, Spanjaard RA: Apoptosis factor EI24/PIG8 is a novel endoplasmic reticulum-localized Bcl-2-binding protein which is associated with suppression of breast cancer invasiveness. Cancer Res. 2005, 65 (6): 2125-2129. 10.1158/0008-5472.CAN-04-3377.
Joseph CG, Hwang H, Jiao Y, Wood LD, Kinde I, Wu J, Mandahl N, Luo J, Hruban RH, Diaz LA, He TC, Vogelstein B, Kinzler KW, Mertens F, Papadopoulos N: Exomic analysis of myxoid liposarcomas, synovial sarcomas, and osteosarcomas. Genes, chromosomes & cancer. 2014, 53 (1): 15-24. 10.1002/gcc.22114.
Firestein R, Cleary ML: Pseudo-phosphatase Sbf1 contains an N-terminal GEF homology domain that modulates its growth regulatory properties. J Cell Sci. 2001, 114 (Pt 16): 2921-2927.
Booth NA, Bennett B, Wijngaards G, Grieve JH: A new life-long hemorrhagic disorder due to excess plasminogen activator. Blood. 1983, 61 (2): 267-275.
Tarpey PS, Behjati S, Cooke SL, Van Loo P, Wedge DC, Pillay N, Marshall J, O’Meara S, Davies H, Nik-Zainal S, Beare D, Butler A, Gamble J, Hardy C, Hinton J, Jia MM, Jayakumar A, Jones D, Latimer C, Maddison M, Martin S, McLaren S, Menzies A, Mudie L, Raine K, Teague JW, Tubio JM, Halai D, Tirabosco R, Amary F, et al: Frequent mutation of the major cartilage collagen gene COL2A1 in chondrosarcoma. Nat Genet. 2013, 45 (8): 923-926. 10.1038/ng.2668.
Ulaner GA, Vu TH, Li T, Hu JF, Yao XM, Yang Y, Gorlick R, Meyers P, Healey J, Ladanyi M, Hoffman AR: Loss of imprinting of IGF2 and H19 in osteosarcoma is accompanied by reciprocal methylation changes of a CTCF-binding site. Hum Mol Genet. 2003, 12 (5): 535-549. 10.1093/hmg/ddg034.
Pant V, Jambhekar NA, Madur B, Shet TM, Agarwal M, Puri A, Gujral S, Banavali M, Arora B: Anaplastic large cell lymphoma (ALCL) presenting as primary bone and soft tissue sarcoma—a study of 12 cases. Ind J Pathology & Microbiology. 2007, 50 (2): 303-307.
Choy E, Hornicek F, MacConaill L, Harmon D, Tariq Z, Garraway L, Duan Z: High-throughput genotyping in osteosarcoma identifies multiple mutations in phosphoinositide-3-kinase and other oncogenes. Cancer. 2012, 118 (11): 2905-2914. 10.1002/cncr.26617.
Tirabosco R, Berisha F, Ye H, Halai D, Amary MF, Flanagan AM: Assessment of MUC4 expression in primary bone tumours. Histopathology. 2013, 63 (1): 142-143. 10.1111/his.12134.
Flores RJ, Li Y, Yu A, Shen J, Rao PH, Lau SS, Vannucci M, Lau CC, Man TK: A systems biology approach reveals common metastatic pathways in osteosarcoma. BMC Syst Biol. 2012, 6: 50-10.1186/1752-0509-6-50.
Ma Y, Ren Y, Han EQ, Li H, Chen D, Jacobs JJ, Gitelis S, O’Keefe RJ, Konttinen YT, Yin G, Li TF: Inhibition of the Wnt-beta-catenin and Notch signalling pathways sensitizes osteosarcoma cells to chemotherapy. Biochem Biophys Res Commun. 2013, 431 (2): 274-279. 10.1016/j.bbrc.2012.12.118.
Leow PC, Tian Q, Ong ZY, Yang Z, Ee PL: Antitumor activity of natural compounds, curcumin and p KF118–310, as Wnt/beta-catenin antagonists against human osteosarcoma cells. Investig New Drugs. 2010, 28 (6): 766-782. 10.1007/s10637-009-9311-z.
Cleton-Jansen AM, Anninga JK, Briaire-de Bruijn IH, Romeo S, Oosting J, Egeler RM, Gelderblom H, Taminiau AH, Hogendoorn PC: Profiling of high-grade central osteosarcoma and its putative progenitor cells identifies tumourigenic pathways. Br J Cancer. 2009, 101 (12): 2064-10.1038/sj.bjc.6605482.
Cai Y, Mohseny AB, Karperien M, Hogendoorn PC, Zhou G, Cleton-Jansen AM: Inactive Wnt/beta-catenin pathway in conventional high-grade osteosarcoma. J Pathol. 2010, 220 (1): 24-33. 10.1002/path.2628.
Van Bezooijen RL, Papapoulos SE, Lowik CW: Effect of interleukin-17 on nitric oxide production and osteoclastic bone resorption: is there dependency on nuclear factor-kappaB and receptor activator of nuclear factor kappaB (RANK)/RANK ligand signalling?. Bone. 2001, 28 (4): 378-386. 10.1016/S8756-3282(00)00457-9.
Avnet S, Longhi A, Salerno M, Halleen JM, Perut F, Granchi D, Ferrari S, Bertoni F, Giunti A, Baldini N: Increased osteoclast activity is associated with aggressiveness of osteosarcoma. Int J Oncol. 2008, 33 (6): 1231-1238.
Vaira S, Alhawagri M, Anwisye I, Kitaura H, Faccio R, Novack DV: RelA/p65 promotes osteoclast differentiation by blocking a RANKL-induced apoptotic JNK pathway in mice. J Clin Invest. 2008, 118 (6): 2088-2097.
Lee J, Pilch PF: The insulin receptor: structure, function, and signalling. Am J Physiol. 1994, 266 (2 Pt 1): C319-C334.
Kim SY, Toretsky JA, Scher D, Helman LJ: The role of IGF-1R in pediatric malignancies. Oncologist. 2009, 14 (1): 83-91. 10.1634/theoncologist.2008-0189.
Canalis E, McCarthy T, Centrella M: Growth factors and the regulation of bone remodeling. J Clin Invest. 1988, 81 (2): 277-281. 10.1172/JCI113318.
Dearth RK, Cui X, Kim HJ, Kuiatse I, Lawrence NA, Zhang X, Divisova J, Britton OL, Mohsin S, Allred DC, Hadsell DL, Lee AV: Mammary tumorigenesis and metastasis caused by overexpression of insulin receptor substrate 1 (IRS-1) or IRS-2. Mol Cell Biol. 2006, 26 (24): 9302-9314. 10.1128/MCB.00260-06.
Sanz-Pamplona R, Garcia-Garcia J, Franco S, Messeguer X, Driouch K, Oliva B, Sierra A: A taxonomy of organ-specific breast cancer metastases based on a protein-protein interaction network. Mol BioSyst. 2012, 8 (8): 2085-2096. 10.1039/c2mb25104c.
Kim MS, Chung NG, Kim MS, Yoo NJ, Lee SH: Somatic mutation of IL7R exon 6 in acute leukemias and solid cancers. Hum Pathol. 2013, 44 (4): 551-555. 10.1016/j.humpath.2012.06.017.
This work was supported by the institutional research funding IUT (IUT20-46) of the Estonian Ministry of Education and Research, by the Centre of Translational Genomics of University of Tartu (SP1GVARENG) and by the European Regional Development Fund (Centre of Translational Medicine, University of Tartu).
The authors declare that they have no competing interests.
ER participated in the research concept and design, data analysis and interpretation, writing of the article and critical revision of the article and performed the statistical analysis. SK involved in the research concept and design, data analysis and interpretation and critical revision of the article. XDH contributed in the research concept and design and carried out the collection and/or assembly of data. KM participated in the research concept and design, collection and/or assembly of data, data analysis and interpretation and critical revision of the article. AM involved in the research concept and design, collection and/or assembly of data, data analysis and interpretation and critical revision of the article. All authors read and approved the final manuscript.