Identification and bioinformatics analysis of differentially expressed milk exosomal microRNAs in milk exosomes of heat-stressed Holstein cows

In summer, heat stress is one of the primary reasons for the compromised health and low milk productivity of dairy cows. Hyperthermia affects milk synthesis and secretion in the mammary glands of dairy cows. As molecules for intercellular communication, milk-derived exosomes carry genetic material, proteins, and lipids, playing a crucial role in mammary tissue growth and milk synthesis in dairy cows. The aim of this study was to explore the milk exosomal miRNA profile of heat-stressed and normal Holstein cows. We isolated and identified milk exosomes to screening for differentially expressed miRNAs using small RNA sequencing. Then, TargetScan and miRanda algorithms were used to predict the putative targets of the differentially expressed miRNAs, whereas GO and KEGG pathway enrichment analyses were performed for the differentially expressed miRNA-target genes. Our results showed that 215 miRNAs were significantly differentially expressed in heat-stressed milk exosomes, of which one was upregulated and 214 were significantly downregulated. GO and KEGG enrichment analyses indicated that differentially expressed miRNAs might play a role in apoptosis, autophagy, and the p38 MAPK pathway. qRT-PCR assay verified that the expression of miRNAs was consistent with the sequencing results, warranting further verification of their specific targets of action. In conclusion, changes in the miRNA expression profile of milk exosomes indicated the role of exosomal miRNAs in regulating heat stress resistance and apoptosis in dairy cows. Our results suggested that milk-derived exosomal miRNAs could increase mammary gland resistance to heat stress, thereby enhancing milk synthesis in dairy cows.


Background
Heat stress (HS) refers to a metabolic disorder in animals caused by an increase in ambient temperature, which causes damage to animal tissues and organs (Cai et al. 2018). In the livestock industry, HS is a key concern because it seriously influences not only reproductive performance but also milk yield (Sakai et al. 2020). Bovine mammary epithelial cells (BMECs) are the secretory cells of the mammary glands, which play an important role in milk synthesis. HS reduces the ability of BMECs to synthesize and secrete milk, further affecting the quality and yield of milk. Milk contains a large number of functional substances, including proteins, lipids, and exosomes, which are vital nutrients needed for a wide array biological activities (Chen et al. 2020b). Milk synthesis in dairy cows is often affected by many environmental, physical, and chemical factors. It has been reported that during an inflammatory response, milk-derived exosomes are mediated by the TLR4/NF-κB and p53 pathways to reduce BMEC apoptosis and protect LPS-induced intestinal epithelial cells against apoptosis (Xie et al. 2019). Porcine milk components, including exosomes, are associated with mammary gland function (Zhang et al. 2018). Zhang et al. showed that exosomal proteins are closely related to the biosynthesis of milk by BMECs (Zhang et al. 2018). It is generally believed that RNA is degraded by RNases in the environment. However, milk exosomes can prevent miRNA degradation (Yu et al. 2017). Exosomes can be directly absorbed in the small intestine upon gastric acid stimulation, conferring a protective effect on the immune system and the health of neonates (Munir et al. 2020). Moreover, exosomes can be used as carriers to deliver various macromolecules into target cells without triggering an immune response (Benmoussa et al. 2016).
Exosomes (Exo) are tiny membranous vesicles secreted by cells. They are composed of a lipid bilayer membrane with a vesicle diameter of 20-150 nm and enriched in proteins, lipids, and nucleic acids; in addition, they act as signaling molecules, regulating many types of cellular functions (Keller et al. 2006). After their entry via endocytosis or membrane receptors, exosomes release proteins, mRNAs, miRNAs, and lncRNAs into the target cells, thereby changing the gene expression profile of these cells (Yu et al. 2019). miRNAs control gene expression by degrading mRNA and inhibiting translation (Li et al. 2020a, b). Exosomes derived from mesenchymal stem cells can protect β cells from hypoxia-induced apoptosis, reduce cellular stress, and inhibit the activation of the p38 MAPK signaling pathway (Chen et al. 2020a). miR-100-5p in human umbilical cord MSC-derived exosomes protects cardiomyocytes from hypoxic injury and inhibits inflammasome activation-mediated pyroptosis (Liang et al. 2020). Qiu et al. have shown that exosomes carrying miR-181a-3p play a critical role in apoptosis by alleviating endoplasmic reticulum stress . In addition, the study showed that bovine milk exosomes enhanced the viability and antioxidant capacity of intestinal cryogenic epithelial cells, and exosome pretreatment increases intracellular miR-146a and miR-155 levels (Wang et al. 2021). When the integrity of BMECs was compromised, the health and milk yield of dairy cows can be highly affected. Previous studies have shown that miRNAs are not only related to apoptosis but involved in the physiological processes of breast development and milk synthesis (Ji et al. 2017). However, whether the differential expression of exosomal miRNA in milk is associated with milk synthesis and quality in Holstein cow is unclear.
In the present study, we aimed to identify differentially expressed miRNAs in exosomes derived from heat-stressed Holstein cow milk using small RNA sequencing technology and to predict the target gene and possible pathways associated with stress and injury. Our results suggested that the expression profile of milk exosomal miRNAs could regulate the development and function of mammary tissue in heat-stressed dairy cows via the MAPK pathway. Our findings serve as a new reference for studies on alleviating HS in cows during summer.

Sample collection
There is continuous monitoring of temperature and humidity index (THI) of cowshed, rectal temperature, and respiratory rate of China Holstein cows under the same management level in summer (August, heat stress group (HS)) and winter (December, control group (Ctr)) in Sihong ranch of Jiangsu province, China. Fifty milliliters of milk from 6 Holstein cows with similar lactation periods and parity in the HS group and Ctr group was collected, transported to the laboratory at low temperature, and stored at − 80℃ for subsequent sequencing. The data of THI, rectal temperature, and respiratory rate are shown in Table 1.

Extraction and identification of exosomes
The frozen milk sample was heated in a water bath at 37℃. After dissolution, the sample was added to a new centrifuge tube and centrifuged at 4℃ for 30 min at 2000 g. The supernatant was transferred to a new centrifuge tube at 10,000 g at 4℃ for 45 min. The supernatant was centrifuged at 100,000 g for 120 min at 4℃ in a new centrifuge tube. The supernatant was discarded and resuspended with 20 mL precooled 1 × PBS, centrifuged at 2000 g for 30 min at 4℃. Supernatant was centrifuged at 4℃ of 2000 g for 30 min. The centrifugation was repeated at 2000 g at 4℃ for 30 min. Then, the supernatant was centrifuged at 100,000 g at 4℃ for 120 min. The supernatant was discarded, and the precipitated pellet was resuspended and precipitated with 1-mL precooled PBS, and temporarily stored at 4℃. Forty percent, 20%, 10%, and 5% iodoxanol were prepared. Iodoxanol of different concentrations was added along the tube wall (3.6 mL for each) according to the high to low concentration. Finally, 1 mL was resuspended, temporarily stored at 4℃, and added to the top layer at 4℃, 100,000 g for 120 min. After centrifugation, it was divided into 12 layers, the liquid in the middle 6-9 layers was taken out, and centrifuged again at 4℃ for 100,000 × g for 120 min. The supernatant was removed, and the exosomes were resuspended with 1 × PBS precooled with 100μL stored at − 80℃.

Transmission electron microscopy
Ten microliters of exosomes was absorbed and added to the copper net precipitation for 1 min, the floating solution was removed by filter paper, and the phosphotungstic acid was added to the copper net precipitation for 1 min. The floating solution was removed by filter paper and dried for several minutes at room temperature. The transmission electron microscopy imaging was performed at 100 kv.

Nanoparticle tracking analysis
The samples were thawed in a water bath at 25℃ and placed on ice. Exosome samples were diluted with 1 × PBS and directly used for NTA detection.

Western blotting
The exosomes were centrifuged at 10,000 g at room temperature for 5 min, then supernatant was discarded, 50μL RIPA lysate was added, and the exosomes were incubated on ice for 30 min. The concentration of exosomes was determined by BCA kit and denaturated in metal bath at 100℃ for 10 min; and then SDS-PAGE was performed. The protein was transferred to PVDF membrane by wet transfer method and was sealed with 5% skim milk for 1 h. Primary antibodies CD81, CD63, TSG101, and Calnexin (1:1000) polyclonal antibodies were added and incubated overnight at 4℃. Washing was done with 1×TBST buffer for 3 times, 5 min each time, then HRP-labeled secondary antibody was added, and incubation at room temperature was done for 90 min. After washing with 1 × TBST buffer for 3 times, 5 min each time, ECL chemiluminescence was developed and scanned in the infrared imaging system.

Small RNA sequencing and analysis
Total RNA was extracted from exosomes, and the Small RNA sequencing library was prepared by Truseq Small RNA Sample Prep Kits (Illumina, San Diego, USA). In simple terms, the total RNA was added with 3′ and 5′ adapter, and the Small RNA connected by the adapter was used for reverse transcription PCR to create the DNA library, and the primer was annealed to connect the tow fringes adapter for PCR, and the Small RNA library was purified by gel. Finally, the High Sensitivity DNA Chip was used to test the sample library. The qualified libraries were sequenced using Illumina Hiseq2000/2500 with a single-ended read length of 150 bp.

Screening of differentially expressed miRNAs
Raw reads were subjected to ACGT101-miR (LC Sciences, Houston, TX, USA) to remove adapter dimers, junk, low complexity, common RNA families (rRNA, tRNA, snRNA, snoRNA), and repeats. The final data is the valid reads, which can be used for subsequent analysis of small RNA data. According to the data of miRNA expression in each sample, T test was used to analyze the differential expression of miRNAs, and according to the significance of expression difference (P < 0.05), differentially conserved miRNAs were screened out.

Target gene prediction and GO and KEGG enrichment analysis
miRNAs bind to target sites mainly through complementary pairing. The 3′ UTR sequence of bovine mRNA was used as the target sequence to predict the target genes of differentially expressed miRNA sequences. We used TargetScan (v5.0) and miRanda (v3.3a) software to predict the target genes of miRNAs with significant differences. The target genes were screened predicted with two software according to the scoring criteria. Target genes with the context score percentile less than 50 were removed in the TargetScan algorithm, and target genes with the maximum free energy (Max Energy) > − 10 were removed in the miranda algorithm (that is, the threshold was targetScan_score ≥ 50, miranda_Energy < − 10). Finally, the intersection of software was selected as the final target gene of the differential miRNAs.
The GO functional significance enrichment analysis first mapped all the differentially expressed genes to each term in the Gene Ontology database (ftp:// ftp. ncbi. nih. gov/ gene/ DATA/ gene2 go. gz), calculated the number of genes in each term, and then applied hypergeometric test to find the GO items significantly enriched in differentially expressed genes compared with the entire genome background. KEGG is the main public database on pathways. Pathway significance enrichment analysis uses KEGG pathway (http:// www. genome. jp/ kegg) as the unit, and hypergeometric tests are used to identify the pathways that are significantly enriched in differentially expressed genes compared with the overall genome background.

Quantitative real-time PCR (qRT-PCR)
Total RNA was extracted with TRIzol™ reagent (Invitrogen, cat: 15,596,026, USA). Before exosome RNA extraction, synthetic cel-miR-39 (GenePharma, Shanghai) was added as an internal control (Title et al. 2015). The density of sample RNA was quantified spectrophotometrically at 260/280 nm. miRNA 1st Strand cDNA Synthesis Kit (by stem-loop) (Vazyme, cat: MR101-01/02, Nanjing) was used to reverse transcribe the total RNA, and mRNA expression was quantified with real-time PCR. Using comparative Ct (2 −ΔΔCt ) value method standardizes the expression levels of all target genes with the endogenous reference gene U6 or Cel-miR-39 (forward primer: ATA TCA TCT CAC CGG GTG TAA ATC; reverse primer: TAT GGT TTT GAC GAC TGT GTGAT). The primer sequences are listed in supplementary material Table S1.

Statistical analysis
Dates are shown as mean values ± SEM. All results were analyzed by T test, using the software GraphPad Prism 8.0.1 (La Jolla, CA, USA). P < 0.05 was considered statistically significant.

Identification of milk exosomes
Exosomes originate from the multi-vesicles in living cells and are secreted as membranous vesicles, with a density of 8.76E + 9 particles/mL and a typical "cup and plate" shape ( Fig. 1A, C-D). TEM revealed that milk exosomes negatively stained with phosphotungstic acid had a diameter of approximately 100 nm. NTA showed that the average particle size of exosomes was 73.82 nm, and the size of exosomes in the range of 30-150 nm accounted for 97.8% of all particles (Fig. 1E-F). CD81, CD63, and TSG101 proteins are the signature proteins of exosomes, and their presence further proves that the vesicles observed using an electron microscope are exosomes. Calnexin is a protein on the endoplasmic reticulum, which acts as a molecular chaperone to participate in the folding and processing of new peptide chains of proteins. Our Western blotting results showed that proteins CD81, CD63, and TSG101, except calnexin, were expressed in exosomes compared to breast epithelial cells, indicating that the extracted exosomes were pure and not contaminated (Fig. 1B).

Sequencing quality analysis and miRNA length distribution
To detect the expression of miRNAs in milk exosomes of cows under HS and normal conditions, the singleended 1 × 50 bp sequencing reading length of Illumina HiSeq2000/2500 platform was used for exosome detection. A total of six exosome samples were detected, including three each from the HS group (HS1, HS2, and HS3) and non-HS group (Ctr1, Ctr2, and Ctr3). The total read length of the original sequencing was 133,884,666, and the average read length of each sample was 22,314,111. Initially, in the original sequencing data, the 3′ connector sequence was removed, and the sequences with base length less than 18 nt were removed. Eighty percent of the sequences were A, C, G, or T; 3 N (not necessarily continuous); only A and C without G or T; or only G and T without A or C, and continuous nucleotide dimers, and trimers were removed (Table 2).
Based on the original sequencing data, we conducted length distribution statistics on the filtered data. Most of the sequences were 20-24 nt in length, indicating the typical characteristics of a Dicer enzyme cleavage ( Fig. 2A). We compared the measured sequences to mRNA, RFam (including rRNA, tRNA, snRNA, and snoRNA), and the Repbase database to yield sequencing data that can be considered data (Table 3). Subsequently, miRNA comparison, identification, and prediction analyses were performed on the valid sequencing data, and a total of 42,159,649 valid reads were identified. Further analysis showed that 80.67% of all the screened sequences were miRNAs, while 6.56%, 10.43%, 0.21%, 0.32%, and 1.80% of the total sequencing data belonged to rRNA, tRNA, snoRNA, snRNA, and other RFam RNAs (Fig. 2B). In addition, we compared the differentially expressed miRNAs on chromosomes, and the circle plot is plotted according to differentially expressed miRNAs, indicating the distribution of differentially expressed miR-NAs on chromosomes (Fig. 2C).

Differential expression analysis of miRNAs
Differential expression of 215 miRNAs was detected in the screened sequencing data (P < 0.05) (Table S2), among which 208 were known and seven were novel miRNAs, compared to the precursors of the selected species in miR-Base (Fig. 3A, Tables S3-S4). Compared to the Ctr groups,  the expression of one miRNA was significantly upregulated, while 214 miRNAs were significantly downregulated (Fig. 3B, Table S2). In addition, significantly different expression of 215 miRNAs was predicted for 12,446 target genes (Table S5). bta-miR-15a had the highest number of target genes. Among them, the most significantly enriched target gene with the most number of differentially miRNAs was UBE2W (66 miRNAs).

GO and KEGG enrichment analyses for differentially expressed miRNAs
We performed GO analysis to identify the biological processes associated with the identified miRNAs (Fig. 4A). These differentially expressed miRNA target genes were significantly enriched in 1432 GO pathways (P < 0.05), among which 910 GO pathways were associated with biological processes, 245 were cellular components, and 277 were molecular functions (Table S6). Moreover, the target genes were significantly   (Table S7). Among the top 20 significantly enriched KEGG pathways, the pathway for cancer (P = 2.02E − 09) was the most significantly enriched (Fig. 4B).

Bioinformatics analysis of differentially expressed miRNAs
Nine of the top fifty miRNAs were differentially expressed with the highest degree of enrichment in Ctr and HS groups. These nine miRNAs were closely related to autophagy, apoptosis, and the MAPK signaling pathways (Fig. 5A). To examine the association between differentially expressed milk exosomal miRNAs and their target genes, miRNAgene interaction networks were generated using Cytoscapses (Fig. 5B, Table S8). Network analysis revealed results consistent with KEGG pathway prediction for the nine miRNAs, including the MAPK signaling pathway.

Discussion
HS not only damages the health of cows but also affects milk synthesis. Milk is rich in various bioactive substances, including lipids, proteins, and exosomes, which are crucial for biological activities. Milk exosomal proteins and miR-NAs play a key role in regulating cell proliferation, inflammatory responses, immune regulation, and cancer cell proliferation (Gonzalez-Rivas et al. 2020;Kim et al. 2020;Zheng et al. 2014). Based on these functions of milk-derived exosomes, we sought to differentially expressed miRNAs in exosomes derived from the mil of heat-stressed Holstein cows using small RNA sequencing technology, and we predicted the target genes and possible enriched pathways involved in stress response and apoptosis.
The mechanism by which mRNA, miRNAs, and lncRNAs in milk-derived exosomes participate in the physiological processes of cells is still unclear. miRNAs play an important role in cellular communication as signal transduction molecules (van Wijnen et al. 2013). Smythies et al. showed that stress can change the content of exosomal miRNA profiles of developing sperms in the epididymis (Smythies et al. 2014).
In this study, we isolated and validated the exosomes derived from the milk of heat-stressed cows. High-throughput sequencing analysis showed that 80.67% of all screened sequences were miRNAs. We detected 215 differentially expressed miRNAs in the screened sequencing data, including one miRNA that was significantly upregulated, and 214 miRNAs, which were significantly downregulated. Therefore, milk exosomal miRNAs have different expression profiles in the cows of the control and HS groups. Differential miRNA expression is often accompanied by changes in cellular function (Fuziwara and Kimura 2017). GO analysis identified the biological processes associated with the differentially expressed miRNA target genes enriched in 1432 pathways (P < 0.05), among which 910 GO pathways were associated with biological processes, 245 were cellular components, and 277 were molecular functions. KEGG pathway Screening of differentially expressed miRNAs and their target genes. A Different miRNA, apoptosis, autophagy, and MAPK signaling pathways of the Venn. B Network analysis of differentially expressed miRNAs and their target genes in milk. Green rhombus shown in network indicate miRNAs whereas red rhombus their target genes enrichment analysis showed that the differentially expressed miRNAs mainly play a role in the MAPK signaling pathway, endotosis, human T-cell leukemia virus 1 infection, the Rap1 signaling pathway, and rheumatoid arthritis, which have been shown to be closely related to apoptosis (Yu et al. 2019;Zhang et al. 2018). Therefore, we speculated that differential expression of milk exosomal miRNAs could play an important regulatory in the injury response of heat-stress dairy cows.
Compared to the Ctr group, the top 9 differentially expressed miRNAs including bta-miR-339a_R-2, bta-miR-339b, chi-miR-340-5p, bta-miR-141_R-1, bta-miR-30b-5p, bta-miR-30f_R-3, bta-miR-27a-3p, and bta-miR-193a-5p in the HS group were significantly downregulated, and bta-miR-193a-5p was significantly upregulated. Previous studies have shown that miR-27a-3p via BTG1 regulates apoptosis in the colon cancer cell line, HCT-116, via the ERK/MEK pathway (Su et al. 2019). Li also demonstrated that miR-27a-3p regulates cellular inflammation and autophagy . Function prediction analysis of the target genes of differentially expressed miRNAs in this study indicated that the miRNAs mainly regulated the internal environment of breast tissues through the MAPK signaling pathway. Therefore, further in-depth investigation of milk-derived exosomes is warranted to elucidate the association between miRNA expression profile and milk synthesis microenvironment under HS conditions.

Conclusions
In conclusion, HS alters the miRNA expression profile of milk exosomes and the differentially expressed miRNAs may play an important role via the MAPK signaling pathway, thereby affecting the survival and lactating function of BMECs in heat-stressed Holstein cows. Therefore, it is crucial to further verify the target genes of milk exosomal miRNAs and associated signaling pathways to provide a theoretical basis for the molecular mechanism by which milk exosomal miRNAs regulate the lactating and physiological cows.