Weixing Zhang

Yan Liu
email: yanliu@gdut.edu.cn

Christelle Hureau

Anne Robert
email: anne.robert@lcc-toulouse.fr

Bernard Meunier

N4-Tetradentate chelators efficiently regulate copper homeostasis and prevent ROS production induced by copper-amyloid-β1-16, even in the presence of an excess of zinc

Keywords: Alzheimer's disease, Copper chelator, Copper homeostasis, Oxidative stress, Zinc

N4-Tetradentate chelators efficiently regulate copper homeostasis and prevent ROS production induced by copper-amyloid-β1-16, even in the presence of an excess of zinc.

brain of patients with Alzheimer's disease (AD) is now considered as a valuable challenge in the chemotherapy of this neurodegenerative condition. 1 However, zinc element is abundant and strictly regulated in the brain to avoid depletion or excess. 2 In addition to zinc-protein complexes, large amounts of chelatable zinc are located in the cerebral cortex and the limbic region, notably in the hippocampal formation, and concentrated within synaptic vesicles of glutamatergic neurons where it acts as a modulator of synaptic transmission. 2c Zinc concentration has been estimated to reach 300 µM in the synaptic cleft system during strong stimulation. 2,[START_REF] Sensi | [END_REF] The mean copper concentration in human frontal lobe and cerebellum has been reported to be in the range of 60-110 µM. 4 The zinc/copper ratio can be as high as 10 to 100, in glutamatergic neurons. 5a Therefore, any attempt to regulate copper homeostasis and to inhibit copper-induced oxidative stress in the brain, must take into account this zinc-rich environment. 5b Therefore, the ligands designed to regulate the homeostasis of redox metal ions in AD brains must be copper specific to avoid the perturbation of zinc. As an example of the importance of the possible copper/zinc competition, clioquinol has been discarded due to its neurotoxicity attributed to zinc chelation. 6 As a matter of fact, the metal selectivity of 8hydroxyquinoline scaffold is very low. [START_REF] Phillips | [END_REF] To regulate copper homeostasis, chelators should be able to transfer copper from its pathological sink, mainly from the copper-amyloid complexes (Cu-Aβ), to glutathione which is a physiological provider of copper to copper enzymes, 8 even in the presence of zinc. In addition, a large excess of zinc should not prevent the ability of these copper chelators to inhibit the production of reactive oxygen species (ROS) generated by copper loaded amyloids. Moreover, such ligands should have logP values expected for membrane crossing, especially the blood brain barrier (BBB).

However, among the copper chelators which are currently under investigation, the selectivity for copper with respect to zinc is usually not (or poorly) documented, 1a except for bis-8-aminoquinoline derivatives that are specific copper chelators. 9 In this view, some of us reported the capacity of a water-soluble negatively charged chelator named L2 (Figure 1) to remove copper from copper-loaded Aβ in the presence of zinc in order to illustrate the importance of the Cu selectivity on ROS formation. 5 This N2O2 chelator, based on a Schiffbase scaffold, was able to inhibit the reductive activation of dioxygen in the presence of 1 molar equivalent (mol equiv) of Zn(II) but failed to stop ROS production in the presence of 5 mol equiv of Zn(II) during the ascorbate consumption. 5a In addition, as pointed out in the article, the ligand L2 is not able to cross the BBB (logP = -3.7).

To enhance the efficiency and drugability of copper chelators, we designed a new series of neutral N4-tetradendate copper ligands named TDMQ, based on an 8-aminoquinoline motif. 10 One of these ligands, TDMQ20 (Figure 1), is a specific ligand for copper compared to zinc. In the present work, we investigated the ability of TDMQ20 to inhibit the Cu-Aβ activation of dioxygen in the presence of a large excess of zinc(II), up to 100 mol equiv. The extraction of copper from Cu-Aβ by TDMQ20, as well as the release of copper from Cu-TDMQ20 in the presence of glutathione, were also evaluated in the presence of zinc, as these two discrete steps are required to put back copper into normal physiological circulation. The obtained data were compared with that of ligand 1, a bis-8-aminoquinoline chelator previously reported to have a very high selectivity for Cu(II) with respect to Zn(II) {log [Kapp Cu-L / Kapp Zn-L] > 12}, and to efficiently transfer copper from Cu-Aβ to glutathione in the absence of zinc. 9 The 8-hydroxyquinoline derivative PBT2 1c,11 was used as a comparator. A detailed experimental section is provided as Supporting Information. The ability of TDMQ20 to extract copper from Cu-Aβ1-16 was evaluated as follows. The Cu(II) complex of Aβ1-16 was prepared by incubation of Aβ1-16 and Cu 2+ at 20 µM (molar ratio = 1/1) in Hepes buffer pH 7.4. One or three mol equiv of Zn 2+ were added into the Cu II -Aβ1-16 solution. One mol equiv of TDMQ20 was then added, and the reaction was monitored by UV-visible spectroscopy. The spectrum of the resulting mixture containing Aβ1-16, Cu(II), Zn(II) and TDMQ20 was recorded (Figure 2, red trace, full line for 1 equiv of Zn, dashed line for 3 equiv of Zn). In both cases, the spectrum exhibited absorbance at 243 nm assigned to Cu-TDMQ20 (bright blue trace), but not at 259 and 342 nm, specific wavelengths of the TDMQ20 free ligand. In fact, the spectra of the mixture Cu II -Aβ1-16 / TDMQ20 in the presence of Zn 2+ (red traces) were very similar to the arithmetic addition of the spectra of free Aβ1-16 and Cu-TDMQ20 (black dashed trace). This result is the same as that obtained in a reaction carried out in the absence of zinc (Figure S1, Supporting Information). It therefore clearly indicates that 1 mol equiv of TDMQ20 was sufficient to fully demetalate the copper loaded amyloid, providing Cu-TDMQ20, in the presence of either 1 or 3 mol equiv of zinc with respect to copper. These data indicate that the presence of 3 mol equiv of zinc does not interfere with the demetalation of Cu II -Aβ1-16 by one equivalent of TDMQ20.

Figure 2. Extraction of copper from Cu II -Aβ1-16 using 1 mol equiv of TDMQ20, in the presence of 1 or 3 mol equiv of Zn 2+ (full red trace, and dashed red trace, respectively) with respect to copper. For detailed experimental conditions, see Supporting Information.

The release of copper from Cu II -TDMQ20 in the presence of GSH was then monitored.

As previously reported, when the Cu(II) complex of the tetradendate ligand 1 was treated by glutathione (GSH), GSH acts both as a reducing agent of Cu(II) and as a competitive ligand for the generated Cu(I). 9c The resulting Cu(I)-glutathione complex is a physiological supplier of copper to the apo form of copper proteins. 8 However, the role of zinc in this process has not been reported yet. Consequently, the influence of zinc on copper circulation in the presence of chelators remained to be established. We decided to monitor, by UV-visible spectrometry, the demetalation of Cu II -TDMQ20 in the presence of both glutathione and Zn 2+ , in Hepes buffer, pH 7.4. Firstly, CuCl2 was incubated with TDMQ20 to generate Cu II -TDMQ20, and 1 mol equiv of ZnCl2 was then added. The spectrum of the reaction mixture exhibited a maximal absorbance at 243 nm, typical of that expected for Cu II -TDMQ20 (Figure 3, red trace). Upon addition of glutathione (by aliquots of 20 mol equiv with respect to Cu-TDMQ), the demetalation of Cu-TDMQ was monitored by formation of the free ligand TDMQ20 detected at 259 and 342 nm (well-defined isosbestic points are detected at 279 and 332 nm; no isosbestic point can be detected below 250 nm, due to the own absorbance of glutathione and Cu(I)-glutathione at these wavelengths) (Figure 3). The plotting of the absorbance at 259 nm as a function of the number of glutathione equivalents (Figure 3, Insert) indicate that 80 equiv of glutathione were required to achieve full demetalation of Cu-TDMQ.

Interestingly, a similar result was observed in the absence of ZnCl2 (Figure S2, Supporting Information), indicating that the presence of zinc did not influence the demetalation of Cu-TDMQ20 by glutathione in these conditions. So, in vitro, TDMQ20 is able to abstract copper from Cu-Aβ1-16 and then to transfer it to glutathione, indicating that this N4-tetradendate ligand has the capacity for regulation of copper homeostasis in the presence of zinc. The in vitro ROS production mediated by Cu-Aβ1-16 in the presence or in the absence of any additive (namely zinc or a putative chelator) can be evaluated based on the kinetics of ascorbate oxidation, easily monitored by UV-visible at 265 nm. 9c, 12,13 We investigated the influence of TDMQ20 (and ligands 1 and PBT2 as comparators, structures in Figure 1) on the oxidation of ascorbate induced by Cu-Aβ1-16 in the presence of Zn 2+ . Firstly, a Zn/Cu molar ratio equal to 3 was investigated since senile plaques contain roughly 3 times more Zn than Cu (25 and 69 µg/g for Cu and Zn, respectively). 14 The results are showed in We also performed similar experiments with PBT2, a bi-or tridentate ligand based on 8hydroxyquinoline scaffold, that has been developed as a potential AD drug. 1c,15 Since the Cu-PBT2 complex has potentially a metal/ligand stoichiometry = 1/2, 16 we therefore used 2.2 mol equiv of PBT2 with respect to Cu-Aβ1-16, to avoid the presence of free Cu 2+ ions in the solution. With 3 mol equiv of Zn with respect to copper, 89% of ascorbate oxidation was observed in 30 min (Figure 4, trace f), close to the 80% value reported in the absence of zinc 9c , indicating that PBT2 did not inhibit ascorbate oxidation under these reaction conditions.

Since both TMDQ20 and compound 1 efficiently inhibited the ascorbate oxidation induced by Cu-Aβ1-16 in the presence of 3 mol equiv of Zn 2+ , we decided to investigate the same reaction with a very large excess of zinc, 100 mol equiv of zinc with respect to copper, a Zn/Cu ratio that has been proposed to be biologically relevant. 5 In these conditions, Cu-Aβ1-16 induced 48% of ascorbate oxidation in 5 min, and 86% in 10 min (Figure 5, trace a). This is a moderate increase of ROS production compared to the result obtained in the absence of zinc (trace b, 39% and 63% in 5 min and 10 min, respectively). A previous indicated that a single zinc equivalent had no effect on Cu-Aβ induced ROS production, even though the coordination of Zn 2+ impacts the Cu(II) and Cu(I) binding sites of Aβ, and may lead to the formation of a hetero-bimetallic Cu-Zn-Aβ complex, without release of Cu. 17 Our result indicates that a 100-fold excess of Zn 2+ resulted in the release of a fraction of Cu 2+ from Cu II -Aβ1-16, which is consistent with a higher ascorbate oxidation. So, the presence of a local huge excess of zinc with respect to Cu-Aβ may significantly exacerbate the redox toxicity of Cu-Aβ in vivo.

In order to evidence the fast action of such N4-tetradendate ligands on the kinetics of the ascorbate oxidation induced by CuCl2/Aβ1-16/ZnCl2 (1/1/100 molar ratio), we added one of these copper chelators, TDMQ20 or 1 (1 mol equiv), after 6.6 min of reaction (Figure 5, traces d and e, respectively). In both cases, the oxidation of ascorbate was immediately and completely stopped, indicating that stoichiometric amounts of TDMQ20 and 1 were able to fully inhibit the ROS production induced by the mixture Cu(II)/Aβ1-16/Zn(II) = 1/1/100. The slopes of traces d and e after 6.6 min (∆Abs265 = 0.002 / min) were not significantly different to that of ascorbate autoxidation (Figure 4, trace a, ∆Abs265 = 0.003 / min). Figure 5. UV-visible (265 nm) kinetic spectra of ascorbate oxidation under air, in the presence of CuCl2/Aβ1-16/ZnCl2 = 1/1/100, then addition of L (1 mol equiv) at 6.6 min (arrow): L = no ligand (trace a), L = TDMQ20 (trace d), L = 1 (trace e). The ascorbate oxidation in the presence of CuCl2/Aβ1-16, (trace b), or CuCl2 alone (trace c) are given for comparison. The increase of absorbance upon addition of TDMQ20 (d) or 1 (e) at 6.6 min, is due to the own absorbance of Cu II -TDMQ20, or Cu II -1 at 265 nm. For detailed experimental conditions, see Supporting Information.

These tetradendate ligands TDMQ20 and 1, that specifically chelate Cu(II) in a N4-square planar fashion, are able to immediately and completely stop the in vitro ongoing aerobic oxidation of ascorbate induced by Cu-Aβ1-16, even in the presence of 100 mol equiv of Zn(II) with respect to copper ions, while the N2O2-tetradentate ligand L2, able to chelate copper with metal/ligand stoichiometry = 1/1, was efficient only in the presence of a stoichiometric amount of zinc. 5 PBT2, whose copper complex has a N3O2 coordination sphere, and a metal/ligand stoichiometry = 1/2, 16 In conclusion, such N4-tetradentate copper-specific ligands, like TDMQ20 and 1, are suitable for a square-planar scaffold with Cu(II) and have the ability to largely stabilize Cu(II) over Cu(I). This structural feature makes them efficient inhibitors of the in vitro ROS production induced by Cu-Aβ, including in a zinc-rich physiological environment. Finally, the calculated logP values are 2.5 and 3.5, for TDMQ20 and 1, 18 respectively, in the optimal range expected for crossing the blood-brain barrier. [START_REF] Hansch | Substituent constant for correlation analysis in chemistry and biology[END_REF] The design of copper ligands able to counteract metal deregulation and deleterious consequences in AD brain requires a good understanding of their coordination properties. Pharmacological studies on TDMQ20 are in progress.

Captions for Figures

Figure 1 .

 1 Figure 1. Structures of the N4-tetradentate ligands, mono-8-aminoquinoline TDMQ20 and bis(8-aminoquinoline) 1, of the N2O2 ligand L2, and of the bi/tridentate 8-hydroxyquinoline PBT2.

Figure 3 .

 3 Figure3. Extraction of copper from Cu II -TDMQ20 by glutathione (GSH) in the presence of 1 mol equiv of Zn 2+ . The spectrum of the free ligand TDMQ20 (dashed black line) is given for comparison. Insert: Variation of the absorbance at 259 nm as a function of the number of GSH mol equiv with respect to Cu-TDMQ / Zn 2+ . For detailed experimental conditions, see Supporting Information.

Figure 4 .

 4 Trace(a) shows the ascorbate autoxidation, without any additive. In the presence of CuCl2, ascorbate was fully oxidized in only 5 min (trace b). Cu-Aβ1-16 induced ca. 40% of ascorbate oxidation in 5 min, and full oxidation occurred in ca. 25 min (trace c). In the presence of 3 mol equiv of Zn 2+ , TDMQ20 or ligand 1 drastically inhibited the ascorbate oxidation induced by Cu-Aβ1-16, the oxidation of ascorbate in 30 min being only 15-16% in both cases (traces d and e, respectively), very close to the ascorbate autoxidation (8%, trace a).

Figure 4 .

 4 Figure 4. UV-visible (265 nm) kinetic spectra of ascorbate oxidation under air, in the presence of a) no additive (ascorbate autoxidation), b) CuCl2, c) CuCl2/Aβ1-16, 1/1.1, d) CuCl2/Aβ1-16/ZnCl2/TDMQ20 = 1/1.1/3/1.1, e) CuCl2/Aβ1-16/ZnCl2/1 = 1/1.1/3/1.1, f) CuCl2/Aβ1-16/ZnCl2/PBT2 = 1/1.1/3/2.2. For detailed experimental conditions, see Supporting Information.

  failed to inhibit oxygen reduction either in the presence or in the absence of zinc. These N4-tetradendate ligands TDMQ20 and 1 have a high selectivity for Cu(II) with respect to Zn(II) {log [Kapp Cu-L / Kapp Zn-L] > 12}, 9a,10 much higher than that of Aβ {log [Kapp Cu-Aβ / Kapp Zn-Aβ] = 4-5}, 5b and thus fulfill the criteria for biologically pertinent copper chelators mentioned in the introduction. The obtained results indicated that both ligands TDMQ20 and 1 have the potential to reduce or inhibit in vivo ROS production induced by Cu-Aβ, even in a very rich zinc physiological environment.

1 ,

 1 Figure 1. Structures of tetradentate mono-8-aminoquinoline TDMQ20 and bis(8aminoquinoline 1, of the N2O2 ligand L2 and of the bi/tridentate 8hydroxyquinoline PBT2.

Figure 2 .

 2 Figure 2.Extraction of copper from Cu II -Aβ1-16 using 1 mol equiv of TDMQ20, in the presence of 1 or 3 mol equiv of Zn 2+ (full red trace, and dashed red trace, respectively) with respect to copper monitored by UV-visible spectrophotometry.

Figure 3 .

 3 Figure 3. Extraction of copper from Cu II -TDMQ20 by glutathione (GSH), in the presence of 1 mol equiv of Zn 2+ , monitored by UV-visible spectrophotometry. The spectrum of the free ligand TDMQ20 (dashed black line) is given for comparison. Insert: Variation of the absorbance at 259 nm as a function of the number of GSH mol equiv with respect to Cu-TDMQ20 / Zn 2+ .

Figure 4 .

 4 Figure 4. UV-visible (265 nm) kinetic spectra of ascorbate oxidation under air, in the presence of a) no additive (ascorbate autoxidation), b) CuCl2, c) CuCl2/Aβ1-16, 1/1.1, d) CuCl2/Aβ1-16/ZnCl2/TDMQ20 = 1/1.1/3/1.1, e) CuCl2/Aβ1-16/ZnCl2/1 = 1/1.1/3/1.1, f) CuCl2/Aβ1-16/ZnCl2/PBT2 = 1/1.1/3/2.2.

Figure 5 .

 5 Figure 5. UV-visible (265 nm) kinetic spectra of ascorbate oxidation under air, in the presence of CuCl2/Aβ1-16/ZnCl2 = 1/1/100, then addition of L (1 mol equiv) at 6.6 min (arrow): L = no ligand (trace a), L = TDMQ20 (trace d), L = 1 (trace e). The ascorbate oxidation in the presence of CuCl2/Aβ1-16, (trace b) CuCl2 (trace c) are given for comparison. The increase of absorbance upon addition of TDMQ20 (d) or 1 (e) at 6.6 min, is due to the own absorbance of Cu II -TDMQ20, or Cu II -1 at 265 nm.

Acknowledgements

This work was supported by CNRS, NSFC (YL, n° 21502023), the Guangdong Province (Program for Innovative Research Teams, BM n° 2050205) and GDUT (BM, grant n° 220418037), CH acknowledges ERC StG-638712 for funding.

Scheme and text for Table of Contents

N4-tetradentate copper-specific ligands generate Cu(II) square planar complexes and stabilize Cu(II) but not Cu(I). This feature makes them able to regulate copper homeostasis and to efficiently inhibit the in vitro oxidative stress induced by Cu-Aβ, including in a zinc-rich physiological environment like Alzheimer's disease brain. The design of ligands able to counteract copper deregulation and deleterious effects in AD brain must take into account their coordination modes.


N C Friggens

† C Duvaux-Ponter

† M P Etienne

T Mary-Huard

P Schmidely

Characterizing individual differences in animal responses to a nutritional challenge: Toward improved robustness measures

Keywords: dairy goat, clustering method, adaptive capacity, resilience, lactation

The aim of this paper was to explore the variation between individuals in the response to and recovery from a nutritional challenge, the repeatability of responses between lactation stages, and the use of shape-clustering methods to classify similar individuals. Sixteen dairy goats were exposed to a 2-d nutritional challenge (underfeeding) at 2 different stages of lactation. Each challenge consisted of a 7-d control period with standard total mixed ration (TMR), 2 d of straw-only feeding, and a 10-d recovery period on the TMR. All feeds were offered ad libitum, as was water. The first challenge was in late lactation on primiparous goats (mean days in milk = 249), and the second challenge was carried out on the same goats early in the following lactation (mean days in milk = 28). The main energetic response traits dry matter intake (DMI), milk yield, body weight, milk fat and protein contents, and plasma glucose, fatty acids (NEFA), β-hydroxybutyrate (BHB), urea, and insulin, were measured daily throughout. A clustering procedure linked to a piecewise mixed model was used to characterize different types of response. As expected, straw feeding caused a large decline in DMI and milk yield, and substantial increases in milk fat and milk protein composition, relative to the prechallenge period on the control TMR. For both DMI and milk yield, the slope of the response, and hence the size of the drop, was strongly related to the prechallenge values, indicating that these 2 measures were tightly constrained by the challenge. Regression slopes between lactation stages for responses to the same nutritional challenge were significant for DMI, milk protein content, plasma BHB and urea, and body weight, indicating that within-animal responses in late and early lactation were repeatable. The clustering procedure generally performed well, classifying both scaling differences and differences in shape. The extent of reranking of cluster designations between late lactation and the following early lactation period was examined. For milk yield, DMI, body weight, and urea, relatively little reranking occurred; the numbers of goats not changing class number were 10, 12, 10, and 13, respectively. In contrast, for milk contents of fat and protein, and also for BHB, no clear association was found between late and early lactation class numbers. For NEFA and glucose, these comparisons were not relevant because either the vast majority of goats were in 1 cluster (NEFA) or because an outlier goat skewed the cluster designation (glucose in late lactation). For insulin, 9 out of 16 goats kept the same rank.

INTRODUCTION

The ability of an animal to respond to, and recover from, environmental challenges is an increasingly important trait. This adaptive capacity is a key component of animal robustness [the ability to maintain life functions in the face of constraining environments [START_REF] Kitano | Biological robustness[END_REF]] in the context of the challenges facing future livestock production. The diversity of environments to which livestock will be exposed is expected to increase as increasing food demand leads to further exploitation of marginal land (Bocquier and GonzálesGarcia, 2010). In addition, climate change is driving an increase in the variability over time in environmental conditions with a higher frequency of extreme conditions [START_REF] Hansen | Perception of climate change[END_REF]. These 2 trends will place increasing demands on the animal's adaptive capacity.

Although there is now general agreement that adaptive capacity is a complex trait conferred by a combination of several underlying components [START_REF] Strandberg | The role of environmental sensitivity and plasticity in breeding for robustness: Lessons from evolutionary genetics[END_REF], there is little agreement on how to quantify it, and on which are the key biological components [START_REF] Friggens | Vers des définitions opérationnelles de la robustesse s'appuyant sur des faits biologiques: l'exemple de la nutrition[END_REF]. This clearly presents a difficulty for characterizing both adaptive capacity and thus robustness. Marked differences are present between breeds and strains in performance responses to environmental challenges [START_REF] Horan | The effect of strain of Holstein-Friesian cow and feeding system on postpartum ovarian function, animal production and conception rate to first service[END_REF][START_REF] Delaby | Effect of different feeding strategies on lactation performance of Holstein and Normande dairy cows[END_REF]. It has also been shown that individual variation in performance, and in the use of body reserves, can be adversely affected in harsh environments [START_REF] Puillet | An individualbased model simulating goat response variability and long-term herd performance[END_REF]. In this context, being able to characterize robustness at the level of the individual animal would be valuable for refining both management and selection strategies with respect to the anticipated increase in diversity of livestock systems.

The present study fits within an overall aim to develop ways to characterize the robustness and adaptive capacity of animals that can then be adapted to become operational methods for quantifying robustness on farm. [Robust animals are those that are able to be healthy and perform well under a wide range of environmental conditions [START_REF] Amer | Turning science on robust cattle into improved genetic selection decisions[END_REF]]. We believe that this is a timely and realistic aim given the increasing panel of automated measures that are becoming available on farm in the context of precision livestock farming [START_REF] Rutten | Invited review: Sensors to support health management on dairy farms[END_REF]. Such measures provide the opportunity to describe the dynamic features of responses of individual animals to environmental perturbations [START_REF] Codrea | Differential smoothing of time-series measurements to identify disturbances in performance and quantify animal response characteristics: An example using milk yield profiles in dairy cows[END_REF]. The advances in on-farm technology also mean that such response dynamics can be measured across multiple performance and physiological traits. Preliminary studies indicate that it is possible to derive operational definitions of components of robustness that are quantifiable. For example, it has been shown that biological entities such as degree of infection can be derived from multivariate measures of the animal's disease response [START_REF] Højsgaard | Quantifying degree of mastitis from common trends in a panel of indicators for mastitis in dairy cows[END_REF].

In the area of nutritional status and nutritional challenges, markers to monitor key metabolites such as BHB have been in use for a considerable time and now exist in automated systems [START_REF] Nielsen | Predicting the risk of ketosis in dairy cows using online measurements of beta-hydroxybutyrate in milk: A biological model[END_REF]. Several indexes are also now available for the purpose of monitoring to identify at-risk animals with respect to health status [START_REF] Bramley | The definition of acidosis in dairy herds predominantly fed on pasture and concentrates[END_REF][START_REF] Moyes | Generation of an index for physiological imbalance and its use as a predictor of primary disease in dairy cows during early lactation[END_REF]. In the present study, we wanted to extend this approach from a monitoring perspective toward phenotyping (i.e., to develop ways to better understand and characterize between animal variation in adaptive capacity using multivariate measures).

Accordingly, an experiment was carried out in which time-series measurements of behavior, performance, and metabolites were made in dairy goats exposed to a nutritional challenge, at 2 different stages of lactation. A major issue of this work relates to the exploration of the resulting data because the focus of this work is not to quantify the overall effect of the nutritional challenge but rather to examine individual differences in responses. Thus, we are interested in the portion of the variation that is usually eliminated (i.e., assigned to the error term) in standard ANOVA approaches. Given this, and the fact that there was no strong a priori rationale for which measures (or aspects of measures) are key components of adaptive capacity, an exploratory data analysis was used in this paper. Therefore, the aim of this paper was to prepare the ground (proof of concept) by providing the methodological basis for quantifying biologically meaningful descriptions of adaptive capacity phenotypes from multivariate measures. More specifically, the aim was to develop a methodology to explore the variation between individuals in the response and recovery to the challenge, the repeatability of responses, and shapeclustering methods to classify similar individuals.

MATERIALS AND METHODS

Animals, Feeds, and Design

Sixteen dairy goats, housed in individual pens (1.2 m by 0.75 m), were exposed to a 2-d nutritional challenge (underfeeding), at 2 different stages of lactation. Each challenge consisted of a 7-d control period with standard TMR, 2 d of straw only feeding, and a 10-d recovery period on the TMR. Prior to the start of each challenge period, the goats had received the standard TMR for at least 15 d. All feeds were offered ad libitum, as was water. The first challenge was in late lactation (mean DIM = 249, SD = 2.8, all goats were primiparous), and the second challenge was carried out on the same goats early in the following lactation (mean DIM = 28, SD = 3.3). The goats were housed in individual pens each with their own feed trough; they were milked twice per day and feed was distributed twice daily shortly after milking. Animals were cared for and handled in accordance with the French legislation on animal experimentation and European Convention for the Protection of Vertebrates Used for Experimental and Other Scientific Purposes (European Directive 86/609). The experiment was carried out between October 12, 2009, andFebruary 14, 2010.

The standard TMR (on a DM basis) consisted of 20% chopped hay, 30% chopped dried alfalfa (Rumiluz, Désialis, Paris, France), 30% sugar beet pulp, and 20% of a commercial dairy concentrate (containing 18% maize, 14% sugar beet pulp, 12% sunflower meal, 10% wheat, 10% soybeans, 9% rapeseed meal, 6% soybean meal, 4% wheat distillers grains, 3.5% linseed, 3% pea seed, 1% rapeseed oil, 3% molasses, and 6.5% mineral and vitamin premix). The standard TMR (DM content 98.3%) had a measured content (on a DM basis) of 12.0% CP, 5.6% starch, 8.9% ash, 37.2% NDF, 20.0% ADF, and 2.5% ADL. The chopped straw (DM content 97.2%) contained (on a DM basis) 2.5% CP, 8.4% ash, 74.1% NDF, 43.1% ADF, and 5.1% ADL.

Measurements

From the beginning of the prechallenge control period until the end of the postchallenge recovery period measurements of feed intake and BW were made daily. Milk yield was recorded and milk samples were taken at each milking for standard milk composition measures (Fossomatic, Hillerød, Denmark). Blood samples were taken daily in heparinized tubes after the morning milking and before feeding (average 20 min pre-feeding), immediately centrifuged for 10 min at 3,000 × g at 4°C, and the plasma was frozen (-20°C) for subsequent analysis, in duplicate, of urea (intra-and interassay CV: 1.7 and 3.2%), glucose (intra-and interassay CV: 0.8 and 1.8%), BHB (intraand interassay CV: 1.4 and 2.0%), and fatty acids (NEFA; intraand interassay CV: 4.1 and 8.2%) using a Cobas Mira-Analyzer (Roche, Mannheim, Germany) with commercial kits for urea (urea/BUN, Roche), glucose (Gluco-quant, Glucose/ HK, Roche), NEFA [NEFA-HR(2), Wako Chemical GmBH, Neuss, Germany], and BHB by the procedure of [START_REF] Barnouin | Micro-dosage automatisé sans déprotéinisation du?hydroxybutyrate plasmatique chez les bovins. (Automated micromethod for determining 3-hydroxybutyrate in bovine plasma without acid precipitation)[END_REF]. Insulin was measured using an ELISA kit (10-1202-01; Mercodia AB, Uppsala, Sweden; intra-and interassay CV: 1.2 and 4.5%).

Feeds were analyzed using the following standard methods: DM estimated from water content [START_REF] Iso | Animal Feedingstuffs. Determination of moisture and other volatile matter content[END_REF], ash [START_REF] Iso | Animal Feedingstuffs. Determination of crude ash[END_REF], and starch [START_REF] Iso | Animal feeding stuffs. Enzymatic determination of total starch content[END_REF]. Total N was determined by the Dumas technique [START_REF] Sweeney | Comparison of Leco-FP-228 nitrogen determinator with AOAC copper catalyst Kjeldahl method for crude protein[END_REF]. Cell wall content was estimated by the NDF method of [START_REF] Van Soest | Use of detergents in the analysis of fibrous feeds. IV. Determination of plant cell-wall constituents[END_REF], modified by [START_REF] Giger | Etude de l'influence préalable de différents traitements amylolytiques sur la teneur en résidu NDF d'aliments du bétail. (Effect of different amylolytic pretreatments on NDF content in feedstuffs)[END_REF] with the use of a heat-stable αamylase but without sodium sulfite and decalin. Lignocellulose (ADF) and lignin (ADL) were obtained using a sequential approach on the NDF residue [START_REF] Giger | Etude de l'influence préalable de différents traitements amylolytiques sur la teneur en résidu NDF d'aliments du bétail. (Effect of different amylolytic pretreatments on NDF content in feedstuffs)[END_REF]. All the cell wall components were expressed on an ash-free basis.

Statistical Analyses

In the present study, we focused on the main energetic response traits (i.e., DMI, milk yield, BW, milk fat and protein contents, and plasma glucose, NEFA, BHB, urea, and insulin). The statistical analysis was carried out on each trait separately and consisted of 2 steps. In the first step, the response profile of each goat was described using a piecewise model with 4 parameters. These parameters that characterize the response recovery profile allowed us to describe (1) the relationship between prechallenge and response to challenge, and (2) the relationship between response to challenge and recovery. In a second step, all response profiles were jointly analyzed to identify cluster of animals exhibiting similar profiles.

Step 1: Individual Analyses. For each trait, and each lactation stage (late vs. early), each of the 16 individual time series were characterized separately using the following piecewise model:

Yt =V1 ×I(t≤0) +V2 ×t×I(0<t≤2) + (V3 ×t +V4 ×t×2)I(2<t≤4) +V5 ×I(t>4 +Et ,
where Yt is the trait measured at time t and the coefficients Vn are estimated by the model, V1 is the prechallenge level (assumed to be constant throughout the prechallenge period), V2 is the linear slope of the response during the 2d challenge, V3 is the linear component of the recovery, V4 is the quadratic component (deceleration) of the recovery, and V5 is the postchallenge level that was assumed to be complete and constant 6 d after the beginning of challenge (confirmed by visual inspection). The error terms Et are assumed to be Gaussian and independent. Indicator variable I{condition} is a dummy variable with value 1 if condition is true, and 0 otherwise. The final postchallenge level was not assumed to be the same as the prechallenge level. In addition, parameters V1 to V5 were estimated subject to the constraint of profile continuity. The model and its associated parameters are summarized in Figure 1. The choice of model, along with its inherent assumptions, was made on the basis of both a visual inspection of the response/recovery profiles and keeping in mind a parsimony principle (keeping an acceptable ratio between the number of parameters and the number of data points in an individual time series). The set of parameters V1 to V5 was then used to characterize the profile of each goat.

The relationship between prechallenge and response to challenge was analyzed using the following model:

V2ij =μ+a×LACTi +b×V1ij +c×LACTi ×V 1ij +Gj +Eij ,
where V2ij is coefficient V2 for goat j and lactation stage i, V1ij is the mean centered coefficient V1 [i.e., V1ij = V1ij -mean(V1i)]. Parameter μ is the intercept for early lactation, LACTi is a dummy variable with values 0 and 1 for early and late, respectively, a is the coefficient for difference in intercepts between lactation stages, b is the slope coefficient measuring the link between V1 and V2, and c is the coefficient for difference in slopes be- tween lactation stages. Last, Gj is a Gaussian random effect with variance σG 2 accounting for possible correlations between coefficients measured on the same goat, and Eij is the residual with variance σ 2 . All random effects are assumed independent. This model was applied to each of the following traits: DMI, milk yield, milk fat content, milk protein content, BW, plasma glucose, plasma NEFA, plasma BHB, plasma urea, and plasma insulin. As across all traits c was found not to be significant, this interaction term was dropped from the model.

The relationship between response to challenge and recovery was analyzed using a model similar to model (1), with the following differences:

• The variable to be explained is now the rate of deceleration in recovery V4ij, • The explanatory variable is now V2ij = V2ijmean(V2i).

Note that V4 is the quadratic component of the recovery and characterizes (half) the rate of deceleration in the recovery. This variable was chosen because the absolute rate of recovery (i.e., the combination of the quadratic and linear components of the recovery) is strongly correlated with the amplitude of the response (r > 0.8 for most traits). One animal in late lactation was excluded as an extreme outlier. As previously, c was found not to be significant across all traits; therefore, this interaction term was dropped from the model. To assess repeatability between lactation stages, within trait, the following linear regression was used:

V22j =a +b×V21j +Ej .
Parameters a and b are the intercept and slope for the relation, respectively. Note that perfect repeatability would result in b = 1. A similar analysis was performed on V4.

Step 2: Clustering Analysis. To identify goats exhibiting a similar response profile to the challenge, a clustering analysis was performed based on a 2-stage model. Each goat is assumed to belong to one of the K clusters, where the number of clusters K is fixed and assumed to be known. We note πk as the proportion associated with class k (i.e., the prior probability to belong to class k). For a goat j that belongs to cluster k, one has Yjt =V1k ×I(t≤0) +V2k ×t×I(0<t≤2) +(V3k ×t +V4k ×t×2)I(2<t≤4) +V5k ×I(t>4) +Ejt , which corresponds to model (1) previously described, except that parameters V1 k ,…, V5 k are common to all animals belonging to class k, and that the variance of Ejt is assumed to be the same for all animals. The goal of the analysis is then to estimate the coefficients (πk, V1 k ,…, V5 k ) but also the posterior probabilities τjk of each animal j to belong to class k. This can be done thanks to the expectation maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. In practice, the number of clusters K being unknown, it was chosen within the range of cluster numbers from 2 to 5 based on the integrated completed likelihood (ICL) criterion [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. Both the EM algorithm and the ICL criterion are described in the Appendix. To allow subsequent comparison of cluster class attributions between traits and lactation stages, the analyses were rerun using a pre-specified fixed number of 3 clusters, which corresponded to the cluster number most frequently found across traits.

Comparison across traits was carried out on the cluster classifications using the adjusted Rand index (ARI; [START_REF] Rand | Objective criteria for the evaluation of clustering methods[END_REF]. This is a pair-counting-based method that calculates across any 2 traits whether the cluster classifications of a given pair of goats in the first trait match the cluster classifications of the same pair of goats in the second trait. When this is done across all possible pair combinations of goats, it provides an index of clustering agreement between traits that makes no a priori assumptions about mapping of one cluster classification onto the other, and that is adjusted for the fact that purely random clustering would still create some matches. Using ARI, a perfect match would give an ARI = 1 and purely random matches would give an ARI = 0. In addition, principal components analysis (PCA) was carried out on the cluster classification membership. The procedure was the following: for each initial variable, the clustering score was obtained (with values 1, 2, and 3 corresponding to a clustering into 3 classes). This clustering score was then replaced by 3 dummy variables: Ind1, Ind2, Ind3, where Indjk = 1 if animal j belongs to cluster k, 0 otherwise; thus, no a priori relation was assumed between cluster classes. A PCA was performed on the set of all dummy variables corresponding to all traits, for each lactation stage.

RESULTS

The average time trends of DMI, milk yield, milk fat, and milk protein contents through the prechallenge, challenge, and postchallenge periods in both late and the following early lactation are shown in Figure 2. The The prechallenge levels were estimated within a piecewise mixed model (see Materials and Methods), the ranges followed by the standard deviations are given in parentheses.

2 One extreme animal excluded.

average prechallenge values for all traits in late and early lactation are given in Table 1.

Response to a Nutritional Challenge (2 d Straw-Only Feeding)

The slopes of response (V2; see Figure 1) to the change from a standard TMR to straw-only feeding in early and late lactation are presented in Table 2. As expected, straw feeding caused a large decline in DMI and milk yield, and substantial increases in milk fat and milk protein composition, relative to the prechallenge period on the control TMR. For both DMI and milk yield, the slope of the response, and hence the size of the drop, was strongly related to the prechallenge values indicating that these 2 traits were tightly constrained by the challenge. The responses in DMI expressed relative to preperiod values (V2/V1) were -0.76 and -0.88 in late and early lactation, respectively. The proportional responses in milk yield were -0.75 and -0.62 in late and early lactation, respectively. In both cases, the differences between lactation stages were significant. For all other traits (milk fat, milk protein, NEFA, BHB, glucose, urea, insulin), the prechallenge values had no significant effect.

For intake, milk protein content, plasma BHB and urea, and BW, significant slope coefficients were observed from regressions relating the rates of response (V2) between lactation stages, albeit with modest R 2 (Table 3). This indicates that within-animal responses in late and early lactation were repeatable (Figure 3). For plasma metabolites other than BHB and urea, the relationship between lactation stages in rates of response (V2) was not significant, indicating that the other metabolic adaptations to nutritional challenge were not repeatable within animal. With the exception of glucose, insulin, and milk protein content, all intercepts of these regressions were significant, indicating The response slopes (V2) and prechallenge levels (V1) were estimated using a piecewise mixed model (see Materials and Methods). Results are shown in equation form: V2 = average early lactation value + late lactation adjustment + coefficient × V1. In the analysis, the prechallenge levels were centered (i.e., expressed as the difference from the mean prechallenge level) and thus the average values of the response slopes in early and late lactation are given directly in the table. Coefficients significantly different from 0 are shown in bold (P < 0.05). As no significant interactions were found between prechallenge level and stage of lactation for any of the measures, this term was excluded from the model. The R 2 was calculated from the residual of the full model (resf) and the residual from the null model (resn) as: (resnresf)/resn. that stage of lactation significantly affected the average size of response (Table 3).

Recovery from a Nutritional Challenge (2 d Straw-Only Feeding)

Given the inherently symmetrical nature of a response recovery profile, the size of the recovery (difference from peak to postchallenge baseline) is highly correlated with the size of response, consequently the stage of lactation effects on size of recovery (not shown) were essentially the same as those for size of response (Table 3). Thus, the quadratic component of the recovery [V4; i.e., (half) the rate of deceleration in the recovery] was focused on. As shown in Table 4, this rate (V4) was significantly greater in early lactation (relative to late lactation) for DMI, milk yield, and BW. Conversely, the rate of deceleration in recovery (V4) was significantly less in early lactation for milk fat content, NEFA, BHB, urea, and insulin. The correlations between the rate of deceleration and the size of the response were above 0.5 for the majority of traits. This reflects the fact that within a similar recovery time, the greater the response, the "sharper" the curve needed to return to a plateau level. Indeed, the slope of response (V2) had a significant effect on the rate of deceleration in recovery (V4) for milk yield, BW, glucose, NEFA, BHB, and urea (Table 4). Thus, the relative rate of recovery was explored using the quadratic component of recovery/ slope of response (V4/V2). This relative rate of recovery In the analysis, the V2 values were centered (i.e., expressed as the difference from the mean value), and thus the average values of the rate of deceleration in recovery in early and late lactation are given directly in the table. Coefficients significantly different from 0 are shown in bold (P < 0.05). As no significant interactions were found between V2 and stage of lactation for any of the measures, this term was excluded from the model. The R 2 was calculated from the residual of the full model (resf) and the residual from the null model (resn) as: (resnresf)/resn.

was significantly greater than zero for DMI, milk yield, milk fat content, BW, plasma glucose, NEFA, and urea (results not shown). This indicates that for these traits the recovery was quicker than would be expected for the size of the initial response. No traits were found for which the recovery was more gradual than expected according to the size of the response. For milk yield, milk fat content, and BW, a significant difference was found between early and late lactation relative recovery rates. However, when examined for repeatability between late and early lactation, no significant effect was found on early lactation values of relative recovery rate in late lactation.

Clustering of Response-Recovery Profiles

The above results focused on different segments of the response-recovery profiles, here the profiles are examined in their entirety (see Figures 1 and4). When optimizing cluster number based on minimizing the ICL criterion, the clustering procedure identified 4 clusters for milk yield, BW, milk protein content, milk fat content in early lactation, BHB in early lactation, and insulin in early lactation. It identified 2 clusters for NEFA in early lactation and urea in late lactation; for the remaining traits, it identified 3 clusters. As can be seen from Figure 4, the clustering procedure generally performed well, classifying both differences that can be described as scaling differences (Figure 4a) and differences in shape (Figure 4b).

However, 2 issues should be noted; the first is that the cluster shapes are limited by the assumptions of the underlying piecewise model (i.e., a plateau value before the challenge and also 6 d postchallenge, a linear slope during the challenge, and a quadratic function for the recovery postchallenge). The second issue is that outlier animals often get assigned a class of their own (e.g., Figure 4a and4c). This had the consequence that when the clustering was done using a fixed number of clusters (to allow comparison of cluster designation between traits the number was fixed to 3), the remaining variation between animals other than the outlier could be forced into only 2 clusters, resulting in a loss of type description. As shown in Figure 5, this was the case for milk protein content in early lactation, and to a lesser extent BHB and insulin also in early lactation. The issue of outlier animals should be further evaluated in studies with larger numbers of animals. Because the clusters were numbered in ascending order of prechallenge level, it was possible to examine the extent of reranking of cluster designations between late lactation and the following early lactation period (whether the relative order of animals in cluster designations 1, 2, and 3 is changed). This was done in 2 ways, by visual inspection using the clustering without constraint on the cluster numbers, and using the ARI comparison with fixed cluster numbers. For milk yield, DMI, BW, and urea, relatively little reranking was observed and the corresponding ARI values were 0.3, 0.02, 0.41, and 0.33, respectively. (The range of ARI values across all comparisons for all traits was from 0 to 0.41). The apparent discrepancy between the ARI value for the late versus early lactation comparison for DMI relative to those of milk yield, BW, and urea was caused by the fact that although little reranking was found for DMI, 8 out of the 16 goats did decrease their class number by 1 (with a further 6 not changing class number). For milk yield, BW and urea, the numbers of goats not changing class number were 10, 10, and 13, respectively. In contrast, for milk contents of fat and protein, and also for BHB, no clear association was found between late and early lactation class numbers (ARI values of 0.02, 0.05, and 0.04, respectively). For NEFA and glucose, these comparisons were not relevant because either the vast majority of goats were in 1 cluster (NEFA) or because an outlier goat skewed the cluster designation (glucose in late lactation). For insulin, 9 out of 16 goats kept the same rank (ARI = 0.1).

When comparing across traits, the highest ARI values in late lactation were for associations between milk yield and milk protein content (0.34), DMI and BW (0.23), DMI and glucose (0.21), milk contents of fat and protein (0.22), milk protein and BHB (0.39), and glucose and urea (0.33). The highest ARI values in early lactation were for associations between milk yield and DMI (0.19), DMI and BW (0.29), BW and insulin (0.28) , and milk fat and BHB (0.19).

To provide a better visualization of the multivariate aspect, we used PCA on the cluster classes for all the andd) in early lactation. For milk protein content using the variable cluster number assignment (a), classes 1, 2, 3, and 4 contained 2, 8, 4, and 2 individuals, respectively. The 2 goats in class 1 combine with those in class 2 to give a new class 1 in the fixed cluster number assignment (c). For BHB, the variable cluster classes 2 and 4 (containing 5 and 3 goats) are merged into class 3 for 5 out of the 8 goats under the fixed cluster classes. traits (Figure 6). In this PCA, within each trait, each cluster class was included as a separate variable [i.e., avoiding any a priori assumptions of relativity between cluster classes (linear or otherwise)]. For example, in Figure 6b it can be seen that DMI cluster classes 1, 2, and 3 fall in 3 distinct quadrants of the plot, which would not have occurred if DMI had been included as one variable with 3 levels (1, 2, and 3). The results of the PCA presented in Figure 6 are for a fixed cluster number across traits but were essentially the same when using variable cluster numbers. From the loading plots for late and early lactation (Figures 6a and6b), some patterns can be discerned. In both late and early lactation, the different cluster classes for milk fat content, milk protein content, and plasma NEFA are each positioned approximately on a linear transect, all running from the lower right to the upper left quadrant. In late lactation, the most closely correlated of these are NEFA and milk fat content, whereas in early lactation milk fat and protein contents are more closely associated. In late lactation, another transect can be discerned for DMI and milk yield, running from lower left to upper right, but in general no strong grouping is present of traits in the loadings plot, or of individual goats in the scores plot (Figures 6a and6c, respectively). However, in early lactation, clearer positioning is present of traits into groups in the loadings plot (Figure 6b). One of these groups (top right, Figure 6b) is essentially associated with only one animal (goat 12, Figure 6d). Within the 2 other groups, contrasts are present between milk yield and milk composition cluster numbers, and these groups lie roughly on the transect lines for plasma NEFA, milk fat, and milk protein.

DISCUSSION

General Response to the Nutritional Challenge

The overall effects of a 2-d nutritional challenge, replacing the standard TMR with straw, reported in this study using goats were broadly similar to other reports in the literature, although these come mainly from dairy cows. The decrease in milk yield in the present study (75% in late lactation and 62% in early lactation) was similar to the 57% decrease reported for black Moroccan goats that were feed deprived for 48 h [START_REF] Hossaini-Hilali | Fluid balance and milk secretion in the fed and feed-deprived black Moroccan goat[END_REF] and was equivalent to that found with a 2-d fast in dairy cows [START_REF] Mcguire | Nutritional modulation of the somatotropin insulin-like growthfactor system-Response to feed deprivation in lactating cows[END_REF][START_REF] Agenäs | Effects of feeding intensity during the dry period. 1. Feed intake, body weight and milk production[END_REF]. Using a more extreme feed restriction, 6 d of total fast, [START_REF] Reid | Fasting and refeeding in lactating dairy cow-Recovery of milk-yield and blood-chemistry following a 6-day fast[END_REF] found that cows needed 49 d of recovery before andd). Within each measure, profile classes are considered as being independent from each other. Profile classes were derived from the combined piecewise fitting and clustering procedure (described in the Materials and Methods) using a fixed number of classes (3) across measures. To aid interpretation, the clusters were numbered, within lactation stage, in ascending order of prechallenge level. In the score plots, each goat is identified by a number, which is the same in both score plots. their milk yield was normal. [START_REF] Toerien | Duration of a severe feed restriction required to reversibly decrease milk production in the highproducing dairy cow[END_REF] found that after 24 h severe restriction cows only needed 18 h of re-feeding before they were back on pretrial levels. In the present study, following 2 d on straw, goats regained prechallenge levels after approximately 3 d. In the recent study of [START_REF] Bjerre-Harpøth | Decreased nutrient density to increase physiological imbalance for dairy cows at different stages of lactation: 1. Milk yield and component responses[END_REF], where cows had a standard TMR diluted with 60% straw for 4 d, milk yield decreased by 40% and was largely recovered in 3 d. The study of [START_REF] Bjerre-Harpøth | Decreased nutrient density to increase physiological imbalance for dairy cows at different stages of lactation: 1. Milk yield and component responses[END_REF] also reported blood and milk metabolite changes for cows grouped according to early, mid, and late lactation. Allowing for the fact that the present challenge in goats was more severe, these 2 studies have general agreement with respect to magnitudes of changes in blood and milk traits, with a couple of exceptions. In the present study, milk protein content increased during the challenge, as did milk fat content, whereas in the study of [START_REF] Bjerre-Harpøth | Decreased nutrient density to increase physiological imbalance for dairy cows at different stages of lactation: 1. Milk yield and component responses[END_REF] only milk fat content increased with milk protein either decreasing (in early lactation) or not changing (in late lactation). [START_REF] Hossaini-Hilali | Fluid balance and milk secretion in the fed and feed-deprived black Moroccan goat[END_REF] also found no increase in milk protein during feed deprivation. Further, in the present study and in the study of [START_REF] Hossaini-Hilali | Fluid balance and milk secretion in the fed and feed-deprived black Moroccan goat[END_REF] urea increased markedly (×2) during the challenge, whereas it decreased in the study of BjerreHarpøth et al. (2012). These differences may be due to the degree of alteration in milk yield relative to the drop in protein supply, and the consequences of this on dilution within milk of the milk constituents. In the study of [START_REF] Bjerre-Harpøth | Decreased nutrient density to increase physiological imbalance for dairy cows at different stages of lactation: 1. Milk yield and component responses[END_REF], the CP content of the straw-diluted diet was adjusted to keep it above 12%, which was not done in the present study. Otherwise, the prechallenge levels and responses in glucose and insulin were lower in early lactation than in late lactation in both studies, and the prechallenge levels and responses in NEFA and BHB were greater in early lactation than in late lactation in both studies. In evaluating these results and the associated individual variation, it is important to remember that these are responses to a short-term challenge, which was useful for testing the methodology proposed in this study but may differ from responses to longer-term challenges.

Exploring Individual Variation in Response-Recovery Profiles

However, the main purpose of the present study was not to quantify the average effects of a nutritional challenge but rather to explore the individual variation in responses and recoveries, and to see if this variation could be used to identify and describe animal types. To do this we primarily used a simultaneous piecewise fitting and clustering procedure. With respect to characterizing individual differences in the profiles of measures prechallenge, during the challenge, and postchallenge, the piecewise mixed model performed well, given the number of parameters to estimate (4) relative to the number of measures in a profile (18). The same was true of the clustering procedure. Statistical analysis of the mixed model parameters revealed some important aspects for predicting responses and for repeatability. With respect to predicting responses, for both DMI and milk yield a significant linear effect was found of the prechallenge level on the slope, and thus the size, of the response (Table 2, Figure 3). Further, the responses in DMI and milk were very similar between early and late lactation when expressed as proportions of the prechallenge levels. This finding is in agreement with the results of [START_REF] Bjerre-Harpøth | Decreased nutrient density to increase physiological imbalance for dairy cows at different stages of lactation: 1. Milk yield and component responses[END_REF] for cows at different stages of lactation . No significant effects were found of prechallenge level on size of response for any of the other traits.

With respect to repeatability, significant coefficients were observed for the regressions between size of response in late lactation and in the following early lactation for DMI, milk protein content, BW, plasma BHB, and urea (and tendencies for milk yield and milk fat content; Table 3, Figure 3). Thus, for these traits, high responders in late lactation were also high responders in early lactation. This finding is highly promising for the characterization of adaptive capacity as it suggests that, despite the substantial homeorhetic changes associated with the transition from one lactation to another, these measures capture some of the innate, genetically driven, variability between individuals. This was not so evident for rate of deceleration in recovery (V4), where only plasma urea had a significant repeatability (data not shown). This quadratic component of the recovery was examined rather than the absolute recovery (i.e., combined linear and quadratic components) because the differences in absolute recoveries would very largely reflect the amplitude of the response (correlation between linear slopes of response and recovery ranged from 0.49 to 0.99 for the different traits). However, examining partial components of recoveries is a relatively crude way of dealing with the correlations among the different phases of the time series (prechallenge, during the challenge, and postchallenge). Another, more elegant way to deal with these autocorrelations is to treat the time series in their entirety and examine their shapes for evidence of clusters that group different individuals. This was done using the methods described earlier permitting visualization of the different shape types, trait by trait (e.g., Figures 4 and5). However, the use of the ARI statistic to examine correlations between cluster types was found not to be particularly informative. In our opinion, this is in part due to the loss of detail when standardizing the number of clusters per trait (to 3), a necessary condition for the ARI comparison. There is an argument for repeating this type of challenge study on a much larger number of animals to increase the statistical power of the cluster comparisons. Another shortcoming of the ARI method is that although it allows pairwise comparison across traits or lactation stages (in this study), it is still difficult to visualize across multiple traits.

Principal components analysis was used for a better visualization of the multivariate aspect. In late lactation a close correlation was observed between NEFA and milk fat content, whereas in early lactation milk fat and protein contents are more closely associated. A strong relationship between NEFA and milk fat content, or milk fat:protein ratio has previously been reported, especially in early lactation [START_REF] Grieve | Relationship between milk composition and some nutritional parameters in early lactation[END_REF][START_REF] Moyes | Generation of an index for physiological imbalance and its use as a predictor of primary disease in dairy cows during early lactation[END_REF]. The clearer positioning of traits in early lactation (Figure 6b) roughly on the transect lines for plasma NEFA, milk fat, and milk protein, and the fact that that the majority of individual goats are positioned along this transect (Figure 6d), is an important multivariate feature. It indicates a possible multivariate characteristic to describe differences between individuals using measures that relate to energy partition and metabolism. Interestingly, neither BHB nor glucose cluster classes exert a strong influence in principal components 1 or 2 in early lactation, although they are more evident in principal component 3 (not shown). The ability to describe differences between individuals in this way, and track these differences across stages of lactation, suggests that this may be a means to examine multivariate homeorhetic relationships.

The use of a simultaneous piecewise fitting and clustering procedure to provide the input variables (within trait cluster classes) to the PCA provides several benefits: the variability that is captured describes the whole shape of the response/recovery, the clustering is independently done for each trait with no prior biological assumptions, and no relativity between cluster classes assumed. The lack of prior grouping hypotheses implicit in this clustering method can be seen as an advantage of this approach relative to some of the indexing methods discussed below. Thus, the results of the PCA provide a simple means to begin to interpret the individual variability in the time-series profiles of response/recovery across multiple traits. Although the PCA visualization is useful for interpretation, this is clearly only one step if we wish to proceed toward methods to index differences in adaptive capacity. A follow-up step could be to derive an index directly from the PCA loadings. An example of this approach has been published for improving quantification of uterine health status from clinical measures [START_REF] Gorzecka | A universal index of uterine discharge symptoms from calving to 6 weeks postpartum[END_REF]. Other approaches are to develop indexes directly from shape profiles without clustering [e.g., using Mahalanobis type statistics [START_REF] Cohen | A novel statistical approach shows evidence for multi-system physiological dysregulation during aging[END_REF]], or to use prior physiological knowledge to identify and combine candidate measures to, for example, index physiological imbalance [START_REF] Moyes | Generation of an index for physiological imbalance and its use as a predictor of primary disease in dairy cows during early lactation[END_REF]. However, not all of these examples may be suited to complex characteristics such as adaptive capacity that are by their nature emergent properties of multiple underlying mechanisms [START_REF] Kitano | Biological robustness[END_REF]. In this context, it is necessary to explicitly recognize the complex nature of the characteristic that we seek to index, and that we cannot directly measure it. What we have measured is the degree of perturbation resulting from a nutritional challenge, which is due to both the adaptive capacity of the animal and to the size of the challenge. The finding that the size of the response is related to the initial level of performance clearly indicates that the true size of challenge, to the internal milieu of the animal, varies from animal to animal even for an apparently constant externally applied challenge. Thus, it seems likely that the next step toward developing indexes for characteristics such as adaptive capacity is to build on multivariate results such as those presented here using models that explicitly invoke latent variables such as degree of perturbation and adaptive capacity. Such models are being developed in the time-series context and in simple form have even been applied to biological problems such as degree of infection [START_REF] Højsgaard | Quantifying degree of mastitis from common trends in a panel of indicators for mastitis in dairy cows[END_REF].

CONCLUSIONS

Responses to nutritional challenge (2 d with access to straw only) can be characterized using piecewise mixed models and shape clustering. Dry matter intake and milk yield responses were found to be predictable from prechallenge levels, and a significant relation was found between responses in late lactation and in the following early lactation for DMI, milk protein content, BW, plasma BHB, and urea. Further, significant variation was observed between individuals in their response/ recovery profiles, which suggests that defining adaptive capacity indexes from response/recovery profiles in multiple traits is possible.

Figure 1 .

 1 Figure 1. Schematic representation of the piecewise model for describing response and recovery time-series profiles relative to a nutritional challenge.

Figure 2 .

 2 Figure 2. The average time trends of (a) DMI, (b) milk yield, (c) milk fat, and (d) milk protein contents through the prechallenge (before d 0) , challenge (d 0 to 2), and postchallenge (d > 2) periods in both late (solid line) and the following early lactation (dashed line ).Table1. Average prechallenge levels in late lactation and in the following early lactation 1

Figure 3 .

 3 Figure 3. The relationship between late lactation and the following early lactation in slope of response to the same nutritional challenge (2 d of straw-only feeding) for (a) DMI and (b) milk yield.

Figure 4 .

 4 Figure 4. Classes of response-recovery profile for (a) milk yield and (b) plasma glucose in early lactation, generated by the combined piecewise mixed model and clustering procedure (see Materials and Methods). The 16 individual time trends from which these profile classes were derived are shown for comparison: (c) milk yield and (d) plasma glucose.

Figure 5 .

 5 Figure 5. Comparison of profile classes based on variable cluster number assignment (a and b) with profile classes for the same measures based on a fixed cluster number (c and d). Two measures are shown: milk protein content (a and c) and plasma BHB (b and d) in early lactation. For milk protein content using the variable cluster number assignment (a), classes 1, 2, 3, and 4 contained 2, 8, 4, and 2 individuals, respectively. The 2 goats in class 1 combine with those in class 2 to give a new class 1 in the fixed cluster number assignment (c). For BHB, the variable cluster classes 2 and 4 (containing 5 and 3 goats) are merged into class 3 for 5 out of the 8 goats under the fixed cluster classes.

Figure 6 .

 6 Figure 6. Loading (a and b) and score plots (c and d) for the first 2 components of a principal components analysis of profile classes across all measures in late lactation (a and c) and in early lactation (band d). Within each measure, profile classes are considered as being independent from each other. Profile classes were derived from the combined piecewise fitting and clustering procedure (described in the Materials and Methods) using a fixed number of classes (3) across measures. To aid interpretation, the clusters were numbered, within lactation stage, in ascending order of prechallenge level. In the score plots, each goat is identified by a number, which is the same in both score plots.

Table 1 .

 1 Average prechallenge levels in late lactation and in the following early lactation 1

	Item	Late lactation 2	Early lactation
	DMI (kg/d)	2.48 (2.02-3.29; 0.36)	2.97 (2.40-3.60; 0.35)
	Milk yield (kg/d)	1.93 (1.28-2.72; 0.46)	3.97 (2.18-5.17; 0.65)
	Milk fat content (g/kg)	35.7 (32.6-52.8; 2.9)	43.7 (36.5-49.4; 3.6)
	Milk protein content (g/kg)	36.5 (30.9-45.8; 3.8)	35.3 (30.6-40.1; 2.4)
	BW (kg)	60 (51-74; 7.4)	62 (54-78; 7.4)
	Plasma glucose (g/L)	0.60 (0.55-0.67; 0.03)	0.57 (0.50-0.63; 0.04)
	Plasma fatty acids (NEFA; μmol/L)	130 (97-149; 16)	493 (200-655 ; 136)

Plasma BHB (mg/dL) 32 (27-41; 3.9) 52 (33-77; 14.2) Plasma urea (g/L) 0.25 (0.18-0.32; 0.04) 0.23 (0.16-0.30; 0.04) Plasma insulin (μg/L) 0.48 (0.28-1.04; 0.19) 0.18 (0.08-0.44; 0.09) 1

Table 2 .

 2 The relationship between the slope of the response during nutritional challenge, stage of lactation (early vs. late), and the prechallenge level of the measure 1

	Response slopes (V2,	Early	Late lactation	Coefficient for
	units per day)	lactation	adjustment	prechallenge level (V1)	Residual	R2
	DMI (kg/d)	-1.31	+0.38	-0.49	0.005	0.93
	Milk yield (kg/d)	-1.25	+0.55	-0.44	0.007	0.95
	Milk fat content (g/kg)	21.6	-3.5	+0.05	6.43	0.46
	Milk protein content (g/kg)	2.34	+2.60	+0.11	3.54	0.49
	BW (kg)	-3.41	+1.85	-0.02	0.19	0.87
	Plasma glucose (g/L)	-0.027	-0.006	-0.150	0.002	0.00
	Plasma fatty acids (NEFA; μmol/L)	845	-410	-0.3	26,600	0.61
	Plasma BHB (mg/dL)	8.34	-8.95	+0.03	57.89	0.38
	Plasma urea (g/L)	0.164	-0.041	-0.283	0.0005	0.66
	Plasma insulin (μg/L)	-0.042	-0.081	-0.111	0.003	0.37
	1

Table 3 .

 3 The relationship between early and late lactation values of the response slope during nutritional challenge derived by linear regression: response slope in early lactation = a + b × (response slope in late lactation) 1

	Item	Coefficient a	Coefficient b	Residual	R 2 adj
	DMI (kg/d)	-0.51*	0.85	0.16	0.40
	Milk yield (kg/d)	-0.73	0.74*	0.29	0.17
	Milk fat content (g/kg)	13.47	0.45*	2.99	0.19
	Milk protein content (g/kg)	0.81	0.31	1.65	0.23
	BW (kg)	-2.34	0.68	0.61	0.38
	Plasma glucose (g/L)	-0.027*	0.017	0.044	0.00
	Plasma fatty acids (NEFA; μmol/L)	1,075	-0.53	211	0.00
	Plasma BHB (mg/dL)	9.92	2.58	9.78	0.37
	Plasma urea (g/L)	0.089	0.608	0.032	0.42
	Plasma insulin (μg/L)	-0.042	0.0011	0.051	0.00

1

Perfect repeatability would result in b = 1. The response slopes (V2) were estimated using a piecewise mixed model (see Materials and Methods). Coefficients significantly different from 0 are shown in bold (P < 0.05). *0.05 < P value < 0.07.

Table 4 .

 4 The relationship between the rate of deceleration in recovery following nutritional challenge, stage of lactation (early vs. late), and the size of the response during the challenge 1The rate of deceleration in recovery (the quadratic term of the recovery; V4) and the size of the response (V2) were estimated using a piecewise mixed model (see Materials and Methods). Results are shown in equation form: V4 = average early lactation value + late lactation adjustment + coefficient × V2.

	Rate of deceleration in recovery (V4,	Early	Late lactation	Coefficient for effect of
	units per day 2 ):	lactation	adjustment	response slope (V2)	Residual	R2
	DMI (kg/d)	-0.241	0.049	0.133	0.0016	0.41
	Milk yield (kg/d)	-0.165	0.102	0.140	0.0020	0.67
	Milk fat content (g/kg)	5.87	-2.95	0.08	1.32	0.62
	Milk protein content (g/kg)	1.02	-0.30	0.07	0.457	0.02
	BW (kg)	-0.804	0.280	0.166	0.033	0.58
	Plasma glucose (g/L)	-0.026	0.007	0.259	0.00009	0.47
	Plasma fatty acids (NEFA; μmol/L)	287	-147	0.27	289.4	0.96
	Plasma BHB (mg/dL)	5.45	-3.25	0.140	1.33	0.76
	Plasma urea (g/L)	0.0756	-0.0139	0.2193	0.00006	0.57
	Plasma insulin (μg/L)	-0.029	-0.044	0.079	0.00062	0.58
	1

ACKNOWLEDGMENTS

We gratefully acknowledge the work of Joseph Tessier, Alexandra Eymard, and the team at the INRA UMR 791 Modélisation Systémique Appliquée aux Ruminants (Paris, France) experimental installation for the care of the animals and their work on the experimental protocol. This project benefited from an INRA PHASE "Crédit Incitatif" grant (INRA, Paris, France).

APPENDIX Description of the Mixture Model Used for the Analysis of Response Recovery Profiles

Let Yj = (Yj1,…,YjT) be the response profile of goat j. The (unknown) membership of goat j to one of the K classes is characterized through a label Zj that takes values in {1,…,K}. The response profile of a goat j that belongs to cluster k (i.e., for which Zj = k), is Yjt =V1k ×I(t≤0) +V2k ×t×I(0<t≤2) +(V3k ×t +V4k ×t×2)I(2<t≤4) +V5k ×I(t>4) +Ejt with the same notations as in the article. Parameters V1 k ,…,V5 k are similar to all goats belonging to class k but different between clusters. Such a model is called a mixture of linear models. The marginal distribution of profile Yj is

where πk is the proportion of goats that belong to cluster k, and where ϕ stands for the (Gaussian) likelihood of profile Yj according to class k.

The mixture model can be used to allocate each goat to a cluster. To do this, the probability of a goat to belong to cluster k given its profile is computed. This probability is called the posterior probability and is defined from equation [A1] by

The mixture parameters (proportions πk, parameters V1 k ,…,V5 k , and the residual variance) are estimated using the EM algorithm. The EM algorithm is dedicated to the class of incomplete data models where the status of the observations is unknown. In this context, the missing data are the information on the cluster for each goat. The algorithm iterates 2 steps:

• In the E step, the parameters are fixed and the posterior probability for each observation to belong to each class is calculated using equation [A2]. • In the M step, the posterior probabilities are fixed the parameters of each class are estimated using a weighted regression, in which the weights are given by the posterior probabilities.

Once convergence is achieved, one can classify an animal into a cluster using the maximum a posteriori rule: one just classifies the animal according to its highest posterior probability.

The whole strategy requires the number of classes to be known, something unrealistic in practice. When K is unknown, the model can be fit for different potential values of K, and the best value is then selected using a penalized criterion. In this article, we considered the classical integrated completed likelihood [START_REF] Biernacki | Assessing a mixture model for clustering with the integrated completed likelihood[END_REF]. This criterion is a penalized likelihood criterion, which tends to select a parsimonious number of relatively homogeneous clusters.