**3. Discussion**

#### *3.1. Similarities between BHBA and NEFA*

Not surprisingly, many of the metabolites identified as having common co-variance with both BHBA and NEFA concentrations are involved in hepatic energy metabolism. These relationships are summarized in Figure 5. Most obvious was the negative relationship between both biomarkers and glucose. Hypoglycaemia has been widely reported in early lactation dairy cows due to the massive demand for glucose for lactogenesis [3,38]. More recently, NMR metabolomics studies have identified serum glucose concentration as being (1) directly correlated to energy balance (r = 0.84) [39], and (2) lower in cows with clinical and subclinical ketosis [25] and fatty liver [27] when compared to healthy controls. Our results offer further evidence of the pivotal role glucose plays in the early lactation metabolic health in dairy cows.

**Figure 5.** Summary of hepatic energy metabolism in early lactation dairy cows. Arrows indicate the direction of the relationship between the metabolites and the reference BHBA (blue) and non-esterified fatty acid (NEFA) (red) concentrations. BHBA = β-hydroxybutyrate; OAA = oxaloacetate; TAG = triglyceride, TCA = tricarboxylic acid, VLDL = very low density lipoprotein.

Lactate and alanine, important gluconeogenic substrates in ruminants [40,41], were also negatively associated with both BHBA and NEFA, as was valine (another gluconeogenic amino acid). Interestingly, Xu et al. [39] found no correlation between calculated energy balance in early lactation dairy cows and the concentrations of any of the branched-chain amino acids or lactate. Conversely, when compared to healthy controls, cows with fatty liver and displaced abomasa have been shown to have lower serum alanine concentrations [27,29], and cows with ketosis have lower lactate and alanine concentrations [25,42]. This suggests that alterations in glucogenic precursors, in particular lactate and alanine, are indicative of a perturbed metabolism, not simply negative energy balance. We previously showed that lactate concentration in pasture-fed dairy cows is heavily influenced by herd-specific managemen<sup>t</sup> factors [43], and as such may not be heavily influenced by genetic factors. Alanine has been shown to be the most important glucogenic amino acid, and the most important gluconeogenic precursor after lactate and propionate, in dairy cows [41]. Therefore, genetic selection for cows with higher serum concentrations of alanine in early lactation may help to increase endogenous glucose supply.

Spectral features attributed to VLDL and LDL were positively correlated with the concentrations of both BHBA and NEFA. These results need to be interpreted with caution as the methanol extraction used in this study removed much of the protein from the samples and may have introduced experimental artefacts. Interestingly, 1H NMR spectroscopy has recently been shown capable of providing high-throughput and accurate quantification of lipoprotein subclasses in human serum and plasma samples [32,44]. It is important to note that these protocols used di fferent pulse sequences and involved the dilution of plasma/serum in a deuterated water/phosphate bu ffer solution without any metabolite extraction, such as the one used in our study. The findings of these studies cannot, therefore, be applied directly to our results. However, lipoprotein metabolism is central to early lactation health in dairy cows, and impaired VLDL production in the liver can result in hepatic triglyceride (TAG) accumulation (Figure 4) and the development of fatty liver [45]. Dyslipoproteinaemia is also an important feature of metabolic syndrome in humans, and the quantification of lipoprotein subclasses is considered critical to the better understanding of this disease [44]. We believe that the investigation of serum lipoproteins using 1H NMR spectroscopy holds grea<sup>t</sup> promise in the research of early lactation metabolic health in dairy cows, and we plan to validate the aforementioned protocols on bovine serum and plasma samples.

The region of the spectrum associated with glycoproteins was also significantly positively correlated with both NEFA and BHBA concentrations. Glycoproteins are acute phase proteins which can be used as indicators of inflammation in cattle [46]. In dairy cattle, increased serum NEFA concentrations in early lactation are associated with uncontrolled inflammation, and this inflammatory dysfunction is hypothesized to be a central link between metabolic and infectious disorders [14,47]. 1H NMR spectroscopy is showing promise for the quantification of glycoprotein A (GlcA) in human research into metabolic diseases such obesity, diabetes mellitus and the metabolic syndrome [33]. Given that these syndromes have much in common with early lactation metabolic disease in dairy cows (e.g., insulin resistance), we believe that further research into GlcA as a biomarker for early lactation health is warranted. Overall, our results o ffer further evidence that inflammation plays an important role in early lactation metabolic health of dairy cows.

Glycine was positively correlated with the concentrations of both BHBA and NEFA. Metabolomics studies comparing healthy and ketotic dairy cows have reported (1) no change in glycine concentrations [25], (2) increased glycine concentrations in cows with sub-clinical ketosis [26], (3) increased glycine concentrations in cows with clinical ketosis [48] and (4) decreased glycine concentrations in cows with clinical ketosis [26] and fatty liver [49]. Glycine concentration has also been shown to increase in response to lipolysis [50]. These di ffering results sugges<sup>t</sup> that changes in glycine concentration may be dependent on the severity of the metabolic disorder (i.e., increased in mild cases, and decreased in more severe cases). Most interesting are the findings of a recent metabolomics study that showed that glycine concentrations in plasma and milk were strongly negatively correlated with energy balance in early lactation dairy cows (r = −0.80 and r = −0.79, respectively) [39]. The authors of this study hypothesized that this relationship was due to an increase in one-carbon or methyl donor metabolism, specifically an increase in the conversion of choline to glycine. Given that all cows in our study were clinically healthy, our results are consistent with glycine being an indicator of negative energy balance, lipolysis, and/or sub-clinical ketosis. Further work is required to better understand the role of glycine metabolism in clinical metabolic disease.

The positive correlations between phosphocholine and both BHBA and NEFA concentrations, and between betaine and BHBA concentration, are consistent with an increase in methyl donor metabolism in cows experiencing negative energy balance. Methyl donor metabolism and nutrition are receiving a grea<sup>t</sup> deal of attention in dairy science due to links with early-lactation cow health (including fatty liver), milk production and immune function [51]. Betaine, phosphocholine and glycine are intermediates in several important one-carbon metabolic pathways including the folate and methionine cycles, and the cytidine diphosphate (CDP)–choline pathway [51] (Figure 6a). The positive correlation between NEFA and phosphocholine may be due to increased breakdown of phosphatidylcholine (Figure 6a). This is consistent with the findings of Imhasly et al. [52] who showed that serum concentrations of lyso-phosphatidylcholines and phosphatidylcholines increase in response to negative energy balance

in post-partum dairy cows. The positive association observed between betaine and BHBA could be due to increased oxidation of choline. A detailed description of these pathways is beyond the scope of this study, however our results sugges<sup>t</sup> that methyl donor metabolism has an important influence on both BHBA and NEFA concentrations in early-lactation dairy cows.

**Figure 6.** Summary of (**a**) phospholipid and one-carbon/methyl donor metabolism [53,54], and (**b**) creatine metabolism in early lactation dairy cows. Arrows indicate the direction of the relationship between the metabolites identified using untargeted 1H NMR metabolomics, and reference BHBA (blue) and non-esterified fatty acids (NEFA) (red) concentrations. ADP = adenosine diphosphate; ATP = adenosine triphosphate; DMG = dimethylglycine; NEB = negative energy balance.

#### *3.2. Di*ff*erences between BHBA and NEFA*

Despite many similarities, we observed some significant differences between the metabolomic fingerprints of BHBA and NEFA. Most obvious was the difference in the direction of correlation between acetate and the two biomarkers. Acetate is a volatile fatty acid produced by microbial fermentation of feedstuffs in the rumen, and is an important energy source [55] (via oxidation or the partial oxidation of acetyl-CoA in the liver) and substrate for de novo milk fat synthesis [56] in cows. The negative relationship we observed between acetate and NEFA is consistent with the findings of Bielak et al. [57], who reported a negative correlation (r = 0.44) between plasma NEFA and acetate concentrations in early lactation dairy cows, possibly due to the down-regulation of the active transport of acetate across the rumen wall. The positive association between acetate and BHBA is consistent with previously discussed metabolomic studies of ketosis and fatty liver [25,27]. These results sugges<sup>t</sup> that differences in acetate metabolism may help to explain the weak correlation between serum BHBA and NEFA concentrations in early lactation dairy cows.

The positive correlation between creatine and BHBA concentration is consistent with previous reports that creatine is a potentially useful biomarker of ketosis and severe energy deficiency in dairy cows [25,26,39]. Creatine is an important intermediate in energy metabolism, and this result may represent increased breakdown of creatine phosphate in skeletal muscle and the release of high-energy phosphate for the conversion of adenosine diphosphate (ADP) to adenosine triphosphate (ATP) (Figure 6b). Interestingly, creatine concentration was negatively correlated with NEFA concentration (albeit weakly and non-significantly (VIP < 1)). That mobilization of energy from skeletal muscle is a feature of the BHBA metabolomic fingerprint, but not that of NEFA, suggests that elevated BHBA concentrations are indicative of a more severe energy deficiency than are elevated NEFA concentrations. However, the ability to rapidly mobilize energy from skeletal muscle may be advantageous to early-lactation dairy cows, and we believe the role of creatine metabolism in transition cow health warrants further investigation. We therefore plan to undertake genome-wide association studies to better understand the genetic relationships between hepatic and skeletal muscle energy metabolism.

The significant negative correlation between DMSO2 and BHBA concentration was an interesting finding of this study. DMSO2 concentration in the milk and rumen fluid of dairy cows has been shown to vary according to feeding system; higher in pasture-fed cows managed outdoors than in cows fed a total mixed ration indoors [58]. Maher et al. [59] showed that the concentrations of DMSO2 in milk and plasma are highly correlated (r = 0.69), so serum DMSO2 may also be an indicator of pasture intake. Given that all animals in this experiment were fed pasture, the negative association we observed between DMSO2 and BHBA concentration may indicate that hyperketonemic cows are consuming less feed.

#### **4. Materials and Methods**

All procedures undertaken in this study were conducted in accordance with the Australian Code of Practice for the Care and Use of Animals for Scientific Purposes (National Health and Medical Research Council, 2013). Approval to proceed was granted by the Agricultural Research and Extension Animal Ethics Committee of the Department of Jobs, Precincts and Regions Animal Ethics Committee (DJPR, 475 Mickleham Road, Attwood, Victoria 3049, Australia), and the Tasmanian Department of Primary Industries, Parks, Water and Environment (DPIPWE Animal Biosecurity and Welfare Branch, 13 St Johns Avenue, New Town, Tasmania 7008, Australia). AEC project approval code 2017-05.

#### *4.1. Animals and Datasets*

A total of 298 Holstein-Friesian cows were used in this experiment. The calibration dataset (Dataset 1) was collected between August and September 2017 from 248 animals located at the Ellinbank Dairy Research Centre, Ellinbank, Victoria, Australia. An independent validation dataset (Dataset 2) was collected in September 2018, from 50 cows located on a commercial dairy farm in Smithton, Tasmania, Australia. All cows were clinically healthy, and had been calved for between 4 and 30 days at the time of sampling. Feeding systems on Australian dairy farms are diverse but can be classified into 5 main feeding systems [60]. Both farms operated under feeding system 2; grazed pasture plus moderate to high level concentrate feeding (>1.0 tonne of concentrate fed in the milking parlour per cow per year).

#### *4.2. Blood Sample Collection and Reference Biomarker Measurements*

A single serum sample was taken from each cow immediately after morning milking (approximately 07:00) according to the protocol described in Luke et al. [43]. Cows were fed their concentrate ration as soon as they entered the milking parlour, meaning that samples were collected approximately 10 min after grain feeding.

An aliquot of each serum sample was shipped on ice to Regional Laboratory Services (Benalla, Victoria, Australia) for BHBA and NEFA analyses. Colorimetric assays were performed using a Kone 20 XT clinical chemistry analyser (Thermo Fisher Scientific, Waltham, MA, USA); an enzymatic kinetic assay for BHBA (McMurray et al., 1984) and enzymatic end point assay for NEFA (Randox Laboratories, Crumlin, UK). The uncertainty of measurement (at a 95% confidence level) was ±0.060 mmol/L at 0.85 mmol/L for BHBA, and ±0.031 mM at 1.45 mM for NEFA. A second aliquot was stored at −20 ◦C until processing for NMR spectroscopy.

#### *4.3. Sample Preparation for NMR Spectroscopy*

Details of the sample preparation and metabolite extraction protocols used in this study can be found in Luke et al. [43]. Briefly, 300 μL of serum was (1) mixed with 600 μL of methanol, (2) vortexed, (3) incubated at −20 ◦C for 20 min, and (4) centrifuged at 11,360× *g* at 21 ◦C for 30 min to pellet proteins. A 600 μL aliquot of the supernatant was then transferred to a clean 2 mL microcentrifuge tube, dried under vacuum at 21 ◦C overnight using a SpeedVac Savant SPD 2010 Concentrator (Thermo Fisher Scientific, Waltham, MA, USA) then reconstituted in a D2O phosphate bu ffer solution (100 mM K2HPO4) containing 0.25 mM DSS-d6 as an internal standard. A 550 μL aliquot of reconstituted extract was transferred to a 5 mm NMR tube for analysis.

#### *4.4. 1H NMR Data Acquisition and Pre-Processing*

One-dimensional proton spectra were acquired using a Bruker Ascend 700 MHz spectrometer equipped with cryoprobe and SampleJet automatic sample changer (Bruker Biospin, Rheinstetten, Germany). A Bruker noesypr1d pulse sequence was used over a −0.76–10.32 ppm spectral range with the following acquisition parameters; (1) a temperature of 298 K, (2) 256 scans after eight dummy scans (3) acquisition time per increment of 2.11 s, and (4) relaxation delay (D1) of 2.00 s. This resulted in 32,768 data points. A line broadening of 0.3 Hz was applied to all spectra prior to Fourier transformation. Spectra were manually phased, baseline corrected and referenced to the internal standard (DSS-d6) at δ 0.00 ppm using the Topspin v.3.6.1 software (Bruker Biospin, Rheinstetten, Germany).

Data pre-processing was performed in MatLab v.R20017b (Mathworks, Natick, MA, USA). Spectra were imported as a matrix of signal intensities using the ProMetab v.1.1 script [61]. Spectral pre-processing involved (1) deletion of the residual water peak region (δ 4.68–5.00 ppm), (2) spectral alignment using the correlation optimized warping algorithm [62], (3) normalization to total signal area (area = 1), (4) deletion of methanol (δ 3.32–3.36 ppm) and DSS-d6 (δ 0.4–0.60 ppm) peak regions, and the non-informative region beyond 9.00 ppm, (5) baseline correction using automatic weighted least squares and (6) mean centering.

## *4.5. Statistical Analysis*

Statistical analysis of experimental metadata was performed in R v3.6.2 [63]. Di fferences between the 2 datasets were analysed using a paired *t*-test or a Wilcoxon signed-rank test depending on the normality of the data.

Multivariate statistical analyses were performed using the PLS Toolbox v. 8.5.2 (Eigenvector Research Inc., Manson, WA, USA). Preliminary data analysis and outlier detection was performed using an unsupervised PCA. Examination of PC1 vs. PC2 scores plot showed 14 samples from Dataset 1 outside the 95% confidence level ellipse (Figure S1). These samples were individually examined, and a single spectrum with poor water suppression and baseline correction was removed from subsequent analyses. The influences of fixed e ffects (DIM, age and herd) on spectra were investigated using ANOVA simultaneous component analysis with 1000 permutations [64]. Untargeted metabolomic fingerprinting was done by regressing reference NEFA and BHBA concentrations against 1H NMR spectra using supervised OPLS regression. Variable importance of projection (VIP) scores for the first latent variable were used to identify the most statistically significant peaks in each model. Peaks of interest were identified using the Chenomx NMR suite software v.8.4 (Chenomx Inc., Edmonton, AB, Canada), comparison to the literature, 2D NMR analysis (COSY, gHMBC and gHSQC), and statistical total correlation spectroscopy [65].

OPLS models were constructed using data from Dataset 1. The robustness of models was assessed using (1) cross-validation using a venetian blind technique (10 sample splits with 1 sample per blind) and (2) external validation using data from Dataset 2. The prediction accuracy of OPLS models was assessed using the coe fficient of determination (R2) and root mean square error (RMSE). Normalized

RMSE (NRMSE) values, calculated as external validation RMSE divided by the interquartile interval of the observed data, were used to compare RMSE estimates for NEFA and BHBA predictions. Permutation testing (50 iterations and statistical significance determined using a Wilcoxon signed-rank test) was performed to ensure that models were not over-fitted.
