Comparative Proteomics Study of Yak Milk from Standard and Naturally Extended Lactation Using iTRAQ Technique

Simple Summary To elucidate the differences in milk protein compositions and mammary gland functions between yaks of standard lactation (TL yaks) and prolonged lactation (HL yaks), iTRAQ technique was used to compare the skim milk proteins in the two yak groups. A total of 202 differentially expressed proteins (DEPs) were revealed, among which 109 proteins were up-regulated and 93 were down-regulated in the milk of HL yaks compared to TL yaks. The bioinformatics analysis revealed that the differences in skim milk protein between the HL yaks and the TL yaks suggests that the mammary gland of the HL yak is at a degeneration stage. Abstract Extended lactation is a common phenomenon in lactating yaks under grazing and natural reproduction conditions. To elucidate differences in milk protein compositions and mammary gland functions between yaks of standard lactation (TL yaks) and prolonged lactation (HL yaks), whole milk samples of TL yaks and HL yaks (n = 15 each) were collected from a yak pasture at the northwest highland of China. The iTRAQ technique was used to compare the skim milk proteins in the two yak groups. A total of 202 differentially expressed proteins (DEPs) were revealed, among which 109 proteins were up-regulated and 93 were down-regulated in the milk of HL yaks compared to TL yaks. Caseins including κ-casein, αs1-casein, αs2-casein, and β-casein were up-regulated in HL yak milk over 1.43-fold. The GO function annotation analysis showed that HL yaks produced milk with characteristics of milk at the degeneration stage, similar to that of dairy cows. KEGG enrichment showed that the metabolic pathways with the most differences are those that involve carbohydrate metabolism and the biosynthesis of amino acids. The present results highlight detailed differences in skim milk proteins produced by HL yaks and TL yaks and suggest that the mammary gland of HL yak is at the degeneration stage.


Introduction
The yak is one of the main unique livestock species in ethnic minority areas of the Qinghai-Tibet Plateau in China [1]. It provides important living materials for local farmers, with milk being one of the major products. There exist some differences in lactation between yaks and dairy cows. The lactation of dairy cows can be extended beyond the standard 305-d lactation through several manipulated strategies. Moreover, these cows are subject to a substantial decline in milk yield and their milk composition changes during the extended lactation stage [2]. Since the feeding and management mode in yaks is natural grazing and reproduction, they are completely dependent on natural grass resources and natural reproduction. Furthermore, considerable proportions of lactating yaks (about 15%) will experience a prolonged lactation stage if they are not pregnant during the calving year [3]. If these yaks are calved in the spring and milked for consumption, can maintain lactation in the following winter by calf suckling only, and can be milked again during the next Milk samples of yaks were collected in August. The experimental yaks were raised on a pasture about 3500 m above sea level in Hongyuan county, Sichuan Province, China. The experimental yaks included TL yaks (n = 15) and HL yaks (n = 15). The TL yaks calved during the spring season (March to May), while the HL yaks calved one year earlier. Their milk secretion was maintained by calf suckling during the winter. All the experimental yaks were 4 to 7 years old and 2 to 4 parities. The experimental yaks and their calves were grazed on the same natural grassland during the daytime and were separated from their calves at night. The lactating yaks were milked by the hands of local farmers in the morning. Approximately 50 mL of whole milk were collected from each yak, which was transferred to the laboratory using dry ice and stored at −80 • C until analysis. All animal care and milking procedures were approved by The Animal Ethics Committee of Southwest Minzu University (No. swun20200138).

iTRAQ Analysis of Skim Milk Proteins of Yaks
The whole milk of each yak was centrifuged at 800× g and 4 • C for 20 min to prepare skim milk. The pooled skim milk samples of TL yaks and HL yaks were prepared by mixing equal volumes of 15 skim milk samples of corresponding yaks, respectively. The two pooled Animals 2022, 12, 391 3 of 15 samples (2 mL each) were transported using dry ice to Shenzhen BGI Technology Co., Ltd. for iTRAQ analysis.
For the iTRAQ assay, the pooled skim milk samples were centrifuged at 25,000× g for 20 min to remove residual fat and cell debris. The supernatant skim milk was removed and mixed with 5 volumes of cold acetone and stored at −20 • C overnight. The mixture was centrifuged again, and the resulting pellet was used for further preparing a protein solution [14]. A total of 100 µg of protein from this solution was digested with Trypsin Gold (protein: trypsin = 20:1) at 37 • C for 12 h. The resulting peptides were labeled using the iTRAQ Reagent 8-plex Kit according to the manufacturer's protocol, followed by fractionation using a Shimadzu LC-20AB HPLC equipped with a 4.6 mm × 250 mm Gemini C18 column (Phenomenex) [15]. The eluted peptides were pooled as 20 fractions and were then vacuum-dried, dissolved, and loaded on an LC-20AD nano HPLC (Shimadzu, Kyoto, Japan) equipped with a 2 cm C18 trap column. Then, the peptides were eluted into an 18 cm analytical C18 column. Mass spectrometry analysis was performed as described in previous studies [13]. Data was acquired using a TripleTOF 5600 System fitted with a Nanospray III source (AB SCIEX, Downtown Redwood City, America).

Bioinformatics Analysis
IQuant software was applied to the quantification of proteins. Proteins with a 1.2-fold change and a Q-value of less than 0.05 were determined as differentially expressed proteins, and they must be defined in at least 1 replicate experiment. All proteins with a false discovery rate (FDR) of less than 1% proceeded with the following analysis, including Gene Ontology (GO), Clusters of Orthologous Groups (COG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway.

Identification of Skim Milk Proteins of Yaks
iTRAQ analysis identified 767 proteins in the skim milk samples of TL and HL yaks based on the 1930 unique peptides obtained (Supplemental Material File S1). The protein coverage was between 0.001 to 0.999, of which 37.1% was identified as using at least two unique peptides. The length of most acquired peptides ranged between 7 and 17 amino acids ( Figure 1A), and the molecular weights of approximately 80% of the identified proteins were less than 100 kDa ( Figure 1B).

GO and COG Annotation of Identified Proteins
All identified proteins in the skim milk samples of TL and HL yaks were subjected to the GO analysis and were categorized into biological processes, cellular components, and molecular functions, based on their GO annotations ( Figure 2, Supplemental Material File S2). The proteins identified were classified by molecular functions and were enriched in 13 different functional terms. The majority were related to binding and catalytic activity (471 and 276 proteins in total, respectively), followed by structural, molecular activity, transporter activity, and enzyme regulator activity. The GO cellular location classifications of the proteins were involved in 16 categories, and the most obvious differences were in the cell, cell part, organelle, and organelle part. There were 23 biological process categories in which the identified proteins were involved, and the largest number of proteins were observed in the cellular process (475), followed by the metabolic process (388), the single-organism process (385), and biological regulation (299).

GO and COG Annotation of Identified Proteins
All identified proteins in the skim milk samples of TL and HL yaks were subjected to the GO analysis and were categorized into biological processes, cellular components, and molecular functions, based on their GO annotations ( Figure 2, Supplemental Material File S2). The proteins identified were classified by molecular functions and were enriched in 13 different functional terms. The majority were related to binding and catalytic activity (471 and 276 proteins in total, respectively), followed by structural, molecular activity, transporter activity, and enzyme regulator activity. The GO cellular location classifications of the proteins were involved in 16 categories, and the most obvious differences were in the cell, cell part, organelle, and organelle part. There were 23 biological process categories in which the identified proteins were involved, and the largest number of proteins were observed in the cellular process (475), followed by the metabolic process (388), the single-organism process (385), and biological regulation (299).
COG annotation of all the identified proteins revealed 24 functional categories ( Figure 3, Supplemental Material File S3). Among these COG categories, the skim milk samples of TL and HL yaks were highly enriched in several major functional COG categories, including post-translational modification, protein turnover, chaperones, energy production and conversion, cytoskeleton, carbohydrate transport, and metabolism. COG annotation of all the identified proteins revealed 24 functional categories (Figure 3, Supplemental Material File S3). Among these COG categories, the skim milk samples of TL and HL yaks were highly enriched in several major functional COG categories, including post-translational modification, protein turnover, chaperones, energy production and conversion, cytoskeleton, carbohydrate transport, and metabolism.

GO Enrichment Analysis of DEPs
GO enrichment analysis of the 202 DEPs between TL and HL yaks demonstrated that 13, 19, and 22 protein categories were highly enriched in the cellular component, molecular function, and biological process categories with p < 0.05 or p < 0.01, respectively ( Figure 4). Among these, the vesicle, cytoplasmic vesicle, and membrane-bounded vesicle were the most abundant categories in the cellular component ( Figure 4A, Supplemental Material File S6). The transporter activity, enzyme regulator activity, and enzyme inhibitor activity were the most abundant categories in the molecular function ( Figure 4B, Supplemental Material File S6). In the biological process, as many as 65 GO terms were enriched with p < 0.05, and the 22 major categories (p < 0.01) were related to oxidation-reduction processes and the generation of precursor metabolites and energy ( Figure 4C, Supplemental Material File S6). 4). Among these, the vesicle, cytoplasmic vesicle, and membrane-bounded vesicle were the most abundant categories in the cellular component ( Figure 4A, Supplemental Material File S6). The transporter activity, enzyme regulator activity, and enzyme inhibitor activity were the most abundant categories in the molecular function ( Figure 4B, Supplemental Material File S6). In the biological process, as many as 65 GO terms were enriched with p < 0.05, and the 22 major categories (p < 0.01) were related to oxidation-reduction processes and the generation of precursor metabolites and energy ( Figure 4C, Supplemental Material File S6).

Pathway Enrichment Analysis of DEPs
The 202 DEPs were used for pathway enrichment analysis, and 25 major pathways were highly enriched via KEGG with p < 0.05 (Supplemental Material File S7). In addition, a scatter plot for the top 20 of KEGG enrichment results is shown in Figure 5. The KEGG pathways that the DEPs mainly participated in were: primary immunodeficiency, staphylococcus aureus infection, cytokine-cytokine receptor interaction, and dilated cardiomyopathy. The metabolic pathways with the most differences involved carbohydrate metabolism, the biosynthesis of amino acids, and fat digestion and absorption.

Pathway Enrichment Analysis of DEPs
The 202 DEPs were used for pathway enrichment analysis, and 25 major pathways were highly enriched via KEGG with p < 0.05 (Supplemental Material File S7). In addition, a scatter plot for the top 20 of KEGG enrichment results is shown in Figure 5. The KEGG pathways that the DEPs mainly participated in were: primary immunodeficiency, staphylococcus aureus infection, cytokine-cytokine receptor interaction, and dilated cardiomyopathy. The metabolic pathways with the most differences involved carbohydrate metabolism, the biosynthesis of amino acids, and fat digestion and absorption. Animals 2022, 12, x FOR PEER REVIEW 12 of 15

Discussion
Our previous research discovered some differences in milk between TL yaks and HL yaks [6]. For example, the contents of protein and fat and the activities of several enzymes, such as γ-glutamyltransferase and alkaline phosphatase, in HL yak milk were significantly higher compared to the milk of TL yaks. It has been reported that the yield and composition of milk are influenced by many factors, including human maternal age, time of delivery and maternal diet, and the stage of lactation, which was the most influential factor [17]. The detailed composition and quality of milk from extended lactation have been studied in cows and humans. One study found higher protein and fat concentrations, an unaffected casein to protein ratio, and protein composition of the bovine milk from the extended lactation [18]. Czosnykowska-Łukacka et al. reported that the concentration of carbohydrates in mother's milk showed a negative correlation with lactation of about two years, while fat and protein concentrations were opposite. Moreover, during prolonged lactation in humans (over 18 months), it was found that the concentration of carbohydrates significantly decreases, and fat and protein concentration significantly increases [19]. The analysis of the composition of prolonged yak milk allows one to assess the nutritional value of milk. In this study, the milk protein profiles of TL yak and HL yak were intensively studied based on the iTRAQ technique. The expressions of some major milk proteins in HL yaks were increased, such as κ-casein, αS1-casein, αS2-casein, and β-casein. In addition, vitamin D-binding protein, retinol-binding protein 4, lactotransferrin, and serotransferrin levels in HL yak milk also increased (Supplemental Material File S4). Proteins are the major nutrients of milk and have many functions except for providing proteins for nutrition purposes. The micronutrients in milk can also affect its function [20]. The present results indicate that HL yaks can provide more nutrients in milk compared to TL yaks.
Based on their GO functional annotations, all identified proteins in the skim milk samples of TL and HL yaks were classified in accordance with molecular function, cellular

Discussion
Our previous research discovered some differences in milk between TL yaks and HL yaks [6]. For example, the contents of protein and fat and the activities of several enzymes, such as γ-glutamyltransferase and alkaline phosphatase, in HL yak milk were significantly higher compared to the milk of TL yaks. It has been reported that the yield and composition of milk are influenced by many factors, including human maternal age, time of delivery and maternal diet, and the stage of lactation, which was the most influential factor [17]. The detailed composition and quality of milk from extended lactation have been studied in cows and humans. One study found higher protein and fat concentrations, an unaffected casein to protein ratio, and protein composition of the bovine milk from the extended lactation [18]. Czosnykowska-Łukacka et al. reported that the concentration of carbohydrates in mother's milk showed a negative correlation with lactation of about two years, while fat and protein concentrations were opposite. Moreover, during prolonged lactation in humans (over 18 months), it was found that the concentration of carbohydrates significantly decreases, and fat and protein concentration significantly increases [19]. The analysis of the composition of prolonged yak milk allows one to assess the nutritional value of milk. In this study, the milk protein profiles of TL yak and HL yak were intensively studied based on the iTRAQ technique. The expressions of some major milk proteins in HL yaks were increased, such as κ-casein, αS1-casein, αS2-casein, and β-casein. In addition, vitamin D-binding protein, retinol-binding protein 4, lactotransferrin, and serotransferrin levels in HL yak milk also increased (Supplemental Material File S4). Proteins are the major nutrients of milk and have many functions except for providing proteins for nutrition purposes. The micronutrients in milk can also affect its function [20]. The present results indicate that HL yaks can provide more nutrients in milk compared to TL yaks.
Based on their GO functional annotations, all identified proteins in the skim milk samples of TL and HL yaks were classified in accordance with molecular function, cellular localization, and biological pathways. It was reported that 66% of the DEPs identified in whey from yak colostrum and mature milk were found to be related to binding activity [21].
Approximately 44% and 22.4% of the identified proteins in bovine colostrum were involved in catalytic activity and binding activity, respectively [22]. Our present results are basically consistent with these reports, and the catalytic activity and binding activity account for the largest proportion of the identified proteins ( Figure 2).
The mammary gland undergoes morphological and functional changes during development. The lactation cycle of cows includes early lactation, middle lactation, late lactation, and dry lactation [23]. Mammary gland degeneration is a key stage in dry lactation, and during this stage, the ability to synthesize milk decreases and cell apoptosis increases [4]. It was shown that the concentration of lactoferrin in mammary secretions during the dry period was significantly increased, which could be used to measure the degree of mammary degeneration in dairy cows [24]. In this study, the expressions of hemoglobin subunit-beta, lactotransferrin, and serotransferrin in HL yak milk increased significantly compared to TL yak milk. Since these components are blood-derived proteins, it suggests that the permeability of mammary gland tight junctions increase in HL yaks, which is consistent with the characteristics of cows in the degeneration stage. In addition, during mammary gland degeneration, the somatic cell count and the apoptosis rate of mammary epithelial cells increase [25]. This study found that HL yaks contained significantly higher levels of several types of keratins in milk than TL yaks (Supplemental Material File S4), which may be a result of increased shedding of mammary epithelial cells in HL yak milk since keratin is a marker of epithelial cells [26]. Moreover, cysteine-rich secretory protein 3 (CRISP-3) precursor, which is a key protein in cell apoptosis [27], and vinculin, which can maintain cell growth and differentiation and promote cell survival [28], were also up-regulated and down-regulated in HL yak milk, respectively. This suggests that the mammary glands in HL yaks are at a state of degeneration, consistent with the characteristics of mammary glands of degenerated cows.
During the dry period, the ability of mammary epithelial cells to synthesize lactose, milk fat, casein, α-lactalbumin, and β-lactoglobulin decreases, while the concentration of lactoferrin in mammary secretions increases significantly [4,7]. GO enrichment analysis of DEPs showed that the vesicle and cytoplasmic vesicle were the most abundant categories in the cellular component ( Figure 4A). In contrast, the transporter activity was the most abundant category in the molecular function ( Figure 4B). Vesicles play a key role in protein transport [29], and the difference in protein transport activity may be related to milk protein components in HL yak milk. The KEGG pathway enrichment analysis, which was based on the DEPs, revealed that the metabolic pathways with the most differences were those that involved carbohydrate metabolism, the biosynthesis of amino acids, and fat digestion and absorption. Previous research found that a high number of proteins in human and ruminant milk serum were related to metabolic processes [30]. In another study, metabolism-related pathways (such as glycolysis/gluconeogenesis and biosynthesis of amino acids) also enriched many differentially expressed whey proteins in human and bovine colostrum and mature milk [31]. The differences in these metabolic pathways were related to the ability of mammary epithelial cells to synthesize lactose, milk fat, and milk protein in HL yak.
In the sample, the pooled skim milk samples of TL yaks and HL yaks were prepared by mixing equal volumes of 15 skim milk samples of corresponding yaks, respectively. This procedure clearly blurs the individual differences between the tested yak females. However, the lactation period between yak groups was similar. The same sample processing method was used in previous studies [32].

Conclusions
The iTRAQ technique was used to compare skim milk proteins in yaks from standard and naturally extended lactation (TL yaks and HL yaks, respectively). A total of 202 differentially expressed proteins were identified, among which 109 proteins were upregulated and 93 were down-regulated in HL yaks compared to TL yaks. The GO function