Abstract
Osteosarcoma is a debilitating bone cancer that affects humans, especially children and adolescents. A homologous form of osteosarcoma spontaneously occurs in dogs, and its differential incidence observed across breeds allows for the investigation of tumor mutations in the context of multiple genetic backgrounds. Using whole-exome sequencing and dogs from three susceptible breeds (22 golden retrievers, 21 Rottweilers, and 23 greyhounds), we found that osteosarcoma tumors show a high frequency of somatic copy-number alterations (SCNA), affecting key oncogenes and tumor-suppressor genes. The across-breed results are similar to what has been observed for human osteosarcoma, but the disease frequency and somatic mutation counts vary in the three breeds. For all breeds, three mutational signatures (one of which has not been previously reported) and 11 significantly mutated genes were identified. TP53 was the most frequently altered gene (83% of dogs have either mutations or SCNA in TP53), recapitulating observations in human osteosarcoma. The second most frequently mutated gene, histone methyltransferase SETD2, has known roles in multiple cancers, but has not previously been strongly implicated in osteosarcoma. This study points to the likely importance of histone modifications in osteosarcoma and highlights the strong genetic similarities between human and dog osteosarcoma, suggesting that canine osteosarcoma may serve as an excellent model for developing treatment strategies in both species.
Significance: Canine osteosarcoma genomics identify SETD2 as a possible oncogenic driver of osteosarcoma, and findings establish the canine model as a useful comparative model for the corresponding human disease. Cancer Res; 78(13); 3421–31. ©2018 AACR.
Introduction
Osteosarcoma is an often fatal form of bone cancer occurring in humans, predominantly in the adolescent age group. It accounts for 2.4% of pediatric cancers and has a greater incidence in males than females (1). In addition, late onset of osteosarcoma is seen in adults over 65 years of age, with a similar male to female ratio, and frequently occurs as a second malignancy comorbid with the osteoclastic Paget disease (2). Osteosarcoma can occur in any bone, though most often it is found in the extremities of long bones near the metaphyseal growth plates (1). The tumor is thought to arise from primitive mesenchymal bone-forming cells that produce malignant osteoid (3). The global 5-year survival rate for osteosarcoma is 80%, but metastatic osteosarcoma is hard to control, with a 5-year survival of only approximately 46% (4).
Early molecular analysis of human osteosarcoma tumors using array-based methodologies highlighted somatic copy-number alterations (SCNA) in cancer genes as one of the significant features of the disease (5). In addition, genome-wide association studies, epigenetic studies, and mRNA and miRNA expression profiles have identified genetic changes affecting the disease (6–8). Recent efforts have used both whole-exome and/or whole-genome sequencing of osteosarcoma tumor/normal pairs to study the pathogenesis of osteosarcoma as well as its somatic progression (9–11). TP53 tumor mutations were observed in 75% of the 59 human patients studied (11), whereas recurrent somatic alterations were found in RB1, ATRX, and DLG2 in addition to TP53 (9). Both these studies demonstrated a predominance of structural variation in the form of SCNAs in osteosarcoma, whereas relatively few somatic point mutations (SPM) and indels were found. The authors also reported complex chains of rearrangements and localized hypermutation across the genome in most samples. Thus, although multiple layers of data exist, the mechanism of progression from germline predisposition to somatic changes in osteosarcoma tumorigenesis is not yet understood (7).
Every year, almost 10 times as many dogs are diagnosed with osteosarcoma as humans (12). In contrast to human osteosarcoma, canine osteosarcoma typically affects middle-aged dogs from large or giant breeds such as golden retriever (Golden), Rottweiler (Rott), and greyhound (GreyH). Sex is not a factor contributing to disease susceptibility among the various breeds, but the size of the dog does matter, with large and giant dogs (>25 kg) accounting for 90% of all cases of osteosarcoma (13). A study of bone tumors among 400,000 dogs in Sweden (14) showed, after correcting for average life span (dog-years at risk, DYAR), that dog breeds have different rates of disease incidence, which is indicative of a genetic predisposition for the disease. For example, Goldens are seen to have a lower incidence when compared with GreyHs and Rotts (incidence rates are 6, 30, and 36 cases per 10,000 DYAR, respectively; ref. 14). Although osteosarcoma has a later onset in dogs than in humans, both have comparable clinical representation, and canine osteosarcoma in predisposed breeds is known to be most similar to human pediatric osteosarcoma (12, 15). In both species, the primary bone tumor consists of malignant stroma (connective tissue cells) and has as a hallmark osteoid formation by the malignant osteoblasts (16). In keeping with common tumor sites between humans and dogs (13), parallel treatment options are offered, including limb-sparing procedures and adjuvant chemotherapy (15).
Earlier canine osteosarcoma studies have provided key insights into disease biology and genetic predisposition. In our previous genome-wide association study of >300 cases and controls from three breeds (GreyH, Irish wolfhound, and Rott), we identified 33 significant loci that explained 50% to 80% of phenotypic variance in these breeds. These loci contain genes such as OTX2 and BMPER that are known to be involved in bone morphogenesis and to affect bone health (17). In an oligonucleotide array comparative genome hybridization (oaCGH) study, it was shown that copy-number alterations in dogs are relatively orthologous to what is found in humans (18). Analyzing 123 dog osteosarcoma tumors (18), the authors found that the genomic imbalance in the two species is similar, and the expression levels of some of the shared genes (e.g., RUNX2, TUSC3, and PTEN) are correlated with the ploidy changes. Nevertheless, the lack of a deep understanding of the genetic causes of tumor variation in canine osteosarcoma limits the use of this naturally occurring model for understanding the human disease. Here, we performed whole-exome sequencing (WES) of matched tumor–normal (T/N) pairs in three breeds (Golden, GreyH, and Rott), all predisposed to osteosarcoma but with different lifetime risks. The goals here were to discover putative candidate driver genes, to evaluate if dogs share mutated genes with humans, and, lastly, to detect potential differences among breeds. The results in all three breeds confirm observations in human osteosarcoma, with a majority of individuals showing high levels of copy-number alterations and TP53 mutations. In addition, we also find multiple mutations in the SETD2 gene, which only once has been described in human osteosarcoma (19).
Materials and Methods
Matched tumor normal sample collection
Tumor and normal tissue samples from a total of 70 matched pairs (23 sample pairs for golden retriever and Rottweiler and 24 greyhound) were collected by coauthors or provided by the Pfizer Canine Comparative Oncology and Genomics Consortium (CCOGC) Biospecimen Repository. Study inclusion was voluntary. All samples (blood samples and leftover tissue from removal of spontaneously acquired tumors) were taken by the dogs' veterinarian. All work was performed in accordance with ethical guidelines and is covered by the ethical approval IACUC 0064-08-15 (Lindblad-Toh, Broad Institute).
Library construction, exome capture, and sequencing
Exome libraries of genomic DNA isolated from osteosarcoma tumors and their matched normal tissue were generated via a Roche NimbleGen custom design canine exome array (120705_CF3_Uppsala_Broad_EZ_HX1) and submitted for 76 base paired-end Illumina sequencing on HiSeq 2000 platform, with a minimum target coverage of 30× for the normal and 60× for the tumor samples.
Exome alignment and somatic variation calling
SPMs and somatic indel mutations (SIM) were identified using a standard workflow methodology (Supplementary Fig. S1). Seventy matched tumor and normal sequence reads from the three breeds were aligned against the dog genome assembly of CanFam3.1 using BWA 0.7.10 with default options. The BWA alignment BAMs, prior to somatic variant calling with the algorithm MuTect2, were refined in accordance with Genome Analysis Toolkit's recommended best practices (Supplementary Fig. S1). Tools and databases used for all of the methods are listed in Supplementary Table S1. SPMs and SIMs were subsequently identified in 66 samples that passed quality threshold criteria with MuTect2. The unfiltered calls were then passed through a set of filters to reduce germline artifacts. The first filter was with Panel of Normal (PON) calls, which was built on a set of normal samples from three dog cancers using SAMTools. Variants present in >2 samples in the PON were filtered out. A second filtering was based on the removal of population variation present in the dog dbSNP database (Build ID: 146), and germline mutations called on 160 dogs across 8 breeds (unpublished data). Further false-positive calls that may arise due to oxidative DNA damage were filtered out in accordance with methods by Costello and colleagues (20). The last filter applied was to remove calls that were outliers based on a boxplot of the read depth values corresponding to the variant calls.
Annotation of the somatic variants using SnpEff
Functional consequences of variants were predicted using the tool SnpEff 3.3 (21). SnpEff annotations of the SPMs and SIMs were categorized as nonsynonymous or missense mutations, synonymous or silent mutations, stop-gain or nonsense mutations, frameshift, and splice site mutations.
Identification of significantly mutated genes with MuSiC
To find putative driver genes, the nonsilent mutations from the SnpEff annotations were run through the program Genome MuSiC 0.4. Default options were used, except for two parameters: the false discovery rate threshold was fixed at 0.1, and the flag to merge concurrent mutations in a given gene was turned on.
Mutational signature discovery with Bayesian NMF
The mutational signature discovery where SPM or SIM counts in each sample, stratified by mutation context, was performed using The Broad Institute's Bayesian nonnegative matrix factorization (NMF) algorithm. The comparison of the signatures discovered in the dog dataset and the standardized COSMIC signatures was performed using standard hierarchical clustering.
Pathway analysis for the significantly mutated genes using ENRICHR/Kyoto Encyclopedia of Genes and Genomes 2016
Functional assessment of the significantly mutated genes (SMG) predicted by MuSiC was performed using ENRICHR with default options. The gene sets were subject to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway (2016) enrichment with a threshold set at a cutoff of P < 0.05 to select significant pathways.
Somatic copy-number detection with VarScan2
SCNAs in the tumor sample that are not present in the normal sample were identified using VarScan2 (http://varscan.sourceforge.net/copy-number-calling.html). The raw calls were subjected to circular binary segmentation to delineate segments by copy number and to categorize major regions with alterations. Finally, a custom Perl script was used to merge similar adjacent events.
Identification of recurrent SCNAs with GISTIC 2.0
To identify recurrent SCNAs, the tool GISTIC 2.0 was used. Default options were used, and 64 of 66 samples passed through the cutoff threshold of 5,000 for the flag max seg, a measure of the copy-number segments and a mandatory parameter for GISTIC 2.0 analysis. Two samples, Golden.04 with 5,436 and GreyH.05 6,352 copy-number segments each, were omitted from the final analysis as these values were > 5,000.
Detection of germline mutations in key osteosarcoma genes
A total of 8 genes previously found associated with human osteosarcoma by genome-wide association study (GWAS; gene-desert on chromosome 2p25.2 not analyzed; refs. 6, 22), the canine top associated locus (17), as well as those implicated in human syndromes featuring a heightened frequency of osteosarcoma (23) were examined for mutations in our germline exome data.
Data access
All data used in this article are available from NCBI BioProject (https://www.ncbi.nlm.nih.gov/bioproject) under the entry PRJNA391455. Raw sequencing reads are available from the NCBI Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra). Identified mutations are displayed in the CamFam3.1 Track Hub Broad Improved Canine Annotation v1 (https://genome.ucsc.edu/cgi-bin/hgHubConnect) under the name osteosarcoma somatic variants.
Results
WES of tumor/normal pairs from three breeds
To characterize somatic mutations in canine osteosarcoma tumors, 70 matched tumor/normal pairs were collected and used to generate whole-exome sequence data (Supplementary Table S2). The samples included 46 male and 24 female dogs, similarly distributed across the three breeds. The age at diagnosis varied from 3 to 12 years, with a median of 8.5 years in GreyHs and 9.0 years in both Rotts and Goldens. The location of tumors included humerus, femur, radius, tibia, and others. The canine exome array contains approximately 57 Mb of sequence representing more than 85% of the coding sequence. Sixty-six pairs of samples (22 Golden, 23 GreyH, and 21 Rott) passed the quality thresholds, with a median sequence depth of 136x (range, 90–176) for the tumor and 95x (range, 31–166) for the normal samples.
SCNAs occur extensively in all three breeds
VarScan2 followed by circular binary segmentation across all three breeds identified a large number of SCNAs. The Goldens had on average more of the genome covered by amplifications (mean ∼230 Mb) compared with Rott (∼200 Mb) and GreyH (∼170 Mb). A similar trend was also observed for deletions (Golden ∼620 Mb, Rott ∼530 Mb, GreyH ∼440 Mb each). Of the 64 samples passing the threshold criteria, 19 broad or large-scale (corresponding to more than half the length of the chromosome in which they occur) and 67 focal (1–11 Mb) amplification and deletion events were found per sample using GISTIC 2.0 (Fig. 1A; Supplementary Table S3). For 27 of 38 autosomes, more than 50% of samples had a focal lesion. For seven of these chromosomes (i.e., 13, 16, 18, 20, 21, 25, and 38), the number of samples with focal lesions surpassed 80% (Supplementary Fig. S2).
Regions of somatic copy-number aberrations in canine osteosarcoma. A, GISTIC identifies regions of recurrent copy-number amplifications (red) as well as deletions (blue) across the autosomes. Cytobands of recurrently altered regions are denoted on the y-axes. The green line indicates significant enrichment of copy-number changes. B, Top eight tumor-suppressor genes and top nine oncogenes with recurrent SCNAs across the breeds. Amp, amplifications; Del, deletions.
Regions of somatic copy-number aberrations in canine osteosarcoma. A, GISTIC identifies regions of recurrent copy-number amplifications (red) as well as deletions (blue) across the autosomes. Cytobands of recurrently altered regions are denoted on the y-axes. The green line indicates significant enrichment of copy-number changes. B, Top eight tumor-suppressor genes and top nine oncogenes with recurrent SCNAs across the breeds. Amp, amplifications; Del, deletions.
We found that the chromosomal regions corresponding to recurrent focal amplifications and deletions contained numerous cancer-related genes (Fig. 1B; Supplementary Tables S4–S6). Downstream analyses of genes affected by SCNAs are limited to focal lesions only, because estimating large-scale events with WES data is at best suggestive. The genes in focal SCNAs that were altered in more than 50% of all samples across all breeds included TP53, CDKN2A/B, HRAS, and PTEN. These are among the top 20 genes previously implicated in human and/or canine osteosarcoma (24, 25). In addition, we identified more cancer genes with no previously known roles in osteosarcoma, including HIPK2 and EDNRB, all of which have focal lesions in more than 80% of all samples across the three breeds. These findings from WES analysis showed extensive correlation with oaCGH profiles of 22 cases from the same cohort, both at the level of large chromosomal segments of DNA copy-number imbalance and for focal lesions involving key cancer genes (Supplementary Fig. S3; Supplementary Material; Thomas and colleagues manuscript in preparation).
Pathway enrichment analysis for the 248 genes classified as recurrently altered by GISTIC was performed using ENRICHR (Supplementary Table S7; ref. 26). The most significant pathway was the cytokine–cytokine receptor interaction pathway. The genes in this pathway include BMPR1B, a receptor involved in bone formation and a molecular target in breast cancer metastases, six interleukin family members that affect inflammation and regulate immune response, and four genes from the TNF receptor superfamily. The second significant pathway was the PI3K–Akt signaling pathway, which is thought to play an essential role in the initiation and development of human osteosarcoma and is considered to be a potential therapeutic target (27).
Variable rates of somatic mutations across breeds
Analysis of the exome sequence data using Mutect 2 revealed a total of 7,900 SPMs and 1,197 SIMs across all individuals. The average numbers of SPMs and SIMs were different among the breeds (one-way ANOVA, P value 0.03; Fig. 2A and B; Supplementary Table S2), and a post hoc analysis showed that the Goldens had significantly higher averages than GreyHs (Tukey test, P value 0.03; Fig. 2A). The Golden group also had a wider range of mutations, showing a higher median of 133 SPMs (129 when excluding two outliers), compared with 113 SPMs in Rotts and 97 SPMs in GreyHs (Table 1). There was little overlap in the SPM and SIM calls among the three breed groups (Supplementary Fig. S4A and S4B). SnpEff, used to annotate the SPMs and SIMs, detected protein-modifying changes in 2,041 genes. Of these, 290 genes were altered in two or more samples irrespective of breed (n = 66), whereas 71 genes had mutations present in three or more samples. Again, the Goldens harbored more nonsilent mutations (mean of 54 per sample) compared with the Rotts (48 per sample) and GreyHs (40 per sample). No correlation was observed across breeds between the somatic mutation burden and sex or tumor location; however, the total number of somatic mutations in all breeds combined was positively correlated with increasing age at diagnosis (Spearman correlation ρ = 0.3, P value 0.03; Supplementary Fig. S5).
Summary of SPMs and SIMs across breeds. A, Distribution of SPMs and SIMs among the three breeds showing median values of 130 in Golden, 79 in GreyH, and 95 in Rotts, with a significant difference in median frequency between Golden and GreyH (*, P < 0.03). B, Mean mutation frequencies normalized per megabase pair (Mb). Color scales: 0 (green) to 4 (pink) for SPMs/Mb, and 0 (gray) to 0.4 (black) for SIMs/Mb.
Summary of SPMs and SIMs across breeds. A, Distribution of SPMs and SIMs among the three breeds showing median values of 130 in Golden, 79 in GreyH, and 95 in Rotts, with a significant difference in median frequency between Golden and GreyH (*, P < 0.03). B, Mean mutation frequencies normalized per megabase pair (Mb). Color scales: 0 (green) to 4 (pink) for SPMs/Mb, and 0 (gray) to 0.4 (black) for SIMs/Mb.
Summary of somatic point and indel mutations across the three breeds
Breed . | Median SPM . | Median SIM . | Mean SPM . | Mean SIM . | Mean SPM/Mbp . | Mean SIM/Mbp . | Mean nonsilent mutations . | Total SMGs . |
---|---|---|---|---|---|---|---|---|
Golden (n = 22) | 133 | 16 | 148 | 21 | 1.4 | 0.19 | 54 | 10 |
GreyH (n = 23) | 97 | 12 | 101 | 14 | 0.9 | 0.13 | 40 | 6 |
Rott (n = 21) | 113 | 22 | 110 | 20 | 1.0 | 0.19 | 48 | 7 |
Breed . | Median SPM . | Median SIM . | Mean SPM . | Mean SIM . | Mean SPM/Mbp . | Mean SIM/Mbp . | Mean nonsilent mutations . | Total SMGs . |
---|---|---|---|---|---|---|---|---|
Golden (n = 22) | 133 | 16 | 148 | 21 | 1.4 | 0.19 | 54 | 10 |
GreyH (n = 23) | 97 | 12 | 101 | 14 | 0.9 | 0.13 | 40 | 6 |
Rott (n = 21) | 113 | 22 | 110 | 20 | 1.0 | 0.19 | 48 | 7 |
SPMs across the three breeds show three distinct mutational signatures including a novel one
To delineate the mutational processes occurring in canine osteosarcoma, we performed a de novo mutational signature analysis across all tumor samples using a Bayesian NMF methodology. Our study identified three distinct mutational processes (Fig. 3A). The first profile, which is the most common, has the hotspot motif C>T in the CpG nucleotide context, similar to what is seen in the COSMIC 1 aging–associated signature (http://cancer.sanger.ac.uk/signatures/Signature-1.png). In this study, the COSMIC 1 signature had a distinct representation of the C>A mutations. The C>A and C>T signals were not separable by NMF and mostly co-occurred across samples. The second signature, COSMIC 17 (http://cancer.sanger.ac.uk/signatures/Signature-17.png), is a motif for which the etiology has not been ascertained in humans. The third signature, characterized by T>G transversions, is novel, having not previously been described in human cancers.
Three mutational signatures differ between breeds. A, Bayesian NMF identifies three distinct mutational signatures including a novel signature: Cosmic 1 signature (brown), Cosmic 17 signature (white), and a novel signature (green). B, The relative proportion of each of the three mutational signatures per sample shows a different distribution across the three breeds.
Three mutational signatures differ between breeds. A, Bayesian NMF identifies three distinct mutational signatures including a novel signature: Cosmic 1 signature (brown), Cosmic 17 signature (white), and a novel signature (green). B, The relative proportion of each of the three mutational signatures per sample shows a different distribution across the three breeds.
Although all three mutational signatures were found in all breed groups, the contribution of each signature per sample indicates that there is a difference in the representation of signatures among the three breeds (Fig. 3B). Tests for enrichment based on the mean contribution of each signature per breed showed that the GreyHs and Rotts have a higher rate of the COSMIC 1–like signature (∼85% each, P value 0.01), whereas the novel signature is more prevalent in the Goldens (24% signature activity in Golden, P value 0.01; Supplementary Fig. S6).
Eleven genes are significantly mutated
Significantly mutated genes (SMG, genes that have a statistically higher rate of mutation than the estimated background rate) that are enriched for nonsilent mutations are likely to be drivers in the development of cancers. Using the MuSiC algorithm, the analysis of all breeds together identified a total of 11 SMGs (default criteria, 10% false discovery rate; Supplementary Table S8; Fig. 4). The top three mutated genes across all samples were, in order, TP53, SETD2, and TANGO. TP53, tumor protein 53, had mutations in roughly 60% of all samples (Goldens 55%, Rotts 67%, and GreyHs 61%). SETD2 was mutated in 21% of samples (Goldens 32%, Rotts 19%, and GreyHs 13%) and showed a variety of mutation types (frameshift, nonsense, splice, and missense mutations). Analyses of nonsilent mutations with MutSigCV (software.broadinstitute.org/cancer/software/genepattern/modules/docs/MutSigCV) and tests that prioritize genes with high-impact mutations, as discerned by Ensembl Variant Effect Predictor (VEP), reiterate TP53 and SETD2 as top putative drivers (Supplementary Tables S9 and S10).
Distribution of mutation types in SMGs across all samples. Note that TP53 and SETD2 have the highest number of mutations in all three breeds. See the inset for the mutation-type color code.
Distribution of mutation types in SMGs across all samples. Note that TP53 and SETD2 have the highest number of mutations in all three breeds. See the inset for the mutation-type color code.
The 11 SMGs were further examined for functional roles in known KEGG pathways using ENRICHR (26). This analysis showed that many cancer-related pathways that include TP53, and one pathway that includes SETD2, were found to be significantly overrepresented (Supplementary Table S7).
Similar to human cancers, many mutations in TP53 are found within protein domains
TP53 was the most frequently altered gene among the three breed groups. Note that 83% of all samples showed TP53 as being affected either by a point mutation and/or by an SCNA (Fig. 5A). TP53 is also the only gene common to both the SMGs and the recurrent SCNA gene lists. A total of 41 nonsilent mutations were distributed along the gene (Fig. 5B). More than 85% of these mutations occurred in protein domains and show a similar pattern with the mutations found in TP53 in human osteosarcoma. VEP predicted 15 of the 41 sites to be “high”-impact mutations, i.e., mutations that are likely to have a disruptive effect on the gene, leading to its loss of function (Supplementary Table S11). Furthermore, there was a significant increase in the amount of SCNAs present in dogs with TP53 loss-of-function mutations (one tailed t test, P value 0.01).
Somatic mutations in TP53 and SETD2, and germline mutations in CDKN2A/B. A, TP53 somatic mutations and SCNAs per individual across breeds. Colored boxes correspond to different mutation types, and the background color represents amplifications (red) and deletions (blue). B, Distribution of canine mutations lifted over to the human TP53 protein shows that the majority of mutations are in the DNA-binding domain. C, Distribution of mutations within the SETD2 protein is similar to that found in other human cancers. D, CDKN2A/B shows recurrent SCNAs and germline mutations in 17 of the samples (missense in 16 samples and frameshift in one sample). Colored boxes correspond to different mutation types, and the background color represents amplifications (red) and deletions (blue).
Somatic mutations in TP53 and SETD2, and germline mutations in CDKN2A/B. A, TP53 somatic mutations and SCNAs per individual across breeds. Colored boxes correspond to different mutation types, and the background color represents amplifications (red) and deletions (blue). B, Distribution of canine mutations lifted over to the human TP53 protein shows that the majority of mutations are in the DNA-binding domain. C, Distribution of mutations within the SETD2 protein is similar to that found in other human cancers. D, CDKN2A/B shows recurrent SCNAs and germline mutations in 17 of the samples (missense in 16 samples and frameshift in one sample). Colored boxes correspond to different mutation types, and the background color represents amplifications (red) and deletions (blue).
SETD2 is mutated in 21% of the samples
The second most frequently mutated gene, SETD2, encodes a histone methyltransferase (HMT) that trimethylates histone H3 lysine 36 (H3K36me3). The gene showed a total of 17 nonsilent mutations (9 frameshift, 2 splice site, 5 nonsense, and 1 missense mutations) across 21% of the tumor samples (Fig. 5C). Of these mutations, 16 are predicted by VEP to be high-impact mutations (Supplementary Table S12). Projecting the mutational changes from the canine SETD2 gene to the orthologous human protein showed that the mutations occurred across the gene, similar to what is observed in other human cancers (COSMIC and cBioPortal databases). As SETD2 is an HMT, we decided to investigate if our nonsilent mutation dataset included variants that could lead to protein changes in other histone-modifying genes. Examining 68 additional genes (excluding SETD2; Supplementary Fig. S7), we found mutations in 14 samples (3 Goldens, 3 GreyHs, and 8 Rotts). Nine samples had mutations in HMTs, four in histone acetyltransferases (HAT), and one in histone deacetyltransferases (HDAC; Supplementary Fig. S7). In total, 36% of samples had a mutation in at least 1 histone-modifying gene. Because these genes were found using a targeted approach, we could not use the resultant gene information for enrichment analysis.
Detection of germline variants in genes previously associated with genetic predisposition to osteosarcoma in humans and dogs
The different lifetime risks for developing osteosarcoma among the breeds in our study suggest the possible presence of inherited risk factors that could confer susceptibility to the disease. We therefore inspected our data for the occurrence of nonsilent coding germline mutations in key genes previously implicated in GWAS for canine (17) and human osteosarcoma (6, 22), as well as genes responsible for familial syndromes with an increased risk for osteosarcoma (Supplementary Fig. S8; refs. 28–31). Overall, 59 of the 66 samples had at least one germline candidate mutation with moderate- or high-impact VEP effects in genes previously reported to increase the risk of osteosarcoma. These included CDKN2A/B (n = 20), GRM4 (n = 12), NFIB (n = 1), TP53 (n = 1), RB1 (n = 1), BLM (n = 23), and WRN (n = 28; Supplementary Fig. S8; Supplementary Table S13). Most prominently, CDKN2A/B previously implicated in canine osteosarcoma predisposition (17) shows germline protein-modifying changes across several samples (Fig. 5D). In the Rotts, 67% (14/21) of the cases harbor a moderate N18D missense mutation (BICF2P440554, chr11:41,226,002). In contrast, only 9% of Goldens (2/22; one frameshift and one missense) and 4% of GreyHs (1/23; one missense) had protein-modifying mutations in CDNK2A. Five missense mutations were observed in the CDKN2B locus (Goldens 3/22, Rotts 2/21). Similarly, germline mutations in BLM are predominantly seen in GreyHs (52%) and Rotts (52%), and mutations in WRN are seen in Goldens (55%) and GreyHs (65%). Thus, the presence of candidate mutations in different genes in different breeds suggests that genetic background shows variability across the breeds and these potentially may also affect the somatic mutational spectrum.
Discussion
For the first time, we perform WES in three canine breeds predisposed to osteosarcoma, and evaluating somatic changes, we recapitulate and extend what has been observed in human osteosarcoma. We find the mutational landscape in dogs reflect what is observed in pediatric osteosarcoma—there are more SCNAs than point mutations, and the average somatic mutation burden across the species is of the same magnitude. At the point mutation level, TP53 is the most frequently altered gene, with most mutations in the DNA-binding domain. In human cancers, inactivating mutations in the TP53 gene may disrupt checkpoint responses to DNA damage, often leading to genomic or chromosomal instability and aberrations and an increased risk for SCNAs (32). In our dog dataset, we also see a correlation for more SCNAs in tumors with TP53 loss of function.
The canine set of 248 recurrently altered genes (Supplementary Table S14) included known tumor suppressors such as TP53, CDKN2A/B, PTEN, and DCL1, as well as genes not previously known to have frequent SCNAs in human osteosarcoma. For instance, the tumor-suppressor gene EDNRB is deleted in more than 80% of the canine samples and is often downregulated in other human cancers (33). Interestingly, EDNRB is a nonspecific receptor for EDN1 that has been implicated in promoting tumorigenesis in canine osteosarcoma (34). Constituents of the endothelin axis (EDN1 and its associated receptors) are known to be impaired/dysregulated in human cancers and have been explored as novel anticancer therapeutic targets (35, 36). HIPK2, homeodomain-interacting protein kinase 2, is also altered in more than 80% of the samples and was shown to increase tumor production upon inactivation (37). HIPK2 also has a role in TP53 regulation and has been proposed as a candidate therapeutic target (38).
The RB1 gene, one of the constituents of the cell-cycle control pathway that includes the CDKN2A/B locus, is often mutated in human osteosarcoma (39). However in the current canine analysis, SCNAs encompassing the RB1 locus were only observed in 6% of tumors. In the absence of RB1 gene mutations, CDKN2A/B is often deleted or inactivated, likely leading to deregulation of the cell-cycle control pathway in high-grade human osteosarcomas (40). We note that CDKN2A/B shows copy-number alterations in 78% of our canine samples and has germline candidate protein-changing mutations in 81% of Rott samples. In addition, nonsilent germline variants in the RB1 gene, which is known to be altered in 30% to 75% of human osteosarcoma (1), only affected 1 of the 66 samples in our study. It therefore can be hypothesized that canine CDKN2A/B mutations could replace the RB1 mutations often seen in human osteosarcoma.
The p14ARF (known as p19ARF in the mouse, encoded by isoform 2 of CDKN2A) protein binds MDM2, thereby preventing p53 degradation (41). Somatic mutations of CDKN2A are very common in human cancers (42), and germline mutations are associated with familial melanoma, glioblastoma, and pancreatic cancer (43). Analysis of the germline risk data and somatic changes showed that 83% of individuals across the three breeds had either a germline risk variant and/or a somatic copy-number aberration affecting CDKN2A, thus strengthening the role of this gene in canine osteosarcoma initiation and progression. Importantly, the germline risk variant common in Rotts only affects CDKN2A isoform 1, thereby only altering p16INK4a. We thus conclude that similarly to the human disease (7), both the p53 and the RB1 pathways are likely key players in canine osteosarcoma.
SETD2, a known tumor-suppressor gene in numerous human cancers (44–46), and only recently found to be mutated in less than 2% of human osteosarcoma samples (19), had previously no recognized roles in canine osteosarcoma. Here, we found SETD2 mutations in 21% of tumors including dogs from all three breeds. Changes in SETD2 in different human cancers have been associated with a variety of clinical consequences. For example, in clear cell renal cell carcinoma, loss of SETD2 through lowered mRNA expression is known to be associated with the aggressive form of the disease, likely leading to poorer survival outcome (47). Similarly, in a breast cancer study, (48), it was shown that the mRNA expression level of SETD2 correlated with the severity of clinical parameters, in that patients with metastasis or recurrence had lower expression levels of the gene compared with patients who were disease free. In a study involving bone and soft-tissue sarcomas in humans, SETD2 was found to be a putative driver in synovial sarcomas (49). The authors also included osteosarcoma as part of their investigation, but reported no mutations in SETD2 for it. Curiously, in the analysis of 24 individuals afflicted by chordoma, a rare bone sarcoma found in the skull, SETD2 was the most frequently altered gene (14/24) either through point mutations or through copy-number losses (50). In our present study, we did not observe any SCNA changes in SETD2, only point or frameshift insertion/deletion events that are likely damaging. The variant allele fraction (VAF; fraction of alleles that are mutated for a given sample) distribution for SETD2 was significantly enriched for high VAF values (Wilcoxon rank test, P value < 0.01) and consequently has an elevated mean VAF of 0.44, compared with 0.20 for the set of all detected mutations. This pattern is similar to that observed for TP53 (Wilcoxon rank test, P value < 0.01, mean 0.54), leading to the hypothesis that the somatic changes in both genes might have been early clonal events that occurred before the onset of carcinogenesis. Thus, it is possible that SETD2 might also be a putative driver for canine osteosarcoma akin to TP53.
Furthermore, when we examined other histone modifiers/methyltransferases that included HATs, HDACs, and HMTs, and taking into account the SETD2 mutations, we found that a total of 36% of samples have mutations in these genes. Histone modifiers are known to be frequently disrupted in cancer causing tumorigenic alterations and tumor proliferation (51). Therefore, they are regarded as potential drug targets (52), and inhibitors of HDACs and HMTs are in clinical trials for a number of human cancers (53). This may suggest that histone modifier inhibitors could be tried in patients with osteosarcoma as well (54).
Patterns of somatic changes in human tumor cells, referred to as “mutational signatures,” have been increasingly characterized in recent years (55). We identified three distinct mutational signatures active in canine osteosarcoma. An aging-associated signature, COSMIC 1–like, has been previously reported in human osteosarcoma (9) and is the most common signature across the three breeds. The second signature, COSMIC 17, though not observed in human osteosarcoma, is a common signature in several other human cancers such as liver, lymphoma, and stomach cancers (55). The third signature is likely a novel signature with frequent T>G transversions not previously reported. Interestingly, in the study of human osteosarcoma in which 31 samples (mostly pediatric cases) were analyzed using WES (10), the dominant mutational signature found was COSMIC 3, a signature that is most often associated with BRCA1/2 mutations in breast cancers. In our dog dataset, we find no evidence for this signature. Although all three signatures are present in all breeds, there is considerable difference in their distribution. The Rotts and GreyHs show a dominance of the aging signature (accounting for ∼85% of mutational profiles for both breeds), whereas COSMIC 17 and the novel signature show enrichments in the Goldens. We hypothesize that the differences might be the result of the genetic background present in the different dog breeds. This along with the mechanisms underlying this signature needs to be further examined.
In addition to finding breed differences in the mutational signatures, we also noted breed differences in the candidate protein-altering germline mutations found in different breeds: the CDKN2A and BLM genes were commonly mutated in Rotts, the WRN and BLM genes had mutations in GreyHs, and the GRM4 and WRN genes were mutated in Goldens. In addition, Goldens, which have a lower disease frequency than the other breeds, appeared to have more SCNAs as well as SPMs and SIMs. (It is likely that GreyHs and Rotts with relatively higher susceptibility to osteosarcoma than Goldens may need fewer mutations for triggering carcinogenesis.) Although it is too early to connect these findings, it will be interesting to try to unravel the effects of genetic background on tumor initiation and progression and study the candidate genes further to elucidate their roles in osteosarcoma.
In summary, by analyzing WES data for canine osteosarcoma tumors across three breeds, we identified both SMGs and SCNAs that are likely to include driver genes and noted multiple similarities between canine and human osteosarcoma, while also identifying several breed differences. Similarly to what is reported for human osteosarcoma, we found a large number of SCNAs and frequent mutations in TP53. In addition, we found that SETD2 is frequently mutated in canine osteosarcoma tumors. Both these genes are known cancer genes, but SETD2 so far has not been explicitly linked with osteosarcoma and along with other histone modifiers that are being altered may offer novel opportunities for targeted therapy. Intriguingly, SETD2 is known to regulate TP53, whereas CDKN2A/B, which shows both germline and somatic mutations, affects both TP53 and RB1. Although we expected multiple differences between the breeds, the recurrent SCNAs and tumor mutations affecting key genes were surprisingly similar; however, germline coding candidate mutations in CDKN2A, GRM4, BLM, and WRN as well as one mutational signature was found more frequently in some breeds than the others. Finally, we believe that the insights learned from this study not only serve to inform dog osteosarcoma biology, but also highlight the usefulness of the canine disease models for enhancing our current understanding of the corresponding human disease.
Disclosure of Potential Conflicts of Interest
C.G. Couto has honoraria from the speakers' bureau of and is a consultant/advisory board member for Idexx Laboratories. No potential conflicts of interest were disclosed by the other authors.
Authors' Contributions
Conception and design: J.F. Modiano, K. Lindblad-Toh
Development of methodology: S. Sakthikumar, J. Turner-Maier, J. Alföldi, M.E. Pettersson, V.D. Marinescu, K. Lindblad-Toh
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): R. Thomas, R. Swofford, J. Johnson, C.G. Couto, W.C. Kisseberth, J.F. Modiano, K. Lindblad-Toh
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S. Sakthikumar, I. Elvers, J. Kim, M.L. Arendt, R. Thomas, S.E. Schumacher, E. Axelsson, M.E. Pettersson, G. Getz, J.R.S. Meadows, M. Kierczak, K. Forsberg-Nilsson, V.D. Marinescu, K. Lindblad-Toh
Writing, review, and/or revision of the manuscript: S. Sakthikumar, I. Elvers, M.L. Arendt, R. Thomas, E. Axelsson, C.G. Couto, W.C. Kisseberth, M.E. Pettersson, J.R.S. Meadows, J.F. Modiano, M. Breen, M. Kierczak, K. Forsberg-Nilsson, V.D. Marinescu, K. Lindblad-Toh
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): S. Sakthikumar, J. Turner-Maier, R. Swofford, J. Johnson, J. Alföldi
Study supervision: I. Elvers, J.R.S. Meadows, M. Kierczak, K. Forsberg-Nilsson, V.D. Marinescu, K. Lindblad-Toh
Acknowledgments
We thank all dogs and their owners for their help and support, the Pfizer Canine Comparative Oncology and Genomics Consortium (CCOGC) and The Ohio State University College of Veterinary Medicine Biospecimen repositories for providing samples, the Broad Institute Genomics platform for sequencing, Uppmax for supplying computational resources, and Dr. Lina Hultin-Rosenberg from Uppsala University for her assistance with statistical analyses. We would also like to thank Leslie Gaffney of the Broad Institute for help with graphics.
Funding is gratefully acknowledged from the NIH (U54 HG003067-08, awarded to the Broad Institute), the European Research Council (ERC), the Göran Gustafsson Foundation, Swedish Cancer Foundation, the American Kennel Club (AKC) Canine Health Foundation, the Golden Retriever Foundation, and the National Cancer Institute (NCI) for the OSU tissue bank (P30 CA016058). I. Elvers is supported by a postdoctoral fellowship from the Swedish Medical Research Council, SSMF. J.R.S. Meadows is supported by the Swedish Research Council, FORMAS. J.F. Modiano is supported in part by the Alvin S. and June Perlman Chair in Animal Oncology at the University of Minnesota. M. Breen is supported in part by the Oscar J. Fletcher Distinguished Professorship in Comparative Oncology Genetics at NC State University. K. Lindblad-Toh is supported in part by the Swedish Research Council and the European Research Council.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.