Gene expression profiling of spontaneously occurring canine mammary tumours: Insight into gene networks and pathways linked to cancer pathogenesis

Spontaneously occurring canine mammary tumours (CMTs) are the most common neoplasms of unspayed female dogs leading to thrice higher mortality rates than human breast cancer. These are also attractive models for human breast cancer studies owing to clinical and molecular similarities. Thus, they are important candidates for biomarker studies and understanding cancer pathobiology. The study was designed to explore underlying molecular networks and pathways in CMTs for deciphering new prognostic factors and therapeutic targets. To gain an insight into various pathways and networks associated with the development and pathogenesis of CMTs, comparative cDNA microarray expression profiling was performed using CMT tissues and healthy mammary gland tissues. Upon analysis, 1700 and 1287 differentially expressed genes (DEGs, P ≤ 0.05) were identified in malignant and benign tissues, respectively. DEGs identified from microarray analysis were further annotated using the Ingenuity Systems Pathway Analysis (IPA) tool for detection of deregulated canonical pathways, upstream regulators, and networks associated with malignant, as well as, benign disease. Top scoring key networks in benign and malignant mammary tumours were having central nodes of VEGF and BUB1B, respectively. Cyclins & cell cycle regulation and TREM1 signalling were amongst the top activated canonical pathways in CMTs. Other cancer related significant pathways like apoptosis signalling, dendritic cell maturation, DNA recombination and repair, Wnt/β-catenin signalling, etc. were also found to be altered. Furthermore, seven proteins (ANXA2, APOCII, CDK6, GATC, GDI2, GNAQ and MYH9) highly up-regulated in malignant tissues were identified by two-dimensional gel electrophoresis (2DE) and MALDI-TOF PMF studies which were in concordance with microarray data. Thus, the study has uncovered ample number of candidate genes associated with CMTs which need to be further validated as therapeutic targets and prognostic markers.

a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 (2DE) and MALDI-TOF PMF studies which were in concordance with microarray data. Thus, the study has uncovered ample number of candidate genes associated with CMTs which need to be further validated as therapeutic targets and prognostic markers.

Introduction
Spontaneously occurring mammary tumours in dog have been demonstrated as useful models for human breast cancer studies owing to their similarities in histological features and disease biology; associated risk factors, clinical progression and response to treatment; biomarkers and molecular targets etc. [1][2]. Although widely used, xenogeneic, as well as, transgenic mouse tumour models fail to mimic various features of human breast cancers like steroid hormone dependency, tumour microenvironment, heterogeneous behaviour etc. [1][2]. With complete sequencing of dog genome revealing close similarities between dog and humans at genetic level, spontaneously occurring canine mammary tumours (CMTs) have emerged as attractive alternatives to artificially induced tumours in mice.
CMTs account for approximately 50% of all tumours of female dogs with approximately 50% cases being malignant in nature [3][4]. Dogs with poorly differentiated tumours have increased risk of recurrent or metastatic disease, with 90% recurrence rate for the most dedifferentiated tumours [4]. In unspayed female dogs, CMTs are the most common neoplasms leading to at least three times higher mortality rates than human breast cancer [5]. Presently, the genes accountable for the aggressive behaviour of mammary cancer are not very clear, and poor prognosis associated with malignant mammary cancers emphasize the necessity to unravel the underlying pathways and genes which could act as targets for therapy. The molecular mechanisms, networks and pathways contributing to the biological behaviour of CMTs are poorly understood and there is a lack of knowledge about reliable tumour markers [6]. The detailed characterization of the dysregulated genes and careful mining of the gene networks and pathways could also be of great help in identification of diagnostic, as well as, prognostic biomarkers.
The study was designed with the objective of identification of new prognostic factors and targets for therapy, as well as, underlying pathological mechanisms and networks. To have an insight into various pathways and networks associated with development and pathogenesis of CMTs, comparative cDNA microarray analysis was performed. The gene expression profiles of malignant and benign CMT tissues were compared with healthy mammary tissues. To gain a deeper knowledge, differentially expressed genes (DEGs), identified from this analysis were further subjected to functional annotation using the Ingenuity Systems Pathway Analysis (IPA) tool. IPA was used to detect upstream regulators, pathways and networks associated with CMTs. The study has uncovered numerous candidate genes involved in pathogenesis of CMTs, which need to be further validated as therapeutic targets and prognostic markers for mammary cancer.

Materials and methods
The overview of the methodology used for gene expression profiling of CMTs is depicted in

Tissue samples
All experimental procedures involving animals were in accordance with Breeding of and Experiments on Animals (Control and Supervision) Amendment Rules, Government of India, 2005 and were approved by Institute Animal Ethics Committee (IAEC) of ICAR-Indian Veterinary Research Institute (ICAR-IVRI), Izatnagar. Cancer tissue samples (n = 10) were collected from clinical cases of CMTs which were referred for surgery to 'Referral Veterinary Polyclinic', ICAR-IVRI, Izatnagar. Reference healthy tissues (n = 3) were collected from healthy dogs with prior written consent from the owner. Animals were given local anaesthetic (1% lignocaine) prior to tissue collection. Tissue sections were collected immediately after surgical resection and stored in RNAlater at -80˚C for RNA isolation. Part of tissue section was also preserved in neutral buffered formalin for histopathology. For proteomic studies, tissues were preserved at -80˚C.

Histopathological classification
For histopathology, formalin-fixed tissues were paraffin-embedded and cut into 5 μm sections. Sections were then mounted on 3-aminopropyl-triethoxy-silane (APTES) coated slides, and air-dried overnight at 37˚C. Prepared slides were deparaffinized in three washes of xylene (10 min each), and rehydrated in graded concentrations of ethanol. Slides were then stained with Eosin-Hematoxylin and studied for identification for tumour type and grade. Histopathological classification of the paraffin embedded tissue sections was done as per Goldschmidt et al. [7]. The tumours of mammary gland used in the current study were classified into malignant (n = 6) and benign (n = 4) nature. The malignant tumours were classified histopathologically as tubuloacinar solid carcinoma, papillary and squamous cell mixed carcinoma, malignant myoepithelioma, squamous cell carcinoma, lipid rich carcinoma and fibrosarcoma. Benign tumours were classified histopathologically as adenomas, papillary adenoma and cystadenomyoepithelioma. Histological grades and types of cancer tissues used in study are provided in supplementary information. Healthy mammary tissue samples (n = 3) were also examined histologically for confirming absence of malignancy and any other pathological condition. Histopathological details of tumour tissues used for the study are provided in S1 Table. concentration of total RNA extracted was checked spectrophotometrically. The quality check of the isolated RNA was performed in Agilent 2100 bioanalyzer as per the manufacturer's recommendations using the Agilent RNA 6000 Nano Kit. RNA samples having RNA integrity number (RIN) greater than 7 were used for further studies.

Microarray hybridization and data analysis
For performing microarray studies, 200ng of RNA was reverse transcribed into cDNA using LowInput QuickAmp labelling kit (Agilent, USA) and further converted into cRNA to cause incorporation of nucleotides containing cyanine 3 (Cy-3-dCTP) as per the manufacturer's instructions. Labelled cRNA was purified using QIAquick purification columns. The quality and concentration of the labelled cRNA was checked spectrophotometrically and hybridization experiments were performed at 65˚C for 17 hours on canine (v2) gene expression 44K array comprising of 43803 canine probes (Agilent Technologies, Cat No. G2519F-021193). Each sample was hybridized to set of duplicate arrays. After generating the microarray scan images, the signal intensities were extracted using feature extraction software (version 10.7.3). The quality of hybridization and overall chip performance was monitored by visual inspection of both internal quality control checks and the raw scanned data. The data generated after feature extraction was exported into GeneSpring software (version 13.0) to identify the differentially expressed genes (DEGs) between malignant v/s healthy and benign v/s healthy. Mean normalised signals from malignant (n = 6), benign (n = 4) and healthy (n = 3) tissues were clubbed as a single experimental group and named as malignant, benign, and healthy respectively and used for comparative expression analysis. Paired t-test using storey with bootstrapping correction (P value � 0.05) was used to identify differentially expressed mRNAs between malignant and benign tumours with respect to healthy mammary tissue. The gene expression fold change between tumour and healthy tissue was calculated by comparison of the normalised signal intensity values. Data was submitted to NCBI GEO database (Accession number GSE104733).

Quantitative real time PCR (qPCR)
To verify the microarray gene expression data, 14 genes were selected for validation using qPCR. The genes selected for validation of microarray results included the top up-regulated genes, some randomly selected genes, as well as, genes having potential role in cancer pathogenesis . These genes included MMP9, CHI3L1, BIRC5, BIRC2, TLR2, CTSS, PGAM1, KIF11,  COL11A1, SRFP2, TOP2A, CPA2KL, CDCA3, RAB31 genes. The primers were designed for qRT-PCR analysis using the Integrated DNA technologies-PrimerQuest Tool. The details of primers sequences used for the study are mentioned in S2 Table. The cDNA was synthesized using Revert Aid First Strand cDNA synthesis kit (Thermofischer Scientific, USA) according to the manufacturer's instructions and qRT-PCR was performed using Applied Biosystems 7500 Fast system using 2X SYBR Green Master mix (Sigma Aldrich, USA). Gene expression in each sample was normalized against the expression of housekeeping gene (β-actin). The relative expression of each sample was calculated using the 2 −ΔΔCT method with healthy mammary tissue as calibrator and log 2 fold change was plotted.

Two-dimensional gel electrophoresis (2DGE) of CMT and healthy mammary tissue
Tissue samples (~20mg) were homogenized in lysis buffer [8M urea, 2M thiourea, 4% CHAPS, 30mM Tris, pH 8.5] and kept on rotator for 15 mins followed by sonicating twice for 20s. The samples were then centrifuged at 13000 rpm for 20 mins at 4˚C. Crude tissue lysate was further subjected to clean up using Ready Prep 2-D clean up kit (Bio-Rad, USA). Post clean up protein was quantified using 2-D Quant Kit (GE Healthcare, USA). Total 400μg of protein sample in 125μl destreak solution (GE healthcare, USA), containing Bio-Lyte 3/10 Ampholyte (40%) (@0.2-2% final concentration), was loaded on 7 cm ReadyStrip IPG Strip [pH 3-10] (Bio-Rad, USA). Protein was loaded onto the IPG strips by a passive rehydration method. After rehydration, isoelectric focusing (IEF) was performed on an Ettan IPGphore III apparatus at 20˚C with increasing voltages 150V-200Vh, 1000V-200Vh, 5000V-4000Vh, 5000-1250Vh. The focused IPG strips were subjected to reduction with 1% w/v DTT in 10 mL of equilibration buffer (6 M urea, 50 mM Tris-HCl pH 8.8, 30% v/v glycerol, and 2% w/v SDS) followed by alkylation with 2.5% w/v iodoacetamide in the same buffer. The strips were then placed on the top of 12% resolving gels and fixed by ReadyPrep Overlay Agarose (BioRad, USA). For spot picking, gels were stained with commassie brilliant blue (R-350) stain. ImageMaster 2D Platinum v7.0 software was used for analysis and spot selection. After an automated spot detection, the protein spots were checked manually in order to eliminate any possible artifacts such as gel background or streaks and selected spots were picked. Proteins at the excised spots were then subjected to MALDI-TOF PMF analysis for their identification.

Functional interpretation of microarray data and network analysis
To understand the functional significance of dysregulated genes in malignant and benign tumours of mammary gland, Ingenuity Pathway Analysis (IPA version 01-07) was used. Using the IPA tool, networks were generated based on an algorithmically generated score. The Zscore, a numerical value was used to rank networks according to how relevant they were to genes presented within the data set. Canonical pathways significant to the input data set were identified from the IPA library of canonical pathways based upon 2 parameters, viz., (1) The ratio of the number of genes within the data set mapping to the pathway divided by the total number of genes mapping to the canonical pathway and (2) a P value (calculated based upon Fischer's exact test) determining the probability that each bio-function assigned to that dataset and the canonical pathway is not due to chance alone. Upstream regulators are defined as the genes that affects the expression of numerous other genes, while canonical pathways are the idealized or generalized pathways that represent common properties of a particular signalling module or pathway.

Differentially expressed genes (DEGs) in canine mammary tumour tissues
To identify differentially expressed candidate genes, gene expression profiles of malignant and benign CMT tissues were comparatively analyzed with healthy mammary tissues by cDNA microarray technique. On analysis, 1700 and 1287 DEGs (P � 0.05) were identified in malignant and benign tumours of mammary gland respectively. Of these, 765 genes were up-regulated and 935 genes were down-regulated in malignant cases. Among the benign tumours, 744 genes were found to be up-regulated and 543 genes were down-regulated. Upon Venn analysis of up-regulated genes (log FC>1), it was observed that 269 genes were commonly up-regulated in all the benign tissues (Fig 2A), whereas 90 common genes were up-regulated in all the malignant tissues. Out of these 32 genes, displayed in Fig 2B, were commonly dysregulated in all the malignant and benign tissues studied. However, there was no such gene which was uniquely up-regulated in all the malignant tissues, as the dysregulated genes which were common in all the malignant tissues were also present in either one or all the benign tissues ( Fig  2C). Among malignant tumours, maximum fold change was observed for COL11A1 gene (Log 2 FC = 4.8), while among benign tumours maximum fold change was observed with MMP3 gene (Log 2 FC = 6.6). Top five up-regulated genes in malignant tumours were COL11A1, SFRP2, LCN2, COL2A1 and H19, while top up-regulated genes in benign tumours were MMP3, MMP1, AREG, PTHLH and SFRP2.

Validation of DEGs by real time PCR
Microarray results were validated by carrying out qPCR analysis for selected genes viz., MMP9, CHI3L1, BIRC5, BIRC2, TLR2, CTSS, PGAM, KIF11, COL11A1, SFRP2, TOP2A, CPA2KL, CDCA3, RAB31 in the same tissues used for microarray gene expression analysis. The genes selected for validation of microarray results included some of the top up-regulated genes, some randomly selected genes, as well as, some genes having potential role in CMT pathogenesis. The expression levels of the selected genes, in both microarray and qPCR technique were compared and results of both the studies revealed dysregulation of the target genes in a similar fashion (Fig 3).

Differentially up-regulated proteins in CMTs
For identification of up-regulated proteins in malignant mammary tumours, total protein content from malignant and healthy mammary tissues was resolved separately using 2DGE. Upon analysis using ImageMaster program, 298 individual spots were identified in the 2-D gel from malignant tissue, whereas, 328 spots protein spots were identified in healthy mammary tissue. Out of these, 178 proteins spots matched in both cancer and healthy mammary tissues. After analysis of 2-D gels, we observed a total of 7 differentially expressed proteins (DEPs) having greater than three times higher expression level in malignant tissue as compared to healthy tissue (Fig 4). Upon MALDI-TOF PMF analysis, these proteins were confirmed as ANXA2, APOCII, CDK6, GATC, GDI2, GNAQ and MYH9. These 7 proteins were also found to be upregulated in malignant tissues by microarray analysis, thus providing further validation of  (1)(2)(3)(4). Analysis revealed that there was no such gene which was uniquely up-regulated in all the malignant tissues, as the dysregulated genes which were common in all the malignant tissues were also present in either one or all the benign tissues.
https://doi.org/10.1371/journal.pone.0208656.g002 microarray data and confirming the overexpression of these genes at protein level in malignant canine mammary cancers.

Relevant functions, pathways and biological networks in CMTs
Significant DEGs (P � 0.05) were used for functional annotation by Ingenuity Pathway Analysis (IPA) tool to explore the relevant biological pathways and functions of candidate genes involved in pathogenesis of CMTs. In IPA, the reference knowledge base was kept as Homo  sapiens and Mus musculus, as most functional studies on cancer has been carried out in these organisms only and there is a severe lack of functional studies in CMTs. Each gene symbol was mapped to their corresponding identifiers in the Ingenuity knowledge base. Of the total 2336 genes uploaded into IPA, 2259 genes were mapped in the IPA knowledge base, whereas 77 genes remained unmapped. Networks were generated based on an algorithmically generated score based on connectivity between genes. A numerical value (Z-score) was used to rank networks according to how relevant they were to genes presented within the data set. Identified networks between genes are presented as a graph, indicating relationships between different genes in the data set. Genes are presented as nodes and relationship between two is indicated as a line.

Canonical pathways deregulated in CMTs
Upon investigating the interactions using the IPA software, 2259 functional pathway eligible genes were found. Canonical pathways, defined as the idealized or generalized pathways representing common properties of a particular signalling module or pathway, were identified by IPA. Top activated canonical pathways in malignant tumours, (z-score > 2) included cyclins and cell cycle regulation, apoptosis, and PPAR signalling pathway etc., and are shown in Table 1 and Fig 5B. Top inhibited pathways in malignant cancers included intrinsic prothrombin activation pathway and neuroprotective role of THOP1 in Alzheimer's disease. For benign tumours, TREM1 signalling, dendritic cell maturation, Integrin linked kinase (ILK) signalling pathway etc., were amongst the top activated canonical pathways, while top inhibited pathways were related to cell cycle regulation, complement system and neuroprotective role of THOP1 in Alzheimer's disease ( Fig 5A and Table 2). The functional relevance of the dysregulated genes in our dataset, related to the top activated pathways in malignant and benign CMTs is depicted in Figs 6-9.
The DNA damage checkpoint, occurs at the end of G2 and prevents DNA damage that might have occurred during replication. Any defects in this pathway results in genome instability and leads to carcinogenesis. Upon comparison of this pathway in malignant and benign tumours, it was seen that, in malignant tissues ATM/ATR kinases having role in controlling BRCA1 phosphorylation in vivo are involved in regulation of G2/M DNA damage checkpoint, while in benign tissues, activation of alternative AURKA-PLK1 pathway may play that role (Fig 10). Other important canonical pathways related to cancer such polo like kinase related pathways, cyclins and cell cycle regulation and Wnt/β-catenin signalling pathways were also compared between malignant and benign CMTs (S1-S5 Figs), highlighting the differential involvement of dysregulated genes contributing to malignant or benign nature.

Upstream regulators identified in CMTs
Upstream regulators have been defined as genes that direct the expression of numerous other genes such that dysregulation of an upstream regulator affects the expression level of downstream gene it affects. In malignant tumours, 24 upstream regulators were identified which were mainly enzymes, transcription regulators, G protein coupled receptors, kinases, etc. (Table 3). The prediction of activation state was based upon the global direction of changes of the modulated genes. The activation Z-score, indicating whether the observed gene responses to upstream regulators agree with expected changes derived from the literature that accrued in the Ingenuity Knowledge Base, was used to predict the activation state. Pathway analysis revealed that secreted phosphoprotein (SPP1) gene was the top activated upstream regulator among malignant tissues and affected expression of various genes including expression levels of breast cancer type 1 susceptibility protein (BRCA1) and cell-division cycle protein 20 (CDC20). Various genes affected by SPP1 are displayed in Fig 11A. Analysis of the malignant data set also revealed F2 (coagulation factor II, Thrombin) as another top activated upstream regulator, which was connected with overexpression of chemokines CCL2, CXCL8, CXCL10, CCL8 and their ligands CXCR4, F2RL2, ITGB1 ( Fig 11B). Other upstream regulators like GCG and PTF1A were found to be downregulated and inhibited the expression of various associated genes as depicted in Fig 12A and 12B. In benign tumours, 25 upstream regulators were identified. Of these 22 were activated, and 3 were inhibited based on Z-score cut-off ( Table 4). The genes affected by the top activated upstream regulators in benign tissues i.e., AREG and TLR2 are displayed in Fig 11C and 11D. Analysis also revealed INSIG1 and NROB2 as the top inhibited upstream regulators among the benign cancers.
Common regulators that were affecting the dysregulation of genes in malignant tumours had a consistency score ranging from 9.9 to -13.584. Further, 47 regulators have a combined consistency score of 9.9 which affected a total of 142 genes in our data set, followed by a group of 27 regulators which affected 173 genes having an overall consistency score of 9.7. These two groups of regulators affected diseases related to female genital neoplasm, hematologic cancer, lymphocytic neoplasm, lymphoid cancer, tumourigenesis of reproductive tract and hematologic cancer, lymphocytic neoplasm, lymphoid cancer. Similarly, common regulators were also identified for benign tumours with consistency score ranging from 3.7 to -12.2. The regulators with the highest score were related to cell movement, invasion of cells and cell movement of tumour cell lines.

Functional relevance of activated pathways and networks in CMTs
For knowing the functional relevance of activated pathways and networks, we investigated the disease and bio-functions related to the top activated, as well as, inhibited pathways using IPA tool. The top most affected diseases with respect to genes in our malignant data set were cancer, organismal injury and abnormalities, etc. The activated bio-functions in malignant tissues were attributed to cell to cell signalling, cellular movement, inflammatory response, immune cell trafficking, DNA replication, recombination, and repair. Genes that were seen to be downregulated in our analysis were related to bio-functions like exocytosis, secretory pathways, bacterial infections, female genital neoplasms etc. (Table 5). Among the benign tumours, top most activated bio-functions were related to immune response of leukocytes, immune response of     (Table 6).
To further understand the inter-relationship between genes that were dysregulated in malignant and benign tumours, networks were identified to present how the genes interact and coordinate their roles in a specific pathway. In malignant tumours, 25 networks were identified that were related to cell cycle, cellular assembly and organization, DNA replication, recombination, and repair. The network with a maximum score was related to cell cycle, cellular assembly and organization, DNA replication, recombination, and repair. Other top networks were related to developmental and hereditary disorders, organismal injury, abnormalities, post-translational modification, DNA replication, recombination and repair, cell cycle etc. Networks presenting interrelationship of genes dysregulated in benign tumours were 25 in total with a maximum score for the network showing genes related to cellular assembly and organization, cell cycle, DNA replication, recombination, and repair. Other notable identified networks were related to cancer, organismal injury, abnormalities and embryonic development, cell injury, DNA replication, cell signalling and cell cycle among others. It was observed that the top scoring network amongst the malignant mammary cancers have a central node of mitotic checkpoint serine/ threonine kinase (BUB1B), while top network in the benign tumours have a central node of vascular endothelial growth factor (VEGF), as shown in Fig 13.

Discussion
Mammary tumour is an important neoplastic disease of dogs, besides its role as a model for human breast cancer studies. Our study has revealed new biomarkers, networks and pathways dysregulated in canine mammary tumours. Previous work on expression profiling of canine mammary tumour has suggested that separate pathological and molecular characterization can complement one another [8]. Previous transcriptomic studies of canine mammary tumours (CMTs) have shown that transcriptomic signatures overlap with human breast cancer profiles and dogs are reliable models for oncogenic and pharmacogenomic studies [9].

Clinical relevance of top up-regulated genes in cancer pathogenesis
The top up-regulated genes among malignant and benign tumours in our study were found to have a significant role in cancer progression and tumour invasiveness, as evidenced by human cancer studies. Top five up-regulated genes in malignant tumours were COL11A1, SFRP2, LCN2, COL2A1 and H19. COL11A1 gene encodes one of the two alpha chains of type XI collagen, a minor fibrillar collagen which has been implicated in tumour progression in humans.
Overexpression of COL11A1 has been reported in mesenchymal derived tumours [10][11]. Wu et al., reported that COL11A1 may promote tumour aggressiveness via the TGF-β1 and MMP3 axis in ovarian cancers [12]. Halsted et al., reported expression of COL11A1 in epithelial cells, stroma, and vessels of normal and cancerous breast tissue [13]. Overexpression of COL11A1 has recently been reported in canine melanomas [14] Secreted frizzled-related proteins (SFRPs) constitute a family of 5 members (SFRP 1-5) that modulate Wnt signalling. Overexpression of secreted frizzled related protein 2 (SFRP2) in CMTs has been shown to induce cancerous transformation in normal mammary epithelial cells [15]. The anti-apoptotic function of SFRP2 is mediated through activation of NF-κB or Janus kinase (JNK) suppression [16]. Further in vitro, as well as, in vivo oncogenic potential of SFRP2 have been demonstrated in renal cancer. Recently, SFRP2 was shown to be related to poor prognosis along with genes associated with epithelial-to-mesenchymal transition [17]. Overexpression of SFRP2 was seen in both malignant and benign CMTs in this study. Lipocalin (LCN2), another top up-regulated gene among malignant tumours in our data set, has been shown to be associated with oestrogen receptor (ER)-negative breast tumours in humans [18]. The gene is significantly increased in the luminal epithelial cells compared with myoepithelial cells, an

Impact analysis of major upstream regulators reveals genes involved in tumour angiogenesis and cancer metastasis
In this study, we found numerous dysregulated upstream regulators in both malignant and benign tumours of mammary gland, which affected expression patterns of various other genes related with tumourigenesis. The top upstream regulators in malignant mammary tumours were secreted phosphoprotein (SPP1) and coagulation factor II (F2). Pathway analysis revealed that SPP1 overexpression affected the expression levels of BRCA1 and CDC20, which have been reported to be associated with several cancers types [50-51]. F2 (thrombin) in our study was found to be dysregulated and its expression has been connected to overexpression of various chemokines and their ligands having role in tumour angiogenesis and cancer metastasis [52][53][54]. Nierodzik and Karpatkin have provided ample evidence to support that thrombin activates tumour cell adhesion to platelets, endothelial cells, and subendothelial matrix proteins. Thrombin also enhances tumour cell growth, metastasis and angiogenesis, apart from its role in preservation of dormant tumour cells in individuals, preventing host eradication. Therefore, it is proposed that tumour malignancy may be regulated by a procoagulant/anticoagulant axis [55]. Top upstream regulator identified for benign mammary tumour was AREG. DNA damage signals caused by radiation and chemo are transmitted by master regulators like NF-kB to generate a powerful, conserved and senescence associated secretory phenotype, and its downstream effectors comprise a large spectrum of extracellular proteins including AREG, SFRP2, HGF, IL8, MMPs. Together these proteins give rise to a pro-angiogenic and pro-inflammatory micro-environmental niche that promotes malignant phenotype [56-57]. Furthermore, AREG affected downstream molecules like BIRC5, CCNA2, CCNB2, TOP2A, MMP9, MMP1, CXCR4, have roles in development and progression of various types of cancers [32, 57-61].
Important upstream regulators in our study, unique to benign tumours were AREG, TLR2, TGF1B, HGF, MAP3K1. TLR2 has been reported in intestinal and breast epithelia Table 6. Diseases and bio-functions associated with benign CMTs and their association with the top signalling pathways.

Categories Diseases or Functions Annotation p-Value Predicted Activation State
Activation zscore  [62]. TLR2 agonists fed to neu transgenic mice significantly inhibits breast cancer growth [63] and leads to inhibition of immune responses by production of IL-10 and regulatory T-cells [64]. Thus, TLR2 stimulation on immune cells may also have opposing immune effects as in the case of PSA and PSK [65].

Role of top activated pathways and networks in cancer pathogenesis
The top activated pathways in malignant tumours revealed a unique cancer landscape wherein induction of certain pathways involved targets associated with cell cycle regulation, cellular proliferation, apoptotic pathways, cellular stress and injury (e.g. pathways of cell cycle regulation, oestrogen mediated S phase entry, granzyme B signalling and apoptosis signalling). The data suggests genes involved in cell cycle regulation, apoptosis and cell signalling as major events in the study. The top scoring network in canine mammary tumour was found to have BUB1B as central node in this study. BUB1B is a protein kinase involved in metaphase to anaphase transition checkpoint [66][67]. Aneuploidy and chromosomal instability (CIN) in human cancers can be attributed to alterations in this checkpoint [68]. BUB1B was seen to be linked with NDC80, an essential protein of kinetochore-associated complex required for chromosome segregation and spindle checkpoint activity [69]. NDC80 complex exists as part of a larger super complex called KMN network previously reported to be overexpressed in CMTs [9]. NDC80 has also been shown to depend on KNL-1 and the CENP/H/I/K complex for kinetochore recruitment [70][71][72]. A subnetwork of centromeric proteins CENP-H/C/U/K/A/N was seen within the network. Overexpression of CENP-H is reported in the development and pathogenesis of human breast, colorectal, oesophageal and oral squamous cell cancers [73][74][75][76].
Among benign tumours, central molecule in the top scoring network was vascular endothelial growth factor (VEGF). VEGF has been identified as a vascular permeability factor, angiogenic cytokine, and a survival factor in mammary tumours. VEGF overexpression and its role in angiogenesis of CMTs have been reported [77]. VEGF can act directly on T lymphocytes and elevated VEGF levels may contribute to the aberrant MMP-9 secretion by mammary tumour bearing T cells [78]. VEGF is known to stimulate IL-18 production, which in turn, promotes cancer cell migration, proliferation, angiogenesis and decreases cancer cell susceptibility to lymphocyte mediated cytotoxicity [79][80]. In this network, overexpression of Vasohibin 2 (VASH2) was also seen. VASH2 is a VASH1 homolog, expressed in mononuclear cells and has been reported as an angiogenesis stimulator in mice [81]. Higher expression of VASH2 induces expression of growth factors and promotes cell proliferation in breast cancer [82]. VASH2 is also involved in the proliferation of hepatic and ovarian cancers [83][84]. Another protein in our dataset lying within VEGF network is maternal embryonic leucine-zipper kinase (MELK), a highly conserved serine/threonine kinase that is essential for cell cycle regulation, organogenesis and stem cell proliferation [85]. Interestingly, MELK is additionally involved in the initiation and propagation of numerous human cancers and correlates with poor prognosis [86][87]. Expression of MELK has been reported in canine prostate carcinoma derived cell lines [88]. In this study we also found increased expression of MELK in CMTs. Further in this network, up-regulation of ST8SIA4 was seen to be directly linked to VEGF. ST8SIA4 is involved in sialylation of proteins associated with cancer progression. ST8SIA4 is significantly up-regulated in breast cancers and its up-regulation is highly correlated with cancer malignancy [89][90]. Dual specificity phosphatase (DUSP6), is another protein within the VEGF network which was found to be overexpressed in CMTs in the present study. DUSP6 expression has been reported to play an oncogenic role in breast cancers [91]. Moreover, DUSP6 is part of a high-risk gene signature for non-small cell lung cancer [92], and its expression is significantly correlated with high extracellular signal-regulated kinase (ERK) 1/2 activity in primary human ovarian cancer cells [93]. Other proteins in the VEGF network, whose overexpression was observed, were cyclin-dependent kinase inhibitor-3(CDKN3) and PROM1 (CD133). CDKN3 is a dual-specificity protein tyrosine phosphatase of the CDC14 group, which is often overexpressed in human cancers [94][95]. Since rapidly growing tumours have more mitotic cells, the high level of CDKN3 in mitotic phase provides the best plausible explanation for the frequent CDKN3 overexpression in human cancers [96]. PROM1 is linked to VEGF directly in this network and is attributed to poor prognosis in triple negative breast cancer due to its nuclear mislocalization [97]. Up-regulation of PROM1 increases the invasive capability, metastasis and drug-resistance of breast cancers [98].

Functional significance of differentially up-regulated proteins in CMTs
Proteomic analysis of malignant CMT and healthy mammary tissues, revealed seven proteins highly up-regulated in CMTs, namely ANXA2, APOCII, CDK6, GATC, GDI2, GNAQ and MYH9. Overexpression of ANXA2 and CDK6 has previously been reported in canine gliomas [99] and CMT cell lines [100]. GNAQ is implicated in canine melanomas and T cell lymphomas by different researchers [101][102]. To the best of our knowledge, none of the other overexpressed proteins has previously been reported in cancers affecting dogs. However, role of these proteins in cancer pathogenesis in humans is well studied. Annexin A2 (ANXA2), a calcium-dependent phospholipid-binding protein localized at the extracellular surface of endothelial cells having role in regulation of cellular growth and in signal transduction pathways [103]. ANXA2 over expression has been reported in cancers of the breast, liver, prostate and pancreas, where it plays role in cancer cell migration, invasion, metastasis and adhesion [104][105][106][107]. Apolipoprotein C2 (APOCII) is a lipid binding protein belonging to the apolipoprotein gene family. This protein activates the enzyme lipoprotein lipase, which hydrolyzes triglycerides and thus provides free fatty acids for cells. APOC2 is one of the biomarkers that have been used in differentiating bladder cancer from hernia [108]. Cyclin-dependent kinase 6 (CDK6) plays a vital role in regulating the progression of the cell cycle. More recently, CDK6 has also been shown to have a transcriptional role in tumour angiogenesis, being a part of a transcription complex that induces the expression of the tumour suppressor p16 INK4a and the proangiogenic factor VEGF-A [109]. Emerging evidence suggests that certain tumour cells require CDK6 for proliferation [110]. Guanine nucleotide-binding protein G(q) subunit alpha (GNAQ) is involved in glutamine metabolism and allows the formation of correctly charged Gln-tRNA(Gln). Researchers have demonstrated that glutamine is a major nutrient involved in multiple aspects of cancer metabolism [111]. Most cancer cells are dependent on glutamine and cannot survive in the absence of exogenous glutamine and targeting glutamine metabolism has been recently looked upon as a promising strategy for cancer treatment [112]. Rab GDP dissociation inhibitor beta (RABGDIB/GDI2) is a member of the GDP dissociation inhibitor family that controls the recycling of Rab GTPases involved in membrane trafficking. GDI2 has been proposed as a tumor suppressor gene and acts as an indicator of tumorigenesis in NSCLC [113]. Increased expression of GDI2 has been reported pancreatic carcinoma [114]. Guanine nucleotide-binding protein G(q) subunit alpha(GNAQ) encodes for alpha subunit of q class of heterotrimeric GTP binding protein (Gq) that mediates signals between G-proteincoupled receptors (GPCRs). Genetic, biochemical and biological analysis has shown that GNAQ behaves as a bona fide human oncogene. Recent studies have shown that GNAQ is a pivotal cancer gene in blue naevi [115]. Myosin Heavy Chain 9, (MYH9) encodes NM II-A protein which exists primarily in the cytoplasm and is involved in cytokinesis, cell motility and maintenance of cell shape. The role of MYH9 in cell migration, invasion, and metastasis has been established [116][117]. The overexpression of MYH9 is related to a poor prognosis in oesophageal, bladder, and gastric cancer [116][117][118][119][120].
Thus, the proteomic studies have identified seven proteins unregulated in CMTs. Overexpression of these proteins has not been reported previously in CMTs, however their role in human cancer pathogenesis is well studied. Thus, these proteins need to be studied further for their role in pathogenesis of different cancers of dog.

Conclusions
This study presents gene expression profile of spontaneously occurring canine mammary tumours. Differential gene expression was analysed in context of gene networks and pathways to have insight into the signature pathways associated with benign and malignant CMTs. The research will also add to the knowledge of human cancer studies as dog mammary cancer is a useful model for human breast cancer studies. Presently, the genes accountable for the aggressive behaviour of mammary cancer are not very clear, and poor prognosis associated with malignant mammary cancers emphasize the necessity to unravel the underlying pathways and genes which could act as targets for therapy. In this study we have identified several genes which play a diverse role in tumour angiogenesis, cancer onset and pathogenesis. We have also identified important canonical pathways and key networks in CMTs. The identified differential genes need to be further validated as therapeutic targets and prognostic markers for mammary cancer. In future, further investigations are required to delineate the gene expression patterns among different types & grades of cancers and their response to specific therapies. Further research for identification of gene signatures associated with prediction of response to specific chemotherapies for developing personalized chemotherapeutic regimens is need of the hour.