Heat stress modifies the lactational performances and the urinary metabolomic profile related to gastrointestinal microbiota of dairy goats

The aim of the study is to identify the candidate biomarkers of heat stress (HS) in the urine of lactating dairy goats through the application of proton Nuclear Magnetic Resonance (1H NMR)-based metabolomic analysis. Dairy does (n = 16) in mid-lactation were submitted to thermal neutral (TN; indoors; 15 to 20°C; 40 to 45% humidity) or HS (climatic chamber; 37°C day, 30°C night; 40% humidity) conditions according to a crossover design (2 periods of 21 days). Thermophysiological traits and lactational performances were recorded and milk composition analyzed during each period. Urine samples were collected at day 15 of each period for 1H NMR spectroscopy analysis. Principal component analysis (PCA) and partial least square—discriminant analysis (PLS-DA) assessment with cross validation were used to identify the goat urinary metabolome from the Human Metabolome Data Base. HS increased rectal temperature (1.2°C), respiratory rate (3.5-fold) and water intake (74%), but decreased feed intake (35%) and body weight (5%) of the lactating does. No differences were detected in milk yield, but HS decreased the milk contents of fat (9%), protein (16%) and lactose (5%). Metabolomics allowed separating TN and HS urinary clusters by PLS-DA. Most discriminating metabolites were hippurate and other phenylalanine (Phe) derivative compounds, which increased in HS vs. TN does. The greater excretion of these gut-derived toxic compounds indicated that HS induced a harmful gastrointestinal microbiota overgrowth, which should have sequestered aromatic amino acids for their metabolism and decreased the synthesis of neurotransmitters and thyroid hormones, with a negative impact on milk yield and composition. In conclusion, HS markedly changed the thermophysiological traits and lactational performances of dairy goats, which were translated into their urinary metabolomic profile through the presence of gut-derived toxic compounds. Hippurate and other Phe-derivative compounds are suggested as urinary biomarkers to detect heat-stressed dairy animals in practice.

Introduction Exposure to high ambient temperature induces several physiological responses in order to maintain body homeostasis. Animals suffer from heat stress (HS) when physiological mechanisms fail to counterbalance an excessive heat load [1]. Exposure of dairy animals to HS results in a decline in their productive [2] and reproductive [3] performances due to a strong metabolic disruption. Dairy animals under HS typically show decreased feed intake, increased water consumption and altered thermophysiological traits, such as respiratory rate and rectal temperature, when compared to thermoneutral (TN) ones. Usually, HS reduces milk yield and impairs milk composition in dairy goats [4]. Although these negative effects on milk production are traditionally attributed to a decline in feed intake, pair-fed TN experiments have shown that intake only accounts for 35 to 50% of milk yield reduction in dairy cows [5,6]. Therefore, there is a specific effect of HS that disrupts body metabolism and milk secretion which remains unknown.
Bio-fluid assessment by Nuclear Magnetic Resonance (NMR) spectroscopy can shed some light on the physiological mechanisms that occur in animals when exposed to HS. Proton ( 1 H) NMR, together with multivariate statistical analysis, has been successfully used as a metabolite profiling method to study the metabolic changes in blood [7], milk [8] and liver [9] of HS dairy cows, as well as in plasma of HS growing pigs [10] and rats [11]. This robust and reliable technique provides vast information on metabolome dynamics and metabolic pathways [12]. The 1 H NMR spectra are derived from thousands of metabolite signals that usually overlap, adding complexity to data processing. Computer-based data reduction and multivariate statistical pattern recognition methods, such as principal component analysis (PCA) and partial least square-discriminant analysis (PLS-DA), have been shown to be beneficial techniques to get the most from the information obtained in the 1 H NMR spectra for classification purposes [13,14].
To our knowledge, no studies have been carried out to evaluate urine metabolomics of dairy goats. The aim of this study is to identify the candidate biomarkers of HS through the application of 1 H NMR-based metabolomic urinalysis of dairy goats.

Animals and treatments
Animal care conditions and management practices of the study were approved by the Ethical Committee of Animal and Human Experimentation (CEEAH Approval No. 09/771) of the Universitat Autonoma of Barcelona (UAB) and agreed the codes of recommendations for livestock wellbeing of the Ministry of Agriculture, Food and Environment of Spain.
Sixteen multiparous Murciano-Granadina dairy does (43.5 ± 1.6 kg body weight), lactating and open, from the herd of the SGCE (Servei de Granges i Camps Experimentals) of the UAB in Bellaterra (Barcelona, Spain), were blocked into 2 balanced groups at mid-lactation (81 ± 3 days-in-milk; 2.00 ± 0.04 L/day). Does were adapted to metabolic cages for 2 weeks before the start of the experiment and the groups randomly allocated to 2 ambient-conditions treatments according to a 2 × 2 (treatment × period) crossover design. There were two 21-day experimental periods (14 days for adaptation, 5 days for measurements, and 2 days for washout) during which both treatments were sequentially applied to each doe. As a result, a total of 16 observations per variable were obtained for each treatment. Treatments were TN (indoor shelter; 15 to 20˚C and 45 ± 5% relative humidity) and HS (climatic chamber 4 m × 6 m × 2.3 m with temperature-humidity control system; Carel Controls Ibérica, Barcelona, Spain; 37 ± 0.5˚C during the day, and 30 ± 0.5˚C during the night; 40 ± 5% humidity and 90 m 3 /h continuous air turnover). Day-night length was set to 12-12 hours. Temperature-humidity index (THI), calculated according to NRC in 1971 [15], resulted in THI TN = 59 to 65 and THI HS = 75 to 83. Experimental conditions were similar to those detailed by Hamzaoui et al. [16].

Sampling and measurements
Thermophysiological traits and lactational performances of the goats. Does were weighed at the start and the end of each period using an electronic scale (True-Test SR2000; Pakuranga, New Zealand; accuracy, 0.2 kg). Rectal temperature (digital clinical thermometer, ICO Technology; Barcelona, Spain; accuracy, 0.1˚C) and respiratory rate (flank movements during 60 s) were recorded daily at 0800, 1200, and 1700 hours throughout the experiment. Milk yield (volume) was recorded daily throughout the experiment, and milk samples were collected weekly for composition (NIRSystems 5000, Foss; Hillerød, Denmark). Feed and water intakes were calculated by weight from the daily refusals and feed samples were collected daily and composited by period for analyses. Feed composition was determined according to analytical standard methods [17].
Urine sampling and preparation. Urine samples from each doe were collected at micturition on the morning of day 15 of each period (n = 32) and stored at −20˚C until 1 H NMR analysis.
Preparation of samples for 1 H NMR spectroscopy was done according to Beckonert et al. [12]. Briefly, a phosphate buffer solution (pH 7.4) was prepared with sodium phosphate dibasic (Na 2 HPO 4 ; 99.95% trace metals basis, anhydrous, Sigma-Aldrich Merck; Darmstadt Germany), sodium phosphate monobasic (NaH 2 PO 4 ; 99.95% trace metals basis, anhydrous, Sigma-Aldrich Merck) and sodium azide (NaN 3 ; Sigma-Aldrich Merck). Deuterium oxide (D 2 O; 99.9 atom % D, Sigma-Aldrich Merck), containing 0.75% 3-(trimethylsilyl) propionic-2,2,3,3-d4 acid (TSP) sodium salt as the NMR reference compound, was added before the flask was filled up to 25 mL with milli-Q water (EMD Millipore; Darmstadt, Germany). The flask was shaken thoroughly and left in a Clifton sonicator (Nickel Electro; Weston-super-Mare, United Kingdom) at 40˚C until the salts were dissolved. The prepared phosphate buffer solution was stored at 4˚C. Urine samples were thawed in a water bath, thoroughly shaken and spun for 5 min at 12,000 × g in a swing-bucket rotor (Hettich; Tuttlingen, Germany) at 4˚C. Then, 400 μL of the urine sample were transferred into Eppendorf tubes and mixed thoroughly with 200 μL of cold phosphate buffer solution. All the tubes were then centrifuged for 5 min at 12,000 × g at 4˚C and 550 μL of the final mixture transferred into 5-mm NMR tubes (VWR International Eurolab; Barcelona, Spain). The prepared NMR tubes were immediately put on ice and sent to the NMR Service of the UAB for 1 H High Resolution NMR Spectroscopy. NMR spectroscopy. 1 H NMR spectra were acquired on a Bruker Avance-III spectrometer (Bruker BioSpin; Rheinstetten, Germany) operating at a 1 H NMR frequency of 600 MHz and a temperature of 298˚K, controlled by a Burner Control Unit-extreme regulator. A 5-mm Triple Resonance Broadband Inverse probe with z-gradients and inverse detection was used and controlled by TopSpin2.1 software (Bruker, Germany). One-dimensional 1 H NMR spectra were obtained using a one-dimensional Nuclear Overhauser Enhancement Spectroscopy (NOESY) pulse sequence. The solvent signal was suppressed by pre-saturation during relaxation and mixing time. A total of 32 scans and 2 dummy scans were performed to produce 32,768 data points for each spectrum using a relaxation delay of 2.0 s with a pulse power level of 54 dB and an acquisition time of 2.6 s. Spectral width (δ) used for all data collected was 12.0 ppm, and 0.3 Hz exponential line-broadening was applied for the Fourier Transform of the raw data. 1 H NMR spectra were phased, baseline corrected, and corrected for chemical shift registration relative to the TSP reference compound previously indicated (δ = 0.0 ppm) in TopSpin 2.1.(data in S1 Dataset).

Statistical analyses
Thermophysiological and performance analysis. Data were analyzed by the PROC MIXED for repeated measurements of SAS v. 9.1.3 (SAS Inst. Inc.; Cary, North Carolina, USA). The statistical mixed model contained the fixed effects of environmental treatment (TN vs. HS), the period (1 and 2) and measuring day (1 to 19), the random effects of the animal (1 to 16), the interactions (treatment × day and treatment × period), and the residual error. Differences between least squares means were determined with the PDIFF option of SAS. Significance was declared as P<0.05.
NMR data pre-processing and analysis. Pre-treatment of raw spectral data is critical for generating reliable and interpretable models using multivariate analysis techniques. Nevertheless, metabolic fingerprinting datasets acquired from 1 H NMR spectrometers suffer from imprecisions in chemical shifts due to temperature, pH, ionic strength and other factors [18]. Therefore, models generated using multivariate analysis may fail to identify separations between classes, and their loadings can be difficult to interpret due to an over-abundance of variables. To mitigate these complications, each spectrum was uniformly divided into 'bins' of 20 signals, and the signal intensities within each bin were integrated to produce a smaller set of variables (i.e., from 0.0003 to 0.007 ppm) using R software v. 3.2.3 [19]. After binning, alignment and normalization of spectra were performed to ensure that all observations were directly comparable. In this sense, urine spectra were normalized to creatinine methyl resonance intensity at δ = 3.05 ppm and then log 2 transformed. Regarding variable selection, raw 1 H NMR spectral data were edited by excluding both the regions outside the chemical shift range of δ = 9.0-0.5 ppm and the residual peak of the imperfect water suppression (δ = 5.5-4.6 ppm). Following the recommendations of Pechlivanis et al. [20], the spectral regions of histidine, 1-metylhistidine, and 3-methylhistidine (δ = 8.17-7.87, δ = 7.15-7.01, and δ = 3.77-3.71 ppm, respectively) were also removed because of the sensitivity to small pH differences among urine samples.
Once 1 H NMR pre-processing data were completed, data were subjected to multivariate statistical analysis. Initially, PCA was performed without considering the class information for samples examination and search for outliers. Then, PLS-DA with leave-one-out cross-validation was also performed on the datasets using the pls package of R software [21]. PLS-DA allowed individual samples to be classified according to the respective class prior to analysis (TN or HS). Model strength was assessed using both R 2 and Q 2 statistical parameters. While R 2 values reported the total amount of variance explained by the model, the Q 2 reported model accuracy as a result of cross-validation. Aside from its theoretical maximum value of 1, for biological models, an empirically inferred acceptable value is � 0.4 [13]. The resulting Q 2 statistic was compared to a null distribution to test model significance (P < 0.05).
Interpretation of multivariate analysis was performed through scores and loadings plots according to their contribution to the separation between groups. For biomarker searches, PLS-DA loadings greater than |0.0005| were selected according to their absolute magnitude values. Consequently, metabolites responsible for the separation between experimental groups were those with the highest values. Moreover, a Volcano plot with paired Student t test analysis between HS over TN cohorts was performed to get a general overview of the data (log 2 fold change thresholds, �1.5 and �1.5; P<0.01) and to identify metabolites with a significant effect. The false discovery rate (FDR) for differentially excreted metabolites was controlled according to the Benjamini-Hochberg procedure [22] with an adjusted P<0.05. Volcano plots are suitable as complementary analysis because both PCA and PLS-DA analysis may be influenced by variable correlations and the intra-and inter-class variance of metabolites may have no significant differences in the univariate analysis [23]. All 1 H NMR data pre-processing, statistical analysis and the generated plots were performed using R.
Metabolite assignment. Chemical shifts linked to the highest loading values found in PLS-DA were annotated for metabolite assignment as HS biomarker candidates. The candidate chemical shifts and corresponding metabolites were assigned using the Human Metabolome Database [24] and queried in KEGG (Kyoto Encyclopedia of Genes and Genomes) database to know in which metabolic pathways they were involved.

Effects of heat stress on thermophysiological and lactational performances of the goats
The effects of the experimental HS conditions on thermophysiological and lactational performances of the dairy goats are summarized in Table 1. Rectal temperature and respiratory rate increased during the day in both groups of does, following the expected circadian rhythm and the daily THI pattern in both TN and HS conditions. The greatest values were observed in the HS does at 1700 hours, the increases being 1.2˚C and 3.5-fold (P<0.001), when compared to TN does. On average, feed intake decreased 35% in HS (P<0.001), when compared to TN does but, in contrast, water consumption increased 74% (P<0.001). Furthermore, HS does lost 115 g/d of body weight, whereas TN goats gained 162 g/d, on average (P<0.001). Results obtained agreed with those reported for the same breed of dairy goats in late-lactation and under similar HS conditions [16].
Reducing feed intake is a way to decrease heat production in warm environments because heat increment of feeding, especially in ruminants, is an important source of heat production [25]. Moreover, increased water consumption under HS conditions is mainly used for boosting latent heat losses by evaporation (e.g., sweating and panting). Despite this, no differences in milk yield were observed, although milk composition markedly worsened. Milk fat, protein and lactose contents varied by −9%, −16% and −5%, respectively (Table 1; P<0.01), which would severely compromise the milk transformation into dairy products [4]. Consequently with the decrease in the content of milk components, fat-corrected milk yield also varied by −14% (P<0.05).
Although our does were less sensitive to HS than were dairy cows, with regard to feed intake and milk yield, the effects of HS on milk fat content and fat-corrected milk were contradictory when compared to cows. So, despite the typical fat depression seen in commercial dairy-cow farms during the summer, Rhoads et al. [5] and Shwartz et al. [26] reported a 9% increase or no change in milk fat content, in the short-or mid-term, respectively, in HS vs. TN dairy cows. On the other hand, the above-indicated negative effect of HS on the milk protein content of our goats (i.e., −16%) was greater than that reported by Rhoads et al. ( [5], −5%) and Shwartz et al. ( [26], −9%) in dairy cows, and also in the same breed of dairy goats in late lactation ( [16], −13%). The negative effects of HS on the lactational performances of dairy ruminants are usually attributed to the decline in feed intake, but pair-fed experiments under TN conditions have shown that feed intake only explains approximately half of the fall in milk yield and body weight in dairy cows [5,6]. Consequently, the other half should be explained by unknown mechanisms induced by HS. Therefore, similar responses were expected in our dairy goats.
As an intermediate conclusion, the thermophysiological and lactational performance responses observed clearly demonstrated that our HS does (kept at THI = 75 to 83) were under severe stress on the days at which the urine samples for 1 H NMR-metabolomics assessment were collected (day 15).

NMR urinary spectroscopy of the goats
A comparison of 1 H NMR urinary mean spectra for the TN and HS lactating does is shown in Fig 1. Resonance assignments reported in the figure were made from the known chemical shifts and coupling patterns of urine spectra previously described in humans [18,27].
At first glance, visible differences in urine metabolites were found between HS and TN groups. The spectral region from δ = 8.0-6.5 ppm showed higher excretion compounds in the HS doe group. On the contrary, all excreted compounds that lay on the δ = 4.5-0.5 ppm spectral region appeared to be at lower concentrations in the HS group. More detailed analyses of metabolic differences between these two thermal conditions were obtained from the multivariate PCA and PLS-DA data analyses and the Volcano plot. First, the Volcano plot (Fig 2) showed that TN does excreted a greater number of urinary metabolites (i.e., a higher number of left-sided spots) than did HS does. Most probably, this was a consequence of the metabolic sparing of nutrients of the HS does, which lost weight as a result of their negative energy balance, to cope with the HS conditions.
Regarding the multivariate analysis, PCA was initially applied to the 1 H NMR spectra. Based on the principle of minimum differentiation, no samples were identified as outliers according to Hotelling's T 2 (95% interval of confidence). Therefore, all samples remained for subsequent PLS-DA in order to identify the metabolic differences between HS and TN dairy does. The PLS-DA scores plot showed a slight distinguishable separation between HS and TN datasets (Fig 3).
The separation along the x-axis (PLS-DA component 1) represents differences related to environmental treatment. All other variations in the NMR data are visualized as separation in the y-axis direction (second component). The cross-validation of urine metabolomics PLS-DA models (first 2 components) gave R x 2 = 0.54, R y 2 = 0.17, and Q 2 = 0.47. The R 2 and Q 2 values in the model were higher than those in the random model (P<0.01). Although the top-ranking metabolites responsible for discriminating HS does were related to gut-derived uremic toxins or mammalian-microbial cometabolites (i.e. hippurate, OH-phenylacetate, OH-phenylacetylglycine, phenylglyoxylate and trimethylamine N-oxide), the thresholds applied for Volcano Plot, allowed to identify a total of 15 metabolites as candidates for urine biomarkers in HS does (Table 2). Thus, by-products of autophagy (i.e. 3-methyladenine) and energy reservoirs for muscle contraction (i.e. phosphocreatine) were also overexcreted under HS conditions. On the other hand, 8 metabolites where detected as underexcreted, some of them related to vitamin metabolism (i.e. cholecalciferol, pyridoxal, β-alanine) and carbohydrate metabolism (i.e. glycogen, galactitol) among others. The increase of gut-derived uremic toxins reflected alterations in the gastrointestinal environment due to the metabolic impact of HS. In fact, it is well known that under HS conditions, mammals redistribute blood to the periphery for heat dissipation purposes, while vasoconstriction occurs in the gastrointestinal tract [28] that leads to tissue hypoxia and oxidative stress [29]. Moreover, lower rumen pH has been reported as a side-effect in HS goats [30] that leads to an abnormal overgrowth of gastrointestinal microbiota, a compromised integrity [2] and hyper-permeability of the gastrointestinal tract barrier [31][32][33]. Therefore, this toxins found in plasma and excreted in urine, cross the cellular and tissue barriers (gastrointestinal epithelium, lymphatic barrier and liver) are absorbed into the blood and mainly cleared by the kidney [33,34].
Hippurate and other Phe-derivative compounds are produced by the aerobic and anaerobic degradations of aromatic amino acids (e.g., Phe and Tyr) and dietary polyphenols by the gastrointestinal microbiota [34][35][36]. Moreover, high levels of gut-derived uremic toxins seem to affect both the cellular protein expression and the activity of the cyclooxygenase-2 (COX-2) enzyme, which plays a major role in the regulation of inflammation through the production of prostaglandins; so, when COX-2 activity is sped up, inflammation increases [37]. Some Phederivatives also produce cytotoxic effects by the inhibition of cell pores opening and the production of reactive oxygen species [38].  Among Phe-derivatives, hippurate has a strong association with diet and gastrointestinal microbiota, and its production requires of both microbial and mammalian metabolisms [39]. Gastrointestinal bacteria produce benzoic acid from dietary aromatic compounds, which is absorbed into the blood. Because of the toxicity of benzoic acid, it is conjugated with glycine in the mitochondrial matrix of the liver and renal cortex to form hippuric acid [31,39], which is later filtered in the kidneys and finally excreted in urine as hippurate [39,40]. The main elimination route for hippurate is the active renal tubular secretion and its disruption results in its accumulation in the blood [39]. Hippurate is a uremic toxin that participates in the correction of metabolic acidosis by stimulating ammoniagenesis, a dominant and adaptive mechanism of proton excretion. Moreover, it interferes with several metabolic processes, such as: inhibition of glucose utilization by the kidney and muscle, modulation of fatty acid metabolism and regulation of the acid-base balance by stimulating the kidneys' ammoniagenesis, among others, as reviewed by Dzúrik et al. [40].
Among these gut-derived metabolic compounds, changing levels of trimethylamine Noxide in plasma and milk were also observed in HS dairy cows [7,8]. Contradictory, these authors pointed out a lower level of this metabolite found in milk and plasma, while we observed an overexcretion of this compound through the urine.
It might also be noted that, in addition to the production of gut-derived uremic toxins from dietary aromatic amino acids by the gastrointestinal microbiota, Phe is known to be an essential amino acid for most animals, including ruminants [41]. It is also the precursor of Tyr, which is essential for the synthesis of thyroid hormones and the levodopa neurotransmitter. Previous studies have shown a strong decrease in plasma thyroid hormones (i.e., TSH, T4 and T3) in different ruminant models [42][43][44], which means that the basal heat production may, in fact, decrease when Phe and Tyr are scarce. Moreover, the rate of milk production is markedly affected by thyroid hormones, which modulate the nutrient partitioning towards milk production [45]. On the other hand, a decrease in the dopaminergic neurons activity was also observed in HS calves [46]. The drop of levodopa synthesis may be the result of the hypersecretion of its antagonist prolactin, as observed in response to HS in goats [47], ewes [48] and cows [49]. Prolactin is not only a hormone related to milk production, but also has a broad variety of biological functions related to thermoregulation and water balance. The increase in plasma prolactin is not reflected in an increase in milk yield, as seen in dairy ruminants under HS conditions [16,50]. Alamer [51] concluded that the mammary gland experiences a downregulation of prolactin-signaling pathways that could partially explain the depressed milk production of dairy cows during HS.
Increased concentration of 3-methyladenine in urine is associated with increased autophagy [52]. Autophagy controls the proteostasis in organisms (reviewed by Dokladny et al. [53]) and HS is an extracellular stressor that alters the folding capacity of a cell leading to the accumulation of misfolded or unfolded proteins [54]. Under stress conditions, eukaryotic cells increase the employ of autophagy to remove misfolded proteins, large protein aggregates, and whole damaged organelles inaccessible to smaller proteolytic systems [55]. Moreover, under negative energy balance, as commonly observed in HS animals, autophagy is an adaptive mechanism that provides biofuel from degraded macromolecules to maintain sufficient ATP production for adaptive macromolecular synthesis to survive stressful conditions [56]. One of the end products of protein catabolism is urea. An increased concentration of urea in blood, milk and urine is commonly observed in HS dairy cows [6,9,26,50] as a result of the strongly up-regulated pathway of nucleotides metabolism during HS [57]. Urea excretion peaks were compared between HS and TN does, but no differences were found in our study (P = 0.48) in agreement with that previously reported by Hamzaoui et al. in the uremia of HS dairy does [16]. Thus, because cows do not have very many active sweat glands, we speculate that a greater portion of urea may be lost in the sweat of goats when compared to cows.
On the other hand, the lower urinary excretion of metabolites related to vitamin metabolism (i.e. cholecalciferol, pyridoxal, β-alanine) may be a reflection of the commonly reported increased vitamin requirements of animals under thermal load [58].

Conclusions
Heat stress caused marked changes in thermophysiological traits and lactational performances of dairy goats, which were translated into their 1 H NMR metabolomic urinary profile. These changes were mainly related to the over-excretion of gut-derived toxic compounds generated by the gastrointestinal microbiota with expected decreases in the bioavailability of aromatic amino acids and impairment of the synthesis of thyroid hormones and neurotransmitters (i.e., levodopa, serotonin), which compromised the milk production of dairy goats. In practice, the use of hippurate and other phenylalanine derivatives are suggested as urinary biomarkers to identify heat-stressed animals.