Goat milk is an important dairy product in its own right as well as in processed forms, especially cheese, and is the most common alternative to cow's milk for those with cow's milk allergy (Wu, Reference Wu2022). There are significant differences in milk production between individual dairy goats kept under the same conditions, mainly influenced by genetic factors, physiological factors, feeding methods and environmental factors (Chang, Reference Chang2014). However, significant protein differences associated with lactation are one of the most important factors affecting the lactation performance in dairy goats. Therefore, determining differentially expressed proteins and key pathways of milk synthesis and secretion in goats with differing milk yield potential is of great production value and practical importance for genetic improvement programmes.
Proteomics is the primary technological method for understanding and analyzing protein components, initially proposed by the Australian scientist Williams (Wasinger et al., Reference Wasinger, Cordwell, Cerpa-Poljak, Yan and Humphery-Smith1995). Proteomics has been applied in cow milk research since the 1990s, and major proteins and their genetic variants have been identified, marking a significant advance (Ren et al., Reference Ren, Yang, Peng, Wang, Zhang, Hui and Jia2021). Quantitative proteomics is a new frontier in proteomics research (Schulze and Usadel, Reference Schulze and Usadel2010). Tandem mass tag (TMT) is an important research technology that utilizes chemical markers to simultaneously identify and quantify proteins in different samples. Shen et al. (Reference Shen, Zhang, Shen, Su, Yu, Cao and Zong2023) used the TMT proteomic method to identify target proteins and molecular pathways effective in the treatment of clinical mastitis in dairy cows using anemoside (an anti-inflammatory saponin). With the improvement of detection methods, milk from various animals at different stages of lactation has been studied to explore the proteins associated with lactation and their potential roles. These animals include humans (Zhu and Dingess, Reference Zhu and Dingess2019; Li et al., Reference Li, Li and Zheng2020), donkeys (Cunsolo et al., Reference Cunsolo, Saletti, Muccilli and Foti2010; Li et al., Reference Li, Yu, Chen, Abdlla, Liu, Song, Zhang, Zhang, Yue and Li2021b), camels (Han et al., Reference Han, Zhang, Luo, Ni, Bansal and Zhou2021) and cows (Coscia et al., Reference Coscia, Orrù, Di Nicola, Giuliani, Varalda, Peila, Fabris, Conti and Bertino2012; Li et al., Reference Li, Li, Song, Liu, Zhang, Zheng and Yue2021a). However, there have been no relevant reports on the milk proteome of Guanzhong dairy goats with different milk yields during the peak lactation period. Therefore, this study identified the key differentially expressed proteins and pathways of Guanzhong dairy goats with different milk production levels during the peak lactation period by TMT and bioinformatics analysis, explored the differential protein regulation mechanisms and provided relevant biomarkers and a theoretical basis for improving milk production and selecting high-yielding dairy goats.
Materials and methods
Ethics statement
All animals in this study were raised in accordance with the provisions of Notice No. 5 of the Ministry of Agriculture (PRC). The animals and experimental protocols were authorized by the Institutional Animal Care and Use Committee of Northwest A&F University (Yangling, Shaanxi, China).
Experimental animals and sample collection
A total of 18 healthy Guanzhong dairy goats in peak lactation were sourced from a goat farm in Gaohanchuan, Longxian County, Baoji City, Shaanxi Province. They were fed under the same conditions, but divided into two groups each of 9 goats with either higher (H) or lower (L) milk yield (average milk production of 3.5 ± 0.17 and 1.2 ± 0.25 kg/d, respectively). The two groups were housed in separate sections of the same building. Feeding and management followed conventional methods, with an average daily feed intake of 3.00 ± 0.5 kg/d, and ad libitum water. A single 2 ml sample of milk was collected from each goat, three of which were pooled to create three replicates milk samples for each group. Sampling was at 60 d in milk for both groups. The samples were then frozen with liquid nitrogen for 3–4 h in preparation for subsequent experiments.
Extraction and digestion of proteins
A total of 200 μl of the frozen sample and 1200 μl of the protein precipitated in acetone were placed in a 1.5 ml centrifuge tube and stored overnight at −40 °C. Samples were centrifuged at 12 000 rpm at 4 °C for 10 min and the supernatant was discarded. The precipitate was dried at room temperature for 5 min, dissolved in the sample cracking solution and dissolved at room temperature for 3 h. Finally, it was centrifuged twice at 12 000 rpm for 10 min. The supernatant was stored at −80 °C and used for further TMT analysis.
A total of 50 μg of protein from each sample was diluted with lysate to the same concentration and volume. Dithiothreitol was added to a final concentration of 5 mmol/l, the solution was mixed and left to reduce at 55 °C for 30 min. Subsequently, the mixture was cooled on ice until it reached room temperature. Then, 50 μg of iodoacetamide was added to a final concentration of 10 mmol/l, the solution was mixed thoroughly and alkylated in the dark at room temperature for 15 min. Next, 300 μg of acetone was added to precipitate the protein, and the sample was left overnight at −20 °C before centrifugation at 8000 g for 10 min at 4 °C, collection of the precipitate and evaporation of the acetone for 2–3 min. Then 100 μl of triethylammonium bicarbonate (TEAB: 200 mmol/l) was added for redissolving and precipitating, before adding 1/50 of the sample mass of 1 mg/ml trypsin – N-p-Tosyl-L-phenylalanine chloromethyl ketone and leaving overnight at 37 °C for digestion. Finally, the enzymatically hydrolyzed sample was freeze-dried and stored at −80 °C.
TMT labeling and reverse-phase chromatography separation
To the freeze dried sample we added 50 μl of TEAB buffer (100 mmol/l, Sigma, Germany), mixed it by vortexing, and labeled it in a 1.5 ml centrifuge tube. After warming to room temperature we added 88 μl of anhydrous acetonitrile, vortexed for 5 min and then centrifuged. 41 μl of TMT reagent was added to the sample, mixed well by vortexing, and placed at room temperature for one hour. Finally 8 μl of 5% hydroxylamine was added and left for 15 min to terminate the reaction, after which the mixture was freeze-dried and stored at −80 °C.
The samples were introduced into a C18 column (Agilent, USA) at a flow rate of 300 μl/min for gradient elution with acetonitrile, and the pH was adjusted to 10 using ammonia. The eluents were collected into centrifuge tubes numbered 1–15 at one-minute intervals, and the samples were collected repeatedly over a total duration of 8–60 min. The collected eluent was subjected to vacuum freeze-drying and stored for mass spectrometry.
Chromatographic and mass spectrometry and trusted protein analysis
The sample was loaded onto a C18 chromatography column at a flow rate of 300 nl/min and then separated by an AccumPepMapRSLC analytical column, 75 μm × 50 cm (RP-C18, ThermoFisher) with the following conditions: Mobile phase A: H2O-FA (99.9:0.1, v/v); Mobile phase B: ACN-H2O-FA (80:19.9:0.1, v/v/v); Gradient elution conditions: 0–50 min, 2–28% B; 50–60 min, 28–42% B; 60–65 min, 42–90% B; 65–75 min, 90% B.
The first-level MS (ThermoFisher, USA) has a quality resolution of 60 000, an automatic gain control value of 3e6, and a maximum injection time of 50 ms. The full-scan mass-to-charge ratio (m/z) range that was collected was 350–1500, and MS/MS scanning was performed on the 20 highest peaks. The activation type of all MS/MS spectra was high-energy collision fragmentation with a collision energy of 32 eV. The resolution of MS/MS was 45 000, the automatic gain control value was 2e5, the maximum ion injection time was 80 ms, and the dynamic exclusion time was 30 s. Data were analyzed using ProteomeDiscover2.4. We utilized a database search to retrieve original data based on ScoreSequestHT > 0 and Uniquepeptide ≥ 1, and removed the blank value to screen trusted proteins. Finally, we performed statistical analysis and visualization of trusted protein results.
Bioinformatics analysis
The identified milk proteins were stratified and clustered using MetaboAnalyst 5.0 software. We used the Gene Ontology (GO) database (http://www.geneontology.org) to annotate differentially expressed proteins and analyzed the annotation function of differentially expressed proteins with the online software OmicShare. We used the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg/) to annotate the pathways of differentially expressed proteins and performed pathway analysis on the identified proteins using OmicShare online software.
Results
Trusted protein analysis
To clarify the grouping and repeatability of H and L groups, principal component analysis (PCA) was performed using the expression levels of trusted proteins. In Fig. 1, the first principal component was 32.82% (PCA1), and the second principal component was 21.60% (PCA2). The trusted proteins in the sample exhibited significant separability within the H and L groups, distinguishable based on the axis of the first principal component. A certain degree of clustering was observed among some samples, indicating a level of similarity in the proteins in goat milk.
Trusted protein analysis. PCA analysis of H and L groups. Each dot in the figure represents a repetition in a group experiment, and different colors represent different groups.
Protein quantification and differential analysis
In order to identify the differentially expressed proteins in groups H and L, TMT technology was used for proteomic analysis and the quantitative data were visualized using a volcano plot with −log10(P) and log2(FC). In Fig. 2, the results of the volcano plot show that a total of 1184 proteins were identified in the H and L groups. Under the condition of FC = 1.2 and P < 0.05, a total of 78 differentially expressed proteins were identified. In our study, 50 proteins were significantly upregulated and 28 proteins were significantly downregulated in the H group compared with the L group. Among them, the expression of HK3, ANXA1, HSPB1 and others was significantly increased, whereas the expression of LALBA, XDH and others was significantly decreased. The 15 most differentially expressed proteins up- and down-regulated according to their P-values are given in online Supplementary Tables S1 and S2.
Differences in protein expression between groups H and L. (a) Differential expression protein volcano plot. Blue represents downregulated differentially expressed proteins, while red represents upregulated differentially expressed proteins. (b) Differential expression protein clustering heat map. Red represents proteins with high expression, while blue represents proteins with low expression.
To verify the rationality and accuracy of the identified differentially expressed proteins, hierarchical clustering analysis was employed to aggregate the proteins based on the trend of expression levels. In Fig. 2b, the expression patterns of three samples from group H and group L were similar, showing two main clusters of differentially expressed proteins and demonstrating good repeatability among the samples. This indicated the reliability of the experimental data.
GO term enrichment analysis of differentially expressed proteins
To define and describe the functions of differentially expressed proteins, GO terms were used to enrich and analyze the differentially expressed proteins in the H and L groups, and they were classified based on biological processes, cellular composition and molecular functions. We selected GO entries with differentially expressed protein numbers greater than one in each of the three categories. Each entry was sorted in descending order by -log10(P), and the top 10 entries for each category were displayed. In Fig. 3, GO terms in biological processes mainly focused on carbohydrate metabolic processes, cell adhesion, complement activation (classical pathway) and positive regulation of angiogenesis. The GO terms in cellular components were mainly concentrated in cytoplasm, extracellular space, extracellular regions, cytoskeleton and membrane taxa. In molecular functions, GO terms mainly focused on calcium ion binding, identical protein binding and serine-type endopeptidase activity.
GO annotation of differentially expressed proteins. The x-axis represents the top 30 enriched GO terms, and the y-axis shows the number and percentage of matched proteins
KEGG pathway enrichment analysis of differentially expressed proteins
To systematically analyze the relationship between differentially expressed proteins and milk production, KEGG pathway enrichment analysis was performed for these differentially expressed proteins. The results showed that the differentially expressed proteins in groups H and L were enriched in 110 pathways. The top 20 KEGG pathways are shown in Fig. 4a. Differentially expressed proteins were mainly enriched in the pathways related to glycolysis/gluconeogenesis, purine metabolism, galactose metabolism and regulation of actin cytoskeleton. The annotation results of the differentially expressed proteins identified by KEGG were classified according to the pathway types in KEGG, as shown in Fig. 4b. Among the top 20 enriched pathways, 2 were related to environmental information processing, 3 to human diseases and cellular processes, 4 to biological processes and 8 to metabolism.
Functional enrichment analysis. (a) KEGG enrichment analysis of differentially expressed proteins in groups H and L. The x-axis represents the enrichment score, while the y-axis indicates the pathway information. The size of the bubbles is positively correlated with the number of differential proteins contained in the entry, and the smaller the P-value, the more significant the difference. (b) Differential expression gene KEGG pathway classification. The x-axis represents the number of differentially expressed proteins annotated to each metabolic pathway, the y-axis represents the pathway name, and different colored columns represent different categories.
Discussion
Low average yields are one of the main problems hindering the high-quality development of the dairy goat industry (Wang and Zhao, Reference Wang and Zhao2021). TMT has become a high-throughput, high-resolution, reproducible and quantitatively accurate method for detecting proteins in different samples and is widely used in research on differential expression analysis of proteins. With the improvement of detection technology, a large number of proteins and their functions related to lactation have been discovered. Differential proteomic analysis of high and low yield Guanzhong goat milk during the peak lactation period can provide a deeper understanding of the regulatory mechanisms behind the differences in goat milk production. Therefore, in this study, the TMT proteomics method was used to determine the key differentially expressed proteins and associated pathways that affect the milk yield of Guanzhong dairy goats during the peak lactation period.
A total of 78 proteins were differentially expressed. Compared with the L group, HK3, HSPB1, and ANXA2 were significantly upregulated in the H group. Glycolysis is a process that converts glucose to pyruvic acid and generates a small amount of ATP and NADH (Chandel, Reference Chandel2021). HK3 plays a role as a key enzyme in the glycolysis process by catalyzing the conversion of glucose to glucose-6-phosphate (Aaron et al., Reference Aaron, Edgard, Sen, Affan, Grant and Michel2020), promoting the glycolysis process and increasing the ATP content of the body. ATP transports chemical energy within the cell for metabolism and provides energy for the milk production process (Malloy et al., Reference Malloy, Sherry and Jeffrey1988). Therefore, it can be speculated that the increased milk production in the high-yield group of dairy goats may be associated with HK3 in the glycolytic pathway. HSPB1 is involved in protein folding as well as degradation of proteasome and autolysosome as a molecular chaperone (Carroll et al., Reference Carroll, Erin, Haven, Tamara, Jessica, Deborah, Damian, Maruf, Tali, James, Alfred, Habil and Karl2023). MAPK is an important regulatory factor that regulates signal transduction cascade responses and plays a crucial role in various aspects of cell physiology, mammary development and milk synthesis (Liu et al., Reference Liu, Hou, Zhang, Seleh-Zo, Wang, Cao and An2020; Li et al., Reference Li, Yuan, Jiang, Zhang, Su, Lv, Sang, Chen, Feng and Chen2022a). P38 MAPK increases mRNA stability by phosphorylating AU-rich element binding proteins, thereby having a positive impact on protein synthesis (Li et al., Reference Li, Cai, Zhou, Yin, Zhang, Loor, Sun, Zhang, Liu and Zhang2016). In our study, HSPB1 was enriched in the MAPK signaling pathway associated with milk protein synthesis. Based on these findings, we speculate that HSPB1 may play a role in promoting the formation of the cytoskeleton and milk secretion in mammary epithelial cells.
It is well known that milk production is related to the activity, proliferation and apoptosis of mammary epithelial cells. Increasing epithelial cell number by improving cell proliferation and decreasing apoptosis contributes to an increase in milk production. In our study, the expression of differential protein HSPB1 was significantly upregulated in the high-yield group. Research has found that HSPB1 is involved in signal transduction and regulation during cell proliferation, differentiation, and apoptosis (Tedesco et al., Reference Tedesco, Cristofani, Ferrari, Cozzi, Rusmini, Casarotto, Chierichetti, Mina, Galbiati, Piccolella, Crippa and Poletti2022). Silencing HSPB1 significantly inhibited cell proliferation and promoted cell apoptosis (Chen et al., Reference Chen, Hou, Xu, Peng, He, Liu, Wang, Huang and Liu2022). HSPB1 is an important member of the heat-shock protein superfamily. It inhibits the caspase cascade response by interfering with the expression of the pro-apoptotic proteins Bax and Bid and separating cytochrome C released by mitochondria, thus playing an anti-apoptotic role at multiple levels (Bruey et al., Reference Bruey, Ducasse, Bonniaud, Ravagnan, Susin, Diaz-Latoud, Gurbuxani, Arrigo, Kroemer and Solary2000). Therefore, we speculate that the increased milk production may be associated with the upregulation of the HSPB1 protein.
ANXA2 is a calcium-dependent phospholipid-binding protein of the annexin family involved in membrane fusion, vesicular transport, cell adhesion, cell proliferation, apoptosis and other activities of exocytosis (Zhang et al., Reference Zhang, Ye, Yang, He, Wang and Zhao2012; Li et al., Reference Li, Wang, Huang, Zhao, Chen and Zhang2022b). ANXA2 is overexpressed in various human malignant tumors. When ANXA2 is knocked down, cell division and proliferation are inhibited and cell apoptosis increases, indicating its anti-apoptotic effect (Lin et al., Reference Lin, Song, He and Dai2016). Research has shown that ANXA2 regulates cell proliferation and apoptosis via the mTOR, SREBP-1c and cyclin D1 signaling pathways, promoting milk synthesis and mammary epithelial cell proliferation (Zhang et al., Reference Zhang, Zhen, Ao, Yuan and Gao2018). Different forms of annexin are differentially regulated in the maintenance of cell number, exocytosis and milk secretion during lactation (Janjanam et al., Reference Janjanam, Singh, Jena, Varshney, Kola, Kumar, Kaushik, Grover, Dang, Mukesh, Prakash and Mohanty2014). ANXA5 is upregulated during the peak and late lactation phases. ANXA3 is downregulated during peak lactation and upregulated in early and late phases. In our study, the expression of annexin A2 was significantly increased in H goats relative to L. therefore, we speculate that ANXA2 may serve as a useful biomarker of milk production in Guanzhong dairy goats.
We observed that proteins such as XDH and LALBA were significantly downregulated in the H group compared with the L group. XDH is the key enzyme responsible for purine degradation in purine metabolism and plays a role in promoting the production of reactive oxygen species (ROS), cell apoptosis and hypoxanthine catabolism (Bortolotti et al., Reference Bortolotti, Polito, Battelli and Bolognesi2021). The accumulation of ROS can lead to an increase in the ratio of Bax/Bcl-2 and promote the activation of caspase-3, leading to cell apoptosis (Wang et al., Reference Wang, Liu, Tian, Du, Lin, Xiao, Wang, Chen, Geng, Subramanian, Niu, Wang and Yue2020). Cancer cells are very sensitive to the rapid increase in intracellular levels of ROS, and XDH activation may play a role in anti-cancer therapy (Xu et al., Reference Xu, Li, Mozziconacci, Zhu, Xu, Tang, Chen, Huang, Holzbeierlein, Schöneich, Huang and Li2019). The XDH protein in the purine metabolism pathway was considerably downregulated in the high-yield group when compared to the low-yield group. We speculate that the increase in milk production may be associated with the downregulation of the XDH protein.
As a mammary-specific protein, LALBA is a subunit of lactose synthase that catalyzes the formation of lactose and thus regulates milk supply (Ostrowska et al., Reference Ostrowska, Zwierzchowski, Brzozowska, Kawecka-Grochocka, Żelazowska and Bagnicka2023), since lactose is an important osmotic regulator of milk volume (Rahmatalla et al., Reference Rahmatalla, Arends, Ahmed, Reissmann and Brockmann2020). For this reason lactose content at the point of secretion is more or less constant, and changes in milk lactose concentration across the course of lactation (typically a decline in late lactation) need to be interpreted with caution since they probably reflect paracellular leakage of lactose into plasma. For example, in their study of lactation-related genes in Holstein cows, Fan et al. (Reference Fan, Arbab, Zhang, Yang and Yang2021) found a cross-lactational positive correlation between the lactose content in milk and milk production as yield and lactose concentration both declined after peak lactation. Interestingly, in our study, LALBA was found to be significantly downregulated in the high-yield group and upregulated in the low-yield group and significantly enriched in the galactose metabolic pathway. This contrasts with expectations based on numerous studies of mechanistic aspects of milk secretion. The straightforward explanation that LALBA may play an important role in reducing milk production in dairy goats is unlikely to be true, and it is more likely that expression is generally in excess of what is needed to maintain optimum lactose synthesis. However, this needs further research.
In conclusion, proteomic analysis of whey protein from high-yielding (H) and low-yielding (L) Guanzhong dairy goats during the peak lactation period was performed using TMT technology. A total of 78 differentially expressed proteins were identified. Compared with L, 50 proteins such as HK3, HSPB1, and ANXA2 were significantly upregulated in H milk, whilst 28 proteins such as LALBA, XDH and RRAS2 were significantly downregulated. Bioinformatics analysis of the differentially expressed proteins showed that galactose metabolism, purine metabolism, glycolysis/gluconeogenesis, MAPK signaling pathway, regulation of actin cytoskeleton and other pathways were closely associated with the milk yield of Guanzhong dairy goats. Proteins such as HK3, HSPB1, ANXA2, LALBA, XDH and RRAS2 were important candidate proteins associated with the milk production characteristics in dairy goats. These data provide relevant biomarkers and a theoretical basis for improving milk production in Guanzhong dairy goats.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0022029924000013
Acknowledgments
We acknowledge the OE Bioinformatics Institute (OE Biotech, ShangHai, China) for TMT-based quantitative proteomics sequencing. This work was financially supported by the Shaanxi Provincial Agricultural Collaborative Innovation and Extension Alliance Project [LMZD202002], the Shaanxi Provincial Agricultural Science and Technology Innovation Drive Project [NYKJ-2021-ST-03] and the Shaanxi Provincial Technology Innovation Guidance Project [2022QFY11-01]
