RESEARCH ARTICLE Heat stress modifies the lactational performances and the urinary metabolomic profile related to gastrointestinal microbiota of dairy goats Alexandra Contreras-Jodar1☯, Nazri Hj. Nayan2☯, Soufiane Hamzaoui1¤, Gerardo Caja ID1*, Ahmed A. K. Salama1☯* 1 Group of Research in Ruminants (G2R), Department of Animal and Food Sciences, Universitat Autonoma a1111111111 de Barcelona, Bellaterra, Barcelona, Spain, 2 Department of Animal Science, Faculty of Agriculture, a1111111111 University Putra Malaysia, Serdang, Selangor, Malaysia a1111111111 a1111111111 ☯ These authors contributed equally to this work. a1111111111 ¤ Current address: Department of Agronomic Science, University of Bouira, Bouira, Algeria * gerardo.caja@uab.es (GC); ahmed.salama@sdstate.edu (AAKS) Abstract OPEN ACCESS Citation: Contreras-Jodar A, Nayan NH., Hamzaoui The aim of the study is to identify the candidate biomarkers of heat stress (HS) in the urine S, Caja G, Salama AAK (2019) Heat stress modifies of lactating dairy goats through the application of proton Nuclear Magnetic Resonance (1H the lactational performances and the urinary NMR)-based metabolomic analysis. Dairy does (n = 16) in mid-lactation were submitted to metabolomic profile related to gastrointestinal microbiota of dairy goats. PLoS ONE 14(2): thermal neutral (TN; indoors; 15 to 20˚C; 40 to 45% humidity) or HS (climatic chamber; 37˚C e0202457. https://doi.org/10.1371/journal. day, 30˚C night; 40% humidity) conditions according to a crossover design (2 periods of 21 pone.0202457 days). Thermophysiological traits and lactational performances were recorded and milk Editor: Juan J Loor, University of Illinois, UNITED composition analyzed during each period. Urine samples were collected at day 15 of each STATES period for 1H NMR spectroscopy analysis. Principal component analysis (PCA) and partial Received: August 1, 2018 least square—discriminant analysis (PLS-DA) assessment with cross validation were Accepted: January 16, 2019 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 Published: February 8, 2019 decreased feed intake (35%) and body weight (5%) of the lactating does. No differences Copyright: © 2019 Contreras-Jodar et al. This is an were detected in milk yield, but HS decreased the milk contents of fat (9%), protein (16%) open access article distributed under the terms of the Creative Commons Attribution License, which and lactose (5%). Metabolomics allowed separating TN and HS urinary clusters by PLS-DA. permits unrestricted use, distribution, and Most discriminating metabolites were hippurate and other phenylalanine (Phe) derivative reproduction in any medium, provided the original compounds, which increased in HS vs. TN does. The greater excretion of these gut-derived author and source are credited. toxic compounds indicated that HS induced a harmful gastrointestinal microbiota over- Data Availability Statement: All relevant data are growth, which should have sequestered aromatic amino acids for their metabolism and within the manuscript and its Supporting decreased the synthesis of neurotransmitters and thyroid hormones, with a negative impact Information files. on milk yield and composition. In conclusion, HS markedly changed the thermophysiological Funding: This study was part of 2 research traits and lactational performances of dairy goats, which were translated into their urinary projects funded by the Ministry of Economy and Finance of Spain (Project AGL2013-44061-R, to metabolomic profile through the presence of gut-derived toxic compounds. Hippurate and AAKS and GC, and Project RTA205-00035-C03-02, other Phe-derivative compounds are suggested as urinary biomarkers to detect heat- to GC) and it was also supported by a research stressed dairy animals in practice. scholarship to ACJ from the same funding source (Reference ES-2012-052602). The funders had no PLOS ONE | https://doi.org/10.1371/journal.pone.0202457 February 8, 2019 1 / 14 Urine metabolome in heat stressed goats role in study design, data collection and analysis, Introduction decision to publish, or preparation of the manuscript. Exposure to high ambient temperature induces several physiological responses in order to maintain body homeostasis. Animals suffer from heat stress (HS) when physiological mecha- Competing interests: The authors have declared that no competing interests exist. nisms 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 meta- bolic 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 produc- tion 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 (1H) 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 1H NMR spectra are derived from thousands of metabolite signals that usually overlap, adding complexity to data processing. Computer-based data reduction and multivariate statis- tical 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 1H 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 1H NMR-based metabolomic urinalysis of dairy goats. Material and methods 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 live- stock 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 experimen- tal 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 observa- tions 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 tem- perature-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 m3/h continuous air PLOS ONE | https://doi.org/10.1371/journal.pone.0202457 February 8, 2019 2 / 14 Urine metabolome in heat stressed goats turnover). Day-night length was set to 12–12 hours. Temperature-humidity index (THI), cal- culated according to NRC in 1971 [15], resulted in THITN = 59 to 65 and THIHS = 75 to 83. Experimental conditions were similar to those detailed by Hamzaoui et al. [16]. Does were milked once a day (0800) with a portable machine (Westfalia-Separator Ibérica; Granollers, Spain) set at 42 kPa, 90 pulses/min and 66% pulsation ratio and provided with vol- umetric recording jars (3 L ± 5%). The milking routine included cluster attachment, machine milking, machine stripping before cluster removal, and teat dipping in an iodine solution (P3-ioshield, Ecolab Hispano-Portuguesa; Barcelona, Spain). Feed was offered ad libitum at 0930 hours (130% feed intake of the previous day) and consisted of a total mixed ration (dry matter, 89.9%; net energy for lactation, 1.40 Mcal/kg; crude protein, 17.5%; organic matter, 87.3%; neutral detergent fiber, 34.4%; acid detergent fiber, 21.8%; on dry-matter basis). Ration ingredients were (as fed): alfalfa hay, 64.2%; ground barley, 9.6%; beet pulp, 9.6%; ground corn, 8%; soybean meal, 3.3%; sunflower meal, 3.2%; molasses, 1%; salt, 0.6%; sodium bicar- bonate, 0.3%; mineral and vitamin complex, 0.2% (Vitafac premix, DSM Nutritional Products; Madrid, Spain). Water was permanently available and offered at room temperature in water bowls connected to individual tanks of 20 L. A sawdust drip tray under each water bowl was used to collect spilled water. 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 micturi- tion on the morning of day 15 of each period (n = 32) and stored at −20˚C until 1H NMR analysis. Preparation of samples for 1H NMR spectroscopy was done according to Beckonert et al. [12]. Briefly, a phosphate buffer solution (pH 7.4) was prepared with sodium phosphate diba- sic (Na2HPO4; 99.95% trace metals basis, anhydrous, Sigma-Aldrich Merck; Darmstadt Ger- many), sodium phosphate monobasic (NaH2PO4; 99.95% trace metals basis, anhydrous, Sigma-Aldrich Merck) and sodium azide (NaN3; Sigma-Aldrich Merck). Deuterium oxide (D2O; 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, Ger- many) 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 cen- trifuged 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 PLOS ONE | https://doi.org/10.1371/journal.pone.0202457 February 8, 2019 3 / 14 Urine metabolome in heat stressed goats immediately put on ice and sent to the NMR Service of the UAB for 1H High Resolution NMR Spectroscopy. NMR spectroscopy. 1H NMR spectra were acquired on a Bruker Avance-III spectrometer (Bruker BioSpin; Rheinstetten, Germany) operating at a 1H 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 con- trolled by TopSpin2.1 software (Bruker, Germany). One-dimensional 1H 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. 1H NMR spectra were phased, baseline corrected, and corrected for chemical shift registration rel- ative 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. Dif- ferences between least squares means were determined with the PDIFF option of SAS. Signifi- cance 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. Neverthe- less, metabolic fingerprinting datasets acquired from 1H 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, align- ment and normalization of spectra were performed to ensure that all observations were directly comparable. In this sense, urine spectra were normalized to creatinine methyl reso- nance intensity at δ = 3.05 ppm and then log2 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 his- tidine, 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 1H NMR pre-processing data were completed, data were subjected to multivariate sta- tistical 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-valida- tion 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 PLOS ONE | https://doi.org/10.1371/journal.pone.0202457 February 8, 2019 4 / 14 Urine metabolome in heat stressed goats (TN or HS). Model strength was assessed using both R2 and Q2 statistical parameters. While R2 values reported the total amount of variance explained by the model, the Q2 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 Q2 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 (log2 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 influ- enced by variable correlations and the intra- and inter-class variance of metabolites may have no significant differences in the univariate analysis [23]. All 1H NMR data pre-processing, sta- tistical 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 candi- date chemical shifts and corresponding metabolites were assigned using the Human Metabo- lome Database [24] and queried in KEGG (Kyoto Encyclopedia of Genes and Genomes) database to know in which metabolic pathways they were involved. Results and discussion Effects of heat stress on thermophysiological and lactational performances of the goats The effects of the experimental HS conditions on thermophysiological and lactational perfor- mances 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 boost- ing 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 PLOS ONE | https://doi.org/10.1371/journal.pone.0202457 February 8, 2019 5 / 14 Urine metabolome in heat stressed goats Table 1. Thermophysiological and lactational performances of dairy goats under thermal neutral (TN) and heat stress (HS) conditions. Values are least square means and standard error of the means (SEM). Item Treatment SEM P value TN HS Rectal temperature, ˚C 0800 hours 38.5 39.1 0.08 <0.001 1200 hours 38.7 39.7 0.07 <0.001 1700 hours 38.7 39.9 0.09 <0.001 Respiratory rate, breaths/min 0800 hours 27 69 4 <0.001 1200 hours 39 131 6 <0.001 1700 hours 37 130 6 <0.001 Performances Dry matter intake, kg/d 2.26 1.47 0.09 <0.001 Water intake, L/d 6.1 10.6 1.0 <0.001 Final body weight, kg 48.6 39.8 1.8 <0.001 Body weight variation, kg 3.5 −2.1 1.0 <0.001 Milk yield, L/d 1.88 1.79 0.11 0.413 FCM1 yield, L/d 2.17 1.86 0.13 0.017 Milk composition, % Fat 3.98 3.64 0.13 0.009 Protein 3.40 2.85 0.10 <0.001 Lactose 4.51 4.30 0.07 0.003 1 Fat-corrected milk at 3.5%; FCM = L × [0.432 + 0.162 × (fat, %)] being L liters of milk. https://doi.org/10.1371/journal.pone.0202457.t001 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 usu- ally 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 1H NMR-metabolomics assessment were collected (day 15). NMR urinary spectroscopy of the goats A comparison of 1H 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 PLOS ONE | https://doi.org/10.1371/journal.pone.0202457 February 8, 2019 6 / 14 Urine metabolome in heat stressed goats Fig 1. One-dimensional 1H NMR spectra at 600 MHz of urine from representative thermoneutral (TN) and heat-stressed (HS) dairy does. Dominant metabolites were: 1, creatinine; 2, creatine; 3, trimethyl-N-oxide; 4, urea; 5, branched-chain amino acids and organic acids; 6, glycine; 7, allantoin; 8, alanine; 9, N-acetyl glycoprotein; 10, glutamate; 11, succinic acid; 12, citric acid; 13, aromatic signals; 14, hippuric acid. https://doi.org/10.1371/journal.pone.0202457.g001 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 multi- variate 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 1H NMR spectra. Based on the principle of minimum differentiation, no samples were identified as outliers according to Hotelling’s T2 (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 Rx2 = 0.54, Ry2 = 0.17, and Q2 = 0.47. The R2 and Q2 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-phenylacetyl- glycine, phenylglyoxylate and trimethylamine N-oxide), the thresholds applied for Volcano PLOS ONE | https://doi.org/10.1371/journal.pone.0202457 February 8, 2019 7 / 14 Urine metabolome in heat stressed goats Fig 2. Volcano plot based on fold change (log2) and P value (−log10) of all spectral bins of 1H NMR urinary spectroscopy of heat-stressed (HS) vs. thermoneutral (TN) lactating dairy does. Red circles indicate the spectral bins that showed significant changes and absolute fold changes greater than 1.5. https://doi.org/10.1371/journal.pone.0202457.g002 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 vita- min 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 envi- ronment 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 vasoconstric- tion 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–33]. Therefore, this toxins found in plasma and excreted in urine, cross the cellular and tissue barriers (gastrointestinal epithe- lium, lymphatic barrier and liver) are absorbed into the blood and mainly cleared by the kid- ney [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 PLOS ONE | https://doi.org/10.1371/journal.pone.0202457 February 8, 2019 8 / 14 Urine metabolome in heat stressed goats Fig 3. PLS-DA scores plot of the first two principal components of 1H NMR urinary spectra of thermoneutral (TN) and heat-stressed (HS) lactating dairy does. https://doi.org/10.1371/journal.pone.0202457.g003 gastrointestinal microbiota [34–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 Phe- derivatives also produce cytotoxic effects by the inhibition of cell pores opening and the pro- duction of reactive oxygen species [38]. Table 2. Selected metabolites contributing to the classification of the urine metabolome of thermoneutral and heat-stressed lactating dairy does. Metabolic pathway Chemical shift (δ, ppm) Metabolite Fold change� P value Phenylalanine (Phe) 7.83, 7.63, 7.54 Hippurate 2.74 <0.001 Tyrosine (Tyr) 7.20 OH-Phenylacetylglycine 2.16 <0.001 Microbial metabolism 3.27 Trimethylamine N-oxide 2.11 <0.001 Nucleotide 3.97 3-Methyladenine 1.82 <0.001 Arginine (Arg) and Proline (Pro) 3.93 Phosphocreatine 1.69 <0.001 Microbial metabolism 7.62 Phenylglyoxylate 1.59 <0.001 Microbial metabolism 7.27 OH-Phenylacetate 1.51 <0.001 Microbial and purine metabolism 5.30 Allantoic acid −1.52 <0.001 Vitamin D 0.56, 0.53 Cholecalciferol −1.57 <0.001 Glucagon signaling 5.40 Glycogen −1.66 <0.001 Galactose 3.69 Galactitol −1.81 <0.001 Neurodegeneration 2.86 6OH-Dopamine −1.87 <0.001 Vitamin B6 5.29, 5.28, 5.27 Pyridoxal −2.03 <0.001 Pantotenate and CoA biosynthesis and pyrimidine metabolism 3.17 β-Alanine (Ala) −2.09 <0.001 Histamine 3.28 Histamine −2.64 <0.001 � Metabolites with positive fold change values mean that they are excreted in greater concentrations under heat-stressed conditions. Those metabolites with negative fold change values are excreted in higher concentrations under thermoneutral conditions. https://doi.org/10.1371/journal.pone.0202457.t002 PLOS ONE | https://doi.org/10.1371/journal.pone.0202457 February 8, 2019 9 / 14 Urine metabolome in heat stressed goats 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 elimi- nation 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 regu- lation 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 N- oxide 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 essen- tial 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–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 hyper- secretion 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 down- regulation of prolactin-signaling pathways that could partially explain the depressed milk pro- duction of dairy cows during HS. Increased concentration of 3-methyladenine in urine is associated with increased autop- hagy [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 aggre- gates, and whole damaged organelles inaccessible to smaller proteolytic systems [55]. More- over, 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 condi- tions [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, PLOS ONE | https://doi.org/10.1371/journal.pone.0202457 February 8, 2019 10 / 14 Urine metabolome in heat stressed goats 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 metabo- lism (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 1H 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. Supporting information S1 Dataset. 1H NMR data matrix of normalized and binned spectral data. HS, heat-stressed lactating dairy does; TN, thermal neutral lactating dairy does. (TXT) Acknowledgments The authors are grateful to the technical team of SGCE (Servei de Granges i Camps Experi- mentals) of the UAB for the care of the animals and to Charles (Chuck) Simmons, native English-speaking University Instructor (Cerdanyola, Barcelona, Spain) for the English lan- guage and style revision of the manuscript. Author Contributions Conceptualization: Gerardo Caja, Ahmed A. K. Salama. Data curation: Alexandra Contreras-Jodar, Nazri Hj. Nayan. Formal analysis: Alexandra Contreras-Jodar. Funding acquisition: Ahmed A. K. Salama. Investigation: Alexandra Contreras-Jodar, Nazri Hj. Nayan, Soufiane Hamzaoui, Gerardo Caja, Ahmed A. K. Salama. Methodology: Nazri Hj. Nayan, Soufiane Hamzaoui, Ahmed A. K. Salama. Project administration: Gerardo Caja, Ahmed A. K. Salama. Resources: Ahmed A. K. Salama. Software: Alexandra Contreras-Jodar. Supervision: Gerardo Caja, Ahmed A. K. Salama. Validation: Nazri Hj. Nayan. Visualization: Alexandra Contreras-Jodar, Nazri Hj. Nayan, Gerardo Caja. Writing – original draft: Alexandra Contreras-Jodar, Nazri Hj. Nayan. PLOS ONE | https://doi.org/10.1371/journal.pone.0202457 February 8, 2019 11 / 14 Urine metabolome in heat stressed goats Writing – review & editing: Gerardo Caja, Ahmed A. K. Salama. References 1. Bernabucci U, Lacetera N, Baumgard LH, Rhoads RP, Ronchi B, Nardone A. Metabolic and hormonal acclimation to heat stress in domesticated ruminants. Anim. 2010; 4:1167–83. https://doi.org/10.1017/ S175173111000090X PMID: 22444615. 2. Baumgard LH, Rhoads RP. Effects of heat stress on postabsorptive metabolism and energetics. Annu Rev Anim Biosci. 2013; 1:311–37. https://doi.org/10.1146/annurev-animal-031412-103644 PMID: 25387022. 3. Hansen PJ. Effects of heat stress on mammalian reproduction. Philos Trans R Soc Lond B Biol Sci. 2009; 364:3341–50. https://doi.org/10.1098/rstb.2009.0131 PMID: 19833646. 4. Salama AAK, Caja G, Hamzaoui S, Badaoui B, Castro-Costa A, Façanha DE et al. Different levels of response to heat stress in dairy goats. Small Rumin. Res. 2014; 121:73–9. 5. Rhoads ML, Rhoads RP, VanBaale MJ, Collier RJ, Sanders SR, Weber WJ et al. Effects of heat stress and plane of nutrition on lactating Holstein cows: I. Production, metabolism, and aspects of circulating somatotropin. J Dairy Sci. 2009; 92:1986–97. https://doi.org/10.3168/jds.2008-1641 PMID: 19389956. 6. Wheelock JB, Rhoads RP, VanBaale MJ, Sanders SR, Baumgard LH. Effects of heat stress on ener- getic metabolism in lactating Holstein cows. J Dairy Sci. 2010; 93:644–55. https://doi.org/10.3168/jds. 2009-2295 PMID: 20105536. 7. Tian H, Wang W, Zheng N, Cheng J, Li S, Zhang Y et al. Identification of diagnostic biomarkers and metabolic pathway shifts of heat-stressed lactating dairy cows. J. Proteomics. 2015; 125:17–28. https:// doi.org/10.1016/j.jprot.2015.04.014 PMID: 25913299. 8. Tian H, Zheng N, Wang W, Cheng J, Li S, Zhang Y et al. Integrated metabolomics study of the milk of heat-stressed lactating dairy cows. Sci. Rep. 2016; 6:24208. https://doi.org/10.1038/srep24208 PMID: 27048914. 9. Fan C, Su D, Tian H, Li X, Li Y, Ran L et al. Liver metabolic perturbations of heat-stressed lactating dairy cows. Asian-Australas. J Anim. Sci. 2018; 31:1244–1251. https://doi.org/10.5713/ajas.17.0576 PMID: 29514433. 10. Dou S, Villa-Vialaneix N, Liaubet L, Billon Y, Giorgi M., Gilbert JL et al. 1HNMR-Based metabolomic pro- filing method to develop plasma biomarkers for sensitivity to chronic heat stress in growing pigs. PLoS ONE 2017; 12:e0188469. https://doi.org/10.1371/journal.pone.0188469 PMID: 29176781. 11. Ippolito DL, Lewis JA, Yu C, Leon LR, Stallings JD. Alteration in circulating metabolites during and after heat stress in the conscious rat: potential biomarkers of exposure and organ-specific injury. BMC Phy- siol. 2014; 14:14. https://doi.org/10.1186/s12899-014-0014-0 PMID: 25623799. 12. Beckonert O, Keun HC, Ebbels TMD, Bundy J, Holmes E, Lindon JC et al. Metabolic profiling, metabo- lomic and metabonomic procedures for NMR spectroscopy of urine, plasma, serum and tissue extracts. Nature Protocols. 2007; 2:2692–703. https://doi.org/10.1038/nprot.2007.376 PMID: 18007604. 13. Westerhuis JA, Hoefsloot HCJ, Smit S, Vis DJ, Smilde AK, van Velzen EJJ et al. Assessment of PLSDA cross-validation. Metabolomics. 2008; 4:81–9. 14. Yamamoto H, Yamaji H, Abe Y, Harada K, Waluyo D, Fukusaki E et al. Dimensionality reduction for metabolome data using PCA, PLS, OPLS, and RFDA with differential penalties to latent variables. Che- mom Intell Lab Syst. 2009; 98:136–42. 15. NRC (National Research Council). A Guide to Environmental Research on Animals. Washington, DC (USA): National Academies Press; 1971. 16. Hamzaoui S, Salama AAK, Albanell E, Such X, Caja G. Physiological responses and lactational perfor- mances of late lactating dairy goats under heat stress conditions. J Dairy Sci. 2013; 96:6355–65. https://doi.org/10.3168/jds.2013-6665 PMID: 23958010. 17. AOAC (Association of Official Analytical Chemists). Official Methods of Analysis of the Association of Official’s Analytical Chemists, 17th Edition. Arlington, Virginia (USA): AOAC; 2003. 18. Worley B, Powers R. Multivariate analysis in metabolomics. Curr Metabolomics. 2013; 1:92–107. https://doi.org/10.2174/2213235X11301010092 PMID: 26078916. 19. R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. 2015. URL https://www.R-project.org/. 20. Pechlivanis A, Kostidis S, Saraslanidis P, Petridou A, Tsalis G, Mougios V et al. H-1 NMR-Based Meta- bonomic Investigation of the Effect of Two Different Exercise Sessions on the Metabolic Fingerprint of Human Urine. J Proteome Res. 2010; 9:6405–16. https://doi.org/10.1021/pr100684t PMID: 20932058. PLOS ONE | https://doi.org/10.1371/journal.pone.0202457 February 8, 2019 12 / 14 Urine metabolome in heat stressed goats 21. Mevik B, Wehrens R, Hovde Liland K. Pls: Partial Least Squares and Principal Component Regression. R package version 2.5–0. 2015. http://CRAN.R-project.org/package=pls. 22. Benjamini Y, Hochberg Y. Controlling the False Discovery Rate: A practical and powerful approach to multiple testing. J R Statist Soc B. 1995; 57:289–300. 23. Ji P, Wei Y, Hua Y, Zhang X, Yao W, Ma Q et al. A novel approach using metabolomics coupled with hematological and biochemical parameters to explain the enriching-blood effect and mechanism of unprocessed Angelica sinensis and its 4 kinds of processed products. J Ethnopharmacol. 2017; 211:101–16. https://doi.org/10.1016/j.jep.2017.09.028 PMID: 28958590. 24. Wishart DS, Knox C, Guo AC. 2009. HMDB: a knowledgebase for the human metabolome. Nucleic Acids Res. 2009; 37 (Database Issue):D603–610. https://doi.org/10.1093/nar/gkn810 PMID: 18953024. 25. Das R, Sailo L, Verma N, Bharti P, Saikia J, Imtiwati Kumar R. Impact of heat stress on health and per- formance of dairy animals: A review. Vet World 2016; 9:260–268. https://doi.org/10.14202/vetworld. 2016.260-268 PMID: 27057109. 26. Shwartz G, Rhoads ML, VanBaale MJ, Rhoads RP, Baumgard LH. Effects of a supplemental yeast cul- ture on heat-stressed lactating Holstein cows. J. Dairy Sci. 2009; 92:935–42. https://doi.org/10.3168/ jds.2008-1496 PMID: 19233786. 27. Bouatra S, Aziat F, Mandal R, Guo AC, Wilson MR, Knox C et al. The human urine metabolome. PLoS One 2013; 8:e73076. https://doi.org/10.1371/journal.pone.0073076 PMID: 24023812. 28. Kregel KC, Wall PT, Gisolfi CV. Peripheral vascular responses to hyperthermia in the rat. J Appl Phy- siol. 1988; 64:2582–8. https://doi.org/10.1152/jappl.1988.64.6.2582 PMID: 3403442. 29. Lambert GP. Stress-induced gastrointestinal barrier dysfunction and its inflammatory effects. J Anim Sci. 2009; 87(14 Suppl):E101–8. https://doi.org/10.2527/jas.2008-1339 PMID: 18791134. 30. Castro-Costa A, Salama AAK, Moll X, Aguiló J, Caja G. Using wireless rumen sensors for evaluating the effects of diet and ambient temperature in nonlactating dairy goats. J Dairy Sci. 2015; 98:4646–58. https://doi.org/10.3168/jds.2014-8819 PMID: 25958285. 31. Nicholson JK, Holmes E, Wilson ID. Gut microorganisms, mammalian metabolism and personalized health care. Nat Rev Microbiol. 2005; 3:431–8. https://doi.org/10.1038/nrmicro1152 PMID: 15821725. 32. Pearce SC, Mani V, Weber TE, Rhoads RP, Patience JF, Baumgard LH, Gabler NK. Heat stress and reduced plane of nutrition decreases intestinal integrity and function in pigs. J Anim Sci. 2013; 91:5183– 93. https://doi.org/10.2527/jas.2013-6759 PMID: 23989867. 33. Ramezani A, Raj DS. The gut microbiome, kidney disease, and targeted interventions. J Am Soc Nephr. 2014; 25: 657–70. https://doi.org/10.1681/ASN.2013080905 PMID: 24231662. 34. Martin AK. The origin of urinary aromatic compounds excreted by ruminants. 1. The metabolism of qui- nic, cyclohexanecarboxylic and non-phenolic aromatic acids to benzoic acid. British J Nutr. 1982; 47:139–54. PMID: 7059566. 35. Macfarlene GT, Macfarlene S. Human colonic microbiota: ecology, physiology and metabolic potential of intestinal bacteria. Scand J Gastroenterol. 1997; 32 (Suppl 222):3–9. https://doi.org/10.1080/ 00365521.1997.11720708 PMID: 9145437. 36. Pallister T, Jackson MA, Martin TC, Zierer J, Jennings A, Mohney RP, MacGregor A, Steves CJ, Cas- sidy A, Spector TD, Menni C. Hippurate as a metabolomic marker of gut microbiome diversity: Modula- tion by diet and relationship to metabolic syndrome. Sci Rep. 2017; 7:13670. https://doi.org/10.1038/ s41598-017-13722-4 PMID: 29057986. 37. Karlsson PC, Huss U, Jenner A, Halliwell B, Bohlin L, Rafter JJ. Human fecal water inhibits COX-2 in colonic HT-29 cells: role of phenolic compounds. J Nutr. 2005; 135:2343–9. https://doi.org/10.1093/jn/ 135.10.2343 PMID: 16177193. 38. Galati G, Lin A, Sultan AM, O’Brian PJ. Cellular and in vivo hepatotoxicity caused by green tea phenolic acids and catechins. Free Radic Biol Med. 2006; 40:570–80. https://doi.org/10.1016/j.freeradbiomed. 2005.09.014 PMID: 16458187. 39. Lees HJ, Swann JR, Wilson ID, Nicholson JK, Holmes E. Hippurate: the natural history of a mamma- lian-microbial cometabolite. J Proteome Res. 2013; 12:1527–46. https://doi.org/10.1021/pr300900b PMID: 23342949. 40. Dzúrik R, Spustová V, Krivosı́ková Z, Gazdı́ková K. Hippurate participates in the correction of metabolic acidosis. Kidney Int. 2001; 78:S278–81. https://doi.org/10.1046/j.1523-1755.2001.59780278.x PMID: 11169026. 41. Patton RA, Hristov AN, Parys C, Lapierre H. Relationships between circulating plasma concentrations and duodenal flows of essential amino acids in lactating dairy cows. J Dairy Sci. 2015; 98: 4707–34. https://doi.org/10.3168/jds.2014-9000 PMID: 25912872. PLOS ONE | https://doi.org/10.1371/journal.pone.0202457 February 8, 2019 13 / 14 Urine metabolome in heat stressed goats 42. Collier RJ, Beede DK, Thatcher WW, Israel LA, Wilcox CJ. Influences of environment and its modifica- tion on dairy animal health and production. J Dairy Sci. 1982; 65:2213–27. https://doi.org/10.3168/jds. S0022-0302(82)82484-3 PMID: 6759540. 43. Kahl S, Elsasser TH, Rhoads RP, Collier RJ, Baumgard LH. Environmental heat stress modulates thy- roid status and its response to repeated endotoxin challenge in steers. Domest Anim Endocrinol. 2015; 52:43–50. https://doi.org/10.1016/j.domaniend.2015.02.001 PMID: 25804834. 44. Hooper HB, Dos Santos Silva P, de Oliveira SA, Merighe GKF, Negrão JA. Acute heat stress induces changes in physiological and cellular responses in Saanen goats. Int J Biometeorol. 2018; 62:2257– 2265. https://doi.org/10.1007/s00484-018-1630-3 PMID: 30368674. 45. Hurley WL. Adrenal and thyroid hormones; 2010 [cited 25 July 2018]. In: Lactation biology website. Uni- versity of Illinois, Urbana-Champaign, IL. http://ansci.illinois.edu/static/ansc438/Lactation/ adrenalthyroid.html 46. Tucker HA, Chapin LT, Lookingland KJ, Moore KE, Dahl GE, Evers JM. Temperature effects on serum prolactin concentrations and activity of dopaminergic neurons in the infundibulum/pituitary stalk of calves. Proc Soc Exp Biol Med. 1991; 197:74–6. PMID: 1902294. 47. Sergent D, Berbigier P, Kann G, Ravault JP. Effect of prolactin inhibition on thermo-physiological parameters and on plasma prolactin and cortisol concentrations in male Creole goats. Reprod Nutr Dev. 1985; 25:629–40. 48. Colthorpe KL, Anderson ST, Martin GB, Curlewis JD. Hypothalamic dopamine D1 receptors are involved in the stimulation of prolactin secretion by high environmental temperature in female sheep. J Neuroendocrinol. 1998; 10:503–9. PMID: 9700677. 49. Ronchi B, Stradaioli G, Supplizi AV, Bernabucci U, Lacetera N. Influence of heat stress or feed restric- tion on plasma progesterone, oestradiol-17b, LH, FSH, prolactin and cortisol in Holstein heifers. Livest. Prod Sci. 2001; 68:231–41. 50. Gao ST, Guo J, Quan SY, Nan XM, Fernández MVS, Baumgard LH, Bu DP. The effects of heat stress on protein metabolism in lactating Holstein cows. J Dairy Sci. 2017; 100:5040–9. https://doi.org/10. 3168/jds.2016-11913 PMID: 28390717. 51. Alamer M. The role of prolactin in thermoregulation and water balance during heat stress in domestic ruminants. Asian J Anim Vet Adv. 2011; 6:1153–69. 52. Wu YT, Tan HL, Shui G, Bauvy C, Huang Q, Wenk MR, Ong CN, Codogno P, Shen HM. Dual role of 3- methyladenine in modulation of autophagy via different temporal patterns of inhibition on class I and III phosphoinositide 3-kinase. J Biol Chem. 2010; 285:10850–61. https://doi.org/10.1074/jbc.M109. 080796 PMID: 20123989. 53. Dokladny K, Myers OB, Moseley PL. Heat shock response and autophagy-cooperation and control. Autophagy. 2015; 11:200–13. https://doi.org/10.1080/15548627.2015.1009776 PMID: 25714619. 54. Tytell M and Hooper PL. Heat shock proteins: new keys to the development of cytoprotective therapies. Expert Opin Ther Targets. 2001; 5(2):267–87. https://doi.org/10.1517/14728222.5.2.267 PMID: 15992180. 55. Klionsky DJ, Codogno P. The mechanism and physiological function of macroautophagy. J Innate immune. 2013; 5:427–33. https://doi.org/10.1159/000351979 PMID: 23774579. 56. Choi AM, Ryter SW, Levine B. Autophagy in human health and disease. The N Engl J Med 2013; 368:1845–6. https://doi.org/10.1056/NEJMc1303158 PMID: 23656658. 57. Contreras-Jodar A, Salama AA, Hamzaoui S, Vailati-Riboni M, Caja G, Loor JJ. Effects of chronic heat stress on lactational performance and the transcriptomic profile of blood cells in lactating dairy goats. J Dairy Res. 2018; 85:423–430. https://doi.org/10.1017/S0022029918000705 PMID: 30236165. 58. Clarkson PM. Nutritional Needs in Hot Environments: Applications for Military Personnel in Field Opera- tions. Institute of Medicine (US) Committee on Military Nutrition Research; Marriott BM, editor. Wash- ington, DC (USA): National Academies Press; 1993. PLOS ONE | https://doi.org/10.1371/journal.pone.0202457 February 8, 2019 14 / 14