Ductal Carcinoma In Situ Progression in Dog Model of Breast Cancer

The mechanisms that drive ductal carcinoma in situ (DCIS) progression to invasive cancer are not clear. Studying DCIS progression in humans is challenging and not ethical, thus necessitating the characterization of an animal model that faithfully resembles human disease. We have characterized a canine model of spontaneous mammary DCIS and invasive cancer that shares histologic, molecular, and diagnostic imaging characteristics with DCIS and invasive cancer in women. The purpose of the study was to identify markers and altered signaling pathways that lead to invasive cancer and shed light on early molecular events in breast cancer progression and development. Transcriptomic studies along the continuum of cancer progression in the mammary gland from healthy, through atypical ductal hyperplasia (ADH), DCIS, and invasive carcinoma were performed using the canine model. Gene expression profiles of preinvasive DCIS lesions closely resemble those of invasive carcinoma. However, certain genes, such as SFRP2, FZD2, STK31, and LALBA, were over-expressed in DCIS compared to invasive cancer. The over-representation of myoepithelial markers, epithelial-mesenchymal transition (EMT), canonical Wnt signaling components, and other pathways induced by Wnt family members distinguishes DCIS from invasive. The information gained may help in stratifying DCIS as well as identify actionable targets for primary and tertiary prevention or targeted therapy.


Introduction
There is evidence that all breast cancer subtypes evolve from a non-invasive ductal carcinoma in situ (DCIS) precursor stage [1]. Seldom encountered before mammographic screening became routine, DCIS now accounts for about 20% of breast cancer diagnoses [2]. Although many cases of DCIS regress or remain indolent, some may progress to invasive cancer, which creates a management dilemma, as there is no reliable molecular marker to identify those cases of DCIS destined to progress. Therefore, the current standard of care for all DCIS patients is surgery, radiation, or hormonal therapy. This management approach has good disease-free outcomes, but entails over-treatment in most cases, subjecting women to unnecessary medical and psychological side effects [3,4]. There is an urgent need to understand the natural history of DCIS and to develop methodologies to distinguish non-progressive versus progressive cases.
with the progression of ADH and DCIS to invasive mammary cancer as compared to normal, we evaluated the differential gene expression results derived from the union of DESeq2 [19] and EdgeR ( [20,21]) software at a threshold of False Discovery Rate (FDR) < 0.05. A total of three atypical ductal hyperplasia (ADH), 2449 ductal carcinoma in situ (DCIS), and 2579 invasive cancer (CARC) differentially expressed genes were determined. A total of 527 differentially expressed genes were exclusively present in invasive cancer, while 397 differentially expressed genes were exclusive to DCIS, and 2050 genes were expressed in both DCIS and invasive cancer ( Figure 1A). The number of differentially expressed genes across samples along with results from sample distribution plots suggests that ADH is closer to normal mammary tissue while substantial transcription changes are acquired in DCIS and invasive cancer ( Figure 1B).

Identification of Genes Associated with DCIS Progression to Invasive Cancer
To identify molecular changes associated with DCIS progression to invasive cancer, the list of 2977 differentially expressed genes was filtered to detect 364 (DCIS log2fold-change ±3) and 28 (DCIS log2fold-change ±4) protein-coding genes ( Table 2). These genes included KLRD1, S100A12, EPHA5, SFRP2, FZD2, STK31, and PCDH15. We selected genes involved in the Wnt singling pathway, which is known to contribute to breast cancer progression Frizzled receptors (FZD2), secreted frizzled-related protein 2 (SFRP2), and genes that only expressed during certain physiological changes such spermatogenesis, serine/threonine kinase 31 (STK31), and lactation, alpha-lactalbumin (LALAB). To verify that the selected genes were over-expressed in canine mammary cancer tissues, we performed immunohistochemistry on normal, DCIS, and invasive cancer in another cohort of dogs. Among the selected genes, RNA and proteins of FZD2, SFRP2, and STK31 were significantly over-expressed in DCIS and invasive cancer (Figure 2A,C). FZD2 and its ligands Wnt5a/b (Wnt signaling pathway) were elevated in many cancers including breast cancer, and their expression correlated with epithelial-mesenchymal transition (EMT), which is a process that allows cancer cells to infiltrate surrounding tissues and metastasize to distant sites [22]. EMT is associated with cancer stem cells, which is a subpopulation of cancer cells with self-renewal and differentiation capacity, which play an essential role in tumor initiation, progression, metastasis, recurrence, and therapy resistance. Surface marker CD4 and intracellular marker Aldehyde dehydrogenase 1 family member A1 (ALDHA1) enzyme identify these cells [23]. To confirm that DCIS transition to invasive cancer associated with EMT in canine, we evaluated the expression of ALDHA1 RNA in these lesions. We found that ALDH1A expression increased in DCIS and invasive cancer ( Figure 2D). Another protein that is upregulated in canine DCIS and invasive cancer is LALBA, although it did not pass the threshold of FDR < 0.05 in our RNAseq data ( Figure 2B,C).

Identification of Genes Associated with DCIS Progression to Invasive Cancer
To identify molecular changes associated with DCIS progression to invasive cancer, the list of 2977 differentially expressed genes was filtered to detect 364 (DCIS log2fold-change ±3) and 28 (DCIS log2fold-change ±4) protein-coding genes ( Table 2). These genes included KLRD1, S100A12, EPHA5, SFRP2, FZD2, STK31, and PCDH15. We selected genes involved in the Wnt singling pathway, which is known to contribute to breast cancer progression Frizzled receptors (FZD2), secreted frizzledrelated protein 2 (SFRP2), and genes that only expressed during certain physiological changes such spermatogenesis, serine/threonine kinase 31 (STK31), and lactation, alpha-lactalbumin (LALAB). To verify that the selected genes were over-expressed in canine mammary cancer tissues, we performed immunohistochemistry on normal, DCIS, and invasive cancer in another cohort of dogs. Among the selected genes, RNA and proteins of FZD2, SFRP2, and STK31 were significantly over-expressed in DCIS and invasive cancer (Figure 2A,C). FZD2 and its ligands Wnt5a/b (Wnt signaling pathway) were elevated in many cancers including breast cancer, and their expression correlated with epithelialmesenchymal transition (EMT), which is a process that allows cancer cells to infiltrate surrounding tissues and metastasize to distant sites [22]. EMT is associated with cancer stem cells, which is a subpopulation of cancer cells with self-renewal and differentiation capacity, which play an essential role in tumor initiation, progression, metastasis, recurrence, and therapy resistance. Surface marker CD4 and intracellular marker Aldehyde dehydrogenase 1 family member A1 (ALDHA1) enzyme identify these cells [23]. To confirm that DCIS transition to invasive cancer associated with EMT in canine, we evaluated the expression of ALDHA1 RNA in these lesions. We found that ALDH1A expression increased in DCIS and invasive cancer ( Figure 2D). Another protein that is upregulated in canine DCIS and invasive cancer is LALBA, although it did not pass the threshold of FDR < 0.05 in our RNAseq data ( Figure 2B,C).   Furthermore, we tested the expression of FZD2, SFRP2, STK31, and LALBA proteins in canine mammary (CF41 and P114) and human breast cancer cells lines (HCC1500, T47D, DCIS.com, MDA231, and MCF7) ( Figure 3A and Figure S1A-E). FZD2 protein is expressed in both canine cell lines (CF41 and P114) and in both triple-negative human cell lines (HCC1500 and MDA231). LALBA was over-expressed in DCIS and invasive cancer compared to normal in human tissues ( Figure 3B). We evaluated the significance of FZD2 and STK31 by matching them against genes' expression of human breast cancer from The Cancer Genome Atlas (TCGA) (https://cancergenome.nih.gov). Kaplan-Meier survival analysis of gene expression (n = 1768 patients) for FDZ2, SFRP2 STK31, and LALBA genes revealed that the high expression of these genes correlated significantly with shorter relapse-free survival of patients ( Figure 4 A-D).
Cancers 2020, 12, x 3 of 20   Kaplan-Meier estimator plots were used to visualize the survival probabilities for the low-risk and high-risk groups, as determined using the median risk score for each dataset.

Genes Expression Similarities between Canine and Humans DCIS and Invasive Cancer
Differentially expressed canine genes had significant overlap with similar DCIS gene signature reported by others in humans [24][25][26][27]. We examined the 74 genes identified by Lee et al., who employed microarray technology for the identification of differentially expressed genes in clinical samples (n = 77 divided across five groups) of human DCIS and cancer. The data showed that the 74 genes classified DCIS versus invasive cancer, and of those, 49 of the identified genes overlapped with one or more of the previous studies above. We found that all 74 genes identified by Lee et al. [27] are expressed in our study ( Figure 5A); however, only 13 (FDR 0.05) or 15 (FDR 0.10) genes were Kaplan-Meier estimator plots were used to visualize the survival probabilities for the low-risk and high-risk groups, as determined using the median risk score for each dataset.

Genes Expression Similarities between Canine and Humans DCIS and Invasive Cancer
Differentially expressed canine genes had significant overlap with similar DCIS gene signature reported by others in humans [24][25][26][27]. We examined the 74 genes identified by Lee et al., who employed microarray technology for the identification of differentially expressed genes in clinical samples (n = 77 divided across five groups) of human DCIS and cancer. The data showed that the 74 genes classified DCIS versus invasive cancer, and of those, 49 of the identified genes overlapped with one or more of the previous studies above. We found that all 74 genes identified by Lee et al. [27] are expressed in our study ( Figure 5A); however, only 13 (FDR 0.05) or 15 (FDR 0.10) genes were significantly differentially Cancers 2020, 12, 418 9 of 18 expressed in our study. These 13 genes were also expressed at different levels in DCIS and invasive cancer by other studies [25,26]. Among the over-expressed genes were INHBA, THBS2, MXRA5, MMP11, COL1A2, COL3A1, COMP, COL8A1, FN1, CDH11, while CCL5, ADIPOQ, FNDC1, and GHR are under-expressed. These findings suggest a substantial similarity between canine and humans.
Cancers 2020, 12, x 5 of 20 significantly differentially expressed in our study. These 13 genes were also expressed at different levels in DCIS and invasive cancer by other studies [25], [26]. Among the over-expressed genes were  INHBA, THBS2, MXRA5, MMP11, COL1A2, COL3A1, COMP, COL8A1, FN1, CDH11, while CCL5,  ADIPOQ, FNDC1, and GHR are under-expressed. These findings suggest a substantial similarity between canine and humans. To investigate the immune escape during DCIS progression to invasive cancer, we examined our gene expression data against data generated by Alcazar et al. [28], who investigated potential mechanisms of immune escape in breast cancer and identified top differentially expressed and high variance immune-related genes in human DCIS. Although immune cells were not profiled in this study, we examined essential human immune genes in canine carcinoma, which are depicted as a heatmap across normal, DCIS, and carcinoma samples ( Figure 5B). Similarly, immune-related genes such as HBEGF, SEMA7A, GLI1, IL10RA, JUNB, PSMD8, PTGES2, ILF3, NF-kB, and TNFRSF1B were over-expressed, while DUOX2 and RORA were under-expressed in DCIS and cancer compared to ADH and normal tissues in our study as well as that of Lee et al [27].
We next compared differentially expressed genes between DCIS and invasive cancer to the most known and widely used the prosigna breast cancer progression gene signature (PAM50) assay used for intrinsic subtyping [29]. Forty of the 50 PAM50 genes were included in our dataset, and seven were significantly differentially expressed ( Figure 5C). To further investigate whether the dog model resembled human breast cancer subtypes, both the PAM50 algorithm and unsupervised hierarchical clustering were used to identify the intrinsic subtypes of our dog mammary tumors. PAM50 assigned most of the canine DCIS examined in this study to the ER-, PR-, and Her-2+ subtypes.
Then, we utilized ingenuity pathway analysis (IPA) (Qiagen, CA, USA) and gene set enrichment analysis (GSEA)-hallmark enrichment [30,31] of the differentially expressed genes identified at each stage. Of the IPA-enriched pathways were the Wnt singling pathways (canonical and non-canonical) To investigate the immune escape during DCIS progression to invasive cancer, we examined our gene expression data against data generated by Alcazar et al. [28], who investigated potential mechanisms of immune escape in breast cancer and identified top differentially expressed and high variance immune-related genes in human DCIS. Although immune cells were not profiled in this study, we examined essential human immune genes in canine carcinoma, which are depicted as a heatmap across normal, DCIS, and carcinoma samples ( Figure 5B). Similarly, immune-related genes such as HBEGF, SEMA7A, GLI1, IL10RA, JUNB, PSMD8, PTGES2, ILF3, NF-kB, and TNFRSF1B were over-expressed, while DUOX2 and RORA were under-expressed in DCIS and cancer compared to ADH and normal tissues in our study as well as that of Lee et al. [27].
We next compared differentially expressed genes between DCIS and invasive cancer to the most known and widely used the prosigna breast cancer progression gene signature (PAM50) assay used for intrinsic subtyping [29]. Forty of the 50 PAM50 genes were included in our dataset, and seven were significantly differentially expressed ( Figure 5C). To further investigate whether the dog model resembled human breast cancer subtypes, both the PAM50 algorithm and unsupervised hierarchical clustering were used to identify the intrinsic subtypes of our dog mammary tumors. PAM50 assigned most of the canine DCIS examined in this study to the ER-, PR-, and Her-2+ subtypes.
Then, we utilized ingenuity pathway analysis (IPA) (Qiagen, CA, USA) and gene set enrichment analysis (GSEA)-hallmark enrichment [30,31] of the differentially expressed genes identified at each stage. Of the IPA-enriched pathways were the Wnt singling pathways (canonical and non-canonical) and signaling pathways regulated by the Wnt pathway such as signaling pathways regulating the pluripotency of stem cells and planar cell polarity (PCP), as well as basal cell carcinoma ( Figure 6). Then, we applied GSEA to test whether a member of each gene set occurs randomly toward the top or the bottom of an entire ranked list of differentially expressed genes. We employed a normalized enrichment score (NES), which reflects the degree of over-representation of a gene set at the top or bottom and takes into account the size of the gene set. In addition, we used a nominal p-value that estimates the statistical significance of NES using a gene-set based permutation procedure and FDR, which controls the proportion of false positives to identify up/down-regulated gene [30]. Among the differentially expressed genes along the progression continuum from normal to ADH to DCIS and finally to invasive cancer, more than 40 gene sets were identified with a positive NES at p-value < 0.05 and FDR < 0.05 ( Figure 7A), and 25 pathways were identified with a negative NES at p-value < 0.05 and FDR < 0.05 ( Figure 7B). Among the significantly enriched gene sets were an epithelial-mesenchymal transition, NF-KB, and MTORC1 ( Figure 7C). Furthermore, pathways that include p53, TNFα via NF-κβMTORC1, and MYC targets were over-expressed along the continuum in ADH, DCIS, and cancer compared to normal ( Figure 7D).
Cancers 2020, 12, x 6 of 20 and signaling pathways regulated by the Wnt pathway such as signaling pathways regulating the pluripotency of stem cells and planar cell polarity (PCP), as well as basal cell carcinoma ( Figure 6). Then, we applied GSEA to test whether a member of each gene set occurs randomly toward the top or the bottom of an entire ranked list of differentially expressed genes. We employed a normalized enrichment score (NES), which reflects the degree of over-representation of a gene set at the top or bottom and takes into account the size of the gene set. In addition, we used a nominal p-value that estimates the statistical significance of NES using a gene-set based permutation procedure and FDR, which controls the proportion of false positives to identify up/down-regulated gene [30]. Among the differentially expressed genes along the progression continuum from normal to ADH to DCIS and finally to invasive cancer, more than 40 gene sets were identified with a positive NES at p-value < 0.05 and FDR < 0.05 ( Figure 7A), and 25 pathways were identified with a negative NES at p-value < 0.05 and FDR < 0.05 ( Figure 7B). Among the significantly enriched gene sets were an epithelialmesenchymal transition, NF-KB, and MTORC1 ( Figure 7C). Furthermore, pathways that include p53, TNFα via NF-κβMTORC1, and MYC targets were over-expressed along the continuum in ADH, DCIS, and cancer compared to normal ( Figure 7D).

Discussion
The natural history of DCIS progression remains mostly unknown because of the surgical removal of DCIS lesions upon diagnosis. The currently accepted view of breast cancer pathogenesis assumes a linear progression in which ADH transitions through DCIS to invasive cancer and metastasis. Although epidemiologic, pathologic, and molecular data in humans and animals support the concept of evolution of invasive cancer from DCIS, in the few retrospective studies of DCIS (initially misdiagnosed as ADH) without curative treatment, only 25-50% of DCIS progressed to invasive cancer within 15-25 years [32]. Therefore, the current standard of care leads to the over-treatment of many patients. However, the factors that drive some DCIS to become invasive are not yet understood, so there is an urgent need for diagnostic markers to stratify DCIS lesions and effective targeted therapy to prevent DCIS transition to invasive cancer.
The molecular analysis and identification of markers to stratify DCIS into progressive and indolent cases ideally requires the study of DCIS as it progresses with comparison before and after the onset of invasiveness. This approach was hindered by the limited availability of human tissues from evolving lesions in the same breast, the time needed for progression (more than a decade), ethical concerns, and the lack of an animal model that faithfully represents human disease. Previously, we described the canine mammary gland and cancer as an excellent model for human disease. We have shown that ADH and DCIS in dog mammary glands resemble those of humans in every aspect, including histopathology, clinical outcomes, and imaging characteristics [14,15]. In addition, we have shown that DCIS is prevalent in asymptomatic dogs as in humans [16]. In this study, we performed transcriptomic studies of mammary lesions along the continuum of cancer progression in the same gland (normal to ADH to DCIS to invasive cancer) in the same dog. We described the molecular similarities of DCIS and invasive cancer between humans and dogs, and identified potentially unique targets and perturbed pathways that relate to DCIS progression to invasive cancer in the dog model.
The differentially expressed genes across samples along with results from sample distribution plots suggest that ADH is similar to healthy mammary tissue, but that substantial transcription changes are acquired in DCIS and invasive cancer, which is in agreement with reports that DCIS is molecularly, histopathologically, clinically, and biologically distinct from ADH [33]. In spite of this distinction, common up-regulated pathways were identified between DCIS and ADH, such as p53 and Myc, which were shared between ADH and DCIS, but not invasive cancer, while epithelial-mesenchymal transition (EMT), TNFα/NF-kβ, and MTORC1 signaling pathways were enriched throughout the cancer progression continuum from ADH to DCIS to invasive cancer. These signaling pathways are involved in cancer progression, invasion, and metastasis [34][35][36] and associated with stem cell-like features [37].
The close similarities in gene over-and under-expression between DCIS and invasive cancer have been reported with many studies unable to consistently differentiate DCIS from paired invasive cancer in humans [25,26,38]. No mutations are unique to invasive cancer in comparison to DCIS. Similarly, the gene expression profile showed a highly similar pattern in DCIS and invasive cancer. Therefore, it was suggested that the genetic program necessary for invasive progression might already exist in the preinvasive stages of breast cancer and that a final series of subtle events in DCIS drives the transition to invasive cancer [25,39]. However, Ma et al.'s study using samples from women found significant global alterations in gene expression at the ADH stage that are maintained in the later stages of DCIS and invasive cancer [25]. In another study by Porter et al., no genes were specific for invasive cancer or DCIS, and the most dramatic and consistent phenotypic changes occurred at the normal-to-in situ carcinoma transition in women [40]. However, other studies that compared the transcriptome of DCIS and invasive cancer using more robust techniques have identified stage-specific markers and a gene expression classifier that differentiate DCIS and invasive cancer [26,41]. This discordance may be because none of the described studies compared DCIS samples before and after invasiveness or with DCIS that did not give rise to invasive disease. Furthermore, when we compared our data to published data by Lee et al. [27] and Knudsen et al. [24], only 13 genes overlapped with their data, which may be due to the reasons mentioned above, the sample size differences between the two studies, and the fact that the comparison of gene expression data obtained from different array platforms results in a marginal overlap [42] and might not suggest species differences.
Recently, growing evidence suggests that the breast microenvironment, including myoepithelial cells, stroma, and extracellular matrix (ECM), plays a crucial role in the transition from DCIS to invasive cancer [43,44]. Tumor-associated myoepithelial cells have a tumor-protective effect [43]. The myoepithelial layer and basement membrane that surrounds DCIS is thought to be lost during progression, allowing tumor cells to invade the stroma and surrounding adipose tissue. We found myoepithelial cells markers, including MMP11, COL1A2, COL3A1, COL8A1, S100A2, and FN1, were over-expressed in DCIS compared to invasive cancer (Table 2) [45,46]. The expression of S100A2 is induced by transforming growth factor-β (TGF-β) in DICS ( Figure 7A), which is a potent inducer of EMT ( Figure 7C) and cancer progression. These findings support a model wherein alterations in myoepithelial cells promote the progression of DCIS to invasive cancer via TGFβ signaling activation. Activation of the TGFβsignaling pathway promoted the EMT, basal-like phenotypes, stemness, and invasiveness of DCIS cells.
Despite the close similarities in gene expression between DCIS and invasive cancer, we identified genes that differed significantly in expression between DCIS and invasive cancer. In our study, the differentially expressed genes in DCIS compared to invasive cancer are FZD2 and SFRP2, which are the core WNT pathway components, STK31, and LALBA. FZD2 is elevated in many cancers, including breast cancer, and its expression correlated with markers of EMT and cell migration through a non-canonical pathway that predicts metastasis and overall survival in patients [22]. Similarly, SFRP2 is shown to be associated with poor prognosis concomitant with the expression of genes associated with EMT [47]. In agreement with the above findings, the pathway analysis of the gene expression data in this study showed enrichment of WNT (both canonical and non-canonical) pathways, the planar cell polarity (PCP) pathway, and signaling pathways regulating the pluripotency of stem cells in DCIS (generally induced by WNT signaling). The Wnt signaling pathway is known to promote cancer progression, metastasis, and is involved in EMT [48]. EMT is essential in tumor invasion because it reduces cell-cell adhesion and increases cell motility. The non-canonical PCP pathway as mentioned above exerts an essential role in cell differentiation by regulating critical components of the cytoskeleton that lead to cell shape and motility changes [49]. These processes are primary events in the multi-step progression from DCIS to invasive cancer.
On the other hand, cancer STK31, a cancer-testis antigen (CTA), plays crucial roles in human cancer through regulation of the cell cycle, and its over-expression increases cell migration and invasiveness, whereas its depletion induces apoptosis [50]. This molecule is expressed in the human germline, but it is not expressed in most adult tissues, where it is restricted to the testis and various malignancies. CTA expression is highly increased in DCIS of the breast, particularly those that are ER-negative [51]. As a result of their restricted expression pattern, these proteins are immunogenic in cancer patients and are considered essential targets of cancer vaccines.
LALBA expression was high in DCIS and invasive cancer, but it did not pass our FDR threshold of <0.05. Our interest in α-lactalbumin stems from the fact that it is expressed in mammary epithelial cells only during lactation [52] but expressed at high levels in the majority of human breast carcinomas [53]. α-Lactalbumin as a vaccine target mediated protection against the development of murine breast cancer in the absence of any detectable autoimmune-induced breast inflammation.
Human breast cancer consists of four subtypes: luminal A and B, triple-negative, and Her-2 positive. Previously, we reported that one-third of canine DCIS and tumors are immunohistochemically negative for hormonal receptors. In this study, PAM50 assigned most of the tumors examined in this study to the HER-2 breast cancer subtype, which is the subtype most commonly associated with ER-negative and HER-2 positive subtype, corroborating the similarities between our dog model and human breast cancer.

Clinical Specimens
The Purdue animal care and use committee approved all animal experiments. For this pilot study, healthy, ADH, DCIS, and invasive cancer samples were collected from the same progressing lesions in mammary glands from the same dogs (n = 3). Tissue blocks were reviewed by a board-certified veterinary pathologist to confirm the diagnosis and define lesions for dissection. Immunohistochemistry with antibodies to estrogen receptor (ER)-α, progesterone receptor (PR), and human epidermal growth factor-2 (HER-2) was used to determine hormonal expression status.

RNA Isolation
Five serial sections of 5-µm thickness were stained; then, normal or abnormal cells were microdissected from each lesion. Total RNA was extracted from each lesion in 300 µL of lysis buffer with 50 µg of proteinase K using MasterPure Complete DNA and RNA Purification Kit (Epicentre Biotechnologies, Madison, WI, USA) according to manufacturer instructions. Genomic DNA contamination was removed by treatment with DNase I (10 U/mL; Ambion, Inc, Austin TX, USA 78744). Then, the RNA was dissolved in RNase-free water and quantified using NanoDrop 1000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) at 260 nm.

Next-Generation Sequencing
Individual RNA samples were submitted to Otogenetics (Norcross, GA, USA) for RNA-Seq analysis. Briefly, the integrity and purity of the RNA samples were assessed using the Agilent Bioanalyzer and OD260/280 (Agilent, Santa Clara, CA, USA). rRNA was depleted from each sample using the Ribo-Zero rRNA Removal Kit (Epicentre, Madison, WI) according to the manufacturer's instructions. cDNA was generated using the Smart PCR cDNA kit (Clontech Laboratories, Mountain View, CA, USA) from 100 ng of total RNA, and adaptors were removed by digestion with Rsa I. This method used a low cycle number PCR to preferentially amplify poly (A+) RNA via a modified oligo (dT) primer. The resultant cDNA was fragmented via sonication (Covaris, Inc., Woburn, MA, USA), profiled via Agilent Bioanalyzer, and subjected to Illumina library preparation using NEBNext reagents (New England Biolabs, Ipswich, MA, USA). The quality, quantity, and size distribution of the Illumina libraries were determined using an Agilent Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA). Then, the libraries were submitted for Illumina HiSeq2000 sequencing as per the manufacturer's recommendations. Paired-end 90 or 100 nucleotide reads were generated, checked for data quality using FASTQC (Babraham Institute, Cambridge, UK), and subjected to data analysis using the DNAnexus platform (DNAnexus, Inc, Mountain View, CA, USA). Transcript-level quantitation was per the DNAnexus White Paper (version 1.1, April 19, 2010, DNAnexus, Mountain View, CA, USA) on RNA-Seq/3SEQ Transcriptome.

Bioinformatics Analysis
Sequence data quality was determined using FastQC (version 0.11.2, Babraham Institute, Cambridge, UK). Quality trimming and filtering were performed using TrimGalore (version 0.4.4). Illumina adapter sequences, any bases with Phred score less than 30, and reads shorter than 50 base pairs were removed during quality control. Quality trimmed reads were mapped against the Canis familiaris reference genome (CanFam3.1) using STAR aligner [54]. Read counts for each gene feature each replicate were determined using the HTSeq tool [55] and the combined count matrix for all replicate was generated using a custom Perl script. For Differential Expression (DE) analysis, genes with zero counts across all replicates were discarded. DE analysis of cancer and healthy samples (ADH versus Normal, DCIS versus Normal, and invasive carcinoma (CARC) versus Normal) was carried out using 'R-Bioconductor' package DESeq2 ( [19]) and EdgeR ( [20]). Since the sequencing of replicates was performed in batches, a batch correction was applied by providing a custom experiment design matrix in DESeq2 and EdgeR analysis. Genes with a False Discovery Rate (FDR) less than 0.05 (i.e., 95% confidence interval) were denoted as DE genes. The Venn diagram denotes a subset of genes that are common across two methods (DESeq2 and EdgeR) and cancer stages (ADH, DCIS, and CARC). A list of genes that are modulated in DCIS and CARC conditions in human models were obtained from the literature survey, and their expression profile in current data is visualized as Heatmaps. Heatmaps are generated by plotting the log2FoldChnage using R package heatmap.

Pathway Enrichment Analysis
Pathway analysis was performed using the Ingenuity pathway Analysis (IPA) and gene set enrichment analysis (GSEA). Enriched pathways and gene sets (p-value < 0.05) associated with cancer stages (ADH, DCIS, and CARC) were determined.

Quantitative Real-Time RT-PCR (qRT-PCR)
Total RNA was isolated from fresh-frozen healthy, benign (ADH, DCIS), and invasive cancer using Trizol reagent (Invitrogen Life Technologies, Carlsbad, CA, USA) according to manufacturer instructions. The RNA concentration was determined using the NanoDrop 1000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Reverse transcription was performed using Transcriptor First Strand cDNA synthesis kit (Roche Diagnostics, Indianapolis, IN, USA) from 1 µg of total RNA and the combination of anchored-oligo (dT) and random hexamer primers according to the manufacturer's instructions. Gene expression analysis was performed with qRT-PCR using the LightCycler 480 SYBR Green I Master (Roche Diagnostics, Indianapolis, IN, USA) and primers. All reactions were performed in triplicate, and the relative expression of target mRNA in each sample was normalized with that of mean GAPDH abundance. The primers were purchased from Integrated DNA Technologies, Inc. (Coralville, IA, USA).

Immunohistochemistry
Paraffin sections, 5-µm thick, were mounted on positively charged Super-frost slides (Fisher Scientific, Chicago, IL, USA), which were then deparaffinized and rehydrated in graded alcohols and subjected to immunohistochemistry according to the manufacturer's protocol (Biocare Medical, Concord, CA, USA). In brief, the hydrated sections were heated in a Biocare Gene retrieval steamer 45 min in 10 mM of citrate buffer, pH 6.0, to unmask the epitopes. Then, 3% hydrogen peroxide was used to block the endogenous peroxidase activity. Nonspecific binding was blocked with background polisher, and slides were incubated for 10 min at room temperature. Then, the slides were incubated for 30 min with primary mouse monoclonal antibody or rabbit polyclonal antibodies as described. Then, the primary antibody complexes were visualized with a biotin-labeled secondary antibody (Universal Goat Link) for 15 min; and finally, with horseradish peroxidase (HRP) for 15 min, reacted for 5 min with 3, 3'-diaminobenzidine tetrahydrochloride, and counterstained with hematoxylin, dehydrated, and mounted. We treated negative control slides with no antibody.

Western Blot Analysis
Proteins were prepared by using Radioimmunoprecipitation assay buffer (RIPA) (0.5 m Tris/HCl, pH 7.4, 1.5 m NaCl, 2.5% deoxycholic acid, 10% NP-40, 10 mm EDTA). Protein concentrations were determined using a Pierce BCA Protein Assay Kit according to the manufacturer's instructions (Thermo Fisher Scientific). Ten micrograms of each sample protein was subjected to SDS/PAGE and transferred to nitrocellulose paper. The blots were reacted sequentially with primary antibodies: HRP-conjugated goat anti-rabbit immunoglobulin (IgG) or goat anti-mouse IgG, and visualized with diaminobenzidine.

Kaplan Meier Srvivial Analysis
To determine the significance of FDZ2, SFRP2, STK31, and LALBA genes that we identified in canine DCIS progression in women with breast cancer survival, we performed Kaplan-Meier survival analysis using TCGA data and the publically available Kaplan-Meier plotter. The cut-off for the analysis was set at the auto select best cut-off. It is 172, 8383, and 9 for FZD2, SFRP2, and STK31, respectively.

Statistical Analysis
All data are presented as means ± standard deviation (SD). Statistical calculations were performed with Microsoft Excel analysis tools. Differences between individual groups were analyzed by paired t test. p-values of < 0.05 were considered statistically significant.

Conclusions
We report gene expression profiles along the cancer progression continuum from ADH to DCIS to invasive cancer in the canine model of breast cancer. In previous studies, we showed that canine DCIS and invasive cancer resemble human DCIS and its invasive progression to histopathology, the expression of many tumor markers, and clinical outcomes, as well as imaging characteristics. The results of this study suggest that alteration of the Wnt signaling pathway may play a role in the progression of DCIS to invasive cancer, similar to data reported in humans. In addition, we have identified gene markers, such as STK31 and α-lactalbumin, that may serve as vaccine targets for the prevention of DCIS progression. With its morphological, diagnostic imaging, and transcriptomic similarities to mammary DCIS and invasive carcinoma in women, the canine model could have a transformative impact on therapeutic and preventative strategies for breast cancer.