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.

KEYWORDS:

• Buffalo milk
• untargeted metabolomics
• label-free proteomics
• metabolic pathways

1. 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⁻¹ of voluminosity and 2.60–2.90 g water g⁻¹ 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 high-nutritive 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 high-performance 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.

2. Materials and methods

2.1. 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.

2.2. 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.

2.3. 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.

2.4. Database searching and metabolite identification

Raw MS/MS spectra were converted into mzXML format using Proteowizard software v3.0.8789. The XCMS package of R v3.1.3 was used to perform peak identification, filtration, and alignment processing to generate a data matrix, including mass-to-charge ratio, retention time, peak area, and precursor molecules in both positive and negative ion modes, which was exported to MS Excel for subsequent analysis (Smith et al., 2006). The metabolite identification was confirmed according to their molecular weight (MW), with an MW error of ≤ 30 ppm. The Human Metabolome Database (HMDB) (http://www.hmdb.ca), METLIN (http://metlin.scripps.edu), BioDeepDB (http://www.biodeep.cn/), MoNA (https://massbank.us/), MetDNA (http://metdna.zhulab.cn/), LipidMaps (http://www.lipidmaps.org), and mzCloud (https://www.mzcloud.org) were used to confirm the annotations according to the MS/MS fragmentation mode.

2.5. 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.

2.6. 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 nano-electrospray 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 peptide-spectrum 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.

2.7. 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).

3. Results and discussion

3.1. 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⁶⁷/His⁶⁷ 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.

Table 1. Milk composition comparison of BBM and DBM.

Milk Source	Fat (%)	Protein (%)	Lactose (%)	Total Solids (%)	A2 β-casein (ng/ml)
BBM	8.23 ± 0.13^a	4.25 ± 0.06^a	5.21 ± 0.06^a	18.22 ± 0.05^a	18.81 ± 0.28^a
DBM	8.16 ± 0.07 ^a	4.18 ± 0.04^b	5.17 ± 0.05^a	18.28 ± 0.05^a	18.68 ± 0.37^a

Different superscript letters between the groups indicate significantly different values (p < .05).

3.2. 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).

Table 2. Determination of major differential metabolite compounds of BBM by UPLC-MS/MS.

Metabolite	Chemical Formula	m/z	rt (min)	Exact Mass	Concentration (ppm)	KEGG	Ion Source	Peak Area
Glycerophosphocholine^a	C8H21NO6P	258.11	4.30	258.11	1.63	C00670	[M]+	4.07E + 10
L-Carnitine^a	C7H15NO3	162.11	2.70	161.11	.07	C00318	[M+H]+	3.30E + 10
Catechol^a	C6H6O2	111.02	1.85	110.04	7.59	C00090	[M+H]+	3.29E + 10
γ-Aminobutyric acid^a	C4H9NO2	104.11	1.74	103.06	7.03	C00334	[M+H]+	1.84E + 10
Acetoacetic acid^a	C4H6O3	85.03	1.42	102.03	1.12	C00164	[M+H – H2O]+	1.62E + 10
2-Furancarboxaldehyde^a	C5H4O2	97.03	1.62	96.02	.30	C14279	[M+H]+	1.51E + 10
Betaine^a	C5H11NO2	118.09	1.97	117.08	4.47	C00719	[M+H]+	9.30E + 09
2-Keto-glutaramic acid^a	C5H7NO4	145.05	2.42	145.04	3.49	C00940	[M]+	9.21E + 09
L-Glutamic acid^a	C5H9NO4	148.06	2.47	147.05	.01	C00025	[M+H]+	6.65E + 09
2-Ketobutyric acid^a	C4H6O3	102.04	1.70	102.03	9.66	C00109	[M]+	6.64E + 09
Quinolinic acid^a	C7H5NO4	167.01	2.78	167.02	.31	C03722	[M]+	5.72E + 09
Butyryl-L-carnitine^a	C11H21NO4	232.15	3.87	231.15	.0007	C02862	[M+H]+	5.05E + 09
D-Fructose^a	C6H12O6	163.06	2.72	180.06	2.51	C00095	[M+H – H2O]+	4.50E + 09
Creatinine^a	C4H7N3O	114.07	1.90	113.06	5.27	C00791	[M+H]+	4.36E + 09
Hippuric acid^a	C9H9NO3	178.05	2.97	179.06	5.54	C01586	[M – H]–	3.91E + 09
L-Malic acid^a	C4H6O5	133.01	2.22	134.02	1.22	C00149	[M – H]–	3.75E + 09
Taurine^a	C2H7NO3S	125.99	2.10	125.01	1.07	C00245	[M+H]+	3.08E + 09
Acetylcholine^a	C7H16NO2	146.12	2.44	146.12	.31	C01996	[M+H]+	2.79E + 09
1,2,3-Trihydroxybenzene^b	C6H6O3	127.04	2.12	126.03	15.21	C01108	[M+H]+	2.78E + 09
cis-Aconitate^b	C6H6O6	173.01	2.88	174.02	3.94	C00417	[M – H]–	2.56E + 09
myo-Inositol^a	C6H12O6	181.02	3.02	180.06	6.01	C00137	[M+H]+	2.43E + 09
N-Methyl-2-pyrrolidinone^b	C5H9NO	100.08	1.67	99.07	1.12	C11118	[M+H]+	2.32E + 09
Phloroglucinol^b	C6H6O3	127.04	2.12	126.03	11.92	C02183	[M+H]+	2.30E + 09
Maltol^a	C6H6O3	127.04	2.12	126.03	13.05	C11918	[M+H]+	1.93E + 09
Fomepizole^c	C4H6N2	83.06	1.38	82.05	16.09	C07837	[M+H]+	1.91E + 09
Pantothenic acid^a	C9H17NO5	220.12	3.67	219.11	2.95	C00864	[M+H]+	1.90E + 09
Cytosine^a	C4H5N3O	112.05	1.87	111.04	.83	C00380	[M+H]+	1.87E + 09
Hydroquinone^a	C6H6O2	110.02	1.83	110.04	6.95	C00530	[M]+	1.82E + 09
Piperidine^c	C5H11N	86.10	1.43	85.09	7.35	C01746	[M+H]+	1.47E + 09
Pipecolic acid^a	C6H11NO2	130.05	2.17	129.08	22.18	C00408	[M+H]+	1.46E + 09
5-Hydroxyindoleacetic acid^a	C10H9NO3	192.06	3.20	191.06	16.04	C05635	[M+H]+	1.43E + 09
Glucose 6-phosphate^b	C6H13O9P	259.02	4.32	260.03	2.22	C00092	[M – H]–	1.41E + 09
L-Proline^a	C5H9NO2	116.07	1.93	115.06	7.20	C00148	[M+H]+	1.05E + 09
Mesaconate^b	C5H6O4	129.02	2.15	130.03	1.77	C01732	[M – H]–	1.04E + 09
Uracil^a	C4H4N2O2	113.03	1.88	112.03	.03	C00106	[M+H]+	9.98E + 08

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.

Metabolite	Chemical Formula	m/z	rt (min)	Exact Mass	Concentration (ppm)	KEGG	Ion Source	Peak Area
Glycerophosphocholine^a	C8H21NO6P	258.11	4.30	258.11	7.10	C00670	[M]+	6.22E + 10
Creatine^a	C4H9N3O2	132.08	2.20	131.07	1.34	C00300	[M+H]+	4.14E + 10
Catechol^a	C6H6O2	111.02	1.85	110.04	8.68	C00090	[M+H]+	3.27E + 10
Phosphorylcholine^a	C5H15NO4P	184.07	3.07	184.07	.19	C00588	[M+H]+	1.51E + 10
Acetoacetic acid^a	C4H6O3	85.03	1.42	102.03	7.79	C00164	[M+H – H2O]+	1.50E + 10
2-Furancarboxaldehyde^a	C5H4O2	97.03	1.62	96.02	3.70	C14279	[M+H]+	1.37E + 10
L-Carnitine^a	C7H15NO3	162.11	2.70	161.11	.001	C00318	[M+H]+	1.20E + 10
γ-Aminobutyric acid^a	C4H9NO2	104.11	1.74	103.06	8.27	C00334	[M+H]+	1.18E + 10
2-Ketobutyric acid^a	C4H6O3	102.03	1.70	102.03	3.43	C00109	[M]+	6.83E + 09
D-Fructose^a	C6H12O6	163.06	2.72	180.06	2.66	C00095	[M+H – H2O]+	4.38E + 09
Cellobiose^a	C12H22O11	341.10	5.68	342.12	29.68	C00185	[M – H]–	3.78E + 09
Butyryl-L-carnitine^a	C11H21NO4	232.15	3.87	231.15	.004	C02862	[M+H]+	3.28E + 09
Taurine^a	C2H7NO3S	125.99	2.10	125.01	.96	C00245	[M+H]+	3.19E + 09
(S)-4-Hydroxymandelate^a	C8H8O4	169.05	2.82	168.04	.48	C03198	[M+H]+	2.89E + 09
Betaine^a	C5H11NO2	118.09	1.97	117.08	.61	C00719	[M+H]+	2.89E + 09
6-Acetyl-D-glucose^d	C8H14O7	221.07	3.68	222.07	5.64	C02655	[M – H]–	2.72E + 09
Creatinine^a	C4H7N3O	114.07	1.90	113.06	5.66	C00791	[M+H]+	2.68E + 09
1,2,3-Trihydroxybenzene^b	C6H6O3	127.04	2.12	126.03	12.22	C01108	[M+H]+	2.62E + 09
Fomepizole^c	C4H6N2	83.06	1.38	82.05	1.04	C07837	[M + 2ACN +2 H]2+	2.27E + 09
Stachyose^a	C24H42O21	667.23	11.12	666.22	3.17	C01613	[M+H]+	1.47E + 09
Hippuric acid^a	C9H9NO3	178.05	2.97	179.06	5.35	C01586	[M – H]–	1.41E + 09
L-Glutamine^a	C5H10N2O3	147.08	2.45	146.07	.002	C00064	[M+H]+	1.39E + 09
11-Dehydrocorticosterone^a	C21H28O4	325.18	5.42	344.20	8.86	C05490	[M – H2O – H]–	1.37E + 09
L-Leucine^a	C6H13NO2	132.10	2.20	131.09	1.15	C00123	[M+H]+	1.15E + 09
Pantothenic acid^a	C9H17NO5	220.12	3.67	219.11	.10	C00864	[M+H]+	1.11E + 09
Quinoline^a	C9H7N	130.16	2.17	129.06	1.20	C06413	[M+H]+	1.02E + 09
Phloroglucinol^b	C6H6O3	127.04	2.12	126.03	13.13	C02183	[M+H]+	1.02E + 09
L-Proline^a	C5H9NO2	116.07	1.93	115.06	2.29	C00148	[M+H]+	1.00E + 09
Gluconic acid^a	C6H12O7	195.05	3.25	196.06	.30	C00257	[M – H]–	9.49E + 08
2-Deoxystreptamine^a	C6H14N2O3	163.12	2.72	162.10	2.86	C02627	[M+H]+	9.24E + 08
Cytosine^a	C4H5N3O	112.05	1.87	111.04	.09	C00380	[M+H]+	9.07E + 08
Phosphocreatine^a	C4H10N3O5P	210.03	3.50	211.04	7.44	C02305	[M – H]–	8.99E + 08
4,5-Dihydroorotic acid^a	C5H6N2O4	158.96	2.65	158.03	7.90	C00337	[M+H]+	8.91E + 08
Arachidic acid^a	C20H40O2	311.17	5.19	312.30	.46	C06425	[M – H]–	8.62E + 08
Methyloxaloacetate^a	C5H6O5	147.03	2.60	146.02	.02	C06030	[M+H]+	8.21E + 08

The superscript letter after the metabolite compound name corresponds to the bioinformatics database used in this identification, where (a): BioDeepDB, (b): mzCloud, (c): MoNA, and (d): MetDNA.

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 zinc-containing 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.,2018, 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.

3.3. 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 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.

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

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.

Furthermore, the analysis of the physical and chemical properties of the identified proteins showed that the iso-electric 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 α-s1-casein (CSN1S1), α-s2-casein (CSN1S2), β-casein (CSN2), κ-casein (CSN3), albumin (ALB), α-lactalbumin (LALBA), lipoprotein lipase (LPL), lactotransferrin (LTF), fatty acid-binding 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., 2021a, 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., 2013, 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²⁺ 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).

3.4. 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, 2020).

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

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).

3.5. KEGG pathways of differential expressed proteins

The KEGG annotations provided a critical reference for investigating the functions and biological systems in buffalo milk. In this study, 52 enriched KEGG maps were obtained using the identified proteins, which were primarily related to the top 20 pathways, and many proteins in those pathways were associated with cholesterol metabolism (bbub04979), Staphylococcus aureus (S. aureus) infection (bbub05150), peroxisome proliferator-activated receptor (PPAR) signaling pathway (bbub03320), coagulation cascades (bbub04610), and Salmonella infection (bbub05132) (Figure 4 and Table S3). Forty-five genes were detected in differential protein expression linked to those functional pathways. KEGG enrichment analysis revealed the highest number of involved genes in the KEGG pathway was exhibited in S. aureus infection with a rich factor of 0.06 (p < .05), including complement C3 (C3), complement component 4A (C4A), keratins (KRT15, KRT18, KRT28), and complement factor B (CFB) genes. While NPC intracellular cholesterol transporter 2 (NPC2), lipoprotein lipase (LPL), CD36 molecule (CD36), apolipoprotein A1 (APOA1), and angiopoietin-like 4 (ANGPTL4) genes were identified in cholesterol metabolism; LPL, CD36, APOA1, ANGPTL4, and fatty acid-binding protein 3 (FABP3) were in the PPAR signaling pathway; C3, C4A, serpin family A member 1 (SERPINA1), clusterin (CLU), and complement factor B (CFB) were in coagulation cascades; and β-actin (ACTB), myeloid cell-specific leucine-rich glycoprotein CD14 (CD14), ras homolog GTPase family member A and B (RhoA, RhoB), and leucocyte-specific cell division cycle 42 homolog precursor (CDC42) genes were in Salmonella infection.

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

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).

3.6. 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⁻⁶). 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.

Figure 5. Ppis of the identified proteins in BBM and DBM. STRING database was used to establish the protein-protein interactions of quantified proteins identified by HPLC-MS/MS. Each node in the PPI network represents an identified protein, and the type of node represents the degree of connection of the protein.

4. Conclusion

This present study provides detailed investigation results for differential metabolites and proteomes of BBM and DBM through an untargeted metabolomics and quantitative label-free 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²⁺ 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.

Acknowledgments

We would like to extend our heartfelt gratitude to the anonymous reviewers and colleagues for their constructive comments and advice on this manuscript. We also thank Professor Huaming Mao in College of Animal Science and Technology, Yunnan Agricultural University, and Suzhou Panomic Biomedical Technology Co., Ltd., China, for experimental support.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/19476337.2023.2264953.

Additional information

Funding

This study was financially supported by the Leading Talent Program of Yun-ling Industrial Technology under Grant number 2014-1782 and the Dairy Processing and Food Safety Program of Yunnan Modern Agricultural Cattle Milk Industrial Technology System under Grant number [2019KJTX0014].