Genome-wide analysis of canine oral malignant melanoma metastasis-associated gene expression

Oral malignant melanoma (OMM) is the most common canine melanocytic neoplasm. Overlap between the somatic mutation profiles of canine OMM and human mucosal melanomas suggest a shared UV-independent molecular aetiology. In common with human mucosal melanomas, most canine OMM metastasise. There is no reliable means of predicting canine OMM metastasis, and systemic therapies for metastatic disease are largely palliative. Herein, we employed exon microarrays for comparative expression profiling of FFPE biopsies of 18 primary canine OMM that metastasised and 10 primary OMM that did not metastasise. Genes displaying metastasis-associated expression may be targets for anti-metastasis treatments, and biomarkers of OMM metastasis. Reduced expression of CXCL12 in the metastasising OMMs implies that the CXCR4/CXCL12 axis may be involved in OMM metastasis. Increased expression of APOBEC3A in the metastasising OMMs may indicate APOBEC3A-induced double-strand DNA breaks and pro-metastatic hypermutation. DNA double strand breakage triggers the DNA damage response network and two Fanconi anaemia DNA repair pathway members showed elevated expression in the metastasising OMMs. Cross-validation was employed to test a Linear Discriminant Analysis classifier based upon the RT-qPCR-measured expression levels of CXCL12, APOBEC3A and RPL29. Classification accuracies of 94% (metastasising OMMs) and 86% (non-metastasising OMMs) were estimated.

lymph nodes is commonly performed during staging, and other lymph nodes may be selected for examination on the basis of lymphangiography. However, where OMM metastasis to regional lymph nodes occurs, it is not always detected by examination of the mandibular and retropharyngeal lymph nodes 16,17 .
Whilst OMM represents the most common type of melanoma in dogs 3,4,18 , oral mucosal melanoma (most often developing in the palate and gingiva 19 ) accounts for only 1-8% of all human melanomas and around 0.5% of all human oral neoplasms 20,21 . However, human OMMs are also aggressive rapidly growing, invasive tumours that display metastatic rates of 66% (regional lymph nodes 22 ), 53% (lung 23 ), 36% (bone 23 ), and 20% (liver and brain 23 ), respectively. The 5-year survival rate is 15-25% 22,24 . Advanced human mucosal melanoma also has a low rate of response to adjuvant chemotherapy 25 . All human ethnic groups are affected by oral mucosal melanoma, although the Japanese appear to have an elevated susceptibility 26,27 .
Similarities between the somatic mutation profiles of human mucosal melanomas and canine OMMs suggests a possible overlapping molecular aetiology for UV-independent tumourigenesis. Putative activating mutations in Kit have been described in 7-16% of human mucosal melanomas 28,29 and in 12% of canine OMM 30 , and mutations in NRAS have been reported in 3.9% of canine OMMs 3 and in 10-22% of human mucosal melanomas [31][32][33] . Around 50% of human cutaneous melanomas have an activating BRAF mutation 34 , but BRAF mutations occur in only 4-9.5% of human mucosal melanomas 35,36 , and have not been found in canine OMMs 3 .
The metastatic cascade 37 comprises a series of steps, which are believed to be at least partially mediated by the acquisition of metastasis-associated genetic and/or epigenetic alterations additional to those that drive tumour development 38 . These somatic changes may affect gene expression, and metastasis-associated gene expression signatures have been identified for many human tumours [39][40][41] . Transcriptional profiling has defined the stages in human cutaneous melanoma development and progression as a series of distinct 'molecular events' , and implicated the involvement of sets of genes in the transition from primary to metastatic melanoma [42][43][44] . Gene expression signatures characteristic of human cutaneous melanoma metastases are detectable in primary cutaneous melanomas. A 1,864 gene expression signature derived from profiling cutaneous melanoma metastases was subsequently shown to delineate primary melanomas into two classes associated with significantly different relapse-free and overall survival 45 . Integration of functional and structural protein interaction data with primary and metastatic melanoma gene expression data enabled derivation of a panel of 6 genes which distinguished human primary and metastatic cutaneous melanoma and predicted melanoma-specific survival 46 . For human oral mucosal melanoma, comparative gene expression analysis of lymph node metastases and paired non-metastatic lymph nodes has elucidated the involvement of long non-coding RNAs in the regulation of metastasis-associated gene expression 47 . Metastasis-associated gene expression signatures identified in primary tumours may predict metastasis, and indeed assay of the expression of 15 genes is the basis of a routine test for human uveal melanoma metastasis 48 .
Gene expression profiling of canine cutaneous melanoma and melanocytoma has implicated the increased expression of genes involved in extracellular matrix-receptor interaction and the phosphoinositide 3-kinase/protein kinase B pathway in metastatic progression 49 , but no equivalent data is available for the more common canine OMM. Although we also recently identified pro-metastatic gene expression in an unrelated canine cancer (mast cell tumours 50 ), the overall low degree of 'overlap' between the metastatic gene expression signatures identified for different human solid tumours 41,51 affords the rationale for the gene expression profiling-based study of metastasis in different canine cancers, mirroring the approach adopted in the study of human tumour metastasis.
The high rate of mortality attributable to (conventional treatment-resistant) canine OMM metastasis is a significant welfare issue. It would be hugely beneficial for clinicians and owners alike to know whether a tumour was going to metastasise, while prevention of OMM metastasis would save the lives of most canine OMM patients. Establishing the role of dysregulated gene expression in canine OMM metastasis is an opportunity for identifying metastasis-associated biomarkers and possible anti-metastasis therapeutic targets.
In this project we performed comparative genome-wide expression profiling of archival biopsies of canine primary OMMs that metastasised and did not metastasise. We sought to identify metastasis-associated gene expression, and assess whether metastasising and non-metastasising OMMs could be delineated by the expression levels of genes associated with metastatic progression.

Results
Tumours subject to gene expression profiling. OMM biopsies from 42 dogs [29 bearing a metastasising (M) tumour and 13 bearing a non-metastasising (NM) tumour] qualified for the study. Through PowerAtlas 52 analysis of human tumour gene expression datasets (from the Gene Expression Omnibus 53 ), it was estimated that a 'Discovery Rate' of 73.1-81.7% at the 0.05 significance level would be afforded using 20 tumour samples in each of two 'outcome groups' .
Several requests for FFPE OMM biopsies from patients fulfilling the NM inclusion criteria were made to all UK veterinary university teaching hospitals and three large multidisciplinary private practices in the UK, but no additional NM OMM biopsies could be recruited because of the regular use of xenogeneic vaccination for patients regardless of metastatic status at the time of presentation. One NM OMM biopsy was subsequently excluded because of a sub-optimal RNA concentration. The integrity of each OMM RNA sample was determined and 20 M OMM and 12 NM OMM sample groups, with similar RNA integrity ranges, assembled (Tables S1 and S2). Tumours included in differential expression analysis. Tumours with exon-level probe set expression profiles that differed from those of the majority of the 32 OMMs were identified by review of associated sample quality metrics 54 , and 2 M and 2 NM OMMs excluded (Tables S2 and S3). The gene-level probe set expression values in 18 M OMMs and 10 NM OMMs, respectively, were compared for 13,422 Transcript clusters ('crosshyb_type' = '1'; 'category' = 'main'), for which the expression (above background) of ≥1 exon probe set was detected in at least 30% of the OMMs in the M and/or M groups.
The details of the dogs that bore the 18 M OMMs and 10 NM OMMs are presented in Table 1. Tumours in the M group were borne by 6 breeds and those in the NM group by 5 breeds. Tumours from Golden Retrievers and Labrador Retrievers, and cross bred dogs, were present in both groups. It is likely that the breed representation of the M and NM OMM groups reflect breed popularity and a predisposition to OMM development 2 ; for example, OMMs from 4 Golden Retrievers were present in both the M and NM tumour groups. The median ages of the dogs with M and NM OMMs were comparable (10.85 and 10.25 years, respectively), whilst 61% and 60% of the M and NM OMM dogs were male or neutered male, respectively. No association between OMM gene expression profile and either gender, or age at diagnosis, could be gleaned by hierarchical clustering of the 28 OMMs according to the expression values of the 20% of Transcript Clusters (2,684) that had the highest variance in expression signal. Biopsies from 5 of the 8 Golden Retrievers were grouped in a single cluster (together with biopsies from a Dachshund and Bullmastiff), whilst the 3 remaining Golden Retriever OMM biopsies were partitioned in a second large sub-cluster with biopsies borne by 6 other breeds (Fig. S1). Although this suggests that further investigation is warranted to assess the impact of 'genetic background' on the somatic molecular profile of canine OMMs, from the perspective of identifying metastasis-associated gene expression it is important that the 8 Golden Retriever OMM biopsies were divided equally between the M and NM OMM groups. Although (after 'outlier array' exclusion) the mean age of a NM OMM FFPE biopsy was 1.6× years higher than the mean age of a M OMM FFPE biopsy there was no correlation between FFPE tumour biopsy age and tumour RNA integrity (Spearman rank correlation coefficient = 0.09, two-sided p-value: 0.61; Table S2). www.nature.com/scientificreports www.nature.com/scientificreports/ Genes displaying differential expression in M and NM OMMs. In total, 331 Transcript clusters displayed a statistically significant difference in expression between M and NM OMMs. Of these, 191 exhibited increased expression in the M OMMs ( Fig. 1). A >1.5-fold difference in expression between the M and NM OMMs was observed for 12 genes ( Table 2). Significant sequence similarity to a mRNA encoded by a single canine gene was established for a Transcript cluster for which gene annotation was unavailable ( Table 2).
Functional annotation enrichment analysis. For Transcript clusters for which an Ensembl Gene ID 55 could be defined, the frequencies of functional annotations assigned to the Transcript clusters differentially expressed (300 of 331) between the M and NM OMMs were compared to those assigned to the Transcript clusters (11,842 of 13,422) for which the expression (above background) of ≥1 exon probe set was detected in at least 30% of the OMMs in the M and/or M groups. Over-represented amongst the genes exhibiting differential expression were 4 Gene Ontology Consortium biological processes and one KEGG pathway (Table 3).
RT-qPCR validation of differential expression. The expression, in the 18 M and 10 NM OMMs included in differential gene expression analysis, of the 12 genes which displayed >1.5-fold differences in expression between M and NM OMMs were assayed by RT-qPCR ( Table 4). The expression of Small nucleolar RNA SNORD61 (which had the lowest median level of expression of the 12 differentially expressed genes) could not be measured reliably in the OMM samples (Cq values >35 were obtained, or amplification was not detected), although its expression could be detected (Geomean Cq = 27.14) in cDNA prepared from 2 µg of a pool of OMM total RNA samples.
For most of the other 11 genes, Cq values ≥35/'undetermined' and/or a Cq standard deviation (triplicate assays) of >35 meant that valid gene expression measurements were recorded for slightly fewer than the 28 OMM biopsies profiled by microarray hybridisation. For 8 genes, the expression levels in individual OMMs measured by microarray (exon-level probe set) and RT-qPCR, respectively, were highly concordant, as were the NM OMM/M OMM fold differences in expression determined by the two techniques (Table 4). For SLC25A51 and SNORA76 there was a negative correlation between the gene expression values measured by microarray and RT-qPCR, and for SLC25A51 there was a difference between the 'direction' of NM: M differential expression as assessed by microarray and RT-qPCR, respectively (Table 4). These results suggest that for both SLC25A51 and SNORA76 the transcript quantified by PCR assay was different to that measured by microarray hybridisation. The differential expression of CXCL12 achieved statistical significance.

Class prediction analysis.
Based on an evaluation of the relative characteristics of the expression values measured for the 13,422 Transcript clusters 'present' in the OMMs, the optimal classification function for prediction of OMM 'metastatic status' (M or NM) was predicted to be Linear Discriminant Analysis (LDA; predicted accuracy = 0.761; lowest predicted accuracy -k-nearest neighbours = 0.408). The ranking of genes for their utility in class prediction may be based on the statistical significance of their difference in expression between classes. However, 'filter' gene selection methods may be based upon other metrics 56 , including fold-change differences in gene expression between classes 57,58 . Consequently, the efficacy of using the 3 genes shown (by RT-qPCR analysis) to exhibit >two-fold differential expression (Table 4, Fig. 2) for class prediction was tested. The relationships between the M and NM OMMs, in the context of the variation in the expression levels of CXCL12, APOBEC3A and RPL29, is effectively visualised by principal component analysis (Fig. 3).
Random sampling cross-validation was initially employed to test the performance of the LDA classifier. Two M OMMs and 1 NM OMM were randomly selected on each of 20 occasions, and the accuracy of their classification (as N or NM) measured after training the classifier using the remaining 15 M and 6 NM OMMs' expression values (Fig. 4). For M OMMs, mean and median classification accuracies of 100% were estimated, whilst mean www.nature.com/scientificreports www.nature.com/scientificreports/ and median classification accuracies of 65% and 100%, respectively, were estimated for the NM OMMs. In a subsequent evaluation of classifier performance by 'leave-one-out cross validation' , 94% of 17 M OMMs and 86% of 7 NM OMMs were correctly classified (Fig. 4).

Discussion
Malignant melanomas are the most common canine melanocytic neoplasm 3,4,18 and the most common oral malignancy in dogs 1 . At the present time there is no means to accurately predict if an individual OMM is one of the significant proportion of tumours that will metastasise, and metastasis is commonly underestimated during tumour staging due to a failure to sample all appropriate lymph nodes, and/or a failure of detection by standard cytology or histopathology. Metastasis is the most common cause of death in dogs treated for OMM, but there is currently no proven effective therapy to delay or prevent metastasis. With no means of accurately predicting OMM metastasis, it is possible that palliative adjuvant treatment is not prescribed for dogs with an unrecognised 'metastasising OMM' .
Molecular genetic and epigenetic 'events' that promote canine OMM metastasis may be both predictive indicators of canine OMM metastasis and the focus for therapeutics intended to prevent metastasis. Through comparative genome-wide gene expression profiling of 18 primary OMMs that metastasised and 10 primary OMMs that did not metastasise, in the present study we aimed to identify dog breed-independent OMM metastasis-associated gene expression.
Increased expression in M OMMs characterised 60% of the genes differentially expressed between M and NM OMMs. Further indication of the potential significance to OMM metastasis of the 'up-regulation' of pro-metastatic gene expression is that a similar proportion of the genes differentially expressed between human oral mucosal melanoma lymph node metastases and paired non-metastatic lymph nodes were expressed at an   Table 2. Genes displaying ≥ 1.5-fold differential expression between 18 metastasising and 10 non-metastasising OMMs as measured by microarray analysis. a Transcript cluster with no gene annotation. The most significant similarity between the sequence (spliced exons) of the Transcript cluster and a canine mRNA is listed. The significance of the sequence similarity is denoted by the E value and the length of the sequence alignment, and the proportion of the Transcript cluster sequence included in the alignment is stated. b Chromosomal location is denoted by the chromosome name and the gene start base co-ordinate 55   www.nature.com/scientificreports www.nature.com/scientificreports/ increased level in the lymph node metastases 47 . An understanding of the mechanisms by which changes in gene expression mediate OMM metastasis is afforded by functional annotation enrichment analysis.
Reversible protein phosphorylation is integral to intracellular signal transduction pathways, and requires the coordinated action of protein kinases and protein phosphatases. Disruption of the balance between phosphorylation and dephosphorylation is associated with carcinogenesis, and protein phosphatases have been recognised as tumour suppressors 59 . Aberrant expression of protein phosphatase inhibitors has been reported in a wide variety of human cancers 60 . Two genes (GPATCH2, PPP1R37) with the 'negative regulation of phosphatase activity' annotation showed decreased expression in the M OMMs and three genes (CASC5, CHP1, CSRNP2) with the same annotation showed increased expression in the M OMMs. PPP1R37, GPATCH2 (also known as PPP1R30), CASC5 (or PPP1R55) and CSRNP2 (or PPP1R72) encode proteins recognised as inhibitory regulatory subunits of phosphoprotein phosphatase 1 (PPP1) 61 , a protein serine-threonine phosphatase that regulates several members of the Transforming growth factor beta signalling pathway 62 , which promotes invasion and metastasis in advanced stages of cancer 63 . Inhibition of PPP1 by the regulatory subunit PPP1R1A in Ewing sarcoma has been shown to promote tumour growth and metastasis 64 . CHP1 (Calcineurin Homologous Protein 1) inhibits the serine/threonine protein phosphatase 2B 65 , which is responsible for the dephosphorylation (and thus activation)    www.nature.com/scientificreports www.nature.com/scientificreports/ of Nuclear factor of activated T cells (NFAT) transcription factors 1-4 66 . NFATs have been shown to have both pro-metastatic 67 and anti-metastastic 68 activities.
The Fanconi anaemia (FA) pathway is a multi-protein DNA repair pathway that resolves DNA interstrand cross-links encountered during DNA replication that would otherwise block replication and transcription, and lead to gross chromosome abnormalities 69 . The pathway forms part of the DNA damage response network 70 which maintains genome integrity. Four genes (FANCB, FANCI, RPA2, TOP3A) with the 'FA pathway' annotation showed elevated expression in M OMMs, whilst a fifth FA pathway gene (FANCC) displayed marginally decreased expression in the M OMMs. FANCB and FANCC encode proteins that form part of the FA core complex, an ubiquitin E3 ligase that monoubiquinates FANCI (and FANCD2), the ID complex, which is subsequently re-localised to the DNA lesion, and (with Group III FA proteins) co-ordinates DNA cross-link repair 71 . The RPA2 and TOP3A proteins perform DNA repair functions associated with the FA pathway. RPA2 is a subunit of the Replication Protein A complex, which is involved in DNA repair in the cellular response to DNA damage 72 , whilst TOP3A encodes a DNA topoisomerase that controls DNA topology during DNA repair 73 . The elevated expression of  On each of 20 occasions, the RT-qPCR-measured expression values of 3 genes (CXCL12, RPL29 and APOBEC3A) which displayed >two-fold differential expression between M and NM OMMs were used to predict the class (M = a square, and NM = a circle) of 3 randomly-selected OMMs (2 M OMMs and 1 NM OMM, which represent 10% of the OMMs and constituted a 'test set'). Prior to class prediction, the LDA classifier was trained using the gene expression data obtained for the remaining 90% of the OMMs (15 M and 6 NM OMMs, which constituted a 'training set'). (B) Leave-one-out cross-validation. The RT-qPCR-measured expression values of 3 genes (CXCL12, RPL29 and APOBEC3A) which displayed >two-fold differential expression between M and NM OMMs were used to predict the class of each of 24 OMMs (17 M = squares, 7 NM = circles). On each of 24 occasions, the class of one OMM was predicted after the LDA classifier had been trained using the gene expression data obtained for the remaining 23 OMMs. Row 1 depicts the actual class of each OMM and row 2 the predicted class of each tumour.
www.nature.com/scientificreports www.nature.com/scientificreports/ DNA repair pathway genes in human primary cutaneous melanoma has been associated with distant metastasis and poor prognosis 74 . Subsequently, the increased expression of FA and DNA damage response pathway genes has been reported in human 'high grade' primary cutaneous melanoma (which are associated with significantly reduced survival), relative to 'low grade' primary tumours 45 , and in human cutaneous melanoma metastases relative to normal skin 75 . As FA pathway gene expression correlates with FA pathway activity 76 , it appears that the FA pathway is activated in metastasising human cutaneous melanoma and canine OMM. FA pathway activation may be a response to increased genome instability in advanced melanomas, and may confer a selective advantage supporting metastasis to distant sites.
Two genes (COX10, MUL1)) with the 'mitochondrial fission' annotation displayed decreased expression in the M OMMs and one gene (MTFR1) with the same annotation showed elevated expression in the M OMMs. Increased expression of mitochondrial fission pathway genes has been associated with the invasiveness and metastasis of some cancers 77,78 . MTFR1 (Mitochondrial Fission Regulator 1) is upregulated in metastatic uveal melanoma 79 , and is a member of a 20-gene panel whose collective high expression is predictive of prostate cancer metastasis 80 . Conversely, suppression of MUL1 (Mitochondrial E3 Ubiquitin Protein Ligase 1) has been associated with the progression of human head and neck cancer 81 . COX10 (Cytochrome C Oxidase Assembly Homolog 10) is a member of a 14-gene classifier for colorectal cancer metastasis, identified by differential gene expression analysis of early and late stage primary colorectal cancer 82 .
Deregulation of the expression of genes involved in cholesterol homeostasis pathways has been associated with cancer development and progression 83 . In melanoma, the increased expression of 7 cholesterol synthesis pathway genes has been correlated with decreased patient survival 84 . Two genes (NR1H3, PLTP) with the 'Positive regulation of cholesterol efflux' and/or '(positive regulation of) cholesterol homeostasis' annotation were expressed at a decreased level in the M OMMs and three genes (APOE, MTTP, ABCA2) with one, or both, of the same annotations, showed elevated expression in the M OMMs. NR1H3 (or Liver X Receptor Alpha isoform) is a Nuclear Receptor superfamily transcription factor which when activated by oxysterol binding drives cholesterol efflux 85 . Agonist activation of the Liver X Receptor Beta isoform (NR1H2) has been shown to suppress the growth and metastasis of melanoma cells by transcriptional induction of apolipoprotein-E 86 . Reduced NR1H3 expression is predictive of decreased recurrence-free survival in muscle-invasive bladder cancer 87 , and associated with reduced overall survival in hepatocellular carcinoma 88 . APOE (apolipoprotein-E) is a lipid transport protein essential for the normal catabolism of triglyceride-rich lipoproteins 89 . Elevated APOE expression is associated with lymph node metastasis of human gastric cancer 90 and lung adenocarcinoma 91 , although it has been identified as a metastasis suppressor in vitro in human cutaneous melanoma 92 . PLTP (Phospholipid Transfer Protein) transfers phospholipids from triglyceride-rich lipoproteins to high density lipoprotein, and is involved in the uptake of cholesterol from peripheral cells and tissues. PLTP expression was increased in Grade IV human glioma relative to low grade glioma, and knockdown in vitro lead to the decreased migration of glioblastoma tumour cells 93 . In concept, the increased expression of PLTP in the NM OMMs observed in the present study may be consistent with the production by the tumours of interleukin 6, which has been shown to inhibit melanoma growth 94 . The increased expression of PLTP in the spontaneous regression phase of canine transmissible venereal tumour has previously been associated with increased IL-6 production 95 . The ATP-binding cassette transporter 2 (ABCA2) is a membrane-associated protein involved in sphingolipid transport. ABCA2 deficiency inhibits prostate tumour metastasis in vivo, potentially through reduction of the intracellular sphingolipid level 96 , whilst ABCA2 expression is increased in ovarian carcinoma metastases relative to primary tumours 97 .
In the present study, RT-qPCR analysis confirmed >two-fold differential expression between M and NM OMMs for 3 genes (CXCL12, APOBEC3A and RPL29). As fold change has been effectively employed to rank genes for their potential efficacy in gene expression level-based classification 57,58 , the 3 genes were selected for use in class prediction. Cross-validation was deployed in a preliminary evaluation to test the accuracy of a Linear Discriminant Analysis-based classifier featuring the 3 genes. The classification accuracies estimated were 94-100% (M OMMs) and 86-100% (NM OMMs), respectively. The LDA classifier performance will need to be validated by further retrospective, and prospective study. If the classification accuracy is confirmed, the cross-breed OMM metastasis-associated 3-gene expression signature would form the basis of an objective and quantitative predictive test for OMM metastasis that could make a significant contribution to the clinical management of canine OMM.
CXCL12 (C-X-C Motif Chemokine Ligand 12, or stromal cell-derived factor-1) is secreted by stromal cells and is a ligand for the G-protein coupled receptors CXCR4 and CXCR7 98 . Binding of CXCL12 to CXCR4 activates four signal transduction pathways that induce cytoskeletal rearrangement, cell growth, angiogenesis, and anti-apoptotic effects 98 . Interaction between CXCL12 and CXCR4 has also been shown to mediate metastasis, and direct metastatic dissemination to organs expressing high levels of CXCL12 99 . Blocking CXCL12 binding to CXCR4 reduced the migration of human uveal melanoma cells in vitro 100 , and pulmonary metastasis of murine B16 cutaneous melanoma cells 101 . Low CXCL12 expression in human primary cutaneous melanomas has been associated with poor prognosis 102 and shown to be predictive of metastasis 103 , and CXCL12 is one of 789 genes displaying reduced expression in 'high grade' human primary cutaneous melanomas that is a member of a 1,864 gene expression signature that delineates two classes of primary cutaneous melanomas with significantly different rates of metastasis 104 .
APOBEC3A (Apolipoprotein B mRNA Editing Enzyme Catalytic Subunit 3A) encodes a cytidine deaminase which preferentially binds to the sequence 5′-(C/T)TCA in RNA or single-stranded DNA and converts cytosine to uracil. Deregulated APOBEC3A expression in cancer is believed to induce double strand breaks in genomic DNA activating DNA damage response pathways 105 . The repair of such breaks triggers the formation of single stranded DNAs which are substrates for APOBEC3A-mediated hypermutation, such that 5′-(C/T)TTA APOBEC3A mutation signatures occur in clusters (on one DNA strand) in multiple human cancers 106 . The extent of APOBEC-associated mutations correlates with APOBEC mRNA expression levels 107 . APOBEC3A-mediated mutagenesis occurs at different stages in different cancers 108 , and is thought to drive tumour evolution, including promoting metastasis 109 . (2019) 9:6511 | https://doi.org/10.1038/s41598-019-42839-x www.nature.com/scientificreports www.nature.com/scientificreports/ RPL29 (Ribosomal protein L29) encodes a component of the 60 S ribosomal subunit. Beyond their role in ribosome assembly and protein translation, differential expression of ribosomal protein genes in cancer has been associated with ribosome-independent regulation of cell growth and proliferation, apoptosis, invasion and metastasis 110 . Ribosome-free ribosomal proteins have been implicated as being both oncogenic and tumour-suppressors 110 . Silencing of RPL29 suppressed the proliferation of human pancreatic tumour cells and enhanced apoptosis 111 suggesting an involvement in cell proliferation. However, RPL29 silencing had no effect on the viability of human metastatic melanoma cells 112 . The expression of a specific ribosomal protein has been shown to be a prognostic indicator for multiple human cancers 113 .
In human medicine, anti-angiogenics and matrix metalloproteinase inhibitors have been licenced for treatment of tumour metastasis 114 . Cellular receptors in signal transduction pathways that control cell to cell and cell to ECM adhesion are targets for anti-metastatics in development 115 . Pro-metastatic gene expression in OMMs is a potential target for anti-metastasis therapeutics. Targeting the interaction between the CXCR4 receptor and CXCL12 has been evaluated as a strategy for inhibiting CXCR4-CXCL12 axis-mediated melanoma metastasis. Small molecule inhibitors of CXCR4 were shown to be effective at disrupting the liver metastasis of uveal melanoma cells in mice 116 , and migration of human cutaneous melanoma cells in vitro 117 . Chemical inhibition of the mutational activity of APOBEC3A (and APOBEC3B) is being evaluated as a cancer therapeutic 118 . Furthermore, demonstration of microRNA post-transcriptional regulation of APOBEC gene expression 119 suggests the use of miR-mimics 120 as a potential means of APOBEC deaminase inactivation 119 . Intriguingly, it has been postulated that APOBEC3A-mediated hypermutation could generate new tumour-specific antigens thereby enhancing the efficacy of immune stimulation therapies 106 . The potential for suppressing melanoma metastasis through activation of Liver X Receptors (LXR) using synthetic agonists has been investigated. Activation of the LXR Beta isoform (NR1H2) was shown to inhibit human and murine cutaneous melanoma cell migration in vitro and murine cutaneous melanoma cell metastasis in a mouse xenograft model 86,121 . If NR1H3 is the predominant LXR Receptor isoform expressed it remains to be seen if LXR agonists would initiate receptor activation with anti-metastatic effect. The up-regulation of FA pathway DNA repair genes in human cutaneous melanomas and canine OMM may be both pro-metastatic by negating the impact of increased genome instability, and contribute to melanoma metastases chemoresistance. Consequently, in concept, the use of FA pathway inhibitors 122 may be an option for treatment of melanoma.
The major limitation of this study was that the number of NM OMMs included was restricted by the need to attempt to ensure that a primary tumour classified as 'non-metastasising' was only classified as such because it did not exhibit pro-metastatic gene expression as opposed to its metastasis potentially being prevented by a systemic adjuvant therapy. The consequence of the widespread use of the xenogenic melanoma vaccine in the UK was the exclusion of OMM biopsies from many dogs that had not developed metastatic disease. Furthermore, due to difficulties associated with collecting fresh canine tumour biopsies, FFPE biopsies of primary canine OMM (surgically removed at specialist veterinary oncology centres) were used in the study. However, data that is both biologically authentic and clinically-relevant has been obtained by Affymetrix microarray-based gene expression profiling of FFPE tissues 123,124 . Validation of the predictive accuracy of the 3-gene LDA classifier will be achieved through further retrospective, and prospective, studies featuring larger numbers of (optimally freshly collected) M and NM OMM biopsies. Ultimately, experimental investigations involving canine OMM cells will be necessary to confirm the functional consequences (e.g. in regard to cell migration) of metastasis-associated differential gene expression that are anticipated given gene function(s) and (in some cases) prior in vitro and in vivo study evidence.
Several of the differences in gene expression observed between primary canine OMMs that metastasised and OMMs that did not metastasise in this study have previously been associated with human cutaneous melanoma metastasis. The genes involved have been targets in proof of principle trials of potential anti-metastatic melanoma therapeutics. Other genes that exhibit differential expression between metastasising and non-metastasising primary canine OMMs may represent potential new targets for both canine and human cancer drugs. The results obtained in the present study suggest that OMMs in dogs may be responsive to anti-human melanoma metastasis therapeutics currently in clinical trials, or being evaluated through preclinical in vitro and in vivo model studies, and may be as likely as human melanomas to be responsive to therapeutics whose efficacy for treatment of human melanoma has yet to be investigated. Quantification of the expression of 3 genes, each of which displays a greater than two-fold differential expression between canine OMMs that do and do not metastasise, may be the potential basis for a test that would accurately predict canine OMM metastasis, and thereby assist a clinician to make an informed decision about the most appropriate treatment for a canine OMM patient.

Materials and Methods ethics statements.
This study was conducted with the approval of the ethics committees of the Animal Health Trust and the University of Liverpool, respectively. All methods were performed in accordance with the relevant guidelines and regulations. An OMM biopsy was only included in the study with the informed, written consent of the owner of the dog who bore the tumour. The treatment that a OMM patient received was unaffected by the inclusion of a biopsy of their tumour in the study. www.nature.com/scientificreports www.nature.com/scientificreports/ diagnostic imaging-confirmed metastasis, and who died or were euthanased because of OMM metastasis <500 days after surgery/biopsy (irrespective of adjuvant chemotherapy, including prednisolone, xenogeneic vaccine, and/or radiotherapy). 'Non-metastasising' (NM) OMM biopsies were from dogs without pathological analysis and imaging-confirmed metastasis, whom did not receive any adjuvant therapy, and who were still alive >540 days post-surgery/biopsy.

RNA sample selection.
Reverse transcription-quantitative PCR (RT-qPCR) measurement of a 126 bp fragment of a 130-150 bp short interspersed nuclear element (SINE) that occurs every 5-8.3 kb in the canine genome 125 was employed to assess the integrity of each FFPE RNA sample as described previously 50  Microarray data analysis. Microarray data was processed using the ' Affymetrix Expression Console Software 1.3' (ThermoFisher Scientific, Paisley, UK). 'Outlier arrays' were first identified by review of exon-level probe set expression values created (using the RMA algorithm 126 ) by quantile normalisation, log 2 transformation and signal summarisation, respectively. Arrays with ≥1 sample quality, labelling quality or hybridisation quality metric value ≥2 standard deviations from the mean for all the arrays 54 were excluded, and the raw probe-level signal intensity data for the remaining arrays re-processed to generate quantile normalised and log 2 -transformed exon and gene-level probe set expression values. Gene-level probe sets ('Transcript clusters') annotated as 'crosshyb_ type' = 1 (unique hybridisation target) and 'category' = 'main' 127 , and for which the expression above background (detection above background p-value < 0.01 128 ) of ≥l exon probe set could be detected in at least 30% of the OMMs in the M and/or NM group, were judged to be expressed in the OMMs and were included in further analyses.
Hierarchical clustering (average linkage; similarity metric = Pearson Correlation Coefficient) was performed using Cluster 129 , and Principal Component Analysis using the R stats package function prcomp 130 , in order to view the relationships between OMMs on the basis of their gene-level expression profiles. A two-tailed t-test for unpaired data was employed to identify genes exhibiting differential expression between M and NM OMMs that was statistically significant, adjusting P-values by permutation testing 131 to correct for false positives arising from multiple testing. BLAST similarity search (against canine and human mRNAs and non-coding RNAs) was employed to attempt to establish the potential identity of Transcript clusters representing 'predicted genes' , or for which gene annotation was unavailable.
Functional annotation analysis. The biological processes and pathways affected by the differences in gene expression observed between M and NM OMMs were identified using DAVID 132,133 . The functional annotations associated with differentially expressed genes were compared with those ascribed to all Transcript clusters ('crosshyb_type' = 1 and 'category' = 'main') for which the expression of ≥1 exon probe set was detected above background in ≥30% of the tumours in the NM and/or M OMM cohort, and over-represented biological processes and pathways identified.
Reverse transcription-quantitative PCR (RT-qPCR). RT-qPCR was employed for validation of differential gene expression. A unique region within the exon probe set that displayed the largest statistically significant fold-difference in expression between M and NM OMMs was the template for design (using Beacon Designer; Premier Biosoft, Palo Alto, USA) of a TaqMan or SYBR Green PCR assay. The expression level of each gene in each OMM sample was measured as the geometric Cq value calculated from triplicate PCR reactions performed using preamplified cDNAs, prepared from cDNAs previously screened for PCR inhibitors. A 71 bp fragment of a SINE 125 , present in the 3′-untranslated region of hundreds of canine mRNAs, was also assayed as a 'reference gene' for target gene expression measure normalisation 134 . Geometric mean Cq ≥35 were considered an unreliable measurement of gene expression and were excluded, as were the results of triplicate PCR assays with a Cq standard deviation of >0.5. Additional information is provided in the Supplementary Information.
RT-qPCR data analysis. Relative quantification of gene expression was performed using a modification of the delta-delta-Ct method which accounts for differences between the amplification efficiencies of target and www.nature.com/scientificreports www.nature.com/scientificreports/ potentially multiple 'reference' genes 135 . Using qbase + (Biogazelle, Gent, Belgium), target gene geometric mean Cq values were converted to relative gene expression measurements ('Normalised Relative Quantity; NRQ 135 ) by application of a canine SINE 125,134 geometric mean C q -derived normalisation factor. A two-tailed t-test for unpaired data was performed on NRQ log 10 transformations to identify genes exhibiting statistically significant differences in expression between M and NM OMMs.
Class prediction analysis. The R package SPreFuGED 136 predicts the performance in class prediction of representatives of 10 classification functions and, through evaluation of the expression values obtained for the Transcript clusters expressed in the OMM, was employed to identify the optimal classification function for prediction of OMM 'metastatic status' (M or NM) on the basis of OMM gene expression. Linear Discriminant Analysis-based class prediction was undertaken using the R Package MASS 130 lda function implemented in the R environment for statistical analyses 137 . The accuracy of class prediction, performed using RT-qPCR-derived gene expression measurements, was estimated by cross-validation. In random sampling cross-validation, on each of 20 occasions the class (M or NM) of two M and one NM OMM, which constituted a 'test data set' , were predicted after the classifier had been trained using the gene expression values ('training data set') obtained for the remaining OMMs (90% of the tumours). In leave-one-out cross-validation, the class of each OMM was predicted after the classifier had been trained using the remaining (n-1) OMM gene expression data.

Data Availability
The microarray gene expression data generated during this study is available from the Gene Expression Omnibus repository (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE129750).