Expression profiling in canine osteosarcoma: identification of biomarkers and pathways associated with outcome

Osteosarcoma (OSA) is the most common malignant primary bone tumor of children and accounts for roughly 5% of all childhood cancers in the United States. Characteristically, OSA is found in the metaphyseal regions of long bones in the appendicular skeleton. More than 15% of patients present with clinically detectable pulmonary metastases and it is estimated that 80% or more have micrometastases at presentation [1]. Advances in treatment such as multi-agent chemotherapy have improved prognosis over the last several decades with five-year survival rates up to 70%. Despite these advances, patients that present with metastases or those whose tumors are refractory to neoadjuvant chemotherapy continue to have a poor prognosis [1]. This suggests that within the same histologic type of tumor, different genetic mechanisms may be operating, altering response to chemotherapy and metastatic capability in some tumors.

Osteosarcoma is also the most common primary bone malignancy in dogs. The majority of these tumors occur in the appendicular skeleton of middle-aged large and giant breeds. Roughly 10,000 new cases of OSA are identified in dogs annually. Standard treatment involves amputation or limb-sparing surgery followed by adjuvant chemotherapy with doxorubicin, a platinum-based drug, or a combination of the two drug types [2]. Median disease-free interval (DFI) following amputation has ranged from 165 to 470 days depending on adjuvant chemotherapy protocol and study size [3–7]. Median survival time in dogs undergoing amputation alone ranges from 134 to 175 days [3–7]. Like their human companions, pulmonary metastases are typically the cause of terminal morbidity. It has been suggested that up to 90% of canine patients may present with microscopic metastases that are undetectable via routine imaging [2]. The high variability in DFI suggests that canine OSA exhibits variable metastatic capability, rate and/or resistance to adjuvant chemotherapy, similar to the disease in humans.

Canine OSA shares many features with human OSA, making dogs a valuable comparative model. Pet dogs develop OSA primarily in the metaphyseal regions of long bones, as do human patients. The lesions are histologically identical [2]. The similarities between the molecular characteristics of human and canine OSA have been established (see [8] for review). Furthermore, Thomas and colleagues recently demonstrated that some of the same genetic aberrations identified in human OSA are also seen in canine OSA with both breed-dependent and independent associations [9]. Among the genetic changes identified, Wilms tumor 1 (WT1), tumor protein p53 (TP53), cyclin-dependent kinase inhibitor 2A (CDKN2A), phosphatase and tensin homolog (PTEN) and retinoblastoma 1 (RB1) tumor suppressors as well as v-myc myelocytomatosis viral oncogene homolog (MYC) and v-kit Hardy-Zuckerman 4 feline sarcoma viral oncogene homolog (KIT) oncogenes were shown to be affected by cytogenetic abnormalities in 76% of their samples [9]. Similarly, comparative analysis of gene expression profiles in human and canine OSA determined that the diseases were indistinguishable by hierarchical clustering [10]. Treatment and chemotherapeutic regimens are also similar with the notable exception that most amputee dogs do not undergo neoadjuvant chemotherapy, so tumors collected at the time of amputation are naïve to drugs. Dogs also provide a valuable model system in that their tumors arise "naturally," they share an environment with humans, and they metabolize drugs at a similar rate. Finally, dogs are more genetically diverse than mouse model systems and share more genetic homology with humans than mice [8]. Thus, genetic prognostic screening in dogs has strong potential applicability to the human disease [11].

In recent years, it has become clear that the tumor microenvironment plays a strong role in metastatic events even if metastatic subclones are only a small proportion of tumor cells [12, 13]. For example, van de Vijver and colleagues demonstrated that gene expression analysis of primary tumors can divide breast cancer patients into "good" and "poor" prognostic groups based on the tumors' intrinsic metastatic ability [14]. Thus, gene expression profiles of primary tumors provide information about metastatic potential and patient prognosis even if distant disease is not detectable or present.

Gene expression analysis of primary tumors can also elucidate novel chemotherapeutic targets by defining individual gene changes and/or whole pathway derangements [15, 16]. Identification of such differences between "good" and "poor" prognostic groups in OSA will allow for more personalized treatment of disease based on an individual's tumor expression profile.

The current study utilized Affymetrix GeneChip Canine Genome 2.0 arrays to explore differences in gene expression between primary OSA tumors taken from client dogs with a DFI of less than 100 days ("poor responders") and those with a DFI greater than 300 days ("good responders") following definitive treatment and chemotherapy. Individual genes with significant changes in expression were validated using qRT-PCR and explored for their ability to correctly classify good and poor responders using four different machine learning schemes. A broader, systems approach was used to identify changes in groups of interacting genes or pathways that may contribute to metastatic progression and resistance to therapy. We have found evidence of altered expression of several genes and pathways and have verified that the Hedgehog signaling pathway is comparatively downregulated in the poor responding group.

Samples were freeze-fractured, homogenized, extracted with Trizol reagent (Invitrogen, Carlsbad, CA, USA) and purified with RNeasy clean up (Qiagen, Valencia, CA, USA) per the manufacturers' protocols. Resultant RNA was quantified via spectrophotometry and assayed for quality on Agilent (Santa Clara, CA, USA) and Bio-Rad (Hercules, CA, USA) bioanalyzers at the Rocky Mountain Regional Center for Excellence (RMRCE) Genomics Core at CSU. Only samples exhibiting minimal degradation as evidenced by RNA Integrity Numbers (RIN) greater than 8 were used for microarrays.

Eight samples from each DFI cohort were selected and array analysis with GeneChip Canine 2.0 Genome Arrays (Affymetrix, Santa Clara, CA, USA) was performed in two batches (batch 1, n = 6; batch 2, n = 10) at CSU's RMRCE Genomics Core per Affymetrix protocols. One sample was removed from analysis after data collection based upon pathologist review and review of hospital records that determined the sample was not OSA but hyperreactive osteoid tissue. Briefly, the One-Cycle Target Labeling and Control Reagents package (Affymetrix, Santa Clara, CA, USA) was used to synthesize cDNA from total RNA spiked with prokaryotic Poly-A RNA as a control. The Sample Cleanup Module (Affymetrix, Santa Clara, CA, USA) was used to purify the cDNA which was then used for synthesis of biotin-labeled cRNA. cRNA was purified, quantified and fragmented before hybridization with the GeneChips. Hybridized chips were washed, stained using streptavidin-conjugated phycoerythrin dye (Invitrogen, Carlsbad, CA, USA) and enhanced with biotinylated goat anti-streptavidin antibody (Vector Laboratories, Burlingame, CA, USA) using an Affymetrix GeneChip Fluidics Station 450 and Genechip Operating Software. The Affymetrix GeneChip scanner 3000 was used to acquire images.

Microarray data was preprocessed using Probe Logarithmic Intensity Error (PLIER) estimation [19] and Robust Multichip Average (RMA) [20] algorithms with log₂ transformations. PLIER and RMA methods were compared as part of the data discovery. A standard unpaired 2-tailed t-test with a false discovery rate correction for multiple comparisons was used. Uncorrected p-values were used to rank probesets. CIMminer was used to generate clustered images of the data with the following parameters: unsupervised clustering on both axes, average linkage and Euclidean distance [21]. Microarray data has been made available through NCBI's Gene Expression Omnibus (GEO) and can be accessed via accession number GSE24251.

The pathway analysis pipeline used in this study has been previously described [28]. Briefly, the University of Pittsburgh Gene Expression Data Analysis suite (caGEDA) [29] with a standard J5 metric, a threshold of 4 and a jackknife of 4 was used to select unique genes for pathway analysis following both PLIER and RMA preprocessing. The DAVID Gene ID conversion tool was used to link canine identifiers to their human counterparts [30, 31] and identifiers absent from the database were hand-annotated by BLAST and BLAT comparisons of the target sequence; GeneGo MetaCore was used to assign functional pathways. Pathways were assigned independently to PLIER and RMA preprocessed data and the resulting pathways were compared.

WEKA software was used to generate classification models to test the analytical value of qRT-PCR-derived expression changes [32]. Classification models were built using a Support Vector Machine (SVM), a J48 decision tree, and logistic regression [33]. Models were generated with the full (n = 20) data set and tested for sensitivity and specificity using stratified tenfold cross-validation. Tenfold cross-validation randomly selects 90% of the data for training the model, and uses the other 10% of the data to test the model. The process is repeated ten times and the ten model error rates are averaged to compute an overall error rate.

In this study, we analyzed gene expression in chemotherapy-naïve primary OSA tumors from 20 dogs with the aim of identifying a gene signature of aggressive metastasis and/or resistance to chemotherapy following definitive treatment with limb amputation and adjuvant therapy with doxorubicin and/or a platinum drug. The purpose of this aim was 3-fold. First, it provides a basis for development of a prognostic screen; such a tool would be of great value to the clinicians diagnosing and treating the more than 8,000 new cases of canine OSA every year. Additionally, pet owners would benefit greatly from a more accurate projected survival time when weighing their dog's quality of life and their own monetary obligations in treatment decisions. Secondly, analysis of gene signatures may allow elucidation of single genes or genetic pathways that may be manipulated for treatment purposes. Finally, dogs are an excellent model for human OSA and identification of markers and pathways leading to disease progression and resistance to therapy in dogs may be translated to the pediatric clinical setting to improve prognosis and treatment of human OSA.

We utilized qRT-PCR to confirm the differential expression of eleven genes between primary OSA from good (DFI > 300 days) and poor (DFI < 100 days) responding dogs (Figs. 2 and 5). Transcriptional profiles of an additional 28 genes selected from fold change analysis of the microarray data were assessed via qRT-PCR and, although differential expression was observed in many, significance criteria for the qRT-PCR analysis were not met (data not shown). Nineteen of these qRT-PCR targets were selected for analysis before pathologist review identified one of the "tumors" as hyperplastic tissue without neoplasia. Removal of that sample and subsequent reprocessing of microarray data removed some of these targets from the RMA and PLIER fold change lists. Despite their failure to achieve significance in our fold-change analysis, the qRT-PCR data for these targets still correlates strongly with the array data on a sample-by-sample basis. Two of these genes, IMP1 and AGTR1, were verified as differentially expressed by qRT-PCR analysis (Figure 2), possibly due to the increased sample number used in each group for this analysis. From the additional 28 genes that failed to show statistically significant differences by qRT-PCR, eleven of these gene targets were selected by the fold change analysis of either RMA or PLIER processed data shown in Figure 1. The failure of these gene targets to reach significance in the qRT-PCR analysis may reflect the variability in microarray preprocessing algorithms as well as differences in expression values based on primer design as primers used in this study were not designed to align with Affymetrix probe locations. In addition, since our qRT-PCR analysis used a larger number of samples than the microarray study, some of the microarray hits may have been false positives that have now been removed from the list of putative biomarkers thanks to the qRT-PCR analysis.

Although these genes were primarily assessed by qRT-PCR for their prognostic potential, they may also have functional roles in metastatic progression and resistance to chemotherapy. IMP1 (6.93 fold up-regulated in the poor-responders), also known as IGF2BP1 and not to be confused with the family of IGF binding proteins, is a member of a family of three oncofetal proteins (IMP1-3) whose function is to bind and regulate mRNA stability in the cytoplasm during development. IMP1 expression is stimulated by Wnt/β-catenin signaling and has many regulatory targets, some of which are implicated in cancer: stabilization of c-myc [34, 35] and CD44 mRNAs [36], translational suppression of IGF2 [37], and localization of β-actin mRNA to sites of actin polymerization [38]. These targets can affect cell growth and survival as well as metastatic mechanisms such as invadopodia formation and cell adhesion [39]. IMP1 over-expression has been associated with poor prognosis in numerous cancer types including human ovarian and colorectal carcinomas [39, 40].

IGF2 (15.4 fold down-regulated in the poor-responders) has been shown to be down-regulated in response to IMP1 as well as to hedgehog pathway inhibition and the observed alterations in these factors/pathways may account for some down-regulation of IGF2 [41]. Additionally, IGF2 expression is modulated by numerous other factors including parathyroid hormone (PTH), cortisol, and transforming growth factor beta (TGF-β) [42]. Finally, our pathway analysis shows reduction in PTH related protein (PTHrP) and subsequent modulation of the PTH pathway suggesting IGF2 may be comparatively under-expressed in poor-responders due to decreased PTHrP expression in that cohort. It is important to note that the observed down-regulation of IGF2 (and all other genes discussed here) is relative between cohorts and that the mRNA was expressed in all samples, but to a lesser degree in poor-responders.

FBP1 (5.94 fold up-regulated in the poor-responders) is involved in gluconeogenesis and is expressed in the liver and, to a lesser extent, most other cell types. Its action opposes that of phosphofructokinase and its expression can lead to increased cellular glutathione and an apoptosis-resistant phenotype [43]. Bigl and colleagues examined FBP1 expression in several types of breast cancer and found it to be up-regulated in invasive lobular carcinoma when compared to normal tissue but down-regulated in other tumor types suggesting a variable role depending on tumor type [44].

ADHFE1 is an iron-activated alcohol dehydrogenase with widely conserved motifs that is found in multiple tissue types. It has been shown to oxidize gamma-hydroxybutyrate and is 3.50 fold down-regulated in the poor-responders [45]. CCDC3, also down-regulated in the poor-responders (7.10 fold), encodes a 270 amino acid protein with a putative coiled-coil domain near the C-terminus. Recent reports indicate that this protein is secreted by both adipocytes and endothelial cells and is under both hormonal and nutritional control [46]. Interestingly, CCDC3 was identified as a factor contributing to ifosfamide resistance in a mouse xenograft model using human OSA cell lines. Bruheim and colleagues reported a 40-fold down-regulation of this gene in resistant tumor cells[47]. As none of the dogs in the current study received ifosfamide, this gene may contribute generally to both metastasis and chemotherapeutic resistance.

Chioni et al. recently elucidated a role for SCN1B in cellular adhesion and migration in breast cancer cell lines. Their mildly metastatic cell line demonstrated increased expression of SCN1B compared to the highly metastatic cell line; furthermore, siRNA-mediated knockdown of SCN1B decreased adhesion and increased migration in the mildly metastatic line [48]. Our observed 3.70 fold down-regulation of SCN1B in the poor-responders indicates that the tumor environment may become more pro-migratory due to reduced expression of this factor.

The putative tumor suppressor gene NDRG2 (4.57 fold down-regulated in the poor-responders) is expressed in an inverse relationship to proliferation in normal tissues and has been observed to be down-regulated in numerous tumor types, especially in response to myc oncogene expression (See [49] for review). Recent cytogenetic analysis of canine OSA revealed breed independent myc amplification in 40% of the cases, suggesting this is a common chromosomal aberration in both canine [9] and human OSA [50]. Tepel and colleagues demonstrated epigenetic promoter modifications as a mechanism for suppression of this gene in glioblastoma [51]. Recent evidence has identified numerous mechanisms by which NDRG2 acts as a tumor suppressor and invasion attenuator: anti-proliferative suppression of AP-1 in colorectal carcinoma [52], anti-invasive suppression of NF-κB in fibrosarcoma and melanoma cell lines [53], pro-apoptotic involvement in the p53 pathway [54], and reduction in invasion and intracellular β-catenin in NDRG2-transfected cell lines [55]. Kim and colleagues demonstrated that NDRG2 expression decreases with increasing tumor stage in colon carcinoma, indicating that this may be an excellent marker for molecular staging [55]. Furthermore, recent studies have shown that the myc oncogene stimulates mitochondrial glutaminolysis resulting in reprogramming of mitochondrial metabolism to depend on glutamine catabolism to sustain cellular viability [56]. In support of this hypothesis, our pathway analysis associated both upregulation of the myc oncogene and alterations in mitochondrial oxidative phosphorylation with poor outcome.

To identify the prognostic potential of these genes, we built several classification models to identify genes with the most predictive power. Of the four models tested, all classified the samples with 100% accuracy when the model was built from all 20 samples. However, when stratified cross-validation was used, the two SVMs and the linear regression model were 90% correct whereas the J48 decision tree was only 75% correct. These stratified cross-validation results are generally thought to more accurately reflect results in subsequent applications of the model. The two SVM models classified with the same success rate regardless of whether built with all fifteen genes or the three most heavily weighted contributors, suggesting that CCDC3, ADHFE1 and FBP1 are highly predictive in this data set and are likely to be robust classifiers in future OSA studies.

While biomarker identification can be successful using traditional fold change methodology, as evidenced by our gene hits above, understanding of the processes of metastasis and chemoresistance can be furthered by all-inclusive pathway analysis. Thus, to eliminate some of the arbitrary nature of traditional fold change analysis, we also examined our microarray data via this methodology. Over 4,000 probesets were selected from microarray data using the J5 metric, annotated and converted to human identifiers using public-access tools including DAVID, and assigned to pathways with the GeneGo MetaCore program. This program assigns pathway significance based upon the number of genes represented within a pathway and the direction of change. The overwhelming benefit to this methodology is that change in a single gene will be ignored unless related genes also demonstrate altered expression. Thus, the downstream impact of chip anomalies, probeset inefficiencies and differences in preprocessing algorithms can be dramatically reduced. This type of analysis allows integration of the typical microarray methodology examining highly expressed genes with the systems biology approach of examining large numbers of genes, some of which may be expressed only at low levels despite their importance to a given pathway [28]. One pertinent example is PTEN deletion which was identified as a chromosomal aberration in 40% of canine OSA [9]. However, loss of PTEN was not detected in our fold change analysis, but was associated with poor outcome by pathway analysis (Additional File 2).

Although discussing each of the modulated pathways is beyond the scope of this study, some notable generalizations can be addressed. Cell adhesion and cytoskeletal remodeling are both strongly represented in pathways we have identified as significantly altered between our two cohorts; this suggests that the aggressiveness of tumor cells with regard to these two elements of metastasis may be just as important as chemoresistance mechanisms in this population. Bone-related developmental and immune-response pathways are also represented, much as one would expect in these osteoblastic/osteolytic tumors. Finally, cAMP/PKA signaling pathways also have strong representation in these analyses. Similar alterations in cAMP/PKA signaling with upregulation of the PKA regulatory subunit 1α have been described in other cancers [57]. However, the differences between good and poor responders are notable and provide evidence for variation in molecular phenotype contributing to aggressiveness within the same histologic subtype of tumor.

The pathway analysis converged with the traditional fold change analysis at the hedgehog signaling pathway. The hedgehog signaling pathway appears to act upstream of Wnt/β-catenin signaling during bone development and aberrant hedgehog signaling has been associated with cancer development and progression [58]. As a result, we decided to examine other genes in the hedgehog pathway with qRT-PCR. Of the eight hedgehog-related genes examined, three were significantly down regulated in the poor responder cohort. These three genes, SMO, PTCH2 and DHH, where not identified on traditional fold change analysis and this is likely due to two factors. First, DHH is not annotated in the canine genome so primers were designed based on a region of the canine genome homologous to the gene in other species. Considering this, the Canine 2.0 microarray does not have a probe for this gene. Probes were present for the SMO and PTCH2 genes, however, in this study, raw array expression values for these genes were very low suggesting that the signal may be nearing the detection limit of the microarrays.

Hedgehog interacting protein (HHIP) was identified by fold change criteria in both RMA and PLIER preprocessed arrays. The up-regulation in the poor responder cohort was also observed as a trend in qRT-PCR but did not meet significance criteria. HHIP antagonizes all three of the hedgehog family of ligands (SHH, DHH and IHH) and has been shown to be down-regulated in numerous epithelial tumor types [59] with the notable exception of basal cell carcinoma where it is upregulated [60]. HHIP is also abundant in endothelial cells but is suppressed during angiogenesis through a VEGF mediated pathway [61]. The up-regulation of HHIP observed in our poor-responders likely has some causative relationship to the down-regulation of DHH and, through feedback loops, SMO and PTCH2 in the same cohort.

Three studies have examined gene expression in primary human OSA to identify chemotherapy-resistance signatures by comparing good and poor responders [62–64]. Among them, they identified over 200 differentially regulated genes but each gene set was unique to each study (i.e. there was no overlap in expression signatures). More recently, Walters and colleagues [65] assayed expression patterns in OSA cell lines with differing aggressiveness and identified 252 differentially regulated genes, four of which overlapped with the Mintz et al. study's gene signature [63]. This lack of similarity in expression patterns is observed in array analyses of various tumor types and is not at all surprising when one contemplates the differences in array preprocessing algorithms. Considering the disparity between the heat maps presented in Figure 1A and 1B, it is plausible that the exact same data processed in two different ways may yield two very different sets of candidate genes. Thus, in addition to traditional fold change analysis of microarray data for biomarker identification, a broader, unbiased systems-biology approach, such as we have done here, is likely to identify biological changes that can be reliably verified in multiple data sets. In fact, this approach was used to analyze multiple independent data sets to show that genes involved in the oxidative phosphorylation pathway were reduced in metastases compared to primary solid carcinomas [28]. Interestingly, in the current comparison of primary sarcomas, increased expression of genes in the oxidative phosphorylation pathway was associated with a poor outcome, suggesting that different metabolic factors may contribute to the initiation of metastasis from a primary tumor and the implantation and successful growth of metastasis at a distant site.

Given the small sample size of the study, we acknowledge that this data serves primarily as a road map for future studies. Our sample size was small primarily due to the stringent selection criteria set forth in Methods limiting our samples to dogs that had appendicular osteosarcoma, undergone amputation and received chemotherapy. Furthermore, we limited our samples to those from dogs with either very low or very high DFIs: the 100 and 300-day cutoffs were intended to straddle our facility's average DFI of 200 days. Selvarajah and colleagues recently studied gene expression in a group of dogs with OSA with good and poor outcome [66]. They utilized a larger sample size (n = 32) but included dogs with axial OSA as well as dogs that did not receive adjuvant chemotherapy. Although they did not observe differences in outcome due to these factors, previous studies have established an effect on DFI [2, 3, 5, 7]. They also based their good and poor responder groups on survival time instead of DFI: this can greatly affect outcome groups in a field of medicine where euthanasia is practiced. Beyond study-design differences, we applied a systems-based model for biomarker/pathway discovery by using the J5 metric to enrich for high to medium expressing genes that are most appropriate for selection as diagnostic/prognostic biomarkers, as opposed to fold-change based input. Hand annotation of many probesets based upon sequence homology allowed us to input a very large and complete data set into the MetaCore pathway analysis. Despite these differences in study design and analysis methods, we identified some pathways with similarity to those they identified in their PANTHER analysis, including hedgehog signaling, WNT signaling and chemokine signaling. Considering the differences in chemotherapy requirements between the two studies, these pathways may be most indicative of increased metastatic potential as opposed to chemotherapeutic resistance.

Work by Paoloni and colleagues provides strong evidence for the validity of dogs as a model for human OSA. They found that canine and human OSA are more similar to each other than to normal tissues from the same species [10]. This, in concert with our growing body of knowledge regarding gene and pathway derangements in canine OSA provides insights into the mechanisms of OSA progression and chemoresistance.

The present study has examined gene expression in primary canine OSA via both traditional fold change analysis and systems-based pathway analysis and found significant differences between dogs that responded poorly to chemotherapy following definitive treatment and dogs that responded well as evidenced by a long disease-free interval. This study has identified candidate biomarkers of aggressive tumors as well as pathways that are deranged in poor responders relative to good responders, opening the door for molecular prognostic screening in canine OSA and further molecular comparison between the human and canine disease. Although further studies, such as protein expression analysis will be necessary to solidify the role of these genes and pathways in OSA, targets identified here provide a strong foundation from which to identify druggable targets and markers of progression in OSA.