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. and ). 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 ), 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 . 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
] 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
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
]. 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 and , 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
]. 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.