Malignant canine mammary epithelial cells shed exosomes containing differentially expressed microRNA that regulate oncogenic networks

Background Breast (mammary) cancers in human (BC) and canine (CMT) patients share clinical, pathological, and molecular similarities that suggest dogs may be a useful translational model. Many cancers, including BC, shed exosomes that contain microRNAs (miRs) into the microenvironment and circulation, and these may represent biomarkers of metastasis and tumor phenotype. Methods Three normal canine mammary epithelial cell (CMEC) cultures and 5 CMT cell lines were grown in serum-free media. Exosomes were isolated from culture media by ultracentrifugation then profiled by transmission electron microscopy, dynamic light scattering, and Western blot. Exosomal small RNA was deep-sequenced on an Illumina HiSeq2500 sequencer and validated by qRT-PCR. In silico bioinformatic analysis was carried out to determine microRNA gene and pathway targets. Results CMEC and CMT cell lines shed round, “cup-shaped” exosomes approximately 150–200 nm, and were immunopositive for exosomal marker CD9. Deep-sequencing averaged ~ 15 million reads/sample. Three hundred thirty-eight unique miRs were detected, with 145 having > ±1.5-fold difference between one or more CMT and CMEC samples. Gene ontology analysis revealed that the upregulated miRs in this exosomal population regulate a number of relevant oncogenic networks. Several miRNAs including miR-18a, miR-19a and miR-181a were predicted in silico to target the canine estrogen receptor (ESR1α). Conclusions CMEC and CMT cells shed exosomes in vitro that contain differentially expressed miRs. CMT exosomal RNA expresses a limited number of miRs that are up-regulated relative to CMEC, and these are predicted to target biologically relevant hormone receptors and oncogenic pathways. These results may inform future studies of circulating exosomes and the utility of miRs as biomarkers of breast cancer in women and dogs. Electronic supplementary material The online version of this article (10.1186/s12885-018-4750-6) contains supplementary material, which is available to authorized users.


Background
Canine mammary tumors (CMT) represent the most frequent tumor in hormonally intact (non-ovariohysterectomizd) female dogs, and CMT have similar incidence and comparable distribution of malignant potential to breast cancer (BC) in women [1][2][3]. Some forms of CMT may represent a useful translational model for human BC as it shares risk factors such as age, hormone exposure, and obesity, although there are some differences in subtypes that limit direct comparison [4][5][6]. Furthermore, CMT and BC share similar genetic alterations, including downregulation of tumor suppressors p16/INK4A, PTEN, BRCA1, and p53, as well as upregulation of oncogenes KRAS, PI3K/AKT, and MAPK [7][8][9]. Recently, CMT has been shown to be classifiable into human molecular subtypes luminal A, luminal B, HER-2, and triple-negative/basal-like according to estrogen receptor alpha (ESR1α), progesterone receptor, and HER-2/ErbB-2 expression by immunohistochemistry in patient tumor tissue and by qRT-PCR in a cohort of well-characterized cell lines [10,11]. In addition, CMT, like human BC, shows a negative correlation between estrogen hormone receptor ESR1 expression and increasing tumor grade [10,12].
Currently, as in women with BC, definitive classification of benign versus malignant CMT, as well as tumor grading, requires histopathology. This is problematic because collecting those samples requires invasive surgery. Current less invasive alternatives such as fine-needle aspirate cytology vary from 67.5-81% accuracy [13,14]. Another significant prognostic factor for CMT is advanced stage, with shortened survival times for dogs with large tumors (> 3 cm) and/or metastasis, highlighting the importance of early detection [15,16]. An accurate, minimally invasive, biomarker for CMT diagnosis and malignant potential could improve outcomes through intervention at a lower stage of disease.
One such potential class of biomarker is microRNAs (miRNAs, miRs), a type of small (18-22 nucleotides), non-coding RNA that are highly conserved across species and play crucial roles in the negative, post-transcriptional regulation of gene expression in both health and disease [17]. Each miRNA recognizes numerous gene targets through hybridization with a complementary "seed sequence" in the 3′ untranslated region (UTR) of mRNA resulting in either degradation of the transcript or inhibition of ribosomal translation [18]. Dysregulation of miR-NAs is particularly prevalent in cancer, where genetic instability of tumors leads to altered miRNA expression profiles that promote oncogenesis [19]. Numerous studies demonstrate miRNA are differentially expressed in women with BC in tissue, exosomes, and serum/plasma, [20][21][22][23].
Multiple miRNAs are already known to be altered in CMT, though the data are complex and sometimes conflicting, particularly depending on the RNA source (cells, exosomes, tumor tissues, serum/plasma, etc), microRNA profiling technique(s), and normalizing strategies. One early study of miRNA expression between CMT and normal mammary gland tissue using a small set of individual qRT-PCR assays found miR-29b and miR-21 were significantly upregulated in neoplastic versus normal tissue, while miR-181b and let-7f were specifically upregulated in tubulopapillary carcinoma [24]. In another, miR-141 specifically was demonstrated to be upregulated in several well-characterized CMT cell-lines and experimentally validated to down-regulate tumor suppressor p16/INK4a [9]. This same study also identified a number of other differentially expressed miRs by qRT-PCR microarray in CMT that were also altered in human BC, including miR-21, miR-155, miR-9, miR-34a, miR-143/ 145, and miR-31 [9]. A separate research group established a CMT line (labeled "SNP") from a primary patient mammary tumor and compared its miRNA expression to normal mammary tissue through miRNA hybridization arrays and found the microRNAs with the greatest increase and decrease were miR-143 and miR-138a, respectively [25]. A separate study of various primary and metastatic canine mammary tumor tissues using qRT-PCR found up-regulated miR-210 in neoplastic versus normal tissue, higher miR-21 in malignant mammary carcinomas (but not benign tumors), and that the metastatic tumors had altered miR-29b, miR-101, mir-125a, miR-143 [26]. Another study evaluated primary versus metastatic mammary carcinomas using RNA hybridization arrays and qRT-PCR, and found a distinct signature of microRNA expression in metastatic canine mammary carcinoma, although the expression of these candidate metastasis markers were not statistically different in peripheral blood [27]. Finally, a recent study evaluating circulating microRNA in blood by qRT-PCR for multiple different types of cancer found miR-214 and miR-126 were significantly up-regulated in serum from dogs with mammary carcinoma (along with numerous other tumor types) [28]. Of note, most of these studies performed no or limited in silico bioinformatics analysis for these miRNA, and only the study of miR-141 and p16/ INK4a experimentally validated the annotated targets [9]. miRNAs make particularly good biomarkers because they can be secreted in biofluids such as serum, urine and breast milk, and protected from endogenous RNases by packaging in exosomes and/or binding to proteins such as Argonaute [17,18]. Exosomes are 30-200 nm in diameter round vesicles with a lipid membrane, and are secreted by cellular organelles called multivesicular bodies. There is some evidence to indicate that exosomes are actively secreted by tumor cells to facilitate cell-to-cell communication to distant cells and tissues [29]. These tumor exosomes and their cargo of miRs, mRNA, and proteins may also modulate the behavior of local stromal and immune cells [30]. One such study provided data that tumor exosomes derived from human patients with lung and pancreatic carcinomas were able to induce myotube apoptosis through miR-21 and TLR7 signaling, recapitulating the cancer cachexia phenotype in an in vitro model [31].
The aims of the current study were to isolate and characterize exosomes shed by normal canine mammary epithelial cells (CMEC) and CMT cells in vitro, analyze the miRNA profile of these exosomes, and to perform in silico bioinformatic annotation of this miRNA signature. We hypothesized that both CMEC and CMT cells grown in serum-free media would shed exosome-like microvesicles containing abundant miRNAs, and that the miRNA signature of the CMT extracellular vesicles would be enriched in a subpopulation of miRs predicted to regulate important molecular targets in CMT.

Cell culture
The following cell lines were used: Three normal canine mammary epithelial cell cultures independently derived from separate canine patients without mammary pathology (CMEC1, CMEC2, CMEC3), and five stable and highly transformed cell lines derived from canine patients with histopathology-confirmed mammary carcinoma including CMT12 (formerly CMT2), CMT27 (formerly CMT4) and CMT28 (formerly CMT5) as well as 2 more recently derived lines including CMT47 (derived from a mammary adenocarcinoma from a pure-bred Miniature Schnauzer) and CMT119 (derived from a mammary carcinoma from a Golden Retriever) [9,11]. The CMT cells used are the product of our laboratory group in collaboration with Dr. Lauren Wolfe (retired). They were all derived/rescued from surplus biopsy specimens recovered following standard of care surgery of canine mammary cancer patients. Each biopsy specimen to be cultured was divided in two at the time of collection: one for epithelial cells to be sorted by flow cytometry and grown in culture, and one placed in formalin to be processed for routine histopathology and reviewed by a board-certified pathologist to identify the cell type and confirm malignancy. All CMT cell lines in this study were confirmed to be derived from mammary carcinoma/adenocarcinoma tumors on histopathology, but tissues were not further classified into mammary tumor histologic subtypes (i.e. simple, complex, micropapillary, etc). All owners of such animals sign a general informed consent that notes that biopsy specimens recovered in this manner may be used for research. No IACUC approval is required for such specimens. All CMEC samples were recovered from normal geriatric female dogs in the breeding colony housed and managed under IACUC PRN 2015-2688. All such specimens were recovered post-mortem following euthanasia performed in the normal management of the colony. In no case were any of these animals euthanized for the current study. All cell lines are routinely analyzed using caninespecific RT-PCR assays for Canine Mammaglobin-A (unpublished data) to ensure the species source.
CMT and CMEC cells were grown in 75 cm 2 flasks in synthetic Xerum-free® media + DMEM media supplemented with 2X penicillin/streptomycin antibiotics at 37°C until 70-80% confluence. Media from the first 24 h of culture was discarded and conditioned media from the second 48 h was harvested on day 3 of growth prior to trypsinization and subculture.

Exosome isolation
Exosomes and exosomal proteins were isolated by progressive centrifugation and ultracentrifugation. Briefly, 5-10 mL conditioned cell culture media were progressively centrifuged at 4°C at 300×g for 10 min, 2000 x g for 10 min, and 10,000 x g for 30 min, each time discarding the pellet and retaining the supernatant, to remove cells and debris. The processed supernatants were then centrifuged at 100,000 x g at 4°C for 70 min, and the resulting supernatant was subjected to another cycle of centrifugation at 100,000 x g at 4°C for 70 min. The final pellet was resuspended in 50 uL PBS [32].

Dynamic light scattering
The size distribution of vesicles in the cell-free conditioned media diluted 1:5 to 1:20 (depending on particle count rate) in 1× DEPC-treated PBS was measured by intensityweighted dynamic light scattering using a Malvern ZetaSizer ZS90 (Malvern instruments, Ltd., Worcestershire, UK) according to the manufacturer's instructions.

Transmission electron microscopy
Cell-free conditioned media from both CMEC and CMT cells, prepared by progressive centrifugation and ultracentrifugation as previously described, was loaded onto copper-formvar grids treated with 1% Alcian blue (to increase hydrophobicity) and negatively stained with 1% uranyl acetate [32]. Grids were loaded into a Zeiss EM10 transmission electron microscope (Carl Zeiss Microscopy, LLC, Thornwood, NY, USA) and imaged at 20,000× to 63,000× magnification with an accelerating voltage of 60 kV (2 s exposure).

Western blot
Protein from ultracentrifuge-precipitated exosomes was quantified by nanospectrophotometry and/or Qubit protein assay. Two micrograms of native exosomal protein from pooled CMEC and CMT samples was mixed with 4× Laemmli buffer heated at 95°C for 15 min. Proteins were resolved by SDS-PAGE on 4-20% precast polyacrylamide gels (Bio-Rad, Hercules, CA, USA) using the Precision Plus Protein Western C standards (Bio-Rad, Hercules, CA, USA) to determine the sizes of the bands, and then transferred to a nitrocellulose blotting membrane (LI-COR Biosciences, Lincoln, NE, USA). After electrophoresis, the fractions were electro-transferred to nitrocellulose membrane (Bio-Rad, Hercules, CA, USA) and blocked for 1 h. with Odyssey Blocking Buffer (LI-COR Biosciences, Lincoln, NE, USA). Membranes were incubated overnight at room temperature with 1:200 primary antibody CD9 Mouse-anti-Human (clone MM2/57, Bio-Rad AbD Serotec Inc., Hercules, CA, USA) in Odyssey Blocking Buffer (LI-COR Biosciences, Lincoln, NE, USA). Next the membranes were washed 3X for 10 min in 1X PBS in 0.1% Tween 20 (Sigma-Aldrich Corp., St. Louis, MO, USA). Next, secondary antibody IRDye Goat-anti-Mouse (LI-COR Biosciences, Lincoln, NE, USA) 1:10,000 dilution was incubated in Blocking Buffer (LI-COR Biosciences, Lincoln, NE, USA) for 1 h at 4°C. Membranes were washed 3 times at room temperature with 1X PBS in 0.1% Tween 20 for 10 min. Fluorescent bands were visualized with Odyssey Near-Infrared Western Blot detection system in Image Studio (LI-COR Biosciences, Lincoln, NE, USA).

RNA extraction & microRNA deep-sequencing
RNA was extracted from 5 mL of cell-free, serum-free conditioned media using the Norgen Biotek Urine/Cell Culture Exosomal RNA Isolation kit (Norgen Biotek, Thorold, ON, Canada) according to manufacturer instructions. After the lysis step, 10 pM final concentration synthetic miRVana cel-miR-39-3p mimic (Thermo Fisher Scientific, Waltham, MA, USA) was spiked into samples as an external control for technical variation. Extracted RNA was stored at − 80°C until being shipped on dry ice to the Genomic Services Laboratory at the Hudson Alpha Institute for Biotechnology. Small RNA libraries were prepared for sequencing from total RNA from each sample using a NEBNext Small RNA Library Prep Set for Illumina (New England BioLabs Inc., Ipswich, MA, USA) according to the manufacturer's protocol. Briefly, 3′ adapters were ligated to total input RNA followed by hybridization of multiplex SR RT primers and ligation of multiplex 5' SR adapters. Reverse transcription (RT) was performed using ProtoScript II RT for 1 h at 50°C. Immediately after the RT reaction, PCR amplification was performed for 15 cycles using LongAmp Taq 2X Master Mix. Illumina indexed primers were added to uniquely barcode each sample. Post-PCR material was purified using a QIAquick PCR purification kit (Qiagen Inc., Valencia, CA, USA). Post-PCR yield and concentration of the prepared libraries were assessed using a Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, California, USA) and DNA 1000 chip on an Agilent 2100 Bioanalyzer (Applied Biosystems, Carlsbad, CA, USA), respectively. Size selection of small RNA was done using 3% dye free agarose gel cassettes on a Pippin Prep instrument (Sage Science Inc., Beverly, MA, USA). Post-size selection yield and concentration of the libraries were assessed using Qubit 2.0 Fluorometer and DNA High sensitivity chip on Agilent 2100 Bioanalyzer, respectively. Accurate quantification for sequencing applications was performed using the qPCR-based KAPA Biosystems Library Quantification kit (Kapa Biosystems, Inc., Woburn, MA, USA). Each library was diluted to a final concentration of 1.25 nM and pooled in equimolar ratios prior to clustering. Single End (SE) sequencing (50 bp) was performed to generate at least 15 million reads per sample on an Illumina HiSeq2500 sequencer (Illumina, Inc., San Diego, CA, USA).
Post processing of the sequencing reads from RNA-seq experiments from each sample was performed as per the Genomic Services Laboratory unique in-house pipeline. Briefly, quality control checks on raw sequence data from each sample was performed using FastQC (Babraham Bioinformatics, London, UK). Raw reads were imported on a commercial data analysis platform (AvadisNGS, Strand Scientifics, CA, USA). Adapter trimming was done to remove ligated adapter from 3′ ends of the sequenced reads with only one mismatch allowed, poorly aligned 3′ ends were also trimmed. Sequences shorter than 15 nucleotides length were excluded from further analysis. Trimmed reads with low qualities (base quality score less than 30, alignment score less than 95, mapping quality less than 40) were also removed. Filtered reads were then used to extract and count the small RNAs which were annotated using microRNAs from the miRBase release 20 database (http://www.mirbase.org/). Samples were subjected to quantification and active region quantification (Ava-disNGS, Strand Scientifics, CA, USA). The quantification operation carries out measurement at both the gene level and at the active region level. Active region quantification considers only reads whose 5′ end matches the 5′ end of the mature miRNA annotation. Samples were then grouped by identifiers and the differential expression of each miRNA was calculated based on the fold change observed between different groups.
qRT-PCR and data analysis microRNA deep-sequencing results were validated by stem-loop quantitative RT-PCR (qPCR) for selected miRNA targets (selection process discussed in detail below). cDNA was created for each miRNA with a unique TaqMan™ stem-loop primer (Thermo Fisher Scientific, Waltham, MA, USA) with 1 ng RNA input using the TaqMan™ MicroRNA Reverse Transcription Kit (Thermo Fisher Scientific, Waltham, MA, USA) according to manufacturer instructions. A 1 μL cDNA product from the RT reaction was used as input for the qPCR reaction with TaqMan Universal Master Mix II no UNG (Thermo Fisher Scientific, Waltham, MA, USA), a specific 20X TaqMan microRNA assay for each target (Thermo Fisher Scientific, Waltham, MA, USA) and run in a BioRad CFX96 thermocycler (Bio-Rad, Hercules, CA, USA) according to the manufacturer protocol. Internal reference target miR-16 and external spike-in cel-miR-39 were used as control genes for normalization. Cq data was normalized using the 2-ΔΔCq Livak method and presented as both log10 relative quantity for individual samples and fold change for the CMT group relative to the CMEC group [33]. miRNA that did not amplify were assigned a Cq number of 40 to allow calculation of normalized Cq values. These qPCR validation experiments were performed in triplicate and results were averaged and presented as mean ± SD.

Gene target predictions
Target prediction for each miRNA was accomplished using the miRDB online resource and analysis platform (http://www.mirdb.org//). This tool was created in 2008 and was comprehensively updated recently when the complete set of miRNA sequences from the miR Base repository were downloaded along with the complete set of 3'UTR sequences contained in the NCBI RefSeq database. Furthermore, the miRDB target prediction algorithm, MirTarget, which was developed using support vector analysis of high throughput expression data, is capable of predicting conserved and non-conserved target genes via weighting target site conservation as a high priority, but not as an absolute requirement. miRDB scores predicted targets in a range from 50 to 100, with a higher score indicating a greater statistical confidence in the prediction. According to the FAQ on the miRDB website, "a predicted target with a score > 80 is most likely to be real." Subsequently, target gene prediction was performed and scores greater than 80 were considered as representing the most confident gene predictions [34]. Gene targets for the complete set of 16 miRNAs were generated by selecting all gene targets having scores greater than 80 for each of the 16 miRNAs. The resulting set of all gene targets was filtered to remove redundant gene targets (i.e. gene targets that were associated with two or more different miRNAs).

Gene ontology and KEGG pathway enrichment
The DAVID database for annotation, visualization and integrated discovery (version 6.8) was used to perform gene ontology enrichment analysis on sets of target genes using the gene symbol produced by the target prediction algorithm in miRDB. Canine gene symbols were uploaded into the DAVID database and the resulting sets of enriched gene ontology terms or KEGG pathways were identified [35].

Statistical analysis
qRT-PCR relative quantity data were assessed for normality by visual inspection and Shapiro-Wilk test. Non-directional, non-parametric Mann-Whitney statistical testing was performed based on data that were not normally distributed. p < 0.05 was considered statistically significant for both RNAseq and qRT-PCR expression comparisons between groups. The DAVID gene ontology software provides both raw and Benjamini-corrected p-values, and a crude threshold of p < 0.06 was selected to screen for potentially relevant pathways [36].

Characterization of extracellular vesicles
Ultrastructural evaluation of exosome-enriched supernatants using transmission electron microscopy confirmed the presence of variable numbers of irregularly round, occasionally cup-shaped vesicles ranging in size from approximately 60-120 nm in diameter (Fig. 1a). These vesicles occasionally formed variably-dense accumulations. Dynamic light scattering of this cell-free fraction showed a broad distribution of particle sizes with an average diameter of approximately 150 to 200 nm (Fig. 1b). The average protein concentrations in these cell-free fractions was 0.13 to 0.6 μg/μL. Putative exosome marker CD9 was detected in cell-free media from both CMEC and CMT cell lines by Western blotting (Fig. 1c). Our findings were consistent with having isolated exosomelike extracellular vesicles.
microRNA profiling by small RNA deep-sequencing and qRT-PCR validation The RNA bioanalyzer profiles were typical of exosomal samples, skewing heavily towards small RNA (~20-200 nucleotides), with minimal detection of rRNA (Fig. 2a). The normalized miRNA in reads per million for all samples are provided in Additional file 1: Table S1. Three hundred thirty-eight unique miRs were detected in the cell-free RNA fractions from CMEC and CMT samples. In a principal component analysis of the miRNA from all eight samples the CMEC and CMT samples clustered into two separate groups, although the CMT group had two significant outliers (Fig. 2b).
Volcano plot analysis illustrates that numerous miRs were significantly over-and under-expressed by CMT exosomes relative to CMEC (Fig. 2c). Using criteria of p < 0.05 and a fold-change ≥ ± 1.5-fold change, there were 170 differentially expressed miRs between CMT and CMEC groups (Table 1). Removing isoform miRs from different chromosomal locations yielded 145 unique differentially expressed miRs, with 118 miRs upregulated and 27 miRs downregulated in CMT as compared to CMEC.
Three of the significantly upregulated miRs (miR-18a, miR-19a, and miR-181a) were selected for qRT-PCR validation. The average relative quantities (log10 2^-ΔCq ) for cfa-miR-18a, cfa-miR-19a, and cfa-miR-181a for CMEC and CMT groups were significantly higher in the CMT group compared to CMEC (Fig. 3a-c). Each of these three miRs showed very similar fold-change between CMT and CMEC for both small RNA deep-sequencing and qRT-PCR assays (Fig. 3d). Table 2 compares foldchange and p-values for RNAseq and qRT-PCR data.
In silico analysis of microRNA targets A set of 16 differentially expressed miRNAs from this data set were selected for in silico analysis of canine predicted miR targets (Fig. 4a). The number of predicted genes identified per miRNA ranged from 124 to 751 genes for total predictions and from 24 to 300 genes for predictions having a score greater than 80 (Table 3). To gain an appreciation for the types of biological processes associated with these 16 miRNAs, the complete set of gene targets for all 16 miRNAs was used for a gene ontology biological process enrichment analysis (Additional file 2: Table S2).
The complete set of all 2323 miRNA target genes with scores greater than 80 were selected. Following removal of redundant target genes in the list, we identified a set of 1849 unique target genes. The filtered list was used in the DAVID database and, of the 1849 gene symbols uploaded, 1819 were successfully mapped to canine genes using the DAVID gene list manager. In contrast, the gene list manager identified 1737 of the 1849 gene symbols as human genes.

Gene ontology enrichment across all gene targets
A total of 145 gene ontology terms were identified as enriched, with 85 having p-values less than or equal to 0.06. Although we considered GO annotations with p-value < 0.06, we focused on those with p-value < 0.05 ( Table 4). The complete set of all identified terms for the combined target genes is included in Additional file 3: Table S3, along with the list of specific genes associated with each enriched GO term. A set of 34 enriched terms are also included (Table 4). These terms provide some insight into the cellular role of the 16 miRNAs.
Among the enriched biological processes are a number implicated in cell division including: G1/S transition of mitotic cell cycle; negative regulation of apoptotic process; intrinsic apoptotic signaling pathway by p53 class mediator; positive regulation of cell proliferation and positive regulation of apoptotic process. Additional terms, associated with specific cellular processes involved in development and differentiation were also identified such as: regulation of establishment of cell polarity; regulation of       3 qRT-PCR validation of RNAseq data. a-c Relative quantification (log10) for selected validation targets miR-18a (a), miR-19a (b), and miR-181a (c). Relative quantification was calculated for each biological replicate according to the eq. 2 -ΔCq , with cel-miR-39 as spike-in exogenous control and miR-16 as an endogenous control; experiments were performed in triplicate from cell culture to RNA extraction, cDNA synthesis and qRT-PCR. Data were not normally distributed and compared by non-parametric, non-directional Mann-Whitney test. p < 0.05 was considered statistically significant. The black horizontal line represents the group mean and the vertical "whiskers" represent ±1 SD. d Comparison of foldchange between microRNA deep-sequencing and manual stem-loop qRT-PCR assays for selected targets miR-18a, miR-19a, and miR-181a. Data were normalized using the 2 -ΔΔCq method. The average group Cq for cfa-miR-16 (endogenous control) and cel-miR-39 (exogenous spike-in control) were used as housekeeping genes for normalization. White bars represent relative fold-change for RNAseq data, black bars represent fold-change for qRT-PCR (3 experimental replicates) H3-K9 trimethylation; histone H2A ubiquitination; and histone H3-K4 trimethylation.

Small subset of miRNAs having potential for clinical biomarker relevance
Based on the initial gene ontology enrichment analysis of all target genes associated with the set of 16 miRNAs, a choice was made to select a relatively small subset that might be suitable candidates for downstream clinical biomarker applications. The overarching goal in selecting the miRNA subset was focused upon (1) maximizing the representation of target genes associated with the set of enriched gene ontology terms while simultaneously (2) minimizing the number of miRNAs selected. Analysis of target gene representation associated with enriched gene ontology terms identified three miRNAs: cfa-miR-18a, cfa-miR-19a, and cfa-miR-181a. Together, these three miRNAs contain target genes for all but one enriched gene ontology term listed in Table 3 (histone H3-K9 trimethylation). Moreover, the overlap of target genes between these miRNAs was less than 7% (Fig. 4b).

Discussion
To the authors' knowledge, this is the first study reporting secretion of exosome-like extracellular vesicles by canine mammary epithelial cells in vitro. Similar to previous reports of canine exosomes found in urine and serum/plasma, the vesicles were irregularly rounded, occasionally "cup-shaped," and immunopositive for the transmembrane tetraspanin protein CD9, known to regulate the progression of cancer [37][38][39]. These findings strongly support these vesicles being exosomes, the small subcellular particles 50-200 nm in diameter that are actively shed from multivesicular bodies of parent cells. Exosomes contain proteins, peptides, mRNA, and microRNA, have been shown to be taken up by distinct cells through endocytosis, and they may play a role in distant cell-to-cell communication, especially in the context of neoplasia [37,38]. As expected, the cell-free conditioned media that contained these exosomes was highly enriched in hundreds of distinct microRNAs. The microRNA profile of normal and malignant canine mammary exosomes was distinct, and yielded a number of significantly up-regulated and down-regulated miRs that may represent putative biomarkers of mammary neoplasia. These findings largely corroborate previous studies on miRNA in canine mammary neoplasia. Several studies of miRNA expression in canine mammary tumor tissue and CMT cells versus normal mammary tissue found significantly increased miR-29b, which was also significantly upregulated in the CMT exosomal RNA in the present study (along with the closely related miR-29c) [24,25]. One of those studies also found miR-181b was significantly upregulated in the tubulopapillary carcinoma subtype [24].
Interestingly, our results differ substantially from those of von Deetzen et al. (2014), although it should be pointed out that in that study, the authors used miR-181b and miR-155 as housekeeping controls for qRT-PCR normalization, and both of those miRs appear to be upregulated in our data and previous studies of CMT [24,26,40]. Our results also diverge from Bulkowska et al. (2017), where several relevant up-regulated miRNA in the present study (such as miR-19a, miR-29b/c, and miR-181a) were shown to be down-regulated in malignant mammary carcinomas with metastasis [27]. This could be explained by the dramatic changes in tumor cell phenotype and gene expression in metastatic lesions compared to their matched primary tumor, such as occurs in the epithelial-to-mesenchymal transition (EMT) that has been documented in metastatic canine mammary carcinoma [41].
miR-126 was previously identified as a circulating biomarker of multiple canine cancers including mammary carcinoma, and it is up-regulated in our CMT exosomal RNA [28]. This would fit with the hypothesis that canine mammary carcinoma cells secrete exosomes containing miR-126 (among other miRNA) into circulation. However, the other putative biomarker in that study, miR-214, was strongly down-regulated in our CMT exosomal RNA. One possibility for this mismatch could include secretion of miR-214 by other cells than mammary carcinoma cells (i.e. cells of the immune system, stroma, or other organs). Another possibility is a mismatch between tumor cell, exosomal, and circulating microRNA profiles. Supporting this hypothesis, previous cell culture work in our lab has shown a complex relationship between exosomal miRNA profiles and the miRNA profile of the parent cell lines. Several miRs, including miR-18a, miR-19a, miR-29c, miR-181a, miR-215, miR-345, miR-371, and miR-1841, are up-regulated in both CMT parent cells and their exosomes [9,40]. However, miR-31 and miR-34c had mixed expression patterns in the parent cells despite being uniformly upregulated in exosomes, and miR-495 was strongly down-regulated in all CMT parent cells while being up-regulated in the exosomal RNA population [9,40]. This preliminary finding may suggest there is an active selection or enrichment process of particular miRNA within exosomes, or a negative feedback loop with the targets they regulate.
(See figure on previous page.) Fig. 4 Overlap of target genes in the miRNA subset. a Supervised absolute expression heat map for 16 canine miRNAs in CMEC vs. CMT exosomal RNA samples. Three biological replicates (CMEC1, CMEC2, and CMEC3) corresponding to normal mammary tissue exhibit relatively low levels of expression for the first eleven miRNAs while the last five miRNAs exhibit considerably higher levels of expression. The pattern is reversed in the mammary tumor samples (CMT12, CMT27, CMT28, CMT47, and CMT119) as visualized in the right side of the heat map. Because of the dichotomous pattern of expression across the control and mammary samples, these miRNAs may represent valuable candidates for clinically relevant biomarkers. b A miRNA subset was selected from among the 16 miRNAs (a). Overlap of target genes is indicated by both numbers of target genes and percentage of total genes (sum of each miRNA's target gene set). Cfa-miR-18a represents an miRNA with relatively low expression in normal mammary tissue and higher expression in mammary tumor samples. Similarly, cfa-miR-19a also exhibits low expression in normal mammary samples and higher expression in mammary tumor samples. In contrast, cfa-miR-181a exhibits considerably higher expression in the normal mammary tissue compared to 18a and 19a. Additionally, 181a expression in mammary tumor samples is greater than any other miRNA among the 16 miRNAs represented in (a) Many of these exosomal miRNAs were predicted to target dozens or hundreds of individual gene targets. Of particular note, miR-18a, miR-18b and miR-22 were highly up-regulated in the CMT exosomal RNA group and predicted to have an extremely high likelihood of targeting to 3' UTR of the estrogen receptor ESR1α mRNA (miRDB score of 99 for miR-18a and miR-18b, and 87 for miR -22). A number of other additional miRNA, including miR-181a, miR-181b, miR-181c, miR-181d, miR-19a, miR-19b, miR-148b, miR-203, miR-323, miR-874, and miR-486 were also predicted to negatively regulate ESR1α, although with a lower probability than the > 80 score threshold (Additional file 4: Table S4). Although targets with score greater than 80 have a greater likelihood of being true targets, some targets with scores below 80 may also be real. While it has long been known that human and canine mammary neoplasms lose ESR1α expression along with increasing grade and stage, this finding may indicate that miRNA such as miR-18a contribute to this loss of hormone receptor activity [12,15]. If this is verified in vivo in clinical patients, it may suggest miR-18a and others represent a non-invasive marker of CMT hormone status and phenotype, and could provide one potential mechanism for the loss of ESR1α with increasing CMT grade. For gene ontology and functional enrichment analysis, the DAVID resource provides both p-values and Benjamini corrected p-values to aid investigators in the analysis process. The higher stringency of Benjamini corrected p-values dramatically reduces the number of significant results. On the one hand, this is a valuable way to reduce false positives from the analysis. On the other hand, a bioinformatics analysis aimed at providing clues as to the biological roles of miRNAs exhibiting altered expression profiles in normal mammary samples versus mammary tumor samples benefits greatly from broader inclusion criteria. Therefore, both the p-value and Benjamini corrected p-value were reported, and inclusion criteria for enriched gene ontology terms was set with a threshold p-value less than or equal to 0.06. This approach parallels the method described in Irizarry et al. (2016) and allows for retrieval of relevant functional annotation occurring at (or near) the boundary of p-value significance [36].
This set of enriched biological processes is particularly interesting in the context of mammary tumorigenesis and progression. The cell cycle-associated processes clearly relate to accelerated proliferation that can contribute to malignant transformation. For example, the transition between G1/S is a critically regulated check point in the cell cycle [42]. Aberrations in the control of G1/S transitions can contribute to aberrant cell division and undesirable cell proliferation [42]. Similarly, positive and negative regulation of apoptosis impacts which cells survive in the cellular population. Altered expression of pro-apoptotic and/or anti-apoptotic gene products may adversely contribute to increased susceptibility to tumorigenesis in mammary tissues.
Equally important are the biological processes associated with cellular differentiation, adhesion and stem cell maintenance. Molecular factors underlying cellular differentiation may contribute to altered genetic programs within the cell. Abnormal expression of these factors may alter cellular programs leading to dysregulation of the cell cycle. Likewise, altered levels of genes implicated in cell adhesion may contribute to metastatic phenotypes that shift the clinical status of tumors from benign to malignant.
Finally, biological processes regulating chromatin have tremendous potential to dramatically alter the long-term genetic programs associated with cells. Modification of histones through methylation, acetylation, deacetylation, and ubiquitination directly modulate which chromatin regions are accessible for gene expression [43]. Silenced regions encoding tumor suppressors may shift cells towards a more oncogenic potential [43].
This study has a number of important limitations. First, the number of biological replicates was relatively small, which was a function of both the high cost of the RNAseq methodology, as well as the difficulty harvesting and maintaining normal canine mammary epithelial cells in culture. However, the RNAseq dataset identified hundreds of miRs, many of which were significantly different, and many of these miRs match other studies in the human and veterinary literature. Another limitation is that there is currently no consensus as to the optimal way to normalize exosomal microRNA qRT-PCR data. This was dealt with through a commonly used approach of relative expression based on ΔΔ-Cq normalization to pooled endogenous (miR-16) and exogenous (cel-miR-39) controls, which yielded very similar fold-change between RNAseq and qRT-PCR (Fig. 3d). Furthermore, the use of specific stem-loop primers and sequence-specific probes, rather than non-specific intercalating dye methods (i.e. SYBR green), increased the specificity and robustness of this data.

Conclusions
These data suggest that as in women with BC, CMT cells shed exosomes enriched in differentially expressed miRNA, especially miR-18a, miR-19a, and miR-181a. Preliminary in silico evidence suggests these miRNAs may modulate biological processes associated with, or contributing to, the balance between normal and neoplastic states. A miRNA population predicted to regulate so many aspects of cellular proliferation and hormone activity suggests that these miRs are not just inert cellular by-products, but may actually play an active part in neoplastic transformation and/or progression, and evidence that they are actively selected for secretion. Furthermore, the identification of these miRs in secreted exosomes raises the possibility that they may be shed into biofluids such as blood, urine, and breast milk, allowing their use as minimally-invasive biomarkers with mechanistic and prognostic relevance, and the similarity between canine and human breast cancer exosomal miRNA profiles may have significance for translational research, and future studies need to experimentally validate that these miRNAs regulate the predicted targets (such as miR-18a and ESR1α).

Additional files
Additional file 1: Table S1. This table contains normalized RNAseq expression data for all miRNA in this study in Reads per Million (RPM). (XLSX 76 kb) Additional file 2: Table S2. This table contains all genes predicted to be targeted and probability of binding their 3' UTR for the set of 16 miRNA in Fig. 4. (XLSX 109 kb) Additional file 3: Table S3. This provides a list of all Gene Ontology (GO) terms enriched in our set of 16 miRNA (Fig. 4). (XLSX 31 kb)