Physiological and Proteomic Responses of Dairy Buffalo to Heat Stress Induced by Different Altitudes

Buffalo are mainly distributed in low-altitude (LA), medium-altitude (MA), and high-altitude (HA) regions characterised by different thermal and oxygen environments in Yunnan province, China. Due to black skin, sparse hair, and the low density of skin sweat glands, buffalo are more sensitive to heat stress. Here, we used data-independent acquisition (DIA) proteomics to reveal a broad spectrum of proteins that play roles in adaptation to the heat stress of buffalo raised at low altitude or hypoxia at high altitude. LA buffalo showed higher body temperatures than MA- and HA buffalo, and HA buffalo had higher levels of GSH and SOD and lower levels of ROS compared to LA and MA buffalo. In 33 samples, 8476 peptides corresponding to 666 high-confidence proteins were detected. The levels of circulating complement proteins in the immune pathways were lower in LA and MA buffalo than in HA buffalo. There were higher levels of alpha-1 acid glycoprotein in LA buffalo than in MA and HA buffalo. Relative to MA buffalo, levels of blood oxygen delivery proteins were higher in LA and HA buffalo. A higher abundance of apolipoproteins was detected in LA and MA buffalo than in HA buffalo. In summary, buffalo adopted similar adaptation strategies to oxidative stress induced by heat stress or hypoxia, including immunological enhancement, high efficiency of blood oxygen delivery, and the inhibition of lipid oxidation.


Introduction
Buffalo are the second-largest source of milk globally. Buffalo milk is characterised by high concentrations of proteins, fat, and total solids compared to cow milk [1]. China has rich buffalo resources and produces 5% of the world's buffalo milk [2]. Yunnan Province is located at altitudes ranging from 76 to 6740 m in Southwestern China and has complex topography. It neighbours Myanmar, Lao, and Vietnam. The geographic distribution of dairy buffalo in Yunnan Province are the states of Dehong, Baoshan, and Dali, which are characterised as low-, medium-, and high-altitude, respectively. There are about 56,000, 8000, and 15,000 buffaloes raised in the three states, respectively. Dehong state, as a typical lowland subtropical area, is characterised by high ambient temperature and humidity in warm seasons, causing chronic heat stress for buffalo in summer [3].
Higher temperatures resulting from global climate change cause low levels of animal production performance, health, and welfare [4,5]. Buffalo are warm-blooded animals, maintaining a constant body temperature of about 38.0 +/−0.5 • C. Due to their black skin, sparse hair, and the low density of skin sweat glands-one-sixth the density of cattle [6]. The optimum temperatures and humidity of buffalo are 13-18 • C and 55-65%, respectively [7], so they are more sensitive to heat stress-induced oxidative stress when the body heat load is greater than heat dissipation. In low-altitude raising regions, high humidity may restrain the efficiency of evaporative heat dissipation in buffalo [8]. Heat stress at low altitudes may adversely affect animal production performance, health, and welfare [9]. In

Experimental Animals and Management
Dairy buffalo farms used in the experiment were situated at a low altitude of 912.0 m (LA), a medium altitude of 1536.0 m (MA), or a high altitude of 1865.0 m (HA). Thirty-three healthy multiparous (parity = 3.2 ± 1.1) Nili-Ravi × Murrah × Local crossbred buffalo (9,12, and 12 LA, MA, and HA buffalo, respectively) at mid-lactation (115 ± 10 days) were group-housed in an open-sided barn and fed two times a day using the same total mixed rations (80% whole-plant corn silage ad libitum, 12.5% concentrate feeding, 7.5% corn protein powder) in June. The nutritional values of the three feed ingredients are shown in Supplemental Table S1.
Temperature and humidity loggers (±0.2 • C; Testo 175H1), placed 2.0 m above the floor of the bedding area, recorded the ambient temperature and relative humidity at 30 min intervals to calculate the temperature-humidity index (THI) according to the following equation [17] in the field trails: THI = (1.8 × T + 32) − (0.55 − 0.0055 × RH) × (1.8 × T − 26) Here, T is the ambient temperature ( • C), and RH is the relative humidity (%). A waterproof micro-temperature sensor (± 0.5 • C; DS1922L, Wdsen Electronic Technology Co., Ltd. Shanghai, China), fixed in the slot of a T-shaped internal controlled drug release device (without progesterone; DEC International, NZ, Ltd. Hamilton, New Zealand), was placed in the vagina of buffalo to record body temperature at 30 min intervals. At the end of the field trial (15 d), the blood samples were collected with vacutainer tubes via jugular venipuncture. Each sample was centrifuged at 1400× g for 10 min to separate the serum. The serum samples were quickly frozen in liquid nitrogen and then stored at −80 • C for later assay.

Sample Preparation and Fractionation for DDA Library Generation
The serum protein concentrations were determined using a BCA protein assay kit. The inter-rater agreement of all the protein samples was evaluated via an SDS-PAGE gel with Coomassie blue staining. About 20-40 µg of protein from each sample was used for data-dependent acquisition (DDA) library construction. The LA, MA, and HA samples were analysed individually.
The FASP procedure was applied for protein digestion. The protein sample (200 µg) was added to urea buffer (8 M) to a final volume of 60 µL. To this was added dithiothreitol (DTT; 1 M, 0.5 µL), and the mixture was incubated at 56 • C for 30 min. Iodacetamide (25 µL) was then added with stirring, and the resulting mixture was incubated at room temperature for 30 min. NH4HCO3 (100 mM, 240 µL) was added to dilute the urea solution to 1.5 M. To each sample tube, trypsin (4 µg) was added at a protein:trypsin ratio of 1:50 and the resulting mixture was incubated in a 37 • C water bath for 20 h. Formic acid (FA; 10%, 20 µL) was added to the protein mixture and then desalted using a Waters SPE column (Oasis MCX SPE 96-well plate, 186,008,304). Finally, the desalted samples were freeze-dried, and the peptide content was estimated via the UV light spectral density at OD280 (Nanodrop 2000C, Thermo Fisher Scientific, Waltham, MA, USA). Pooled peptides (200 µg) were fractionated into 10 fractions using a High-pH Reversed-Phase Peptide Fractionation Kit (Thermo Scientific™ Pierce™, 848,868).

Data-Dependent Acquisition Mass Spectrometry Assay
Data-dependent acquisitions were used to generate a spectral library to form the query database for the DIA mass spectra in later data analysis. Serum protein samples were fractionated using high pH reversed-phase separation to increase proteomic depth. The gradient elution of the peptides was performed on a homemade column (100 µm × 20 mm, 5 µm-C18). The linear gradient in the column was as follows: 10-30% solution B at 0-97 min, 30-100% solution B at 97-110 min, 100% solution B at 110-120 min and 100% solution B. The peptides were separated using a linear gradient of buffer B (84% acetonitrile and 0.1% FA) at a flow rate of 250 nL min −1 .
The DDA runs were performed on a Q-Exactive HF mass spectrometer (Thermo Scientific) coupled with an Easy nLC-1200 system (Thermo Scientific) for DDA analysis. Q-Exactive HF mass spectrometry was performed in the data-dependent mode to switch automatically between MS and MS/MS acquisition. The full-scan MS survey was in the positive ion mode, and the scan range was 300-1800 m/z. The mass resolution for MS1 was 60,000 at m/z 200 (3e6 AGC, 200 ms maximum injection time). The resolution for twenty sequential MS2 scans was 30,000 at m/z 200 (3e6 AGC, 120 ms maximum injection time). For MS/MS, the normalised collision energy was set at 27 eV. The DDA spectra were analysed using MaxQuant and were filtered for a false discovery rate (FDR) of 1% at the peptide and protein levels.

Mass Spectrometry Assay for Data-Independent Acquisition
The peptides in each sample were analysed by LC-MS/MS operating in the DIA mode by Shanghai Applied Protein Technology Co., Ltd., Shanghai, China. After the construction of the spectral library for database query, peptide (2 µg) mixed with iRT peptide was injected for DIA mass spectrometry. The Easy nLC-1200 system was used for chromatographic separation using a binary solvent system consisting of solvent A (0.1% FA) and solvent B (84% acetonitrile, 0.1% FA). Solvent A (95%) was used to equilibrate the analytical column (homemade tip column, 75 µm × 250 mm, 1.9 µm-C18). The peptide mixture was separated via gradient elution with a 25 cm tip column at a flow rate of 250 nL.min −1 using a linear gradient in the trap column (100 µm × 20 mm, 5 µm-C18; 10-30% solution B at 0-97 min, 30-100% solution B at 97-110 min and 100% solution B at 110-120 min).
HF mass spectrometry (Q-Exactive) was performed in the DIA mode to switch automatically between MS and MS/MS acquisition. The run time was 120 min in the positive ion mode with a linear gradient of buffer B (80% acetonitrile and 0.1% FA) at a flow rate of 250 nL min −1 . The quality control samples were injected in the DIA mode to monitor the MS performance. DIA scans covered a mass range of 350-1650 m/z with the following settings: SIM full scan resolution 120,000 at 200 m/z, AGC 3e6, 50 ms maximum injection time). Thirty sequential DIA MS2 scans were run at a resolution of 30,000 at 200 m/z, AGC target 3e6, automatic maximum injection time and normalised collision energy of 30 eV.
The MS search parameters were as follows: cutting enzyme, trypsin; maximum missed cleavages, 2; fixed modification, carbamidomethyl (C); dynamic modifications, oxidation (M) and acetyl (Protein N-term). The reported data were based on 99% confidence for protein identification by false discovery rate (≤1%).
The spectral library was constructed by importing the original raw files, and DDA search results were imported into Spectronaut Pulsar XTM_12.0.20491.4 (Biognosys AG, Schlieren, Switzerland). Spectronaut was used to analyse the DIA data with the above constructed spectral library. The main software parameters were set as follows: retention time prediction type, dynamic iRT; interference on MS2 level correction, enabled; cross run normalisation, enabled. All of the results were filtered based on a Q value cut-off of 0.01 (equivalent to FDR < 1%).

Principal Component Analysis, GO and KEGG Pathway Annotation
The relationship between the serum samples obtained from different altitude buffalo was visualised using a principal component analysis of the quantified proteins using R packages. Significantly differential serum proteins were identified at a fold-change of greater than 1.5 (<0.67) and a p-value of less than 0.05. Gene ontology (GO) terms and the Blast2GO software (http://www.blast2go.com/b2ghome, 18 September 2018) were used for the annotation of protein sequences, and the GO annotation results were plotted using R scripts. The differential protein sequences were blasted against the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://geneontology.org/, 30 September 2018) to identify KEGG Orthology. Differential proteins were mapped to KEGG pathways.

LC-PRM/MS Quantitative Validation of Targeted Proteins
The levels of serum proteins that are associated with heat stress or hypoxia in buffalo were quantified using parallel reaction monitoring (PRM). Each sample containing peptide (4 µg) was mixed with standard peptide ELGQSGVDTYLQTKSAAGAFGPELSR (200 fmol) after protein digestion. The gradient elution was performed on a C18 HPLC column. Solvent A (95%) was used to equilibrate the analytical column (Thermo Scientific EASY). The peptides were separated using a linear gradient of solvent B (84% acetonitrile and 0.1% FA) at a flow rate of 300 nL min −1 . The linear gradient in the column was as follows: 5-10% solution B at 0-2 min, 10-30% solution B at 2-45 min, 30-100% solution B at 45-55 min and 100% solution B at 55-60 min.
HF mass spectrometry (Q-Exactive) was performed in the PRM/MS acquisition mode. Full-scan MS surveys (m/z 300-1800) were acquired with a mass resolution (R) of 60,000 at m/z 200 (3e6 AGC 200 ms maximum injection time). This was followed by twenty sequential MS2 scans for higher-energy collisional dissociation (HCD) that conformed to the following parameters: isolation window, 1.6 Th; R, 30,000 at m/z 200; AGC, 3e6; maximum IT, 120 ms; MS2 activation type, HCD. For MS/MS, the normalised collision energy was 27. Skyline v3.5.0 was used to analyse the PRM raw data. The peptide settings used for the Skyline import were consistent with the MaxQuant search parameters (i.e., enzyme set to trypsin, max missed cleavages set to 2). Three consecutive high-intensity peptides were selected to import into Skyline for each set of Skyline analyses. The number of targeted proteins, the sequences of the targeted peptides, the charges of the parent ions, and the peak areas were determined. The total peak areas were calibrated against internal peptide standards using heavy isotopes, and the target peptides were quantified based on the total peak areas.

Physiological Parameters of Buffalo Raised at Different Altitudes
The THI values and buffalo body temperatures are shown in Figure 1. The THI values of LA and MA buffalo (66.3 and 67.2, respectively) were higher than for HA buffalo (62.6). LA buffalo had significantly higher body temperature (38.8 • C) than MA and HA buffalo (38.0 and 38.1 • C, respectively). No significant differences were observed in the serum levels of CRF, PRL, T3, ACTH, COR, CP, HSF, ROS, and HIF-1α in buffalo raised at different altitudes. In contrast, HA buffalo had higher levels of IGF-1, GSH and SOD and lower levels of ROS as compared to LA and MA buffalo ( Figure 2).

Spectral Library Construction of DIA Proteomics
The SDS-PAGE gel assay with Coomassie blue staining indicated that the proteins had good consistency in the serum samples ( Figure 3). The BCA protein assay showed differences in the amount of total protein (10.36, 10.44 and 11.82 µg µL −1 for LA, MA, and HA buffalo, respectively), so the same amount of protein was used for further analyses of all the samples (Supplemental Table S2). The samples were digested with the in-solution digestion method and identified using the DIA proteomic technique. The DDA library from Maxquant 1.5.3.17 had 509 protein groups and 3990 peptides, while the DIA library from Spectronaut Pulsar 12.0.20491.4 had 535 protein groups and 3790 peptides. When these libraries were merged, the resulting spectral library had 749 protein groups and 5388 peptides. The peptides were separated using a linear gradient of solvent B (84% acetonitrile and 0.1% FA) at a flow rate of 300 nL•min − 1. The linear gradient in the column was as follows: 5-10% solution B at 0-2 min, 10-30% solution B at 2-45 min, 30-100% solution B at 45-55 min and 100% solution B at 55-60 min. HF mass spectrometry (Q-Exactive) was performed in the PRM/MS acquisition mode. Full-scan MS surveys (m/z 300-1800) were acquired with a mass resolution (R) of 60,000 at m/z 200 (3e6 AGC 200 ms maximum injection time). This was followed by twenty sequential MS2 scans for higher-energy collisional dissociation (HCD) that conformed to the following parameters: isolation window, 1.6 Th; R, 30,000 at m/z 200; AGC, 3e6; maximum IT, 120 ms; MS2 activation type, HCD. For MS/MS, the normalised collision energy was 27. Skyline v3.5.0 was used to analyse the PRM raw data. The peptide settings used for the Skyline import were consistent with the MaxQuant search parameters (i.e., enzyme set to trypsin, max missed cleavages set to 2). Three consecutive high-intensity peptides were selected to import into Skyline for each set of Skyline analyses. The number of targeted proteins, the sequences of the targeted peptides, the charges of the parent ions, and the peak areas were determined. The total peak areas were calibrated against internal peptide standards using heavy isotopes, and the target peptides were quantified based on the total peak areas.

Physiological Parameters of Buffalo Raised at Different Altitudes
The THI values and buffalo body temperatures are shown in Figure 1. The THI values of LA and MA buffalo (66.3 and 67.2, respectively) were higher than for HA buffalo (62.6). LA buffalo had significantly higher body temperature (38.8 °C) than MA and HA buffalo (38.0 and 38.1 °C, respectively). No significant differences were observed in the serum levels of CRF, PRL, T3, ACTH, COR, CP, HSF, ROS, and HIF-1α in buffalo raised at different altitudes. In contrast, HA buffalo had higher levels of IGF-1, GSH and SOD and lower levels of ROS as compared to LA and MA buffalo ( Figure 2).

Spectral Library Construction of DIA Proteomics
The SDS-PAGE gel assay with Coomassie blue staining indicated that the proteins had good consistency in the serum samples ( Figure 3). The BCA protein assay showed differences in the amount of total protein (10.36, 10.44 and 11.82 μg μL -1 for LA, MA, and HA buffalo, respectively), so the same amount of protein was used for further analyses of all the samples (supplemental Table S2). The samples were digested with the in-solution digestion method and identified using the DIA proteomic technique. The DDA library from Maxquant 1.5.3.17 had 509 protein groups and 3990 peptides, while the DIA library from Spectronaut Pulsar 12.0.20491.4 had 535 protein groups and 3790 peptides. When these libraries were merged, the resulting spectral library had 749 protein groups and 5388 peptides.

Data Reliability and Principal Component Analysis
Spectronaut™ Pulsar X was used to evaluate the correlations in the quantitative data from the protein expression analyses. During data acquisition, there was an average of 7.4 data points per chromatographic peak, which was adequate for accurate quantification. The average peak capacity was 503.2, reflecting the highly efficient chromatographic sep-

Spectral Library Construction of DIA Proteomics
The SDS-PAGE gel assay with Coomassie blue staining indicated that the proteins had good consistency in the serum samples ( Figure 3). The BCA protein assay showed differences in the amount of total protein (10.36, 10.44 and 11.82 μg μL -1 for LA, MA, and HA buffalo, respectively), so the same amount of protein was used for further analyses of all the samples (supplemental Table S2). The samples were digested with the in-solution digestion method and identified using the DIA proteomic technique. The DDA library from Maxquant 1.5.3.17 had 509 protein groups and 3990 peptides, while the DIA library from Spectronaut Pulsar 12.0.20491.4 had 535 protein groups and 3790 peptides. When these libraries were merged, the resulting spectral library had 749 protein groups and 5388 peptides.

Data Reliability and Principal Component Analysis
Spectronaut™ Pulsar X was used to evaluate the correlations in the quantitative data from the protein expression analyses. During data acquisition, there was an average of 7.4 data points per chromatographic peak, which was adequate for accurate quantification. The average peak capacity was 503.2, reflecting the highly efficient chromatographic sep-

Data Reliability and Principal Component Analysis
Spectronaut™ Pulsar X was used to evaluate the correlations in the quantitative data from the protein expression analyses. During data acquisition, there was an average of 7.4 data points per chromatographic peak, which was adequate for accurate quantification. The average peak capacity was 503.2, reflecting the highly efficient chromatographic separation achieved by the DIA analysis. The elution times of the 11 standard peptides indicated the high stability of the chromatographic analysis. False-positive proteins were filtered at a 1% protein FDR threshold so that only proteins with high credibility were kept and all of the average CScore values were greater than 1. The principal component analysis of the quantified proteins provided a visualisation of the relationships between the samples in this study. The biological replicates clustered appropriately and could be distinguished to display the dataset with reasonable reproducibility (Figure 4). cated the high stability of the chromatographic analysis. False-positive proteins were filtered at a 1% protein FDR threshold so that only proteins with high credibility were kept and all of the average CScore values were greater than 1. The principal component analysis of the quantified proteins provided a visualisation of the relationships between the samples in this study. The biological replicates clustered appropriately and could be distinguished to display the dataset with reasonable reproducibility (Figure 4).

Differential Protein Identification
There were 8476 peptides identified in the serum, which corresponded to 666 highconfidence proteins in the Spectronaut Pulsar X search engine. The results from the protein quantitative analyses and differential expression analyses are shown in supplemental  Table S3. Analysis of variance revealed 60 proteins that were relevant to adaptation to environmental stresses that were significantly differentially abundant in the serum of buffalo raised at different altitudes (Table 1). These differentially expressed proteins in the serum were selected for further analysis based on the fold-change (≥1.50 or <0.67) and pvalue < 0.05 (supplemental Table S3). The serum proteins were displayed in a volcano plot (supplemental Figure S1), and hierarchical clustering showed they could be grouped according to their differential abundance (supplemental Figure S2). Functional analysis showed that these differentially expressed proteins were involved in immunoregulation, lipoprotein metabolism, the regulation of blood coagulation, oxygen binding and delivery, and antioxidant stresses ( Table 2).

Differential Protein Identification
There were 8476 peptides identified in the serum, which corresponded to 666 highconfidence proteins in the Spectronaut Pulsar X search engine. The results from the protein quantitative analyses and differential expression analyses are shown in Supplemental  Table S3. Analysis of variance revealed 60 proteins that were relevant to adaptation to environmental stresses that were significantly differentially abundant in the serum of buffalo raised at different altitudes (Table 1). These differentially expressed proteins in the serum were selected for further analysis based on the fold-change (≥1.50 or <0.67) and p-value < 0.05 (Supplemental Table S3). The serum proteins were displayed in a volcano plot (Supplemental Figure S1), and hierarchical clustering showed they could be grouped according to their differential abundance (Supplemental Figure S2). Functional analysis showed that these differentially expressed proteins were involved in immunoregulation, lipoprotein metabolism, the regulation of blood coagulation, oxygen binding and delivery, and antioxidant stresses ( Table 2).

Gene Ontology Analysis
The GO terms that were highly enriched in the 60 differentially expressed proteins were identified by DIA proteomics (Table 2). In a pairwise comparison of LA and MA buffalo, 10 GO terms were mainly related to major histocompatibility complex (MHC) protein, oxygen binding and delivery. Eight GO terms were mainly and significantly related to oxygen binding and delivery in a pairwise comparison of MA and HA buffalo. In another pairwise comparison of LA and HA buffalo, 20 GO terms were significantly related to the regulation of immune cells, lipoprotein lipid oxidation, blood coagulation and protein ubiquitination. In the three pairwise comparisons (LA vs. MA, MA vs. HA, and LA vs. HA), buffalo shared two GO terms: GO 0002376-immune system process and GO 0050896-response to environmental stress (Supplemental Figure S3, Supplemental Table S4).

Quantitative Analysis of Targeted Proteins with Parallel Reaction Monitoring
Twenty-five of the 60 differentially expressed proteins from the DIA analyses relating to immunoregulation, lipoprotein metabolism, and the regulation of blood coagulation were selected for PRM analysis. Based on the identification of peptides for which FDR < 30%, 12 target peptides of seven proteins were selected for quantitative PRM analysis. As shown in Table 3, there were 2.06-and 2.48-fold changes in the levels of alpha-1 acid glycoprotein (D2U6V0) in LA buffalo compared to MA and HA buffalo, respectively (p < 0.05). No significant differences were detected in the levels of lipoprotein lipase (E1ACW2) and beta-casein (B7VGH4) between LA and HA buffalo, but there were significantly higher levels of both proteins in LA and HA buffalo compared to MA buffalo (p < 0.05). The levels of apolipoprotein C-III (Fragment; L8IHQ3) and kappa-casein (A0A1L6BP00) were significantly higher in LA buffalo than in MA buffalo (p < 0.01), but there were no differences in the two pairwise comparisons of LA and HA buffalo or MA and HA buffalo.

Blood Physiological Parameters of Buffalo
The temperature-humidity index is widely used to evaluate the level of heat stress [18], with 72 being the threshold for dairy cows [3]. The threshold of heat stress, however, may be lower than 72 for dairy buffalo due to black skin, sparser hair, and fewer sweat glands compared to dairy cows [6,19]. LA buffalo may suffer from chronic heat stress, causing oxidative stress with an imbalance between oxidants and antioxidants. To reduce heat load from the ambient environment with the decrease in temperature gradient between internal organs and external environment [20], LA buffalo had significantly higher body temperature than MA and HA buffalo.
ROS at low levels play important roles in immune function and redox regulation [21], but excessive ROS that are not scavenged under long-term heat or oxidative stress induce lipid peroxidation and oxidative damage to protein [22][23][24]. SOD and GSH are important cellular antioxidants for scavenging ROS. In our study, increased production of ROS was detected, which suppressed the levels of serum antioxidants [25]. Lower levels of GSH and SOD were present in LA buffalo as compared to MA and HA buffalo. HA buffalo live in a comfortable thermal environment (17.5 • C and THI 62.6) and had higher levels of GSH than the LA and MA buffalo. HA buffalo also showed higher serum levels of HIF-1α, but this increase was not significantly different compared to LA and MA buffalo.

Blood Immunoglobulins of Buffalo
Heat stress suppresses the immune function and antioxidant capacity of farm animals [26,27], and increased susceptibility to diseases was observed in heat-stressed animals [28]. Complement pathways and blood coagulation are important natural antibacterial systems [29], and complement proteins play important roles in the elimination of pathogens with complement pathways [30]. We detected that the serum levels of immunoglobulins, including complement factor H, complement component C8 alpha chain, complement factor I, complement factor H-related protein 3 (Fragment), complement component C6, complement factor B, complement C5, and properdin were lower in LA and MA buffalo than in HA buffalo (Table 1). Complement factor H is a member of the complement regulator and plays an important role in complement activation pathways [31].
Most of the serum immunoglobulins were at lower levels in LA buffalo than in MA and HA buffalo, which indicated that buffalo raised at low altitudes suffered from heat stress. Alpha-1 acid glycoprotein, an acute-phase protein produced by the liver and peripheral tissues, which plays a role in immune system regulation [32], was at higher levels in LA buffalo than in MA and HA buffalo. Beta-casein, which regulates the innate and adaptive immune responses in ruminants [33], showed higher levels in LA buffalo than in MA and HA buffalo. Selenoprotein P plays an important role in the transportation of selenium and is one of the main blood antioxidants [34], and glutathione peroxidase is also an important antioxidant [35][36][37]. Both selenoprotein P and glutathione peroxidase were at higher levels in LA and HA buffalo than in MA buffalo. It was inferred that the impact of oxidative stress resulting from heat stress at low altitude or hypoxia at high altitude on the immune responses of buffaloes could be mediated with high levels of serum antioxidants.

Blood Oxygen-Delivery Proteins of Buffalo
As a metalloprotein and the most abundant haem-containing protein in red blood cells, haemoglobin is responsible for blood oxygen delivery. Haemoglobin subunit alpha-1, haemoglobin subunit alpha-2, haemoglobin foetal subunit beta, haemoglobin subunit beta-A, haemoglobin subunit alpha-I/II and haemoglobin beta were present at higher levels in LA and HA buffalo than in MA buffalo. Unexpectedly, there were no differentially abundant proteins or significantly enriched GO terms related to blood oxygen delivery when comparing LA and HA buffalo ( Table 2). Buffalo raised at low altitudes with high ambient temperature and relative humidity in summer may suffer from chronic heat stress, and HA buffalo adapt to hypoxia with highly efficient blood oxygen delivery. Heat stress induces local tissue hypoxia due to blood flow redistribution for increased skin heat dissipation [38][39][40], and it was reasonable that no significantly enriched GO terms related to blood oxygen delivery were detected between LA and HA buffalo.
HA buffalo showed higher levels of haemogen, haemopexin, lactoferrin, and serotransferrin, which possibly increased the amount of blood oxygen to meet the physiological demands. HA buffalo also had higher levels of coagulation factors V and IX compared to LA buffalo. Kallikrein and plasminogen were at higher levels in HA buffalo than in LA buffalo, which have roles in the renin-angiotensin system by converting prorenin into renin and in blood coagulation [41].
Vitamin K-dependent protein S, which circulates in the blood, functions as a cofactor and is important for the inactivation of coagulation factors V and VIII on membrane surfaces [42]. Buffalo raised at high altitudes had higher levels of vitamin K-dependent protein S with anticoagulant properties [43,44], and this higher abundance of vitamin K-dependent protein S may play a role in increasing blood flow to meet the blood oxygen supply of dairy buffalo raised at high altitudes. Heat stress induces oxidative stress and hypoxia [39]. To satisfy the requirements of oxygen delivery for internal organs, an increase in blood flow is necessary for dairy buffalo under chronic heat stress. A higher abundance of vitamin K-dependent protein Z was detected in LA and HA buffalo compared to MA buffalo. As a cofactor, vitamin K-dependent protein Z inhibits coagulation factor Xa on phospholipid surfaces to avoid thrombosis and the coagulation response, thereby promoting smooth blood flow [45]. This indicated that buffalo adapt to hypoxia at high altitudes and heat-stress-induced hypoxia at low altitudes through a high blood oxygen delivery capacity.

Serum Apolipoproteins of Buffalo
Apolipoproteins are involved in the combination and transportation of lipids [46]. A lower abundance of apolipoprotein C-IV and apolipoprotein F (Fragment) was detected in LA and HA buffalo than in MA buffalo. The trend for phospholipid transfer protein and proactivator polypeptide, which are involved in lipoprotein metabolism and the regulation of lipid metabolic processes, was similar [47]. However, HA buffalo were higher levels of apolipoprotein D, which regulates the protection from oxidative stress and increases stress resistance [46,48]. The HA buffalo also had higher levels of vitamin D binding protein than LA and MA buffalo. Vitamin D binding proteins are the major carriers of vitamin D [49,50], which is involved in the regulation of fatty acid β-oxidation and energy metabolism [51], inhibiting the activation of stress-activated protein kinases [52].

Conclusions
Our study indicated that low-altitude buffalo adapted to heat stress-induced oxidative stress with high levels of alpha-1 acid glycoprotein and apolipoproteins, whereas highaltitude buffalo adapted to hypoxia-induced oxidative stress with high complement protein levels for immunological enhancement. In summary, buffalo adopted similar adaptation strategies to oxidative stress induced by heat stress or hypoxia with high levels of serum antioxidants, immunological enhancement, the inhibition of lipid oxidation, and high efficiency of blood oxygen delivery. The research findings may be beneficial to improve the health and welfare of buffaloes with feed additives that enhance immune function and cutaneous circulation.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/metabo12100909/s1, Figure S1: Volcano plot of the proteins identified as differentially abundant in the three pairwise comparisons of serum; Figure S2: Hierarchical clustering of serum proteins identified as differentially abundant; Figure S3: Pairwise comparisons of the gene ontologies (GO) of the serum proteins between the three farms; Table S1: Nutrient content of feed samples (%);