- Research article
- Open Access
MicroRNA expression profiles of bovine milk exosomes in response to Staphylococcus aureus infection
BMC Genomics volume 16, Article number: 806 (2015)
Milk exosomes are a rich source of microRNAs (miRNAs) that are protected from degradation. Ingestion of milk and subsequent absorption of miRNAs into recipient cells by endocytosis may play a role in the regulation of neonatal innate and adaptive immunity. In contrast, the miRNA content of milk exosomes may also be indicative of a lactating animal’s health; whereby, the presence or absence of specific miRNAs could serve as biomarkers for early detection of bacterial infection that can lead to mastitis. In the present study, we therefore analyzed and compared miRNA expression profiles of milk exosomes from four Holstein cows obtained during mid-lactation prior to and after infection (48 h) of the mammary gland with Staphylococcus aureus.
Milk exosomes, purified from control and S. aureus infected cows, were extracted for RNA. Following preparation indexed libraries from both groups the samples were subjected to next generation sequencing.
Next generation sequencing of eight, unpooled small RNA libraries derived from milk exosomes produced about 60.5 million high-quality, bovine-specific sequence reads for comparison of miRNA expression between treatments. Sequence identity analysis showed the miRNAs make up about 13 % of the average RNA content of these exosomes. Although 417 known bovine miRNAs were identified, miRNAs represented the least diverse class of RNA accounting for only 1 % of all unique sequences. The 20 most prevalent unique sequences within this class accounted for about 90 % of the total miRNA-associated reads across samples. Non-annotated, unique reads provided evidence for another 303 previously unknown bovine miRNAs. Expression analyses found 14 known bovine microRNAs significantly differed in frequency between exosomes from infected and control animals.
Our survey of miRNA expression from uninfected milk exosomes and those produced in response to infection provides new and comprehensive information supporting a role for delivery into milk of specific miRNAs involved in immune response. In particular, bta-miR-142-5p, and −223 are potential biomarkers for early detection of bacterial infection of the mammary gland. Additionally, 22 mammary-expressed genes involved in regulation of host immune processes and response to inflammation were identified as potential binding targets of the differentially expressed miRNAs.
MicroRNAs (miRNAs) are a class of endogenous non-coding small RNA molecules, typically 22 nucleotides in length, which bind primarily to the 3'UTR of target mRNAs to repress translation and/or accelerate the decay  of up to 30 % of all expressed transcripts . Numerous studies have found roles for miRNAs in regulation of gene expression for important biological processes including cellular proliferation and differentiation, tissue development, and immune response. Of interest for our study, miRNAs help regulate the development of immune cells and modulate the innate and adaptive immune responses . Some specific examples include the role of miR-150 in inhibiting synthesis of the transcription factor c-Myb to help regulate B-cell differentiation , and miR-126 targeting of insulin regulatory subunit-1 transcripts to positively control the fate of B-cells . Similarly, miR-146a functions as a negative regulator of TNF receptor-associated factor 6 and Interleukin-1 receptor-associated kinase 1 transcripts during and/or after the innate immune system responds to a bacterial infection . These examples demonstrate how microRNAs can exert some regulatory control of the immune system within a cell. More recent studies have shown that extracellular miRNAs stably present in body fluids through encapsulation into exosomes or microvesicles; may also play a significant role in intercellular communication during an immune response [7–9]. Therefore, the presence of specific miRNAs in exosomes may be indicative of various pathological conditions and provide biomarkers for detection of certain disease conditions .
Despite the recognized importance of miRNAs in host immune response to disease in humans and mice, their comparative roles in regulating the immune response of livestock are still in preliminary stages of study. Expression surveys of bovine miRNAs have been reported for embryos , muscle , adipose [12, 13], mammary gland , testicular and ovarian tissues [14–16], and various tissues involved in immune response ; but none of these reports have interrogated miRNA expression changes in response to infection. More recently, Lawless and colleagues reported two studies surveying miRNA expression changes in response to bacterial infection of a bovine cell line  and circulating monocytes from blood and milk . In the latter report, bovine specific miRNAs were identified as amplifiers of monocyte inflammatory response networks and repressors of several metabolic pathways in response to Staphylococcus uberis (S. uberis) infection of the mammary gland.
Based on these findings and the hypothesis that exosomes may contain specific miRNAs indicative of infection means that characterizing the miRNA profile of bovine milk exosomes could provide further understanding of host immune response signaling during lactation. It is well known that human milk is the most essential nutritional source for infant health during postnatal development , containing nutrients that closely match infant requirements for brain development, growth, and a healthy immune system [21, 22]. In the past decade, many immune-related substances like secretory immunoglobulins, leukocytes, and antimicrobial factors such as lysozyme, lactoferrin, and oligosaccharides have been detected in milk . More specific analyses of the exosome component of human milk found abundant immune-related proteins  as well as miRNAs [8, 9]. The presence of miRNAs has also been confirmed from cursory surveys of bovine [25–27] and porcine milk . Taken together, the immune-related factors in milk exosomes could impart some pathogen resistance to the newborns lacking a fully developed immune system . Some miRNAs may play a role in these effects [25, 26], since there is evidence that exosome miRNAs are readily absorbed within the neonate digestive tract after ingestion of milk .
Thus, the goal of our study was to comprehensively survey by deep sequencing the miRNA content of bovine milk exosomes derived from both healthy and infected mammary glands, and then attempt to identify if host response to infection significantly changes the content of specific miRNAs in exosomes. We chose to harvest exosomes from lactating Holstein cows challenged with Staphylococcus aureus (S. aureus); because Holsteins are the primary source of milk-related products consumed by humans globally, and S. aureus is a leading causative agent of bovine mastitis. Furthermore, identification of differentially expressed exosome miRNAs in this infection model serve as potential molecular targets for development of biomarkers assays to provide early detection of sub-clinical mastitis.
Results and discussion
Sequence analysis of milk exosome miRNAs
Next generation sequence analysis of the eight small RNA libraries derived from milk exosomes prior to (N = 4 control) and 48 h post S. aureus infection (N = 4 infected) yielded more than 76.4 million 36 nt sequences. After filtering, approximately 17.7 and 54.3 million high-quality sequence reads from the control and infected libraries were available for reference genome alignment, respectively (Additional file 1). The unbalanced read results between treatments was due to an intentional increase of sequencing coverage of the infection sample 2 library, which showed the rate of unique sequence discovery was reduced with additional depth of sequencing but still comparable to the other libraries with fewer sequence reads (Additional file 1 – last row). In total, all reads represented about 1.8 million unique sequences, and SOAP alignment  to the bosTau6 reference genome assembly  resulted filtered this set of unique sequences down to 747,169 bovine specific sequences corresponding to about 60.5 million total sequences (84 % of total) for further analyses. Consistent with a previous report , the relatively high sequence depth (avg. of 80 reads per unique sequence) and multiple, unpooled samples (N = 8) provided enough replication for bona fide discovery of potentially novel miRNAs expressed from the bovine genome.
Sequence identities were assigned to the 747,169 unique sequences by alignments to Rfam and bovine-specific sequences within mirBAse, Repeat and Reference mRNA databases. The average RNA content of milk exosomes based on alignments of mapped, annotated reads was approximately: 1) 58 % other non-coding RNA (mostly tRNA with some rRNA, snRNA and snoRNA – Additional file 1B), 2) 13 % known bovine miRNA, 3) 23 % repetitive sequences, 4) 1 % mRNA (unspliced and mature), and 5) 4 % non-annotated sequence (Additional file 1). Comparing the ratio of total sequence reads to unique sequences for each identity category revealed miRNA associated sequences were the least diverse class of RNA present in milk exosomes with mRNA and other non-coding RNA associated sequences having 500- and 50-fold more observed diversity, respectively. Because the purpose of this study was to characterize miRNA-associated sequences and attempts to identify bovine snoRNA and piRNA identities were hampered by lack of annotation for the bovine reference genome, none of the other categories of sequence were analyzed further. However, non-miRNA sequence data is available for future analyses through GEO accession GSE55144.
In order to examine the potential miRNA content of milk exosomes more deeply, the 4732 miRNA-associated unique sequences (<1 % of all unique RNA sequences) were clustered into homolog groups based on annotation of best match to miRBAse. Using SOAP and mirDeep2  methods, identities of unique sequences could be assigned to 337 miRNA homolog groups corresponding to 376 miRNA loci and 411 miRNA homologs corresponding to 452 miRNA loci, respectively (Additional file 2A). The union of these two results yielded 417 known miRNAs corresponding to over 53 % of known bovine miRNAs in miRBase Release 19.
The reduction of 4732 unique sequences to a set of 417 homologs is most likely attributable to the presence of isomiRs, which are heterogeneous variants of canonical miRNA species . This finding was supported by variation in sequence read length (1–5 nt) for some unique sequence alignments within a homolog groups (Additional file 3). However, miRNA length distribution analysis (Fig. 1) showed the highest count of miRNA-associated reads was 22 nt in length (40.14 %), and the majority of unique sequences reads ranged within the expected length of 21–23 nt. Interestingly, a significant portion of unique reads was also present in the length ranging from 18–20 nt. This result requires further investigation, but the length variation may be attributed to enzyme modifications and imprecise processing of primary or precursor miRNAs by Drosha and Dicer enzymes, such as RNA editing in miRNA-mediated gene silencing , 3'-editing , and degradation of microRNAs by a family of exonucleases .
Looking at diversity within the 417 miRNA homolog groups, the number of total sequences representing isomiRs ranged from one to hundreds. As expected, the most highly expressed isomiR usually corresponded to the bovine reference sequence in miRBase (Additional file 3). A total of 201 miRNAs were found to have at least one isomiR expressed at a level of ≥10 sequence reads and more than 1432 different isomiRs were identified in total. The group designated bta-mir-2904 had the highest diversity of isomiRs at a level of ≥10 sequence reads with counts of 106 variants at the 5' end and 102 variants at the 3' end. The presence of such miRNAs variants was previously suggested to be cell type specific, have functional differences, and vary in their response to biological stimuli . In our dataset, about 38 % of the isomiRs were variant at the 5' end, which should change the seed sequence (2–7 nt at the 5' end), and thereby effect changes in target mRNA and possibly miRNA function . Overall, these results fully illustrated that miRNA diversity in exosomes is greatly extended through the presence of multiple isomiRs, which is a result that cannot be elucidated by real-time PCR or microarray expression analyses.
After accounting for isomiRs, we evaluated the overall distribution of miRNA diversity based on homolog groupings. The 15 most prevalent miRNAs accounted for between 76-88 % of the total normalized sequence reads depending on analysis method, while the members of this list were nearly conserved between methods (Additional file 2B). Furthermore, eight of these miRNAs (miR-30a-5p, −148a, 22-3p, −141, −186, −26a, −181a and -27b) shared common ranking as top 10 expressed miRNAs when parsing the sequence data based on treatment (control and infected exosomes), suggesting these miRNAs could be important nutritional components of milk. Recently, similar miRNAs were found by sequencing of lactating caprine mammary gland [40, 41], bovine mammary epithelial cells , and porcine  and human  milk. In the latter two studies; miR-30a, −148a, −141 and -27b were also represented as the top 10 miRNAs. Previous miRNA expression studies using bovine mammary tissue [42, 43] or milk  also found miR-148a and -181a as significantly elevated during lactation, respectively. However, most of the miRNAs reported as elevated in mammary tissue during ovine and bovine lactation do not match the most prevalent miRNAs from exosomes . Some of these differences could be because previous studies only profiled expression for a pre-selected subset of miRNAs, which did not overlap with our Top 10 list. However, miR-148a-3p, which accounted for ~18.6 % of the total miRNA associated sequence reads, has previously been suggested as a nutritional biomarker corresponding to the protein content of various bovine-derived milk products .
Interestingly, feeding studies that have shown exogenous plant  or milk  miRNAs can be found in the sera and tissues and influence regulation of target genes in recipient animals. Combining the findings of these reports with the pleiotropic roles of the 10 most prevalent miRNAs in our study suggests the enrichment of specific miRNAs derived from mammary secretory cells may also influence development of the immune system of neonates. Specifically, the miR-30a-5p and -30d (members of the miR-30 family) are known to be involved in regulation of autophagy in cancer progression and treatment by suppressing the expression of beclin 1  and also cellular invasion and immunosuppression by targeting GalNAc transferase GALNT7 to increase synthesis of the immunosuppressive cytokine interleukin-10 . Also widely studied for possible involvement in tumor progression, miR-148-3p has been shown to directly target the TGIF2 gene in ovarian cancer , the drug metabolizing PXR gene in human cancer metastasis , and CAND1 expression in human prostate cancer to promote prostate cell growth . In the case of miR-27b, mRNA stability of PPARgamma is destabilized, which is often associated with chronic inflammatory diseases provoked by an immune response . MiR-27b has also been found to be degraded by a viral transcript in lytic murine cytomegalovirus (MCMV) infection , further highlighting its role in immunity. Another prevalent milk miRNA, miR-181a, regulates T-cell selection by altering sensitivity to peptide antigens, which is partly achieved through the down regulation of multiple phosphatases that act as negative regulators of T cell receptor signaling . Additional supporting evidence includes the roles of plasma miR-141 as a biomarker for colon cancers , miR-186 as a tumor-suppressor for the development and progression of lung adenocarcinoma , miR-26a in miRNA biogenesis to target Lin28B and Zcchc11 as a suppression mechanism for tumor growth and metastasis , and miR-22 targeting of estrogen receptor α mRNA to inhibit estrogen signaling associated with some forms of breast cancer . Finally, miR-191, a well characterized oncogenesis-related miRNA, is a biomarker of colorectal cancer , primary effusion lymphoma , and hepatocellular carcinoma . Even though there is compelling evidence that the most prevalent miRNAs in our study could potentially exert an influence on immune response, the specific functional roles of each miRNA needs further detailed investigations to obtain a thorough understanding of the specific targets and mechanistic effects of consumption of miRNA-loaded bovine milk exosomes by a recipient animal. This is especially relevant considering that effects of microRNA consumption are still not well validated .
Novel miRNA discovery
After assignment of sequence identities by various alignments (Additional file 1), nearly 38 % of all unique sequences had no identity or were considered non-annotated. A mireap analysis  of these sequences identified 562 potentially novel bovine miRNAs corresponding to 593 genomic loci (Additional file 4A). Interestingly, the total read counts for each of these potential miRNAs were much lower on average than for those matching known bovine miRs already present in miRBase. For instance, there were only 57 unique sequences with a total read count >100, and many of our initial novel miRNAs (N = 259) had read counts <10. Because known miRNAs use as little as 6–8 nucleotides on the 5' end to recognize target mRNAs [63, 64], the sequence of this seed region for the potentially novel miRNAs was compared to those seed regions found among 25,141 miRNAs in miRBase. A total 88.6 % of novel miRNAs had a seed sequence identical to those reported for known mature miRNAs, and this sequence conservation supports functionality for many of the miRNAs predicted by mireap analysis of non-annotated reads.
After removing the 259 unique sequences with <10 total reads counts, the remaining 303 unique sequences were also aligned to miRBase using BLAST. Only three unique sequences had 100 % identity with known miRNAs found exclusively in other species (eca-miR-1842, dre-miR-1388-5p, mmu-miR-219c-3p). In addition, 24 and 23 unique non-annotated sequences had significant homology with the bta-mir-2285 and bta-mir-2284 families, respectively (Additional file 4B). These results were consistent with a previous report that mir-2285 and −2284 families had possibly over 40 members spanning the entire bovine genome . Although we are confident these analyses identified previously undiscovered bovine miRNAs, none of these sequences were used for the differential expression analyses between exosomes from control and infected milk.
Differentially expressed miRNAs in response to S. aureus infection
Compared to previous NGS studies investigating the milk miRNAs in bovine , porcine  and human , we generated sequence reads from unpooled miRNA libraries corresponding to individual animal replicates under different treatment (pre- and post-infection). This allowed exploration of statistically significant changes in miRNA content of milk exosomes in early response to S. aureus infection of the mammary gland. A separation of control (pre-infection) from infected samples was confirmed by principal component analysis (PCA) of the normalized sequence reads representing the 417 miRNA homolog clusters (Fig. 2). PC1 explained 55.17 % of the overall miRNA expression variability, whereas PC2 explained 28.51 % of the variability.
The primary differential expression (DE) analysis was done using EdgeR and compared expression values of the unique miRNAs as determined by miRDeep2 normalization (Additional file 5). This approach detected a total of 14 DE miRNAs (Table 1). Heatmap analysis of these 14 DE microRNAs shows the clustering of infected from control samples (Fig. 3), ageing with the infection challenge statistics of these animals (Additional file 6).
As a comparison to the EdgeR analysis, a differential expression analysis was also done by SOAP (Additional file 7), and six common miRNAs (bta-miR-101, −142-5p, −183, −2285 g-3p, −223 and -99a-5p) with differential presence in milk exosomes between uninfected and infected animals were identified (Table 1).
Comparing our DE miRNA results from milk to previously reported changes in bovine mammary-derived miRNA expression revealed observations in support of our findings. For example, we detected increases of bta-miR-142-5p and −223 in milk exosomes 48 h post-infection, while both lactating mammary epithelium  and bovine monocytes  under infection challenge with S. uberis also were found to increase levels of miR-223. Other reports surveying bovine milk [26, 27] found higher levels of miR-223 in colostrum possibly in support of increased immunity for early neonates or the mammary gland during a period of higher susceptibility (post-partum) to bacterial infection. Lawless and colleagues also found significant increases in miR-142-5p in mammary monocytes post-infection . The cellular source of increases in our results, mammary epithelium or monocytes or both, remains to be determined; however, the lower levels of relative expression in comparison to the most prevalent miRNAs in exosomes suggests only a subset of cells responding to infection are contributing these RNAs to milk. The decreases in -15b, and -193a-5p found only using SOAP analysis are supported by similar findings in mammary epithelium  and monocytes , respectively. In contrast, we did not find any expression of bta-miR-2898, which was previously found to be upregulated in response to bacterial infection of mammary gland, where a potential role as a modulator of immune response was supported through differential regulation of A2M transcript binding .
Regarding the roles of these miRNAs in immune response, a previous study demonstrated elevated miR-142 in systemic lupus erythematosus patients compared to healthy controls . MiR-142 has also been shown to be abundant in T cells , which implied its role as an immune-relevant miRNA while its function remains unclear. On the other hand, miR-223 has a potential role in balancing metabolism and immune response during infection , as previous bioinformatic analysis provides support for its importance in down-regulating lipid metabolism ; a process very important during milk production. Beyond its role in response to mammary infection in cattle, miR-223 was also recently reported to regulate granulopoiesis . In spite of the limited knowledge of milk enriched miRNAs in response to bacterial infection, our results present significant information about highly expressed miRNAs in bovine milk exosomes from healthy and infected animals, and provides specific miRNAs that can be targeted for biomarker development as a potential means for early diagnostic detection of sub-clinical mastitis.
Prediction and analysis of target genes
A total of 7961 unique genes were predicted by RNAhybrid to be targeted based by the six DE miRNAs found in common between both analysis methods (bta-miR-101, −142-5p, −183, −2285 g-3p, −223 and -99a-5p). Of these, 4511 genes were predicted to be targeted by up-regulated miRNAs; 1500 were targeted by down-regulated miRNAs; and 1950 were targeted by both up and down-regulated miRNAs (Additional file 8). Because many of these target mRNAs may not actually be expressed during lactation, we filtered our candidate list of 7961 genes for mammary specific genes by only selecting those encoding proteins identified in previous reports [71, 72]. From this reduced list of 2350 genes encoding mammary derived protein found in exosomes, a total of 219 overlapping target genes were identified (Additional file 9A). DAVID analysis results showed 11 Biological Processes GO terms for 22 genes, which were significantly related to host immune processes and inflammation (Additional file 9B). We therefore hypothesized that differentially expressed miRNAs can target sequences in these 22 genes to regulate the immune responses with S. aureus infection. This new knowledge of milk miRNA expression between healthy cows and cows with mastitis will provide information important for the immune function in the mammary gland.
In conclusion, we have comprehensively analyzed the RNA content of milk exosomes from infected and uninfected animals using next generation sequencing to identify 417 known bovine miRNAs and 303 novel candidate miRNAs. Expression analyses identified six miRNAs that were significantly differentially present in exosomes in response to bacterial infection of the mammary gland, and provided two promising targets in bta-miR-142a and −223 for biomarker development of bacterial infection. Overall, this study expands the repertoire of bovine miRNAs and provides some specificity for the most prevalent miRNAs in milk, which may have effects on downstream gene expression through ingestion.
Milk exosomes collection
The United States Department of Agriculture (USDA), Agricultural Research Service (ARS), National Animal Disease Center (NADC) animal care and use committee approved all experimental animal procedures used in this study. No human subjects or data were used in this study. Milk exosomes were derived from Holstein cows (N = 4) housed at the NADC. Prior to infection all cows were healthy, and milk was confirmed free of bacteria. The exosome isolation methods and S. aureus challenge parameters were described previously for a proteome study of milk exosomes . Briefly, control milk was aseptically collected from each quarter during mid-lactation (90 days in milk) and infected samples were collected 48 h post-infection from the matching quarter infused with approximately 400 colony forming units (CFU) of S. aureus (Newbolt strain) suspended in 10 ml of PBS. The crude exosomes were purified by sucrose gradient centrifugation and filtration . The S. aureus infections were confirmed as previously described  and are summarized in Additional file 6.
Small RNA library sequence analysis
Total RNA was extracted from pelleted exosomes using the manufacturer’s protocol for TRIzol® Reagent (Life Technologies, Carlsbad, CA, USA). Presence of post-extraction RNA was assessed by using a RNA 6000 pico assay run an Agilent Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA). Eight Illumina-indexed libraries (four representing milk exosomes prior to infection and four from S. aureus infected samples 48 h post-infection) were constructed with an Illumina TruSeq small RNA library kit according to the manufacturers’ protocol (Illumina, Inc., San Diego, CA), which included size selection of the final library amplicons. All small RNA libraries were subjected to quality and purity check using Agilent’s 7500 DNA Bioanalyzer assay (Agilent Technologies, Santa Clara, CA, USA) and Qubit Fluorometric Quantitation (Life Technologies, Carlsbad, CA, USA). Quantitated libraries were diluted to 1 nM concentration and equal volumes were mixed for multiplex sequence analysis, and diluted to a final of 5 uM concentration for clustering on a Truseq flow-cell prior to analysis by single end sequencing (40 nt) on Illumina Genome Analyzer II.
Preliminary quality control analysis of the eight Fastq files generated using the FastQC package v0.10.0  aimed to filter low quality short reads, remove poly A, trim 3'/5' adapter, and format into non-redundant Fasta files. Cutadapt v1.3  was implemented to trim 3'/5' adaptor sequences, and reads less than 18 nt were removed from the dataset for analysis. The occurrence of each unique sequence read was tallied (Additional file 1). All unique sequence tags that passed through the above filtering criteria were mapped onto the bosTau6 version of the bovine reference genome assembly  using the SOAP algorithm  allowing up to 3 bp mismatches to detect isomiRs. Uniquely mapped tag sequences were aligned against Rfam v11.0, bovine miRNAs in miRBase v19.0, repeat database produced by RepeatMasker  and the coding genes of the reference genome . Based on the annotation of the best matches, these unique sequences were classified as other non-coding RNA, known miRNA, genomic repeats or degradation products of mRNA. Sequences that mapped to the reference genome, but did not fit any of the above annotation categories were classified as ‘non-annotated’ reads.
A second method of data processing, named “miRDeep2-method” , was also used with some improvements. Under this method, Trimmed reads were then further filtered using Sickle  to remove low quality sequences with quality values <20. Reads that successfully passed filtering were aligned to the bovine miRBase v19.0 using miRDeep v126.96.36.199 with a Perl script “quantifier.pl”. All miRNAseq data was submitted to the NCBI Gene Expression Omnibus (GEO) database with the experiment series accession number GSE55144.
miRNA expression analysis
To compare differentially expressed miRNAs between multiple samples through the miRDeep2-method, read count of each identified miRNA was normalized between libraries using upper-quartile normalization . The R v3.0.2 Bioconductor package EdgeR v2.4.6  was used to identify statistically differentially expressed miRNAs. Those miRNAs expressed at least once in control or infected samples were considered for further analysis. Differentially expressed miRNAs were defined as having a Benjamini and Hochberg corrected P value <0.05 . The heat map of the differentially expressed miRNAs was generated using Cimminer  under the Canberra and Ward settings for Distance Method and Cluster Algorithm settings, respectively.
As a secondary validation analysis using SOAP , count data was individually normalized across 8 libraries using the formula: normalized expression (NE) = actual miRNA count/total count of clean reads × 106 using SOAP. P-values were calculated from the normalized expression  in each replicate group compared between control and infected samples by:
where x and y represent normalized expression levels, and the N1 and N2 represented total sequence read count for a given miRNA from the infected and control small RNA libraries, respectively. A specific miRNA was deemed to be significantly differentially expressed when the P-value was <0.01, and there was at least a 2-fold change in normalized sequence counts across all four replicates. The significant differentially expressed miRNAs identified by both SOAP and miRDeep2 approaches were selected for further target prediction analysis.
Novel miRNA prediction
To investigate discovery of potentially novel miRNAs, all unmapped reads from the “SOAP-method” were run through Mireap , a prediction software that explores secondary structure characteristic, the Dicer cleavage site and the minimum free energy of input sequence reads. The following criteria  were used for identifying potentially novel miRNAs from unique sequence reads: 1) aligned to non-annotated regions of a bovine reference genome, 2) folded into hairpin secondary structures (with flanking sequence added) and unique sequence read was present in one arm of miRNA precursor, 3) demonstrated a 2 nt 3'-overhang between the sequence read (presumed miRNA mature strand) and its complementary strand (miRNA*), 4) folded hairpin miRNA precursors lacked large internal loops or bulges, 5) had a free energy of hybridization ≤ −18 kcal/mol as a folded hairpin structure, 6) carried a higher minimal free energy indexes (MFEIs) as a predicted precursor miRNA than other different types of RNAs.
Differentially expressed miRNAs from the exosomes were further analyzed to predict potential regulatory target transcripts. Bos taurus mRNA sequences were downloaded from the NCBI database and input into RNAhybrid software  based on a series of rules: 1) no mismatch present between base pairs 1–9 on the 5' end of the miRNA/target RNA duplex, 2) G:U pairing was permitted for up to 3 bases. In addition, we also used Reinhardt et al.’s data  gained by analysis of bovine milk exosome proteome infected with S. aureus to narrow down the pool of potential targets predicted by computational analysis. Only target sites which were identified by the two approaches were further considered. Using DAVID Bioinformatics Resources 6.7  only predicted targets which were identified by both computational and experimental analysis from infected exosomes were submitted to David Gene Ontology (GO) terms.
Bartel DP. MicroRNAs: genomics, biogenesis, mechanism, and function. Cell. 2004;116:281–97.
Lewis BP, Shih IH, Jones-Rhoades MW, Bartel DP, Burge CB. Prediction of mammalian microRNA targets. Cell. 2003;115:787–98.
Xiao C, Rajewsky K. MicroRNA control in the immune system: basic principles. Cell. 2009;136:26–36.
Xiao C, Calado DP, Galler G, Thai TH, Patterson HC, Wang J, et al. MiR-150 controls B cell differentiation by targeting the transcription factor c-Myb. Cell. 2007;131:146–59.
Okuyama K, Ikawa T, Gentner B, Hozumi K, Harnprasopwat R, Lu J, et al. MicroRNA-126–mediated control of cell fate in B-cell myeloid progenitors as a potential alternative to transcriptional factors. Proc Natl Acad Sci. 2013;110:13410–5.
Williams AE, Perry MM, Moschos SA, Larner-Svensson HM, Lindsay MA. Role of miRNA-146a in the regulation of the innate immune response and cancer. Biochem Soc Trans. 2008;36:1211–5.
Weber JA, Baxter DH, Zhang S, Huang DY, Huang KH, Lee MJ, et al. The microRNA spectrum in 12 body fluids. Clin Chem. 2010;56:1733–41.
Kosaka N, Izumi H, Sekine K, Ochiya T. MicroRNA as a new immune-regulatory agent in breast milk. Silence. 2010;1:7.
Zhou Q, Li M, Wang X, Li Q, Wang T, Zhu Q, et al. Immune-related microRNAs are abundant in breast milk exosomes. Int J Biol Sci. 2012;8:118–23.
Mondou E, Dufort I, Gohin M, Fournier E, Sirard MA. Analysis of microRNAs and their precursors in bovine early embryonic development. Mol Hum Repro. 2012;18:425–34.
Sun J, Li M, Li Z, Xue J, Lan X, Zhang C, et al. Identification and profiling of conserved and novel microRNAs from Chinese Qinchuan bovine longissimus thoracis. BMC Genomics. 2013;14:42.
Jin W, Dodson MV, Moore SS, Basarab JA, Guan LL. Characterization of microRNA expression in bovine adipose tissues: a potential regulatory mechanism of subcutaneous adipose tissue development. BMC Mol Biol. 2010;11:29.
Gu Z, Eleswarapu S, Jiang H. Identification and characterization of microRNAs from the bovine adipose tissue and mammary gland. FEBS Lett. 2007;581:981–8.
Huang J, Ju Z, Li Q, Hou Q, Wang C, Li J, et al. Solexa sequencing of novel and differentially expressed microRNAs in testicular and ovarian tissues in Holstein cattle. Int J Biol Sci. 2011;7:1016–26.
Hossain MM, Ghanem N, Hoelker M, Rings F, Phatsara C, Tholen E, et al. Identification and characterization of miRNAs expressed in the bovine ovary. BMC Genomics. 2009;10:443.
Tripurani SK, Xiao C, Salem M, Yao J. Cloning and analysis of fetal ovary microRNAs in cattle. Anim Repro Sci. 2010;120:16–22.
Coutinho LL, Matukumalli LK, Sonstegard TS, Van Tassell CP, Gasbarre LC, Capuco AV, et al. Discovery and profiling of bovine microRNAs from immune-related and embryonic tissues. Physiol Genomics. 2007;29:35–43.
Lawless N, Foroushani AB, McCabe MS, O’Farrelly C, Lynn DJ. Next Generation Sequencing Reveals the Expression of a Unique miRNA Profile in Response to a Gram-Positive Bacterial Infection. PLoS One. 2013;8:e57543.
Lawless N, Reinhardt TA, Bryan K, Baker M, Pesch B, Zimmerman D, et al. MicroRNA regulation of bovine monocyte inflammatory and metabolic networks in an in vivo infection model. G3. 2014;4:957–71.
Hoddinott P, Tappin D, Wright C. Breast feeding. BMJ. 2008;336:881–997.
Lawrence RM, Pane CA. Human breast milk: current concepts of immunology and infectious diseases. Curr Probl Pediatr Adolesc Health Care. 2007;37:7–36.
Paramasivam K, Michie C, Opara E, Fewell A. Human breast milk immunology: a review. Int J Fertil Womens Med. 2006;51:208–17.
Petherick A. Development: mother’s milk: a rich opportunity. Nature. 2010;468:S5–7.
Admyre C, Johansson SM, Qazi KR, Filén JJ, Lahesmaa R, Norman M, et al. Exosomes with immune modulatory features are present in human breast milk. J Immunol. 2007;179:1969–78.
Hata T, Murakami K, Nakatani H, Yamamoto Y, Matsuda T, Aoki N. Isolation of bovine milk-derived microvesicles carrying mRNAs and microRNAs. Biochem Biophys Res Commun. 2010;396:528–33.
Chen X, Gao C, Li H, Huang L, Sun Q, Dong Y, et al. Identification and characterization of microRNAs in raw milk during different periods of lactation, commercial fluid, and powdered milk products. Cell Res. 2010;20:1128–37.
Izumi H, Kosaka N, Shimizu T, Sekine K, Ochiya T, Takase M. Bovine milk contains microRNA and messenger RNA that are stable under degradative conditions. J Dairy Sci. 2012;95:4831–41.
Gu Y, Li M, Wang T, Liang Y, Zhong Z, Wang X, et al. Lactation-related microRNA expression profiles of porcine breast milk exosomes. PLoS One. 2012;7:e43691.
Goldman AS. The immune system in human milk and the developing infant. Breastfeed Med. 2007;2:195–204.
Baier SR, Nguyen C, Xie F, Wood JR, Zempleni J. MicroRNAs Are Absorbed in Biologically Meaningful Amounts from Nutritionally Relevant Doses of Cow Milk and Affect Gene Expression in Peripheral Blood Mononuclear Cells, HEK-293 Kidney Cell Cultures, and Mouse Livers. J Nutr. 2014;144:1495–500.
Li R, Li Y, Kristiansen K, Wang J. SOAP: short oligonucleotide alignment program. Bioinformatics. 2008;24:713–4.
Bos Taurus reference genome assembly. [http://hgdownload.cse.ucsc.edu/goldenPath/bosTau6/bigZips/refMrna.fa.gz].
Friedländer MR, Mackowiak SD, Li N, Chen W, Rajewsky N. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 2012;40:37–52.
Cloonan N, Wani S, Xu Q, Gu J, Lea K, Heater S, et al. MicroRNAs and their isomiRs function cooperatively to target common biological pathways. Genome Biol. 2011;12:R126.
Kawahara Y, Zinshteyn B, Sethupathy P, Iizasa H, Hatzigeorgiou AG, Nishikura K. Redirection of silencing targets by adenosine-to-inosine editing of miRNAs. Science. 2007;315:1137–40.
Landgraf P, Rusu M, Sheridan R, Sewer A, Iovino N, Aravin A, et al. A mammalian microRNA expression atlas based on small RNA library sequencing. Cell. 2007;129:1401–14.
Ramachandran V, Chen X. Degradation of microRNAs by a family of exoribonucleases in Arabidopsis. Science. 2008;321:1490–2.
Neilsen CT, Goodall GJ, Bracken CP. IsomiRs–the overlooked repertoire in the dynamic microRNAome. Trends Genet. 2012;28:544–9.
Grimson A, Farh KK, Johnston WK, Garrett-Engele P, Lim LP, Bartel DP. MicroRNA targeting specificity in mammals: determinants beyond seed pairing. Mol Cell. 2007;27:91–105.
Ji Z, Wang G, Xie Z, Wang J, Zhang C, Dong F, et al. Identification of Novel and Differentially Expressed MicroRNAs of Dairy Goat Mammary Gland Tissues Using Solexa Sequencing and Bioinformatics. PLoS One. 2012;7:e49463.
Li Z, Lan X, Guo W, Sun J, Huang Y, Wang J, et al. Comparative transcriptome profiling of dairy goat microRNAs from dry period and peak lactation mammary gland tissues. PLoS One. 2012;7:e52388.
Wang M, Moisá S, Khan MJ, Wang J, Bu D, Loor JJ. MicroRNA expression patterns in the bovine mammary gland are affected by stage of lactation. J Dairy Sci. 2012;9:6529–35.
Li Z, Liu H, Jin X, Lo L, Liu J. Expression profiles of microRNAs from lactating and non-lactating bovine mammary glands and identification of miRNA related to lactation. BMC Genomics. 2012;13:731.
Gigli I, Maizon DO. microRNAs and the mammary gland: A new understanding of gene expression. Genet Mol Biol. 2013;36:465–74.
Zhang L, Hou D, Chen X, Li D, Zhu L, Zhang Y, et al. Exogenous plant MIR168a specifically targets mammalian LDLRAP1: evidence of cross-kingdom regulation by microRNA. Cell Res. 2011;22:107–26.
Zhu H, Wu H, Liu X, Li B, Chen Y, Ren X, et al. Regulation of autophagy by a beclin 1-targeted microRNA, miR-30a, in cancer cells. Autophagy. 2009;5:816–23.
Gaziel-Sovran A, Segura MF, Di Micco R, Collins MK, Hanniford D, Vega-Saenz de Miera E, et al. miR-30b/30d Regulation of GalNAc Transferases Enhances Invasion and Immunosuppression during Metastasis. Cancer Cell. 2011;20:104–18.
Lujambio A, Calin GA, Villanueva A, Ropero S, Sánchez-Céspedes M, Blanco D, et al. A microRNA DNA methylation signature for human cancer metastasis. Proc Natl Acad Sci. 2008;105:13556–61.
Takagi S, Nakajima M, Mohri T, Yokoi T. Post-transcriptional regulation of human pregnane X receptor by micro-RNA affects the expression of cytochrome P450 3A4. J Biol Chem. 2008;283:9674–80.
Murata T, Takayama K, Katayama S, Urano T, Horie-Inoue K, Ikeda K, et al. miR-148a is an androgen-responsive microRNA that promotes LNCaP prostate cell growth by repressing its target CAND1 expression. Prostate Cancer Prostatic Dis. 2010;13:356–61.
Jennewein C, von Knethen A, Schmid T, Brüne B. MicroRNA-27b contributes to lipopolysaccharide-mediated peroxisome proliferator-activated receptor γ (PPARγ) mRNA destabilization. J Biol Chem. 2010;285:11846–53.
Marcinowski L, Tanguy M, Krmpotic A, Rädle B, Lisnić VJ, Tuddenham L, et al. Degradation of cellular mir-27 by a novel, highly abundant viral transcript is important for efficient virus replication in vivo. PLoS Pathog. 2012;8:e1002510.
Li QJ, Chau J, Ebert PJ, Sylvester G, Min H, Liu G, et al. miR-181a is an intrinsic modulator of T cell sensitivity and selection. Cell. 2007;129:147–61.
Cheng H, Zhang L, Cogdell DE, Zheng H, Schetter AJ, Nykter M, et al. Circulating plasma MiR-141 is a novel biomarker for metastatic colon cancer and predicts poor prognosis. PLoS One. 2011;6:e17745.
Cai J, Wu J, Zhang H, Fang L, Huang Y, Yang Y, et al. miR-186 downregulation correlates with poor survival in lung adenocarcinoma, where it interferes with cell-cycle regulation. Cancer Res. 2013;73:756–66.
Fu X, Meng Z, Liang W, Tian Y, Wang X, Han W, et al. miR-26a enhances miRNA biogenesis by targeting Lin28B and Zcchc11 to suppress tumor growth and metastasis. Oncogene. 2014;33:4296–306.
Pandey DP, Picard D. miR-22 inhibits estrogen signaling by directly targeting the estrogen receptor α mRNA. Mol Cellular Biol. 2009;29:3783–90.
Xi Y, Formentini A, Chien M, Weir DB, Russo JJ, Ju J, et al. Prognostic values of microRNAs in colorectal cancer. Biomark Insights. 2006;2:113–21.
O’Hara AJ, Vahrson W, Dittmer DP. Gene alteration and precursor and mature microRNA transcription changes contribute to the miRNA signature of primary effusion lymphoma. Blood. 2008;111:2347–53.
Elyakim E, Sitbon E, Faerman A, Tabak S, Montia E, Belanis L, et al. hsa-miR-191 is a candidate oncogene target for hepatocellular carcinoma therapy. Cancer Res. 2010;70:8077–87.
Wagner AE, Piegholdt S, Ferraro M, Pallauf K, Rimbach G. Food derived microRNAs. Food Funct. 2015;6:714–8.
mirreap - MicroRNA Discovery By Deep Sequencing. [http://sourceforge.net/projects/mireap/].
Bartel DP. MicroRNAs: target recognition and regulatory functions. Cell. 2009;136:215–33.
Lewis BP, Burge CB, Bartel DP. Conserved seed pairing, often flanked by adenosines, indicates that thousands of human genes are microRNA targets. Cell. 2005;120:15–20.
Guduric-Fuchs J, O’Connor A, Cullen A, Harwood L, Medina RJ, O’Neill CL, et al. Deep sequencing reveals predominant expression of miR-21 amongst the small non-coding RNAs in retinal microvascular endothelial cells. J Cell Biochem. 2012;113:2098–111.
Naeem A, Zhong K, Moisá SJ, Drackley JK, Moyes KM, Loor JJ. Bioinformatics analysis of microRNA and putative target genes in bovine mammary tissue infected with Streptococcus uberis. J Dairy Sci. 2012;95:6397–408.
Wang XG, Huang JM, Feng MY, Ju ZH, Wang CF, Yang GW, et al. Regulatory mutations in the A2M gene are involved in the mastitis susceptibility in dairy cows. Anim Genet. 2014;45:28–37.
Dai Y, Huang YS, Tang M, Lv TY, Hu CX, Tan YH, et al. Microarray analysis of microRNA expression in peripheral blood cells of systemic lupus erythematosus patients. Lupus. 2007;16:939–46.
Wu H, Neilson JR, Kumar P, Manocha M, Shankar P, Sharp PA, et al. miRNA profiling of naïve, effector and memory CD8 T cells. PLoS One. 2007;2:e1020.
Fukao T, Fukuda Y, Kiga K, Sharif J, Hino K, Enomoto Y, et al. An evolutionarily conserved mechanism for microRNA-223 expression revealed by microRNA gene profiling. Cell. 2007;129:617–31.
Reinhardt TA, Lippolis JD, Nonnecke BJ, Sacco RE. Bovine milk exosome proteome. J Proteomics. 2012;75(5):1486–92.
Reinhardt TA, Sacco RE, Nonnecke BJ, Lippolis JD. Bovine milk proteome: Quantitative changes in normal milk exosomes, milk fat globule membranes and whey proteomes resulting from Staphylococcus aureus mastitis. J Proteomics. 2013;82:141–54.
Fast QC Software. [http://www.bioinformatics.babraham.ac.uk/projects/fastqc/].
Cutadapt Algorithm. [https://cutadapt.readthedocs.org/en/stable/].
RepeatMasker Database. [http://www.repeatmasker.org].
UCSC Genome Browser - Coding Genes of the Bovine Reference Genome. [http://hgdownload.cse.ucsc.edu/goldenPath/bosTau6/bigZips/refMrna.fa.gz].
Sickle Algorithm. [https://github.com/najoshi/sickle/].
Bullard JH, Purdom E, Hansen KD, Dudoit S. Evaluation of statistical methods for normalization and differential expression in mRNA-Seq experiments. BMC Bioinform. 2011;11:94.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinform. 2010;26:139–40.
Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J Royal Stat Soc Series B (Meth). 1995;57:289–300.
Cimminer Heatmap Algorithm. [http://discover.nci.nih.gov/cimminer/].
Audic S, Claverie JM. The significance of digital gene expression profiles. Genome Res. 1997;7:986–95.
Xie FL, Huang SQ, Guo K, Xiang AL, Zhu YY, Nie L, et al. Computational identification of novel microRNAs and targets in Brassica napus. FEBS Letter. 2008;581:1464–74.
Rehmsmeier M, Steffen P, Höchsmann M, Giegerich R. Fast and effective prediction of microRNA/target duplexes. RNA. 2004;10:1507–17.
DAVID Bioinformatics Resources 6.7. [http://david.abcc.ncifcrf.gov/home.jsp].
Mention of trade names or commercial products in this article is solely for the purpose of providing specific information and does not imply recommendation or endorsement by the US Department of Agriculture. The USDA is an equal opportunity provider and employer. The contribution by scientists in the Animal Genomics and Improvement Laboratory and the National Animal Disease Center are supported by appropriated projects 1245-31000-104-00D and 5030-32000-102-00, respectively. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
The authors declare that they have no competing interests.
Research idea and experimental design: SGS, TAR, GDL, TSS; Experiments performed: GDL, TAR, TSS and KA; Data analysis: SGS, JJS and KA; Manuscript preparation: JJS, KA, and TSS. All the authors have read and approved of this manuscript.
Jiajie Sun and Kshama Aswath contributed equally to this work.
Summary (A) of Unique and Total Sequence Reads Counts and Summary (B) of Sequence Identities for Matches to Rfam. (XLSX 17 kb)
Total (A) and Normalized (B) Sequence Counts per Known Bovine miRNA Homolog Group. (XLSX 73 kb)
Analysis Results for Potential Bovine IsomiRs from Milk Exosomes. (XLSX 259 kb)
Identification of Novel Bovine miRNAs Based on Seed Sequence Alignments (A) and Sequence Identity of Novel Bovine miRNAs from BLAST Analysis of miRBase 19 (B). (XLSX 112 kb)
Full Results of EdgeR Differential Expression Analysis of Normalized Sequence Reads Corresponding to Known Bovine miRNAs by miRDeep. (XLSX 33 kb)
Animal Infection and Milk Exosome Data from Staph infected quarters. Control exosomes came from same quarter prior to infection. (XLSX 9 kb)
Full Results of SOAP-based Differential Expression Analysis of Normalized Sequence Reads Corresponding to Known Bovine miRNAs (Animal/Replicate 1–4). (XLSX 60 kb)
Predicted mRNA Targets of Differentially Expressed miRNAs in Milk Exosomes. (XLSX 2115 kb)
A) List of 219 Genes Encoding Proteins Differentially Expressed (DE) in Milk Exosomes and Predicted Targets of DE miRNAs and B) Gene Ontology Results of DAVID Analysis of 219 mRNA Targets Expressed in Mammary Gland. (XLSX 35 kb)
About this article
Cite this article
Sun, J., Aswath, K., Schroeder, S.G. et al. MicroRNA expression profiles of bovine milk exosomes in response to Staphylococcus aureus infection. BMC Genomics 16, 806 (2015). https://doi.org/10.1186/s12864-015-2044-9