In 27 sequencing libraries, there were 558 unique miRNA candidates sequenced across all species, and only 395 miRNAs which were identified at least in four libraries were considered further (Table S1A). Of these, the top 10 highly expressed miRNA in four buffalo libraries, which accounted for 75.45±1.18% of all aligned reads, were evolutionarily conserved across multiple species. These miRNAs represented four different miRNA families; miR-let-7 (bta-let-7a-5p, bta-let-7b, bta-let-7c, bta-let-7e and bta-let-7f), miR-30 (bta-miR-30a-5p, bta-miR-30d and bta-miR-30e-5p), miR-148 (bta-miR-148a), and miR-26 (bta-miR-26a). In addition, four of these miRNAs (miR-148a, miR-30a-5p, bta-miR-30d, bta-let-7c) shared common ranking as top 10 expressed miRNAs when parsing the sequence data based on different species (buffalo, cow, pig, human and panda milks) (Table S1A), suggesting these miRNAs could be important nutritional components of milk [16]. In details, the most abundant was bta-miR-148a accounting for 38.88±3.12% in each buffalo library, which was consistent with previous reports in cow [17], pig [12] and human [9] milks. Recent characterizations of the human salivary exosome have also highlighted the miR-148a expression at high levels that expected to modulate oral cavity defence [18]. To date, although the underlying regulatory mechanism of the most predominant miRNA (miR-148a) packaged into milk exosomes was still unclear, it was intriguing to suggest that miR-148a was a potential biomarker for the quality control of mammalian milk [19].The miR-30a was highly enriched in exosomes from the serum of acute myocardial infarction patients in vivo, and inhibition of miR-30a or exosome release contributed to maintaining of autophagy after hypoxia [20]. In addition to miR-30a, exosome-enriched miR-30d in endometrial fluid was proposed to be taken up by the pre-implantation embryo, thereby resulting in the observed modifications to the transcriptome and embryo adhesion [21]. The enrichment of let-7 miRNA family in the extracellular fractions, particularly, in the exosomes from gastric cancer cells reflected their oncogenic characteristics including tumorigenesis and metastasis [22]. These data implied the roles of highly expressed miRNAs with a specified function in the milk exosomes. To further identify the potential function of the top expressed miRNAs, target prediction was performed using miRanda software. A total of 400 putative target sites for the top 10 expressed miRNAs were identified (Table S1B), and then all of these target candidates were submitted for homology and functional annotation using an online version of the DAVID program. These target genes belonged to 5 KEGG annotated categories that significantly related to a wide variety of disease resistance, immune responsiveness and basic metabolism events, including microRNAs in cancer, lysosome, glycerophospholipid metabolism, proteoglycans in cancer, galactose metabolism and neurotrophin signaling pathway (Table S1C).
For validation and identification of species-specific miRNAs in buffalo milk, we compared the sequencing libraries between buffalo and non-buffalo milks and found that 27 indexed libraries were obviously divided into five groups based on breed differences (Fig. 1). To comparison with miRNA expression in non-buffalo milks, a total of 32 miRNAs in buffalo milk were significantly up-regulated (Fig. 2A), while 16 were significantly down-regulated (Fig. 2B) using the R Bioconductor package EdgeR analysis with TMM-normalized algorithm. In addition, target prediction revealed a total of 4372 miRNA-target interactions between significantly expressed miRNAs and NCBI RefSeq genes. This resulted in 1759 unique genes targeted by up-regulated miRNAs (Table S1D), 1761 genes targeted by down-regulated miRNAs (Table S1F), and 477 genes targeted by both up and down-regulated miRNAs. Especially, pathway analysis revealed that 39 predicted targets of up-regulated miRNAs were significantly enriched in categories annotated as a role in endocytosis (Table S1E). The endocytosis pathway, for instance, has been highlighted as an important step in the delivery of bovine exosomes and their cargo to human vascular endothelial cells and peripheral tissues by endocytosis [23]. Additionally, we also revealed many of the other pathways which have been previously implicated in series basic metabolism processes in animals. These pathways included metabolic pathways enriched with 127 targets, glycerophospholipid metabolism enriched with 18 targets, protein processing in endoplasmic reticulum enriched with 26 targets, glycosaminoglycan biosynthesis enriched with 6 targets, biosynthesis of amino acids enriched with 13 targets, etc (Table S1E). Taken together, these analyses strongly suggested that up-regulated miRNAs that were differentially expressed between buffalo and non-buffalo milks were key regulators of the host metabolism processes shuttling within exosomes by endocytosis. On the contrary, the targets of down-regulated miRNAs were annotated with a total of 31 significantly KEGG pathways (Table S1G). In addition to metabolic pathways shared with the target annotation of up-regulated miRNAs, computational analysis of the down-regulated miRNA targets has identified the statistical over-representation of several terms previously implicated in response to immunity events, such as T cell receptor signaling pathway, NF-kappa B signaling pathway, AMPK signaling pathway, HTLV-I infection, lysosome, inflammatory mediator regulation of TRP channels, etc. The evidences indicated that the down-expressed miRNAs identified in this study were likely central regulators of the immune response and thus represented potential targets or novel biomarkers of infection and inflammation.