Untargeted metabolomics and quantitative label-free proteomics analysis of whole milk protein from Chinese native buffaloes

ABSTRACT This study aimed to analyze metabolite abundances and proteome differences between buffalo milk from Binglangjiang (BBM) and Dehong (DBM). Untargeted UPLC-MS/MS, label-free proteomics approach, and bioinformatics analyses including GO, KEGG, and PPI were performed. The results revealed that catechol, L-proline, γ-aminobutyric acid, hippuric acid, creatinine, and betaine were detected in BBM and DBM as the major metabolites. One hundred thirteen proteins, mainly within the molecular weight range of 20–80 kDa, were identified, and the differential protein numbers of BBM were higher than those of DBM. These proteins were primarily involved in cholesterol metabolism, Staphylococcus aureus and Salmonella infections, the peroxisome proliferator-activated receptor signaling pathway, and coagulation cascades. Relative abundances of α-s1-casein, α-s2-casein, and β-casein were plentiful in BBM, while κ-casein, albumin, and α-lactalbumin were in DBM. PPI network analysis indicated that the identified proteins were highly interconnected. These findings promote a better understanding of quality biomarkers in dairy product improvements.


Introduction
Milk and its by-products are one of the most important protein sources produced by the mammary glands of mammals, providing complete dietary nutrition, bioavailable amino acid sources, and essential micronutrients, as well as readily digestible food for the human diet (Arslan et al., 2021;Davoodi et al., 2016;Leischner et al., 2021;Smith et al., 2022;Vargas-Ramella et al., 2021).At present, the global production of unprocessed milk is dominated by cattle (81%), followed by buffalo (15%), and a combination of goats, sheep, and camels (4%), totaling 902 million tons of solid milk and predicted to grow at 2.2% per annum by 2030 (OECD/ FAO, 2021;van Heerden, 2021).Nonetheless, animal-based milk sources produced by minor dairy animals like buffalo, goat, sheep, camel, yak, horse, and donkey are also nutritious and economically impactful in some countries (Alichanidis et al., 2016;Faccia et al., 2020).In Asia continent, buffalo milk accounts for 35.30% (124,958,493 tons) of the total unprocessed milk production, where the majority is supplied by the five most prominent producers, including India, Pakistan, China, Egypt, and Nepal (Minervino et al., 2020).Across Yunnan province in the Southwestern region of China, recognized as the second biggest population center of buffaloes after Guangxi province, Binglangjiang (river type) and Dehong (swamp type) buffaloes are the most famous native dual-purpose breeds to produce meat and milk (Fan et al., 2019;Qiu et al., 2021;Yang et al., 2007;Zhou et al., 2020).
Previous studies reported that buffalo milk has been widely investigated for a long time in association with the particular use of the food industry because of its technological suitability, the differences in functional properties, and higher energy per unit volume (Arora et al., 2022;Jiang et al., 2022;Mejares et al., 2022;Tsuji et al., 2022).Buffalo milk contains high nutritive values, rich flavor profiles, and is more appropriate for processing to make various by-products, such as condensed milk, cheese, cream, yogurt, butter, and probiotics, as compared to cow milk (Jiang et al., 2022;Stobiecka et al., 2022;Vargas-Ramella et al., 2021).Buffalo milk poses a noticeable difference in physicochemical features relative to other domesticated ruminant species, such as cows, sheep, and goats (Bittante et al., 2022;Verduci et al., 2019;Yuan et al., 2022).It has specific nutrition properties with a higher level of proteins (e.g.α-s1-casein and κ-casein), fatty acids, lactose, vitamin A, calcium, a lower percentage of water, dry matter, vitamin E, riboflavin, larger casein micelles (ranging between 110 and 170 nm, with 2.68-3.72 mL g −1 of voluminosity and 2.60-2.90g water g −1 of solvation), an absence of carotene, and reduced cholesterol in comparison with other ruminant milk (Bertoni et al., 2021;D'Ambrosio et al., 2008;El-Salam & El-Shibiny, 2011;Figliola et al., 2021;Garau et al., 2021;Guo & Hendricks, 2010).The gross chemical composition and particular biological characteristics of buffalo milk may be strongly affected by the differences in metabolic mechanisms due to geographical location, rearing system, feeding regime, and breed (Garau et al., 2021;Han et al., 2007;Sun et al., 2014;Yuan et al., 2022).Therefore, owing to the highnutritive values and metabolites in buffalo milk being beneficial as potential ingredients in various food industry applications, there is a desire to study these biological differences among breeds in the light of new available emerging technologies.
Recent metabolomics and proteomics approaches are recognized as robust methods for the comprehensive investigation of milk and its by-products.Nuclear magnetic resonance (NMR) spectroscopy, gas chromatography-mass spectrometry (GC-MS), and liquid chromatography-mass spectrometry (LC-MS) are considered the most extensively used analytical methods to determine the metabolomes and proteomes in milk.Among those methods, LC-MS has quickly become the preferred analytical method due to its high reproducibility and sensitivity for analyzing an extensive range of chemical compounds in a single detection (Jiang et al., 2022;Lv et al., 2020;Shi et al., 2021;Trobbiani et al., 2017).In this study, we used an untargeted ultraperformance liquid chromatography -tandem mass spectrometry (UPLC-MS/MS) and quantitative label-free highperformance liquid chromatography -tandem mass spectrometry (HPLC-MS/MS)-based proteomics approach and performed gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and protein-protein interaction (PPI) analyses to identify the differentially expressed proteins in BBM and DBM in order to expand the molecular understanding of their metabolomics and proteomics profiles.To the best of our knowledge, this is the first study using label-free mass spectrometry with high mass accuracy in MS/MS scans for buffalo milk from Yunnan native purebreds.

Milk sample collection
Sixty-five buffalo raw-milk samples were obtained from two native Yunnan breeds located on two different dairy farms in Dehong prefecture (Dehong buffalo breed, n = 35) and Tengchong county (Binglangjiang buffalo breed, n = 30), Yunnan province, China.The selection criteria for the buffalos were based on the clinical records, where the body temperature range, the absence of systemic signs, and mastitis status were normal.Fresh milk samples were transferred to 50 mL sterile centrifuge tubes, kept on ice, and immediately transported to the laboratory.All milk samples were separated into three parts and stored at − 80°C for further analysis.The use of all animals in this research was approved according to the Institutional Animal Care and Use Committee of Yunnan Agricultural University guidelines.On the first part of milk samples, routine tests were performed, including chemical composition analysis (fat, protein, lactose, and total solids) using a MilkoScan FT2 analyzer (FOSS Electric A/S Analytics, Hillerd, Denmark); protein quantification using a bicinchoninic acid (BCA) assay kit (Beijing Solarbio Life Science Co., Ltd., China); and A2 β-casein protein isoform concentration measurement using a sandwich enzyme-linked immunosorbent assay (ELISA) kit (Biosensis Pty., Ltd., Australia), following the manufacturer's instructions.Previously, the buffalo DNA samples were genotyped using a sequencing approach, which revealed that all samples carried the A2A2 β-casein genotype and that the milk samples were A2 milk (Prabakusuma et al., 2022).The Mann-Whitney U test was performed in SPSS software v26.0 (IBM Corp., Armonk, NY, U.S.A.) for statistical analysis.The second and third parts of milk samples were prepared for untargeted metabolomics and label-free proteomics approaches, respectively.

Milk sample preparation for untargeted UPLC-MS/ MS analysis
Before conducting untargeted UPLC-MS/MS analysis, the skim milk samples obtained from each breed of buffalo had been thawed at a temperature of 4°C.These samples from each breed of buffalo were then mixed to make a single sample of BBM and DBM.Those mixed samples were prepared with three biological replications for further analysis.An aliquot of 150 µl of each sample was transferred into 2 mL centrifuge tubes.A total of 200 µl methanol and 200 µl methyl tertiary-butyl ether were added into each centrifuge tube and thoroughly vortexed for 60 s.All samples were centrifuged at 4°C for 10 min at 12,000 rpm to separate large biomolecules.Before being subjected to autosampler vials, the supernatant was filtered through a 0.22-μm hydrophobic membrane to obtain the samples for LC-MS analysis.The metabolite from buffalo milk samples was extracted in accordance with several previous studies, with slight modifications (Shi et al., 2021;Tian et al., 2016;Villaseñor et al., 2014).The obtained filtrate in the centrifuge tube was kept at − 80°C until further analysis.To check the quality of the milk metabolite extract, the same amount of BBM and DBM samples was taken.

Untargeted UPLC-MS/MS analysis
An Acquity UPLC separation module with a column heater system (Waters Corp., Milford, MA, U.S.A.) was occupied.Chromatographic separation of milk metabolites was accomplished on a Thermo Ultimate 3000 system equipped with a UPLC ® HSS T3 (150 × 2.1 mm, 1.8 µm) coupled with a Thermo Q Exactive HF-X mass spectrometer, which was installed with an electrospray ionization source (ESI) with a spray voltage of 3.5 kV (ESI + ) and − 2.5 kV (ESI − ).Ten µL of each milk sample was injected into the column at a flow rate of 0.45 mL/min.The data from mass spectrometry were obtained in both positive and negative ion modes.The analyzer scanned over a range of mass-to-charge ratios of 81-1,000 m/z for a full-scan at a mass resolution of R = 60,000.Data-dependent acquisition (DDA) MS/MS was performed with a high collision dissociation scan, and the normalized MS collision energy was 30 eV for precursor ions.

Milk sample preparation for quantitative label-free proteomics analysis
All milk samples were thawed at 4°C, 100 μL of each sample was taken, and 400 μL of 2X sodium deoxycholate (DOC) lysis buffer was added.Sample solutions were further sonicated on ice for 8 min with a 15% amplitude.Protein enzymatic hydrolysis was then performed according to the quantitative results of electrophoresis.Twenty μL of protein lysate for each sample was prepared for digestion.Thirty μL of 50 mM ammonium bicarbonate (ABC) pH 8.0 and 0.5 μL of 1 M dithiothreitol were added to each tube, mixed with a pipette, heated at 100°C for 5 minutes, centrifuged, and cooled to room temperature for denaturation.The protein lysate was transferred into a 10-kDa filter, centrifuged at 14,000 rpm for 10 min at 4°C, and 100 μL of 50 mM ABC was added and centrifuged again to obtain the filtrate.Fifty μL of 50 mM ABC and 0.5 μg of pancreatin enzyme were added to the filtrate in a tube and initially digested at 37°C for 4 h.The digested filtrate was digested again overnight at 37°C with the addition of 1 μg of pancreatin enzyme.The digested peptides were collected the next day, washed twice with 50 mM ABC, mixed with 20 μL of 10% formic acid (FA), precipitated with DOC, and centrifuged at 14,000 rpm for 10 min, and the peptide supernatant was taken for desalting.

HPLC-MS/MS analysis
HPLC-MS/MS analysis was performed with an Easy-nLC 1200 system coupled to an Orbitrap Fusion Lumos mass spectrometer (Thermo Fisher Scientific, Waltham, MA, U.S.A.) to generate peak lists and detect peptide signals from MS data acquisition.In this mass spectrometric analysis, a nanoelectrospray ionization source (nano-ESI) with a spray voltage of 2,200 V and an ion transmission capillary temperature of 320°C was used.The prepared sample, as outlined in section 2.5, was dissolved in 100 μL of mobile phase A, consisting of 100% deionized water with a 0.1% FA aqueous solution for liquid chromatography separation.Data from mass spectrometry experiments used DDA in the positive ion mode method for acquisition.The total ion current chromatogram was obtained after scanning the mass spectrum signal (Figure S1).After first-round peptide-spectrum matching was completed, quantitative analysis was performed using the Sequest tool in the Proteome Discoverer software v2.4 (Thermo Fisher Scientific, Waltham, MA, U.S.A.).A ≤ 1% false discovery rate (FDR) was used at equal peptide and protein levels to increase the confidence and remove false-positive identification.The second round of peptidespectrum filtering eliminated proteins supported by fewer than two distinct peptide identifications, and only identical proteins were conceded and grouped.The Proteome Discoverer software filtered the credible peptide spectrum matches (PSMs) from the search results, which contained at least one PSM with a credibility ≥ 99%.The UniProt database (Taxon ID: 89462) with protein records of 46,754 was used to annotate further matched spectra.

Bioinformatics analysis
The enrichment analysis was conducted using the Bubalus bubalis organism database in UniProtKB.The obtained protein list was exported to Microsoft Excel and subjected to the BioDeep bioinformatics resource to find GO annotations.The GO categories were processed and analyzed using the DAVID bioinformatics tool (https://david.ncifcrf.gov/),a free online resource that provides functional-related genes' interpretation of enriched biological themes derived from genomic studies, and only the genes with a p-adjusted value of FDR of less than 0.05 were involved, and subsequent GO terms were plotted (Huang et al., 2009).GO enrichments were classified using Uniprot based on functional annotation involvement, including biological processes, molecular functions, and cellular components (Götz et al., 2008).The metabolic pathway analysis involved in differentially expressed genes was annotated using the KEGG database (https:// www.genome.jp/kegg/)(Kanehisa et al., 2017).The clustering heat map was generated using the pheatmap package v2.0.3 (https://cran.r-project.org/web/packages/pheatmap/index.htmL).The in silico PPI analysis was performed using the STRING database platform v10.5 (https://string-db.org/),which allows for the critical integration and evaluation of the interactions among proteins, including physical direct and functional indirect associations (Szklarczyk et al., 2021).The results of the PPI analysis were then downloaded in XGMML format and input into Cytoscape56 software v3.9.0 (http:// www.cytoscape.org/) to make a network diagram of the functional PPIs (Otasek et al., 2019).

Routine analysis of buffalo milk composition
Before subjecting BBM and DBM to an untargeted metabolic analysis and a label-free proteomics approach, we analyzed the difference between the two milk types in the main nutritional composition, and the results are shown in Table 1.The concentrations of milk fat, lactose, total solids, and A2 β-casein protein isoform in BBM and DBM were not significantly different.However, only the protein content of BBM (4.25% ± 0.06) was significantly higher than that of DBM (4.18% ± 0.04) (p < .05)and slightly lower than that of the Murrah and Nili-Ravi buffalo breeds studied in Guangxi province (4.92% and 4.54%, respectively) (Zhou et al., 2018).Using a genomic sequencing approach for single nucleotide polymorphism (SNP) genotyping, all buffalo DNA samples that corresponded with the milk samples observed in this study were determined to carry the A2A2 β-casein genotype only and all milk samples were categorized as A2 β-casein milk variant.It indicated that the mutation Pro 67 /His 67 in the amino acid sequence located at exon 7 of chromosome 6 in the CNS2 gene encoded β-casein milk phenotype was absent in the buffaloes.This secondary result was also consistent with the findings reported by Pineda et al. (2019) and de Oliveira et al. ( 2021), who genotyped several buffalo breeds of water and swamp types in the Philippines and Brazil, respectively.Our present study has demonstrated the reliability of quantifying A2 β-casein protein isoform concentration in fresh buffalo milk samples through the use of the sandwich ELISA.This method used a chicken anti-bovine A2 β-casein detection antibody, a pre-coated rabbit anti-bovine A2 β-casein polyclonal capture antibody, and a horseradish peroxidase-conjugated secondary donkey anti-chicken IgY antibody for signal amplification.The results revealed that the average concentration of A2 β-casein protein isoform in DBM was higher than that in BBM (Table 1).The optical density values of the two milk types detected by a sandwich ELISA assay were associated with the positive interaction between the A2 β-casein protein isoform and the enzyme-labelled antibodies.

Buffalo milk metabolite profiles
This present study is the first investigation to simultaneously elucidate differential metabolites and proteomes of water and swamp buffalo milk from two famous native breeds in Yunnan province.To investigate the difference in metabolite profiles of BBM and DBM, an Acquity UPLC-MS/MS coupled to a Thermo Q Exactive HF-X high-resolution mass spectrometry was used, and 500 and 319 metabolites were identified in ESI+ and ESI -modes, respectively (Table S1).The information on the determination of major differential milk metabolite compounds showed an apparent propensity for separation between the two buffaloes (p ≤ .05),and the information on the determination of those milk metabolite compounds is presented in Tables 2 and 3. A continuous scan of average m/z was controlled between 80 and 790 at different retention times within 1.38-13.07,indicating that the most differential metabolite compounds were detected within 1-4 min.The UPLC-MS/MS analysis results indicated that pipecolic acid, fomepizole, 5-hydroxyindoleacetic acid, 1,2,3-trihydroxybenzene, and maltol metabolite compounds detected in [M+H]+ were more relatively abundant with a higher peak area than other metabolites in BBM, which were 22.18 ppm, 16.09 ppm, 16.04 ppm, 15.21 ppm, and 13.05 ppm, respectively (Table 2).
In bovine animals, particularly in buffalo, the greater concentrations of L-pipecolic acid in milk may be linked to and involved in the lysine degradation pathway (bbub00310), especially to form L-2-aminoadipate 6-semialdehyde (Jorge-Smeding et al., 2021).In addition, fomepizole compounds play an essential role in blocking zinccontaining alcohol dehydrogenase (ADH) to promote acute alcoholic liver damage in buffalo (bbub04936) and other mammals (Tvrdý et al., 2021).5-hydroxyindoleacetic acid (5HIAA), a biogenic amine related to the C3 carbon atom of an indole that is a primary metabolite of serotonin, metanephrines, or catecholamines (bbub_M00042), is catalyzed by the mitochondrial enzyme aldehyde dehydrogenase and commonly used as a neuroendocrine tumor secretory biomarker (Corcuff et al., 2017).Related to bovine milk processing, Kokkinidou and Peterson (2014) reported that the addition of phenolic compound 1,2,3-trihydroxybenzene prior to ultrahigh temperature pasteurization significantly decreased the key transient of Maillard reaction (MR) intermediates and linked off-flavor compounds.Jo et al. found that during pasteurization, changes in the amino groups and lactose milk proteins in MR were caused by the heat treatment process.These changes caused the amino groups and lactose milk proteins to unfold, which made sweet aromatic volatile compounds like maltol (Jo et al., 2017).
In DBM, the most relatively abundant metabolites were cellobiose (29.68 ppm), 2-deoxystreptamine (20.86 ppm), phloroglucinol (13.13 ppm), 1,2,3-trihydroxybenzene (12.22 ppm), and creatine (10.34 ppm) (Table 3).Cellobiose was identified in [M -H]-and is involved in several pathways, including starch and sucrose metabolism (ko00500), metabolic pathways (ko01100), ATP-binding cassette transporters (ko02010), and the phosphotransferase system (ko02060).The pyranosyl groups in cellobiose and lactose are connected to the C-4 atom of d-glucose by a β-(1,4') glycosidic linkage (Ouellette & Rawn, 2019).The high concentration of 2-deoxystreptamine found in DBM, identified in ESI + [M+H]+, was strongly associated with biosynthesis of secondary metabolites (ko01110) and had an essential function in antimicrobial activity against protozoal and bacterial infections (ko00524).Through balancing metabolic hydrogen in the bovine rumen, the phenolic compound phloroglucinol contributes to stimulating ruminant methane reduction (Martinez-Fernandez et al., 2017;Sarwono et al., 2019).In DBM, the 1,2,3-trihydroxybenzene differential metabolite compound was also found in a higher concentration, although lower than in BBM.Additionally, creatine plays a vital role in arginine and proline metabolism (bbub00330) and any metabolic pathway (bbub01100) in buffalo.Creatine, along with important metabolic sources like glycine, glutathione, and serine plays an important role in a number of biological processes, such as the activity of lactating mammary epithelial cells, the balance of energy in bovine animals, and their response to heat stress (López et al., 2021;Sun et al., 2017).
Catechol, L-proline, γ-aminobutyric acid (GABA), hippuric acid, creatinine, betaine, D-fructose, taurine, cytosine, L-carnitine, and butyryl-L-carnitine were also detected in both buffalo milk groups as the major metabolite compounds with a high concentration.The mass-to-charge in mzXML format of those metabolites converted from MS/MS spectra is shown in Figure S2.Previous studies have reported that the differential metabolite and proteome profiles of buffalo milk vary by breed, location, and condition of the rearing region, and it was estimated that technological properties in making dairy by-products also change (de Nicola et al., 2020;Garau et al., 2021;Li et al.,2018Li et al., , 2021;;Pu et al., 2021;Salzano et al., 2020;Shi et al., 2021).These results were also in line with our current research, which found that the metabolomes and proteomes of buffalo milk were different in light of breed, especially between the BBM and DBM.

Buffalo milk proteome profiles
One hundred thirteen proteins were identified and quantified through analyzing the BBM and DBM datasets resulting from the label-free HPLC-MS/MS-based proteomics approach.The specific number of identified proteins in BBM and DBM is demonstrated in Figure 1.Of these identified proteins, there were 96 quantitative proteins identified in BBM, which is higher than those in DBM, with 64 proteins detected, and 47 proteins were found in both buffalo milk samples (Table S4).For quality control of the distributions of all the identified protein quantitative values in the The superscript letter after the metabolite compound name corresponds to the bioinformatics database used in this identification, where (a): BioDeepDB, ( corresponding samples, analyses of the relative standard deviation (RSD) and correlation coefficients between two samples were performed and are showcased in Figure 2. The measurements showed that the coefficient of variation (CV), which is the ratio of the standard deviation to the mean value, was greatly reduced after data correction.This led to a lower CV, which showed good repeatability and reproducibility in identifying the same analyte under different conditions and with certain adjustments.Furthermore, the analysis of the physical and chemical properties of the identified proteins showed that the isoelectric point (pI) value ranged between 4.5 and 5.5, the MW was mainly observed within 20-80 kDa (Figure S3a), the mass error was centered at 0 and concentrated in the range below 10 ppm (Figure S3b), and most of the peptide lengths were distributed between 7 and 23 amino acid residues (Figure S4).The results indicated that the pI, MW, mass error, and peptide length met the requirements and conformed to the rule of trypsin digestion of peptides.The cumulative distribution of quantitative protein value and the degree of dispersion of the data distribution in the corresponding samples were also evaluated and shown in Figure S5, where it was concluded that the protein quantitative value distribution between two buffalo milk samples was not completely the same.Both BBM and DBM suggested containing the same common proteomes in large numbers, which also included four types of phosphoprotein casein family; these were α-s1casein (CSN1S1), α-s2-casein (CSN1S2), β-casein (CSN2), κcasein (CSN3), albumin (ALB), α-lactalbumin (LALBA), lipoprotein lipase (LPL), lactotransferrin (LTF), fatty acidbinding protein (FABP3), and progestogen associated endometrial protein (PAEP) (Table S4).Three out of four caseins were detected in higher quantities in BBM than those in DBM; those were α-s1-casein, α-s2-casein, βcasein (Table S4).Among the casein family, α-s2-casein is the most hydrophilic, with more stable calcium-phosphate micelles, non-thermostable, non-acidic peptides, and has a lower efficiency in translation than other caseins (Cosenza et al., 2011;Rehman et al., 2021).Previous studies reported that α-s1-casein, α-s2-casein, β-casein, and β-lactoglobulin were found to be among the most abundant proteins in bovine milk based on qualitative similarity profiling analysis and were utilized further as the key protein ingredients in infant formula production as well (Raikos & Dassios, 2014;Smolenski et al., 2014;Xu et al., 2022).
The α-s1-casein protein (Uniprot ID O62823) is highly related to compositional aspects in buffalo milk, where a higher content of α-s1-casein exerts higher protein and total solids, increasing appropriateness in cheese-making properties (Kashwa, 2016).In humans, α-s1-casein contributes to the immune system's development, delivers caseins from the endoplasmic reticulum to the Golgi apparatus, and provides a source of amino acids for newborns (Stocco et al., 2018;Zhang et al., 2021aZhang et al., , 2021b)).The α-s2-casein protein-coding gene (O62825) plays an essential role in the capable transportation of calcium phosphate (Cosenza et al., 2021;Rehman et al., 2021).Lactotransferrin (O77698), also known as lactoferrin, is a major iron-binding transport protein that regulates antimicrobial activity in terms of the prevention of microbiological infection diseases, which is related to its ability to sequester free iron and to release lipopolysaccharides from the bacterial outer membrane (Giacinti et al., 2013).In relation to human health, lactoferrin is responsible for various critical functional actions in preventing necrotizing enterocolitis and sepsis and reducing infection risks associated with respiratory and gastrointestinal pathogens in neonates, infants, and young children (Björmsjö et al.,  2022; He et al., 2018;Manzoni, 2016;Manzoni et al., 2018).Moreover, albumin (P02769) has a ligand-binding capacity in the transportation of a diverse range of molecules, nutrients, metals, and metabolites and controls the growth of enterobacteriaceae, such as Escherichia coli and Salmonella enterica, through its ability to bind the siderophore enterobactin and inhibit enterobactin-mediated iron uptake of bacteria (Behnsen et al., 2021;Majorek et al., 2012).
Bovine milk casein contributes to multifunctional functions, including stimulating cellular immune functions, providing amino acids and calcium, and stabilizing chemotactic responses (Kawasaki et al., 2011;Vordenbaumen et al., 2013Vordenbaumen et al., , 2016;;Wada & Lonnerdal, 2014;Yang et al., 2016).In dairy by-product making, such as cheese manufacturing, casein structure and composition are technically influential.The understanding of buffalo milk protein expressions expands the insights about the behavior of curd formation and the aggregating mechanisms of Ca 2+ sensitive proteins via calcium-phosphate bridges in casein micelle, in particular strongly related to κ-casein content and reversible protein post-translational modifications (PTMs) (D'Ambrosio et al., 2008;Fan et al., 2019;Gai et al., 2021;Li et al., 2016;Nguyen et al., 2017).However, the low relative abundance of proteins presented in BBM and DBM also plays an essential function in several signaling pathways in buffalo milk production, including membrane/protein trafficking (Ras-related protein Rap-1b, Ras-related protein Rab-18, ADP-ribosylation factor 2), cell signaling (CIDE-N domain-containing protein), defense/immunity (HHIP like 2), and fat transport/ metabolism (diazepam binding inhibitor, acyl-CoA binding protein).

GO enrichments of identified proteins
All identified proteins in BBM and DBM were subjected to GO enrichment analysis to obtain insights into three distinctive functional classifications and characteristics, particularly in exploring the possible molecular mechanisms influencing buffalo milk production.GO analysis demonstrated that 73 identified proteins were enriched in a broad range of complex biological processes, molecular functions, and cellular components, which were classified into 22, 6, and 16 functional categories, respectively (p < .05)(Figure 3).The GO enrichments for mapping upregulated differentially expressed proteins in BBM and DBM are listed in Table S2 accordingly.Among those proteins, the biological processes are mainly associated with cellular processes (11%), biological regulation (9%), regulation of biological processes (9%), metabolic processes (8%), response to stimulus (8%), and multicellular organismal processes (6%) (Figure 3a).Expressed proteins related to cell, biological regulation, and metabolic processes, as well as homeostasis, are well known in relation to buffalo milk production and composition (Du et al., 2019(Du et al., , 2020)).
The cellular process terms (GO:0009987) with the largest number of expressed proteins (38) in the biological processes and networks play an important role in physiological changes in the buffalo mammary gland and ovarian follicle to produce milk, which is mainly related to cell differentiation, proliferation, apoptosis, and intracellular transportation (Buaban et al., 2022;Ye et al., 2022).In terms of the molecular functions, those expressed proteins included binding (53%), catalytic activity (20%), molecular function regulator (14%), structural molecule activity (8%), and antioxidant activity (3%) (Figure 3b).The binding terms (GO:0005488) with the highest number of expressed proteins (34) in molecular function describe the non-covalent, selective, often stoichiometric, interaction of a molecule with one or more molecules or certain sites.While the primary cellular component terms involved in cells, cell parts, organelles, extracellular regions, and membranes were 17%, 16%, 14%, 10%, and 8%, respectively.The cell terms (GO:0005623) with the most expressed proteins (38) in cellular components are in charge of basic structural and functional units in buffalo (Figure 3c).
We also found eight significantly upregulated differentially expressed protein-coding genes in BBM and DBM in the top 20 enriched pathways, including family keratin of epithelial cell markers (KRT15, KRT18, KRT28); ACTB, determining cell motility in relaying chemical signals; C3, regulating the innate immune system; RhoA, RhoB, transducing extracellular signals to the actin cytoskeleton for modifying axon growth and outgrowth cone motility; and CDC42, involved in epithelial cell polarization processes (Table S3).Together with lymphocytes, neutrophils, and macrophages, exfoliated epithelial cells comprise milk somatic cells (Lu et al., 2016;Robenek et al., 2006).Mammary epithelial cells can develop apocrine secretion to produce milk fat, referred to as milk fat globule proteins, after being formed in the endoplasmic reticulum membranes and delivered to the apical plasma membrane (Walter et al., 2020).Milk fat globule proteins are linked to a number of biological functions and are a key indicator of the nutritional value of bovine milk (Zhang et al., 2021a).

PPI networks
Protein expression data were converted to co-expression networks in this experiment for a network-level comparison of the BBM and DBM proteomes, visualizing protein-protein interactions in both breeds.The interaction score is set to 0.4, and the density of observed high interactions was significant according to the statistical calculation available in STRING (p < 10 −6 ).One hundred fifteen experimentally determined interactions and 359 nodes, with approximately 90% of the proteins having at least two or more connections, were observed.The combined score ranged from 0.401 to 0.999, with the highest score of the co-expression protein network being 0.925 (Table S5).The preponderance of multitudinous interactions among proteins in the PPI network demonstrated a strongly interconnected network where most proteins influenced the expression of other proteins related to buffalo milk synthesis (Figure 5).The PPI network for BBM and DBM proteomes also exhibited a dense core with highly interconnected proteins.

Conclusion
This present study provides detailed investigation results for differential metabolites and proteomes of BBM and DBM through an untargeted metabolomics and quantitative labelfree proteomics approach.Main polypeptide fractions, including α-s1-casein, β-s2-casein, β-casein, κ-casein, and αlactalbumin, were identified in higher protein abundances.Differentially expressed proteins were primarily involved in regulating cholesterol metabolism, the infection of S. aureus and Salmonella, the peroxisome proliferator-activated receptor signaling pathway, and coagulation cascades.By elucidating some of the primary differential metabolites and proteomes influencing cheese-making properties and organoleptic qualities, together with total protein/fat content, Ca 2+ concentration, and pH level, this study demonstrates a fundamental basis for the coagulation mechanisms applied in dairy production.Additionally, the number of minor protein fractions detected in this study also emphasizes the multiple functions of buffalo milk, particularly in providing nutrition for the newborn through its significant metabolites and protein components, supporting neonate development, and promoting the protection of the gastrointestinal tract from pathogens.Our findings provided a more in-depth reference for both human health and dairy food manufacturers based on buffalo milk processing.

Figure 1 .
Figure 1.Venn diagram of total identified protein numbers in BBM and DBM.

Figure 2 .
Figure 2. (a) RSD coefficient distribution before and after data correction of all the protein quantitative values in the corresponding sample.A lower CV coefficient indicated better repeatability of all samples.(b) Pearson correlation coefficient of all quantitative protein expressions between the two corresponding samples.

Figure 3 .
Figure 3.The distribution of enriched GO terms in BBM and DBM on three functional categories.

Figure 4 .
Figure 4.The distribution of enriched KEGG pathway in BBM and DBM.

Table 1 .
Milk composition comparison of BBM and DBM.
a Different superscript letters between the groups indicate significantly different values (p < .05).

Table 2 .
Determination of major differential metabolite compounds of BBM by UPLC-MS/MS.
The superscript letter after the metabolite compound name corresponds to the bioinformatics database used in this identification, where (a): BioDeepDB (b): mzCloud, and (c): MoNA; rt: retention time; KEGG refers to metabolic pathway mapping in KEGG mapper; ion source is associated to the sensitive atmospheric pressure ESI, both in ESI + and ESI -mode; and peak area represents the intensity measure of the metabolite compound.

Table 3 .
Determination of major differential metabolite compounds of DBM by UPLC-MS/MS.