The Role of Milk Exosomal microRNAs in Balancing the Nutrition between Buffalo Breeds with Low and High Milk Production

Milk is produced by the mammary glands of mammals and composed of multiple bioactive components, such as proteins, peptides, fatty acids, oligosaccharides, beneficial microorganisms. It is well recognized that milk comprises the primary source of nutrition for infants (Van Winckel et al., 2011). Breastfeeding is the dominant means to feed the milk to infants in almost all mammals. Accumulating studies have demonstrated that breast milk plays an irreplaceable role in compensating for developmental delays in immune function in the neonate and it can also reduce the permeability of the intestine to prepare it for extrauterine life (Goldman, 2000). Besides, breast milk can affect infants’ growth, development and self-regulation of food intake since it is enriched with hormones, neuropeptides and growth factors (Savino and Liguori, 2008). As a result, researchers have suggested that breast milk serves as the best nutrition, immune protection and regulation of growth, development and metabolism for infants of mammals (Goldman, 2012).
       
Interestingly, milk composition and production (normalized by body weight) vary from species to species, even different breeds of the same species (Aqib et al., 2019). For example, the milk yield of the yak cow is much less than that of dairy cattle. Nevertheless, their respective infants don’t exhibit significant differences in health conditions. This raises an interesting scientific question: Why could a species (or breed, individual) with lower milk production raise their infants as healthy as those with higher milk production?.
       
Effective food digestion is essential to guarantee the health of infants. It is well known that milk contains multiple enzymes (amylase and lipase) which aid digestion (Hamosh, 1996). Recently, some research reported that exosomes, a type of extracellular vesicle containing constituents (protein, DNA and RNA) of the cells that secrete them (Kalluri and LeBleu, 2020), can be isolated from both colostrum and mature milk (Admyre et al., 2007). Importantly, it has been reported that milk exosome microRNAs can be taken up by intestinal epithelial cells and stimulate gastric/pancreatic digestion (Liao et al., 2017). Moreover, Kahn et al. found that milk exosomal microRNAs can directly stimulate infant gut digestion (Kahn et al., 2018). These studies highlight that, apart from digestive enzymes, milk exosomes serve an essential role in facilitating food digestion in infants. We thus hypothesize that the property of milk exosomes in individuals with low milk production may be different from those with high milk production, which may contribute to the comparable health conditions of respective infants.
       
To test this hypothesis, we extracted the exosomes from milk of two buffalo breeds (Murrah and Nili-Ravi buffalos), both of which have individuals with low and high milk production. The microRNAs were profiled to unravel their differential expression (and their target genes) between low and high productive groups, underlying the potential roles of microRNAs in the two groups. This study supplies genetic data to helps improve buffalo molecular breeding.

All animal procedures (milking) were conducted in accordance with the Animal Care and Use Committee of the Guangxi Buffalo Research Institute, Chinese Academy of Agricultural Sciences. The experimental animals were kept on the farm of Guangxi Buffalo Research Institute, with proper feeding conditions. We selected six healthy Murrah female buffalos (three with low milk production M6629, M333, M4578 and the other three with high milk production M991, M3102, M1190) and six Nili-Ravi buffalos (three with low milk production M4556, M6815, M1975 and the other three with high milk production M6137, M4336, M4516) and milked them by hand. For pasteurization, the milk was heated at 71.6°C for 15 seconds. We chose fat-free pasteurized milk rather than ultra-heat treated milk, as the latter causes a loss of milk microRNAs and presumably exosomes (Kirchner et al., 2016). Exosomes were isolated by ultracentrifugation as described previously (Izumi et al., 2015; Wolf et al., 2015). The combined use of fat-free milk and exosome purification by ultracentrifugation minimized sample contamination with fat globules, which may also contain microRNAs (Alsaweed et al., 2016). Exosomes were authenticated by transmission electron microscopy, nanoparticle tracking analysis (NanoSight NS300, Malvern, Inc.) and western blot analysis following the guidelines provided by the International Society for Extracellular Vesicles (Kusuma et al., 2016; Lötvall et al., 2014). For each exosome, RNA was isolated from milk and milk exosomes following the protocol of treated with TRIzol reagent (ThermoFisher) and RNA with length 18-30nt was recovered by PAGE electrophoresis. cDNA synthesis was performed by using the miScript Reverse Transcription Kit following the manufacturer’s instructions (Qiagen), with the products subjected for library construction, followed by sequencing on BGISEQ-500 sequencing platform.
       
We used SOAPnuke (1.5.0) (Chen et al., 2018) to filter the raw sequencing data. Briefly, we removed the reads: 1) with more than four bases with base quality less than ten or more than six bases with base quality less than 13; 2) contaminated with 5' adaptor sequence; 3) without 3' adaptor sequence; 4) without inserts; 5) with consecutive ‘A’ bases longer than 60% of the total length; with length less than 15 nucleotides. Then, we used Bowtie2 (2.2.9) (Langmead et al., 2009) to align the clean data to the reference genome of buffalo with parameters: -q -L 16 -p 6 -phred64 -rdg 1,10 -rfg 1,10. Only the reads with at most one mismatch were retrieved for further analysis.
       
To carry out microRNA annotation on the retrieved sequences of a sample, we first aligned the mature miRNA sequences of all animals in the miRBase (V22) database (Kozomara et al., 2019) onto the buffalo genome by Bowtie2 (2.2.9) (Langmead and Salzberg, 2012). For those miRNAs that were completely mapped onto the genome, we extracted the aligned regions along with 150 bp flanking fragments and fed them into RNAfold (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi) to predict hairpin structure. The ones with a hairpin structure were used as miRNA database of buffalo. Next, we mapped the retrieved reads onto the miRNA database using Bowtie2, with the ones exactly matched considered as miRNAs. As for those retrieved reads that fail to map onto the miRNA database, we further annotated them based on Rfam (V12.1, https://rfam.xfam.org/) database using cmsearch (Nawrocki and Eddy, 2013). Finally, the annotation results were integrated as the final annotated miRNA for the studied sample.
       
To determine the expression of each miRNA, we calculated the TPM (UMI count of this miRNA × 1000000/ total UMI count) for each miRNA in each sample. We further applied DESeq2 (Love et al., 2014) to conduct differential expression analysis between low milk production and high milk production groups, with corrected p-value (FDR) ≤0.5 being considered to be significant different.
       
Next, we used miRanda (3.3a) (John et al., 2004) to predict the target genes for the identified miRNA with parameters set as: -en -20 -strict. We also used RNAhybrid (Kruger and Rehmsmeier, 2006) (parameters: - b 100 -c -f 2,8 -m 100000 -v 3 -u 3 -e -20 -p 0.8) to perform this prediction. The target genes supported by either miRanda or RNAhybrid were kept for further analysis. Finally, we performed functional enrichment analysis based on GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) on these target genes using Hypergeometric test, with the term having a corrected p-value (FDR) ≤0.05 defined to be significant different.

We generated approximately 25 million clean tags per sample, with over 98% having a Phred score ≥20, suggesting the high quality of our sequencing data (Table 1). Annotation of these clean tags identified multiple known miRNA and novel miRNA in these samples, ranging from 437 to 808 and from 44 to 999, respectively. RNAhybrid and miRanda identified a total of 20, 201 target genes for all of these miRNAs. Differential expression analysis on these miRNAs identified 388 miRNA up-regulated and 208 down-regulated in low milk production group relative to high milk production group for Murrah buffalos and 128 up-regulated and 1,106 down-regulated in low milk production group for Nili-Ravi buffalos.
 

       
To find out the common features in groups with low milk production, we first extracted the genes significantly up-regulated by miRNA in low milk production groups for both Murrah and Nili-Ravi buffalos. GO enrichment analysis identified ten significantly enriched GO terms, including intrinsic to membrane, integral to membrane, proteinaceous extracellular matrix, plasma membrane, membrane part, cell periphery, membrane, extracellular matrix part, extracellular space, intrinsic to plasma membrane (Table 2). KEGG enrichment analysis detected 12 significantly enriched pathways, including protein digestion and absorption, amoebiasis, AGE-RAGE signaling pathway in diabetic complications, ECM-receptor interaction, PI3K-Akt signaling pathway, platelet activation, focal adhesion, cytokine-cytokine receptor interaction, other types of O-glycan biosynthesis, cocaine addiction, vibrio cholerae infection and small cell lung cancer (Table 3). Protein digestion and absorption pathway exhibited the most significant difference among all pathways (FDR =1.75E-24).
 

 

       
Next, we analyzed he genes significantly down-regulated by miRNA in low milk production groups for both buffalo breeds. These genes were significantly enriched on only two GO terms (plasma membrane, cell periphery) and eight KEGG pathways (Protein digestion and absorption, Amoebiasis, ECM-receptor interaction, Notch signaling pathway, Focal adhesion, AGE-RAGE signaling pathway in diabetic complications, Platelet activation, Insulin secretion) (Table 4). Of particular note is that the pathway showing the most significant difference is still protein digestion and absorption pathway (FDR=9.38E-14), the same as the findings mentioned above (Table 3 and Table 4).
 

       
We were interested in looking at the protein digestion and absorption pathway in detail. Protein digestion and absorption pathway mainly describe the processes of protein hydrolysis in stomach and absorption of multiple amino acids in the small intestine and distal colon. When mapping all of the genes targeted by down-regulated miRNAs in low milk production groups, we found that these genes cover over 80% of the genes in this pathway (Fig 1). On the other hand, analysis of the genes targeted by up-regulated miRNAs revealed that less than 50% of the genes in the studied pathway were covered (Fig 1). Particularly, all of these genes overlapped with the ones targeted by down-regulated miRNAs (Fig 1). This observation could be explained as these down-regulated or up-regulated miRNAs targeted on different isoforms in the same gene locus. Under such circumstance, the effect on such gene expression was difficult to determine, which warrants further investigation in the future (Fig 1). Since genes with down-regulated miRNAs have a higher coverage on the Protein digestion and absorption pathway compared with those with up- regulated miRNAs and considering the unknown effects on genes targeted by both down- and up-regulated miRNAs, we focused on the genes targeted by exclusively down-regulated miRNAs in this pathway. In this study, we identified nine such genes in this pathway (Fig 1).
 

       
Given that miRNA usually represses the expression of its target gene (Cannell et al., 2008), the down-regulation of miRNA implies that the suppression effects on the expression of the target gene are relieved. The expression of the gene regulated by down-regulated miRNA should be higher than that regulated by up-regulated miRNA. One of the genes regulated exclusively by down-regulated miRNA is a gene encoding Pepsin, expressed in the gastric chief cells of the stomach lining. Pepsin is the principal acid protease of the stomach, serving to break down proteins into smaller peptides (Tang and Lin, 2004). Another noteworthy gene down-regulated by miRNA is SLC15A1 (Solute Carrier Family 15 Member 1), highly expressed in the small intestine and encoding a protein that mediates the uptake of small peptides from the lumen into the enterocytes (Frazier et al., 2008). It has been reported that exosomal microRNAs in milk from mothers can be taken up by infants’ (e.g. intestinal cells) (Liao et al., 2017), implying that exosomal microRNA can potentially regulate the expression of genes from infants. In this respect, individuals with low milk production may down-regulated these milk exosomal microRNAs, which can relieve the suppression on the expression of pepsin-encoding gene and SLC15A1 in their infants to enhance the efficiency of protein digestion and peptide absorption.
       
Apart from the two abovementioned genes, SLC1A1 (Solute Carrier Family 1 Member 1), which encodes a member of the high-affinity glutamate transporters that play an essential role in transporting glutamate across plasma membranes (Arnold et al., 2006), was found to be exclusively regulated by down-regulated microRNAs. Glutamate is an important neurotransmitter present in over 90% of all brain synapses in the central nervous system, thus representing the most abundant excitatory neurotransmitter in the vertebrate nervous system (Meldrum, 2000). Moreover, glutamate is also involved in cognitive functions such as learning and memory in the brain (McEntee and Crook, 1993). From this perspective, the down-regulated microRNAs targeting SLC1A1 could increase the expression level of SLC1A1 to increase the glutamate transport efficiency, thus maintaining the development of the nervous system in infants.
       
Interestingly, we found that genes encoding transporters of neutral and basic amino acids were included in the gene set in which the genes were regulated by microRNAs down-regulated in low productive milk. For example, the SLC1A5 (Solute Carrier Family 1 Member 5) gene encodes a sodium-dependent amino acids transporter that accepts all neutral amino acids, including glutamine, asparagine branched-chain and aromatic amino acids as substrates (Garaeva et al., 2018; Kekuda et al., 1996). The SLC38A2 (Solute Carrier Family 38 Member 2) gene, encoding a sodium-dependent amino acid transporter, mediates Na⁺-coupled cellular uptake of small neutral α-amino acids (AAs) and is extensively regulated in response to humoral and nutritional cues with a stoichiometry of 1:1 (Hoffmann et al., 2018). It may function in the transport of amino acids at the blood-brain barrier. The SLC7A8 (Solute Carrier Family 7 Member 8) gene encodes a sodium-independent, high-affinity transport, with small and large neutral amino acids such as alanine, serine, threonine, cysteine, phenylalanine, tyrosine, leucine, arginine and tryptophan as substrates (Zaragozá, 2020). The SLC3A1 (Solute Carrier Family 3 Member 1) gene encodes a type II membrane glycoprotein, one of the components of the renal amino acid transporter, transporting neutral and basic amino acids in the renal tubule and intestinal tract (Kahya et al., 2021). Conclusively, individuals with low milk production may down-regulate the expression of microRNAs targeting genes encoding transporters for neutral and basic amino acids for better absorption of amino acids in the offspring.

In conclusion, buffalo individuals compensate low milk production with upregulation genes involving in the protein digestion and absorption pathway through downregulation miRNA production and consequently increase the efficiency of protein digestion and absorption of amino acids in infants.

The study was funded by the National Natural Science Foundation of China (No.31660649), Guangxi Science and Technology Major Project (GKAA16450002) and Guangxi Natural Science Fund program (2017GXNSFAA198108).

None.