Plasma concentrations of branched-chain amino acids differ with Holstein genetic strain in pasture-based dairy systems

In pasture-based systems, there are nutritional and climatic challenges exacerbated across lactation; thus, dairy cows require an enhanced adaptive capacity compared with cows in confined systems. We aimed to evaluate the effect of lactation stage (21 vs. 180 days in milk, DIM) and Holstein genetic strain (North American Holstein, NAH, n = 8; New Zealand Holstein, NZH, n = 8) on metabolic adaptations of grazing dairy cows through plasma metabolomic profiling and its association with classical metabolites. Although 67 metabolites were affected (FDR < 0.05) by DIM, no metabolite was observed to differ between genetic strains while only alanine was affected (FDR = 0.02) by the interaction between genetic strain and DIM. However, complementary tools for time-series analysis (ASCA analysis, MEBA ranking) indicated that alanine and the branched-chain amino acids (BCAA) differed between genetic strains in a lactation-stage dependent manner. Indeed, NZH cows had lower (P-Tukey < 0.05) plasma concentrations of leucine, isoleucine and valine than NAH cows at 21 DIM, probably signaling for greater insulin sensitivity. Metabolic pathway analysis also revealed that, independently of genetic strains, AA metabolism might be structurally involved in homeorhetic changes as 40% (19/46) of metabolic pathways differentially expressed (FDR < 0.05) between 21 and 180 DIM belonged to AA metabolism.

Pasture-based dairy systems are gaining scientific interest due to the increasing demand for dairy products associated with lower environmental impact and better animal welfare 1,2 . However, these systems impose restrictive conditions which determine that high genetic merit cows are not able to express their genetic potential for milk yield as their nutritional requirements are not fulfilled due to limited dry matter intake 3 . Furthermore, pasturebased dairy systems are associated with greater nutritional and metabolic challenges for dairy cows because of climatic conditions and variability of pasture allowance and quality along the year 4 . Consequently, adaptive responses of dairy cows are of main importance when they are managed in these systems 5 .
In the dairy cow, at the beginning of lactation, dry matter intake is not sufficient to sustain the high rate of milk synthesis, thus a negative energy balance is established, and adipose and skeletal muscle are mobilized 6,7 . In consequence, it has been widely demonstrated that early lactating cows are characterized by a catabolic state of peripheral tissues associated with the uncoupling of somatotropic axis and a state of insulin resistance 8,9 . In this regard, grazing dairy cows are confronted with more challenging environments than confined dairy cows which determine that these homeorhetic changes are even more dramatic and last for longer during early lactation 10 . However, metabolic adaptations-particularly in early lactation-are affected by the dairy genotype associated with different priorities in nutrients partitioning between milk synthesis vs. body reserves 11,12 . When compared at grazing, NAH (selected for individual milk yields in confined systems) vs. NZH (selected for milk solids and reproductive efficiency in grazing systems) cows had a stronger uncoupling of somatotropic axis and an increased insulin resistance 8,9,13 which resulted in greater milk yield and a more pronounced loss of body reserves. In addition, differences in homeorhetic strategies between these Holstein strains would also include AA and protein metabolic pathways. Preliminary results obtained by our research group showed greater milk Metabolomic profiling and cluster analysis. The final metabolomic dataset consisted of 172, out of 200 metabolites assayed by the GCToF/MS method, which were effectively identified above the limit of detection (i.e. metabolite's identifier ion with peak intensity greater than the limit of detection in at least 80 % of the plasma samples) (Supplementary Table S2 online). Days in milk had a stronger effect than the genetic strain on principal component analysis (PCA) clustering (Fig. 1a). In addition, PCA showed that it was easier to discriminate early vs. mid-lactating NZH than NAH cows. Principal components 1 and 2 accounted for 33.7 % of the variation. The partial least square discriminant analysis (PLS-DA) confirmed that the animals could be better discriminated by stage of lactation for NZH than NAH cows (Fig. 1b). The best classification PLS-DA model, based on the first 3 components, had an acceptable predictive ability Q 2 ( R 2 = 0.94, Q 2 = 0.54). Fifty-nine, out of the 172 metabolites, had a value in importance projection (VIP) > 1.0. Among them, seven metabolites including oxoproline, 5-hydroxynorvaline, isoleucine, p-tolyl glucuronide, leucine, α-aminoadipic acid and L-pipecolate had a VIP > 2.0. For all of these metabolites , the lowest concentration was observed for NZH at 21 DIM (Fig. 1c).
Univariate analysis and correlations. According to ANOVA analysis, 67 of the quantified metabolites varied with DIM (FDR < 0.05; Table 2) whereas there was no metabolite affected by genetic strain after P-value Table 1. Least square means, standard error of the mean (SEM) and ANOVA fixed effects of productive variables and metabolite and endocrine concentrations for North American (NAH, n = 8) and New Zealand (NZH, n = 8) Holstein genetic strains (GS) at 21 and 180 days in milk (DIM) on grazing systems. 1: Parameters abbreviations: BW: body weight; BCS: body condition score; NEFA: non-esterified fatty acids, BHB: β -hydroxybutyrate. FPCM: fat and protein corrected milk, presented as average milk yield of +/− 2.5 days the date in which the blood plasma samples were collected. FPCM was estimated according to Østergaard et al. 63 as follows: FPCM=Milk yield[(0.0383×fat% + 0.0242×protein% + 0.7832)/3.14].  Table 2). Despite none of the quantified metabolites having an FDR ≤ 0.05 when comparing genetic strains, 16 metabolites had a raw-P ≤ 0.05, most of them being AA or AA-related compounds (n = 12) and fatty acids (n = 2) with greater concentrations for NAH than NZH cows (Table 2). Additionally, only alanine was affected (FDR = 0.02) by the interaction between genetic strain and DIM as it was greater at 21 DIM but lower at 180 DIM for NAH vs. NZH dairy cows ( Table 2, Fig. 2c). However, changing patterns of metabolites between genetic strains in a time-dependent manner were also assessed by both, ANOVA simultaneous component analysis (ASCA) and the Bayes statistical time-series analysis (MEBA) ranking. According to the ASCA analysis, a valid model was obtained for the interaction effect (genetic strain × DIM, P = 0.04, 4/100 permutations) (Supplementary Fig. 1 online). Four metabolites were well modelled (leverage > 0.040, SPE < 1.5 × 10 −30 ) by this interaction model (Fig. 2b,c). Among these 4 metabolites, dehydroabietic acid was decreased and 3-aminoisobutyric acid was increased for NAH vs. NZH cows only at 21 DIM (Fig. 2c). In addition, serine increased from 21 to 180 DIM only for NZH cows, while uracil was lower for NAH than NZH cows at 180 DIM but no significant differences were detected after Tukey's test for multiple comparisons (Fig. 2c). In addition, the top-10 metabolites (MEBA ranking top-10, Hotelling's T 2 > 10.0 ) with the greatest differences in temporal profiles between genetic strains were alanine, BCAA (valine, isoleucine, leucine), arachidic acid, oxoproline, squalene, p-tolyl glucuronide, erythritol and serine (Supplementary Table S3 online). Specifically, concentrations of BCAA were greater for NAH than NZH cows only at 21 DIM (Fig. 2c). In addition, p-tolyl glucuronide increased from 21 to 180 DIM for NAH cows reaching greater concentrations than their NZH counterparts (Fig. 2c). Except for uracil (raw-P = 0.054), arachidic acid (raw-P = 0.08), erythritol , and oxoproline (P > 0.10), all of these metabolites identified by both ASCA and/or MEBA ranking had a raw-P ≤ 0.05 for the interaction effect (genetic strain × DIM) according to ANOVA results.
Other metabolic pathways being affected (FDR < 0.05) by the lactation stage included 9 metabolic pathways that participate in lipid metabolism (Table 3; Fig. 3). These metabolic pathways were mostly related with greater concentrations of several lipids and metabolic intermediates for 21 vs. 180 DIM. These metabolites included myristic, palmitic, arachidonic acids , and glycerol, among others ( Fig. 4; Supl. Figure 3D). Finally, metabolic pathways belonging to nitrogen metabolism as well as redox metabolism, and vitamins and coenzymes metabolism were also observed to be affected by DIM.
Effect of genetic strain on metabolic pathways. Metabolic pathways differing (FDR < 0.05) between genetic strains were only detected at 21 DIM. The BCAA biosynthesis and degradation, and lysine degradation were shift-regulated (FDR = 0.022) between genetic strains as they were associated with greater concentrations of several metabolites for NAH vs. NZH cows (Table 3, Figs. 3 and 5). Differences in BCAA pathways were associated with greater concentrations (FDR < 0.05) of the three BCAA (valine, leucine, isoleucine) for NAH vs. NZH cows, as well as intermediary metabolites (e.g.: 2-ketoisocaproic, 2-ketobutyric) that were greater but did not significantly differ (Fig. 6). Lysine degradation was also linked with greater concentrations of L-pipecolic acid and α-aminoadipic acid for NAH vs. NZH cows, but no significant differences were detected for lysine. In addition, seleno compounds metabolism and aminoacyl-tRNA biosynthesis tended to be affected (FDR < 0.07) when comparing genetic strains as most of its involved metabolites (16/19) were greater for NAH vs. NZH cows in early lactation (Figs. 3, 5).

Discussion
Modulating effects of the genetic background on energy and lipid homeorhetic adaptations of dairy cows facing challenging conditions imposed by pasture-based dairy systems have been widely reported 10 . However, the role of genetic strain on AA metabolism in relationship with energy and lipid metabolism has not yet been reported. Our results demonstrated that homeorhetic changes of AA metabolism in early lactation are likely affected by the genetic strain associated with differences in energy metabolism. The greater milk yield and lower body reserves of NAH vs. NZH cows associated with a greater degree of insulin resistance and uncoupling of the somatotropic axis in NAH cows 8,9 , could be also linked with chronic activation of the mTOR pathway through increased BCAA degradation 20 in the NAH cows. These results are further discussed around two main aspects: a) differential effects of the Holstein genetic strain (NAH vs. NZH) on AA and energy metabolism during early lactation, b) effects of early vs. mid-lactation (21 vs. 180 DIM) on AA, energy, and lipid metabolic changes in grazing cows from a systemic point of view.
Our results indicated that differential effects of the Holstein genetic strain (NAH vs. NZH) on energy and AA metabolism, were evident only in early lactation. This suggests that the effect of the genetic strain on metabolic    www.nature.com/scientificreports/ adaptive responses is exacerbated in highly challenging conditions such as the onset of lactation. Differences were mostly related to plasma concentrations of AA or AA-related compounds as most of these metabolites were increased for NAH vs NZH cows in early lactation. Indeed, increased plasma BCAA could be related to the greater insulin resistance previously demonstrated for NAH vs NZH cows under grazing conditions 9 (Fig. 6). Circulating BCAA, which are poorly catabolized in the liver, can act as signaling molecules sensing the nutritional state and activating cellular signaling cascades 21 . However, a relationship between insulin resistance and BCAA has been recently suggested in ruminants 22 , and the exact mechanisms still remain to be known. In this sense, it is likely that greater BCAA would have upregulated the mTOR activity 23 in NAH cows during early lactation (Fig. 6). Despite upstream activation pathways of the mTOR differ between insulin and BCAA, insulin fails to stimulate the mTOR downstream cascade in the absence of AA 24 . Moreover, in humans, it has been proposed that high BCAA concentrations would indirectly impair insulin sensitivity through the chronic activation of the mTOR pathway leading to over-phosphorylation of the insulin receptor substrate 1 20,25 which ultimately determines a decreased expression and activity of glucose transporter 4 26 . However, the linking mechanism between insulin resistance and BCAA is not completely understood and whether high concentrations of BCAA are a cause or a consequence of insulin resistance development still remains to be known. Additionally, greater BCAA concentrations for NAH than NZH cows could also be the consequence of increased muscle protein mobilization in early lactation. Hypothetically increased activity of mTOR for NAH cows would enhance protein synthesis 27 . Therefore, it is possible that NAH cows had not only greater muscle mass breakdown but also protein synthesis, leading to increased protein turnover when compared with NZH cows. Indeed, Ghaffari et al. 18 recently suggested that lower insulin sensitivity in obese dairy cows would be associated with increased protein turnover when compared to lean cows.
Plasma concentrations of alanine and α-ketoisocaproic were greater for NAH than NZH cows suggesting increased catabolism of the BCAA in NAH animals 28 (Holecek et al., 2018). It has been observed that low intracellular glucose concentrations, due to reduced glucose uptake by peripheral tissues associated with reduced insulin sensitivity 29 , would enhance BCAA degradation through the glucose-Krüppel-like factor 15 (KLF15)-BCAA axis 30 . These latter authors demonstrated that low glucose levels down-regulate the negative feedback of the KLF15 transcriptional regulator factor at the BCAA decarboxylation step which in turn is the limiting   www.nature.com/scientificreports/ rate step of BCAA catabolism. In addition, alanine as well as the degradation products of BCAA replenish the oxalacetate pool or enter into the TCA cycle in order to produce glucose or energy, reflecting the interconnection between protein and energy metabolism 26,31 . Actually, a greater oxalacetate redirection from the TCA cycle towards gluconeogenesis would determine a down-regulation state of this cycle (cataplerosis) leading to the reduced citrate, aconitate and α-ketoglutarate concentrations observed for NAH vs. NZH cows (Fig. 6) during early lactation 32,33 . In addition, the expected greater DMI of NAH cows (due to its greater BW, + 8%), in addition to increased valine catabolism leading to greater propionyl-CoA, could explain the increased plasma concentrations of succinate and fumarate observed for NAH vs. NZH cows 34 . Taken together our results point at a greater BCAA degradation rate and decreased TCA cycle activity for NAH vs. NZH cows.
Interestingly, other metabolites such as α-aminoadipic, which were increased for NAH than NZH cows, are probably also related with greater insulin resistance. Indeed, α-aminoadipic acid, a product of lysine degradation, has been previously reported as a pre-diabetic biomarker in rodents and humans 35 as well as in dairy cows 18 . In addition, both, α-aminoadipic acid and L-pipecolate, another metabolite of the lysine degradation, were reported to be elevated in blood samples of patients with liver injury and peroxisomal disorders in humans 36 . The greater concentrations of these metabolites observed for NAH vs. NZH cows would be in agreement with the reduced liver mitochondrial functionality observed by our team for NAH cows 37 .
Independently of their genetic strain, all cows showed several metabolic changes across lactation stages which confirmed several links widely reported between AA, and energy and lipid metabolism during negative energy balance in early lactation. The general trend for decreased plasma concentrations of AA in early lactation was in agreement with previous reports 38 reflecting a state of negative protein balance at 21 DIM 39,40 because of increased AA requirements for milk protein synthesis and gluconeogenesis together with a decreased dry matter intake after calving 39,41,42 . In contrast to the general trend for decreased plasma AA, the tendency of greater plasma glycine and greater concentrations of trans-4-hydroxyproline in early than mid-lactation suggest an increased muscle protein breakdown coupled with increased AA de novo synthesis at 21 vs. 180 DIM 38,43 . Indeed, the greater concentrations of trans-4-hydroxyproline in early lactation could reflect increased muscle protein breakdown at this time as its plasma concentrations would probably be indicative of connective tissue amounts of musculoskeletal system 44,45 .
In addition to decreased AA plasma concentrations, our results indicate that the urea cycle was down-regulated at early lactation probably linked with the increased demand for energy production 17 in agreement with studies that reported a decreased activity of urea cycle enzymes around parturition 46,47 . Indeed, all measured compounds related to the urea cycle (e.g.: urea, ornithine, citrulline, aspartate), except fumarate, were lower for early than mid-lactation. In this sense, our results reinforce the hypothesis formulated by Kuhla et al. 43 which states that the enhanced muscle breakdown during early lactation provides substrates for milk synthesis in parallel with a slowed-down activity of the urea cycle. Consequently, an enhanced nitrogen use efficiency would be expected to occur during early lactation, which is consistent with the lower plasma and milk concentrations of urea observed in our study 48,49 . In agreement with recent ly reported data by Luo et al. 17 , greater fumarate together with decreased glutamate, which can act as a shuttle between the urea and the TCA cycles, suggest that the slowed-down activity of the urea cycle could happen to sustain increased TCA activity through the synthesis of fumarate. Indeed, the metabolic pathway analysis indicated that the TCA cycle was enhanced in early vs. www.nature.com/scientificreports/ mid-lactation as previously reported for confined dairy cows 50 . The increased plasma concentrations of citrate and isocitrate agreed with increased activities of citrate synthase and isocitrate dehydrogenase previously reported during negative energy balance of early lactation 51,52 . Moreover, greater TCA activity during early lactation is concordant with increased NEFA and saturated free fatty acids (e.g.: palmitic and stearic acids) concentrations due to adipose tissue mobilization after calving as they are main components of adipose tissue 53,54 . Additionally, greater concentrations of unsaturated fatty acids such as linoleic, linolenic, and arachidonic acids suggest that early lactation was associated with oxidative stress and a proinflammatory state 55 . In fact, polyunsaturated fatty acids (PUFA), which are well-known precursors for oxylipid synthesis 56 through its metabolization by the cyclooxygenase/lipoxygenase pathway, provide a link between lipid mobilization, and oxidative stress and inflammation during the transition period 55,57 .
Lastly, decreased concentrations of ketogenic AA such as lysine, tyrosine, and phenylalanine in early lactation might have accounted for increased ketone body synthesis 58 , which is further depicted by increased plasma concentrations of BHB. Indeed, metabolic pathway analysis revealed that the butanoate metabolism was differentially expressed between stages of lactation. The decreased concentrations of α-ketoglutarate in early lactation could be due, at least in part, to increased activity of 2-hydroxyglutarate dehydrogenase leading to the observed greater concentrations of 2-hydroxyglutarate 59 , which is further used in the synthesis of acetoacetyl-CoA. In turn, acetoacetyl-CoA is metabolized via acetoacetyl-CoA and β-hydroxy-β -methylglutaryl-CoA (HMG-CoA) in the liver mitochondria 59 .
Regarding the overall effects of the lactation stage on energy and lipid metabolism, our results confirm data widely reported in both confined and pasture-based dairy systems. However, pasture-based dairy systems impose great variation in the quality and availability of pastures which determines often changes in AA intake along the year 60 ; therefore, further research is needed to better know the homeorhetic interactions between AA metabolism, and energy and lipid metabolism when grazing dairy cows are subjected to changing nutritional environments.

Conclusions
In agreement with reported data for confined systems, our results confirmed decreased plasma AA and increased fatty acid concentrations concomitant with a down-regulation of the urea cycle, and an up-regulation of the TCA cycle in early vs. mid-lactating grazing cows. However, more studies are needed for a better understanding of homeorhetic interactions between specific AA and energy metabolism in cows managed in pasture-based dairy systems. In addition, BCAA possibly act as signaling molecules behind differences in adaptive metabolic responses of NAH vs. NZH genetic strains when managed under grazing conditions.

Methods
Experimental design and management. The experiment was carried out at the Experimental Research Station "La Estanzuela" (34°20' S, 57°40' W) belonging to the Instituto Nacional de Investigación Agropecuaria (INIA) of Uruguay. This study was carried out in compliance with the ARRIVE guidelines and all procedures were approved by the Ethic Committee on Animal Experimentation of INIA (form INIA_2017.2).
Cows were randomly selected from a larger experiment designed to study the interaction between genetic strains and grazing-based feeding strategies previously reported by Talmon et al. 61 . Multiparous dairy cows ( 3.1 ± 0.9 lactations; fall calving 5/15/2018 ± 12 days, mean ± SD) of NAH (n = 8) and NZH (n = 8) genetic strains were used. Cows were fed a strategy that maximized pasture grazing. Previous to calving, cows had a LW of 593 ± 17 and 560 ± 17, and a BCS of 2.94 ± 0.06 and 3.28 ± 0.06 for NAH and NZH cows, respectively (mean ± SD). At least 87.5% of each cow's ancestors (three generations) had an American (USA or Canada) or New Zealand proved origin for NAH and NZH genetics strains, respectively (Mejoramiento y Control Lechero Uruguayo; https:// www. genet icale chera. com. uy/). The 305-days expected milk yield was 7500 and 5500 kg and the economic and productive selection index was 104 and 130 on average for NAH and NZH cows, respectively. The NAH cows had an expected progeny difference of +70.4 kg, +0.01 % fat and -0.01 % protein for milk yield, milk fat and milk protein content, respectively, when compared to the national herd. The NZH cows had an expected progeny difference of -185.7 kg, +0.23 % fat and +6.19 % protein for milk yield, milk fat and milk protein content, respectively, when compared to the national herd.
The feeding strategy was designed to maximize pasture grazing according to weekly pasture growing rate, and concentrate was offered twice a day individually (at a rate of 33% of the predicted daily dry matter intake) at the milking parlor. Cows were outdoors all year-round and throughout the experiment cows grazed orchard grass + lucerne (Dactylis glomerata + Medicago sativa, 75% of the grazing time) or fescue (Festuca arundinacea, 25% of the grazing time) on a rotational grazing system with strips assigned after milking and free access to freshwater. In order to avoid ingestion behavior interferences or dominance between genetic strains due to body size differences, each group was offered its own daily strip of pasture to keep a similar herbage allowance relative to their BW. Greater details on pasture and grazing management are available on Supplementary Table S1 online as well as in the work published by Talmon et al. 61 . Conserved forage was offered in a feeding parlor immediately before the afternoon milking. On average, the diet was comprised of 64% of grazed pastures, 5% of concentrate, and 31% of conserved forage (corn silage and pasture haylage mix; 73:27 ± 8% on dry matter basis, respectively). Further details of diet offered are available in Supplementary Table S1 online.
Animal measurements and sample collection. Cows were milked twice a day at 0400 and 1400 h.
Milk yield was recorded daily and weekly milk samples were collected for composition analysis (Combi FOSS FT+, Foss Electric, Hillerhød, Denmark). Cow body weight and BCS (scale from 1 to 5, Edmonson et al. 62 ) were recorded once every two weeks. Fat and protein corrected milk (FPCM) yield was calculated according to Østergaard et al. 63 Blood samples were collected at early and mid-lactation ( www.nature.com/scientificreports/ (mean ± SE) by caudal venipuncture using 10 mL heparinized Vacutest ®tubes (Vacutest Kima, Arzergrande, Italia). Plasma samples were harvested by centrifugation (4000 × g, 12 min) and immediately stored at -80°C until analysis. Plasma glucose, NEFA, β-hydroxybutyrate (BHB), urea, and insulin concentrations were determined as previously reported 64 .
Targeted metabolomic analysis and annotation. Metabolomic analysis and ion annotation procedures were performed by using gas chromatography time-of-flight mass spectrometry (GCToF/MS) according to Fiehn et al. 65 at the core lab West Metabolomic Center (UC Davis Genome Center, Davis, CA, USA) and 200 metabolites were analyzed. A column of 30 m length × 0.25 mm internal diameter with 0.25 µ m film made of 95% dimethyl/5%diphenylpolysiloxane (Rtx-5Sil MS, Restek®Corporation,) was used for chromatography analyses. A Leco®Pegasus IV mass spectrometer (St Joseph, MI, USA) was used with unit mass resolution at 17 spectra/s from 80 to 500 Da and ionization energy set in -70 eV and equipped with an 1800 V detector voltage, 230 °C transfer line, and a 250°C ion source. Chromatographic data were pre-processed without smoothing and peak width was set at 3 s, the baseline subtraction was done just above the noise level, and automatic mass spectral deconvolution and peak detection at signal/noise levels of 5:1 were automatically performed throughout the chromatogram. Annotation of derivatized ions was performed using the BinBase algorithm based on deconvoluted spectra and peak metadata (retention index, unique ion, spectral similarity, signal/noise ratio, peak purity) from the LecoChromaTOF software using a multi-tiered filtering system with stringent thresholds, specifically developed for GCToF/MS data annotation 66 . Data is reported as peak height normalized by the sum of all peak heights (mTIC) of each sample as previously reported by Fiehn et al. 65 .
Statistical analysis. Productive performance and metabolite and endocrine concentrations were analyzed as repeated measures using the MIXED procedure (SAS University Edition, SAS Institute Inc., Cary, NC, USA). The model included genetic strain, lactation stage (21 vs. 180 DIM) and their interaction as fixed effect, and cow as a random effect. Metabolomic data pre-processing and statistical analysis were performed using MetaboAnalyst v4 (https:// www. metab oanal yst. ca/). Prior to statistical analysis, data were normalized by the median, cube root transformed, and auto-scaled according to Chong et al. 67 . Data quality was assessed through multivariate analysis comparing individual samples' data against the pooled samples. Data were subjected to multivariate analysis: clustering was assessed by PCA and classification models were assessed by PLS-DA 68 .
Data were also subjected to ANOVA analysis considering a time-series model in which genetic strain (NAH vs. NZH), lactation stage (21 vs 180 DIM), and its interaction were considered as fixed effects, while the cow was considered as a random effect. Raw-P values were adjusted for multiple hypothesis testing 69 at a false discovery rate of 5 % (FDR ≤ 0.05). The time-trend of metabolites was complementarily assessed both by ANOVA simultaneous component analysis (ASCA) 70 and Bayes statistical time-series analysis (MEBA) 71 . The ASCA analysis was performed on the basis of a model for each fixed effect, this is genetic strain, lactation stage, and interaction models. Each model was validated according to a permutation test being the significance threshold P < 0.05. If a model was significant, then the metabolites were evaluated according to leverage and SPE. Therefore, a given metabolite was considered to be affected by the genetic strain in a time-dependent manner if it was well-modeled by the ASCA interaction model, this is when the metabolite had high leverage (> 0.04) and low SPE (< 1.5 × 10 −30 ) values. Additionally, MEBA was performed based on Hotteling's T 2 , which is a Student's t-statistic generalization for multivariate analysis. It was used to rank the metabolites according to their difference in temporal profile when comparing genetic strain 71 . The higher position in the ranking, the higher differences in temporal profiles among genetic strains.
Finally, metabolic pathways analysis was based on the Bos taurus KEGG database, combining the Globaltest, which calculates the association between the metabolite sets and the phenotype without referring to a background 72 , and a topologic analysis based in the betweenness centrality, which is a measure of the importance of a compound within a given metabolic pathway 67 . Significant enrichment of metabolic pathway was set at FDR ≤ 0.05, and only metabolic pathways with at least two metabolites quantified in the current data set were further considered for discussion purposes. The metabolic pathway analysis was performed based on three comparisons: www.nature.com/scientificreports/ conducted the field experimental platform. All authors were involved in reviewing the original draft and approved the manuscript.