Identication and Bioinformatics Analysis of Differentially Expressed Milk Exosomal Micrornas in Milk Exosomes of Heat-Stressed Holstein Cows

and identied 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 signicantly differentially expressed in heat-stressed milk exosomes, of which one was up-regulated and 214 were signicantly downregulated. GO and KEGG enrichment analyses indicated that differentially expressed miRNAs might play a in apoptosis, autophagy, and the p38 MAPK pathway. qRT-PCR assay veried that the expression of was consistent with the sequencing results, warranting further verication of their of


Background
Heat stress (HS) refers to a metabolic disorders in animals caused by an increase in ambient temperature, which causes damage to animal tissues and organs [1]. In the livestock industry, HS is a key concern because it seriously in uences not only reproductive performance but also milk yield [2]. 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. Our previous study showed that HS triggers apoptosis in BMECs and causes changes in their morphology, structure, and function [3][4][5]. Although the mechanism of HSinduced injury of BMECs has been lucidated, the role of exosomal miRNAs in the HS-regulating mechanism of BMECs is still unclear. Therefore, this study aimed to investigate the miRNA pro les of milk exosomes of heat-stressed Holstein cows.
Milk contains a large number of functional substances, including proteins, lipids, and exosomes, which are vital nutrients needed for a wide array biological activities [6]. Milk synthesis in dairy cows is often affected by many environmental, physical, and chemical factors. It has been reported that during an in ammatory 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 [7]. Porcine milk components, including exosomes, are associated with mammary gland function [8]. Zhang et al. showed that exosomal proteins are closely related to the biosynthesis of milk by BMECs [9]. It is generally believed that RNA is degraded by RNases in the environment. However, milk exosomes can prevent miRNA degradation [10]. 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 [11].
Moreover, exosomes can be used as carriers to deliver various macromolecules into target cells without triggering an immune response [12].
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 [13]. After their entry via endocytosis or membrane receptors, exosomes release proteins, mRNAs, miRNAs, and lncRNAs into the target cells, thereby changing the gene expression pro le of these cells [14]. miRNAs control gene expression by degrading mRNA and inhibiting translation [15; 16]. 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 [17]. miR-100-5p in human umbilical cord MSC-derived exosomes protects cardiomyocytes from hypoxic injury and inhibits in ammasome activation-mediated pyroptosis [18]. Qiu et al. have shown that exosomes carrying miR-181a-3p plays a critical role in apoptosis by alleviating endoplasmic reticulum stress [19]. 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 [20]. 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 [21]. In line with this, Hou et al. reported that 5, 10methylenetrahydrofolate reductase modulates goat milk synthesis by regulating the binding activity of miR-1266 and miR-616 [22].
In the present study, we aimed to identify signi cantly 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 pro le of milk exosomal miRNAs could regulate the development and function of mammary tissue in heat-stressed dairy cows via the MAPK pathway. Our ndings serve as a new reference for studies on alleviating HS in cows during summer.

Sample collection
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. 50mL milk from 6 Holstein cows with similar lactation periods and parity in the HS group and Ctr group were 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 identi cation 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 30min at 2,000g. The supernatant was transferred to a new centrifuge tube at 10,000g at 4℃ for 45min. The supernatant was centrifuged at 100,000g for 120min at 4℃ in a new centrifuge tube. The supernatant was discarded and resuspended with 20mL precooled 1×PBS, centrifuged at 2,000g for 30min at 4℃. Supernatant was centrifuged at 4℃ of 2,000g for 30min. The centrifugation was repeated at 2,000g at 4℃ for 30min. Then the supernatant was centrifuged at 100,000g at 4℃ for 120min. The supernatant was discarded, resuspended and precipitated with 1mL precooled PBS, and temporarily stored at 4℃. 40%, 20%, 10% and 5% iodoxanol were prepared. Iodoxanol of different concentrations was added along the tube wall (3.6mL for each) according to the high to low concentration. Finally, 1mL resuspend temporarily stored at 4℃ was added to the top layer at 4℃, 100,000g 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 resuspend with 1×PBS precooled with 100μL stored at -80℃.
Transmission electron microscopy 10μL of exosomes was absorbed and added to the copper net precipitation for 1 min, the oating solution was removed by lter paper, and the phosphotungstic acid was added to the copper net precipitation for 1min. The oating solution was removed by lter 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,000g at room temperature for 5min, then supernatant was discarded, 50μL RIPA lysate was added, and the exosomes were incubated on ice for 30min. The concentration of exosomes was determined by BCA kit, and denaturated in metal bath at 100℃ for 10min; 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 1h. Primary antibodies CD81, CD63, TSG101 and Calnexin (1:1000) polyclonal antibodies were added and incubated overnight at 4℃. Washing with 1ÍTBST buffer for 3 times, 5min each time, then add HRP-labeled secondary antibody and incubate at room temperature for 90min. After washing with 1×TBST buffer for 3 times, 5min 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 puri ed by gel. Finally, the High Sensitivity DNA Chip was used to test the sample library. The quali ed libraries were sequenced using Illumina Hiseq2000/2500 with a single-ended read length of 1Í50bp.

Screening of differentially expressed miRNAs
Raw reads were subjected to ACGT101-miR (LC Sciences, Houston, Texas, USA) to remove adapter dimers, junk, low complexity, common RNA families (rRNA, tRNA, snRNA, snoRNA) and repeats. The nal data is the valid reads, which can be used for subsequent analysis of small RNA data. According to the data of miRNAs expression in each sample, T test was used to analyze the differential expression of miRNAs, and according to the signi cance 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 miRNAs sequences. We used TargetScan (v5.0) and miRanda (v3.3a) software to predict the target genes of miRNAs with signi cant 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 nal target gene of the differential miRNAs.
The GO functional signi cance enrichment analysis rst mapped all the differentially expressed genes to each term in the Gene Ontology database (ftp://ftp.ncbi.nih.gov/gene/DATA/gene2go.gz), calculated the number of genes in each term, and then applied hypergeometric test to nd the GO items signi cantly enriched in differentially expressed genes compared with the entire genome background. KEGG is the main public database on Pathways. Pathway signi cance enrichment analysis uses KEGG Pathway (http://www.genome.jp/kegg) as the unit, and hypergeometric tests are used to identify the pathways that are signi cantly enriched in differentially expressed genes compared with the overall genome background.

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 signi cant.

Identi cation 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-150nm accounted for 97.8% of the 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 miRNAs length distribution
To detect the expression of miRNAs in milk exosomes of cows under HS and normal conditions, the single-ended 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 133884666, and the average read length of each sample was 22314111. 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; 3N (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 ltered 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, identi cation, and prediction analyses were performed on the valid sequencing data, and a total of 42159649 valid reads were identi ed. 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 generated circle polt showed that miRNAs were evenly distributed on each chromosome (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 miRBase (Fig. 3A, Tables S3-S4). Compared to the Ctr groups, the expression of one miRNA was signi cantly upregulated, while that 214 miRNAs were signi cantly downregulated (Fig. 3B, Table S2). In addition, signi cantly 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 signi cantly 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 identi ed miRNAs (Fig. 4A). These differentially expressed miRNA target genes were signi cantly enriched in 1,432 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 signi cantly enriched in 188 KEGG pathways (Table S7). Among the top 20 signi cantly enriched KEGG pathways, the pathway for cancer (P = 2.02E-09) was the most signi cantly enriched (Fig. 4B).

Bioinformatics analysis of differentially expressed miRNAs
Nine of the top fty miRNAs were differentially expressed with the highest degree of enrichment in in the two 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, miRNA-gene 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 miRNAs play a key role in regulating cell proliferation, in ammatory responses, immune regulation, and cancer cell proliferation [24][25][26]. Based on these functions of milk-derived exosomes, we sought to differentially expressed miRNAs in exosomes derived from the mil of heatstressed 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 [27]. Smythies et al. showed that stress can change the content of exosomal miRNA pro les of developing sperms in the epididymis [28]. 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 signi cantly upregulated, and 214 miRNAs, which were signi cantly downregulated. Therefore, milk exosomal miRNAs have different expression pro les in the cows of the control and HS groups. Differential miRNA expression is often accompanied by changes in cellular function [29]. GO analysis identi ed the biological processes associated with the differentially expressed miRNA target genes enriched in 1,432 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 enrichment analysis 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 [30; 31]. 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 10 differentially expressed of miRNAs included 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 signi cantly downregulated, and bta-miR-193a-5p, which was signi cantly upregulated. Previous studies have shown that miR-27a-3p via BTG1 to regulate apoptosis in the colon cancer cell line, HCT-116, via the ERK/MEK pathway [32]. Li also demonstrated that miR-27a-3p regulates cellular in ammation and autophagy [33]. 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 pro le and milk synthesis microenvironment under HS conditions.

Conclusions
In conclusion, HS alters the miRNA expression pro le 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 exosomel 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.

Con ict of Interest
The authors declare no con ict of interest. Authors' contributions Y.W. conceptualized the study, participated in its design and research. J.F. and H.F.Z carried out the molecular studies and sample collection. K.L.C and H.X.L analyzed data and drafted the manuscript. All authors read and approved the nal manuscript for publication.

Group
Average Rectal temperature /℃  Figure 1 Exosomes identi cation. A, Electron micrograph shows the morphology of exosomes from milk. B, Western blotting were used to analyze the expression of CD63, CD81, TSG101, Calnexin. C and D, The concentration of exosomes from milk. E and F, The particle size of exosomes from milk.  Identi cation of signi cantly differentially expressed miRNAs in exosomes. ND means not detected.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.