Advertisement

Climate sensitivity of milk production traits and milk fatty acids in genotyped Holstein dairy cows

Open ArchivePublished:March 10, 2021DOI:https://doi.org/10.3168/jds.2020-19411

      ABSTRACT

      The aim of this study was the evaluation of climate sensitivity via genomic reaction norm models [i.e., to infer cow milk production and milk fatty acid (FA) responses on temperature-humidity index (THI) alterations]. Test-day milk traits were recorded between 2010 and 2016 from 5,257 first-lactation genotyped Holstein dairy cows. The cows were kept in 16 large-scale cooperator herds, being daughters of 344 genotyped sires. The longitudinal data consisted of 47,789 test-day records for the production traits milk yield (MY), fat yield (FY), and protein yield (PY), and of 20,742 test-day records for 6 FA including C16:0, C18:0, saturated fatty acids (SFA), unsaturated fatty acids (UFA), monounsaturated fatty acids (MUFA), and polyunsaturated fatty acids (PUFA). After quality control of the genotypic data, 41,057 SNP markers remained for genomic analyses. Meteorological data from the weather station in closest herd distance were used for the calculation of maximum hourly daily THI. Genomic reaction norm models were applied to estimate genetic parameters in a single-step approach for production traits and FA in dependency of THI at different lactation stages, and to evaluate the model stability. In a first evaluation strategy (New_sire), all phenotypic records from daughters of genotyped sires born after 2010 were masked, to mimic a validation population. In the second strategy (New_env), only daughter records of the new sires recorded in the most extreme THI classes were masked, aiming at predicting sire genomic estimated breeding values (GEBV) under heat stress conditions. Model stability was the correlation between GEBV of the new sires in the reduced data set with respective GEBV estimated from all phenotypic data. Among all test-day production traits, PY responded as the most sensitive to heat stress. As observed for the remaining production traits, genetic variances were quite stable across THI, but genetic correlations between PY from temperate climates with PY from extreme THI classes dropped to 0.68. Genetic variances in dependency of THI were very similar for C16:0 and SFA, indicating marginal climatic sensitivity. In the early lactation stage, genetic variances for C18:0, MUFA, PUFA, and UFA were significantly larger in the extreme THI classes compared with the estimates under thermoneutral conditions. For C18:0 and MUFA, PUFA, and UFA in the middle THI classes, genetic correlations in same traits from the early and the later lactation stages were lower than 0.50, indicating strong days in milk influence. Interestingly, within lactation stages, genetic correlations for C18:0 and UFA recorded at low and high THI were quite large, indicating similar genetic mechanisms under stress conditions. The model stability was improved when applying the New_env instead of New_sire strategy, especially for FA in the first stage of lactation. Results indicate moderately accurate genomic predictions for milk traits in extreme THI classes when considering phenotypic data from a broad range of remaining THI. Phenotypically, thermal stress conditions contributed to an increase of UFA, suggesting value as a heat stress biomarker. Furthermore, the quite large genetic variances for UFA at high THI suggest the consideration of UFA in selection strategies for improved heat stress resistance.

      Key words

      INTRODUCTION

      Phenotypic expressions of a genotype can differ with changing environmental gradients. Such a phenomenon is known as environmental sensitivity or reaction norm (
      • Falconer D.S.
      Selection in different environments—Effects on environmental sensitivity (reaction norm) and on mean performance.
      ). Environmental sensitivity has a genetic component, expressed in scaling or re-ranking effects (
      • Lynch M.
      • Walsh B.
      Genetics and Analysis of Quantitative Traits.
      ;
      • Ravagnolo O.
      • Misztal I.
      Genetic component of heat stress in dairy cattle, parameter estimation.
      ). The scaling effect is associated with changes in genetic variances for the trait of interest in dependency of environmental alterations. The animal re-ranking effect is due to genetic correlations lower than one between the same traits recorded in different environments (
      • Falconer D.S.
      • Mackay T.F.C.
      Introduction to Quantitative Genetics.
      ). Especially the re-ranking effect has, if ignored, a significant effect on the efficiency of cattle breeding programs and on the accuracy of selection (
      • Bohmanova J.
      • Misztal I.
      • Tsuruta S.
      • Norman H.
      • Lawlor T.
      Genotype by environment interaction due to heat stress.
      ;
      • Hammami H.
      • Rekik B.
      • Bastin C.
      • Soyeurt H.
      • Bormann J.
      • Stoll J.
      • Gengler N.
      Environmental sensitivity for milk yield in Luxembourg and Tunisian Holsteins by herd management level.
      ;
      • Cheruiyot E.K.
      • Nguyen T.T.T.
      • Haile-Mariam M.
      • Cocks B.G.
      • Abdelsayed M.
      • Pryce J.E.
      Genotype-by-environment (temperature-humidity) interaction of milk production traits in Australian Holstein cattle.
      ).
      There is increasing interest studying environmental sensitivity of dairy sires because phenotypic daughter records are distributed along a broad pattern of environments (e.g., in climatic conditions beyond the thermoneutral zone reflecting heat stress). Consequences of heat stress were a decline in feed intake and an increased respiration rate (
      • Umphrey J.
      • Moss B.
      • Wilcox C.
      • Van Horn H.
      Interrelationships in lactating Holsteins of rectal and skin temperatures, milk yield and composition, dry matter intake, body weight, and feed efficiency in summer in Alabama.
      ;
      • Herbut P.
      • Angrecka S.
      • Walczak J.
      Environmental parameters to assessing of heat stress in dairy cattle—A review.
      ), ultimately contributing to inferior milk production (
      • Brügemann K.
      • Gernand E.
      • Von Borstel U.
      • König S.
      Genetic analyses of protein yield in dairy cows applying random regression models with time-dependent and temperature x humidity-dependent covariates.
      ;
      • Hammami H.
      • Bormann J.
      • M'hamdi N.
      • Montaldo H.H.
      • Gengler N.
      Evaluation of heat stress effects on production traits and somatic cell score of Holsteins in a temperate environment.
      ), as well as impaired milk quality (
      • Hammami H.
      • Vandenplas J.
      • Vanrobays M.-L.
      • Rekik B.
      • Bastin C.
      • Gengler N.
      Genetic analysis of heat stress effects on yield traits, udder health, and fatty acids of Walloon Holstein cows.
      ) and female fertility (
      • Ravagnolo O.
      • Misztal I.
      Effect of heat stress on nonreturn rate in Holstein cows: Genetic analyses.
      ;
      • Haile-Mariam M.
      • Carrick M.
      • Goddard M.
      Genotype by environment interaction for fertility, survival, and milk production traits in Australian dairy cattle.
      ).
      • Ravagnolo O.
      • Misztal I.
      Genetic component of heat stress in dairy cattle, parameter estimation.
      and
      • Hill D.L.
      • Wall E.
      Dairy cattle in a temperate climate: The effects of weather on milk yield and composition depend on management.
      identified environmental sensitivity of test-day production traits (i.e., in response to changes in air temperature in combination with relative humidity). A comprehensive evaluation of various temperature-humidity indices (THI) for dairy cattle breeding purposes was presented by
      • Bohmanova J.
      • Misztal I.
      • Cole J.
      Temperature-humidity indices as indicators of milk production losses due to heat stress.
      . In genetic analyses, most of the introduced THI were accurate predictors for the decline in dairy cow productivity (
      • Bohmanova J.
      • Misztal I.
      • Cole J.
      Temperature-humidity indices as indicators of milk production losses due to heat stress.
      ;
      • Hammami H.
      • Bormann J.
      • M'hamdi N.
      • Montaldo H.H.
      • Gengler N.
      Evaluation of heat stress effects on production traits and somatic cell score of Holsteins in a temperate environment.
      ). From a genetic-statistical modeling perspective, reaction norm models (RNM) with a regression on a THI function (
      • Hammami H.
      • Vandenplas J.
      • Vanrobays M.-L.
      • Rekik B.
      • Bastin C.
      • Gengler N.
      Genetic analysis of heat stress effects on yield traits, udder health, and fatty acids of Walloon Holstein cows.
      ;
      • Bohlouli M.
      • Alijani S.
      • Naderi S.
      • Yin T.
      • König S.
      Prediction accuracies and genetic parameters for test-day traits from genomic and pedigree-based random regression models with or without heat stress interactions.
      ) have been applied, to estimate genetic components of heat stress. Such modeling approach allows more flexibility in statistical modeling of environmental sensitivity rather than defining a specific fixed THI threshold (
      • Brügemann K.
      • Gernand E.
      • Von Borstel U.
      • König S.
      Genetic analyses of protein yield in dairy cows applying random regression models with time-dependent and temperature x humidity-dependent covariates.
      ). The RNM provides genetic parameters and breeding values of animals across THI, enabling selective improvement of thermotolerance (
      • Ravagnolo O.
      • Misztal I.
      Genetic component of heat stress in dairy cattle, parameter estimation.
      ;
      • Cheruiyot E.K.
      • Nguyen T.T.T.
      • Haile-Mariam M.
      • Cocks B.G.
      • Abdelsayed M.
      • Pryce J.E.
      Genotype-by-environment (temperature-humidity) interaction of milk production traits in Australian Holstein cattle.
      ).
      • Bohlouli M.
      • Alijani S.
      • Naderi S.
      • Yin T.
      • König S.
      Prediction accuracies and genetic parameters for test-day traits from genomic and pedigree-based random regression models with or without heat stress interactions.
      and
      • Cheruiyot E.K.
      • Nguyen T.T.T.
      • Haile-Mariam M.
      • Cocks B.G.
      • Abdelsayed M.
      • Pryce J.E.
      Genotype-by-environment (temperature-humidity) interaction of milk production traits in Australian Holstein cattle.
      identified increased prediction accuracies for dairy cow performances in distinct THI classes when considering RNM with a component reflecting environmental sensitivity (e.g., an additional genotype by environment interaction term).
      Environmental sensitivity analyses in dairy cows focused on milk production and reproduction traits (
      • Ravagnolo O.
      • Misztal I.
      Effect of heat stress on nonreturn rate in Holstein cows: Genetic analyses.
      ;
      • Brügemann K.
      • Gernand E.
      • Von Borstel U.
      • König S.
      Genetic analyses of protein yield in dairy cows applying random regression models with time-dependent and temperature x humidity-dependent covariates.
      ;
      • Brügemann K.
      • Gernand E.
      • König von Borstel U.
      • König S.
      Defining and evaluating heat stress thresholds in different dairy cow production systems.
      ). Regarding fine milk composition traits, quite large genetic variances and heritabilities were estimated for fatty acids (FA) across lactations (
      • Bastin C.
      • Soyeurt H.
      • Gengler N.
      Genetic parameters of milk production traits and fatty acid contents in milk for Holstein cows in parity 1–3.
      ) and THI (
      • Hammami H.
      • Vandenplas J.
      • Vanrobays M.-L.
      • Rekik B.
      • Bastin C.
      • Gengler N.
      Genetic analysis of heat stress effects on yield traits, udder health, and fatty acids of Walloon Holstein cows.
      ). Both studies suggested FA, especially UFA, as predictors for cow energy balances and methane emissions. Recently, FA are available on a large scale through mid-infrared spectroscopy of milk samples (
      • Soyeurt H.
      • Dehareng F.
      • Gengler N.
      • McParland S.
      • Wall E.
      • Berry D.P.
      • Coffey M.
      • Dardenne P.
      Mid-infrared prediction of bovine milk fatty acids across multiple breeds, production systems, and countries.
      ). Phenotypically, oleic acid and mid-infrared predictors (UFA, MUFA, PUFA, and long-chain FA) responded highly sensitive to heat stress (
      • Hammami H.
      • Vandenplas J.
      • Vanrobays M.-L.
      • Rekik B.
      • Bastin C.
      • Gengler N.
      Genetic analysis of heat stress effects on yield traits, udder health, and fatty acids of Walloon Holstein cows.
      ) and were suggested as heat stress candidate biomarkers in milk (
      • Tian H.
      • Zheng N.
      • Wang W.
      • Cheng J.B.
      • Li S.L.
      • Zhang Y.D.
      • Wang J.Q.
      Integrated metabolomics study of the milk of heat-stressed lactating dairy cows.
      ).
      To our knowledge, no previous study focused on environmental sensitivity of FA, simultaneously considering high-throughput SNP marker from cow reference populations combined with dense meteorological data. Consequently, the objectives of this study were (1) to construct a genomic RNM suitable for large-scale genetic analyses of thermal stress, (2) to evaluate environmental sensitivity for milk production and milk FA considering a broad grid pattern for lactation stages combined with THI, and (3) to examine prediction accuracies of genomic estimated breeding values (GEBV) for sires (without daughter records in heat stress environments; i.e., the model stability of genomic RNM).

      MATERIALS AND METHODS

      Cow Traits

      The phenotypic data consisted of 47,789 test-day records for milk production traits and 20,742 test-day records for FA from 5,257 first-lactation genotyped Holstein dairy cows. The cows were kept in 16 large-scale cooperator herds, located in the region of former East Germany and comprising the federal states of Berlin-Brandenburg and Mecklenburg-West Pomerania. The 16 large-scale herds reflect the cow genotype basis as implemented for national genomic evaluations considering cow training sets (as described by
      • Yin T.
      • König S.
      Heritabilities and genetic correlations in the same traits across different strata of herds created according to continuous genomic, genetic, and phenotypic descriptors.
      ). The test-day production traits included milk yield (MY), fat yield (FY), and protein yield (PY). The FA (measured in g/100 g of MY) were C16:0, C18:0, MUFA, PUFA, SFA, and UFA. The FA were determined by mid-infrared spectra using pre-installed prediction equations from the Foss MilkoScan FT6000 spectrometer, installed at the milk recording center in Güstrow, Germany. Trait recording comprised the years 2010 to 2016, considering cows with an age at first calving between 20 to 39 mo. Test-day records within first lactation spanned the period from DIM 5 to 305. Cows had at least 4 test-day records. The herd-test-day size comprised at least 10 records.

      Genotypes and Pedigree

      A total of 1,618 animals (1,445 cows and 173 sires) were genotyped with the Illumina BovineSNP50 BeadChip (Illumina), and 4,612 animals (3,812 cows and 800 sires) were genotyped with the Illumina Bovine Eurogenomics 10K low-density SNP chip. The animals genotyped with the 10K SNP chip were imputed to 50K SNP panel as done for official national genetic evaluations by project partner VIT Verden (
      • Segelke D.
      • Chen J.
      • Liu Z.
      • Reinhardt F.
      • Thaller G.
      • Reents R.
      Reliability of genomic prediction for German Holsteins using imputed genotypes from low-density chips.
      ). After imputation, the genotype data consisted of 45,613 SNP from 6,230 genotyped animals (5,257 cows and 973 sires). In the quality control process, 866 SNP located on the X and Y chromosomes and 3,688 SNP with a minor allele frequency lower than 0.05 were discarded. Consideration of Hardy-Weinberg equilibrium excluded 2 SNP because the difference between observed and expected heterozygous frequencies was larger than 0.15 (according to
      • Wiggans G.
      • Sonstegard T.
      • VanRaden P.
      • Matukumalli L.
      • Schnabel R.
      • Taylor J.
      • Schenkel F.
      • Van Tassell C.
      Selection of single-nucleotide polymorphisms and quality of genotypes used in genomic evaluation of dairy cattle in the United States and Canada.
      ). Finally, 41,057 SNP remained in the genotype data set for all genomic analyses. Quality controls of the SNP data were performed using the preGSf90 program (
      • Aguilar I.
      • Misztal I.
      • Tsuruta S.
      • Legarra A.
      • Wang H.
      PREGSF90 – POSTGSF90: Computational tools for the implementation of single-step genomic selection and genome-wide association with ungenotyped individuals in BLUPF90 program. Proc. 10th WCGALP, Vancouver, Canada.
      ). Pedigrees of phenotyped cows could be traced back to at least 4 generations, and included 26,180 animals with genetic relationships to the cows with phenotypic records. The phenotyped and genotyped cows were daughters of 344 genotyped sires, implying at average 15 genotyped daughters per sire. Table 1 gives an overview for the genotype and phenotype data sets as used in the present study.
      Table 1Number of animals and number of records as considered in the present study
      Information sourceNumber
      Pedigree
       Cows26,180
       Sires2,457
       Dams18,786
      Genotyped
       Cows5,257
       Sires973
      With phenotyped daughters344
       Dams5,257
      With phenotyped daughters316
      Phenotyped and genotyped cows5,257
       Records for test-day production traits47,789
       Records for test-day fatty acids20,742

      Meteorological Data

      The weather station in nearest farm distance was identified considering the respective longitude and latitude, and using the GEOSPHERE software package in R (
      • Hijmans R.J.
      • Williams E.
      • Vennes C.
      • Hijmans M.R.J.
      Package 'geosphere’.
      ). The maximum distance between a herd and the weather station was 26 km, and the minimal distance was 6 km. Hourly THI were calculated considering hourly temperature (T) and relative humidity (RH) as follows (
      • NRC
      A Guide to Environmental Research on Animals.
      ):
      THI=(1.8×T+32)-(0.55-0.0055×RH)×(1.8×T-26).


      Afterward, we identified the maximum hourly daily THI. Maximum hourly daily THI were merged with the respective test-date for the milk traits. For ongoing association studies, maximum hourly THI from 4 d (from the test-date and the previous 3 d) were averaged (according to
      • Bohmanova J.
      • Misztal I.
      • Tsuruta S.
      • Norman H.
      • Lawlor T.
      Genotype by environment interaction due to heat stress.
      ). Small proportions of phenotypic records were available at both extreme ends of the THI scale. A balanced data distribution across THI is important because a low proportion of phenotypic records for specific THI might result in unexpected trajectories of genetic parameters (
      • Misztal I.
      • Strabel T.
      • Jamrozik J.
      • Mäntysaari E.A.
      • Meuwissen T.H.E.
      Strategies for estimating the parameters needed for different test-day models.
      ;
      • Cheruiyot E.K.
      • Nguyen T.T.T.
      • Haile-Mariam M.
      • Cocks B.G.
      • Abdelsayed M.
      • Pryce J.E.
      Genotype-by-environment (temperature-humidity) interaction of milk production traits in Australian Holstein cattle.
      ).
      Therefore, to avoid possible artifacts of RNM on (co)variance component estimates, averaged maximum THI were grouped into 10 THI classes (Figure 1). Cows had phenotypic records in at least 3 different THI classes; 34% and 37% of cows had records in 6 and 7 different THI classes, respectively. Nevertheless, the THI grouping strategy does not necessarily reflect physiological mechanisms, which may be different for the different traits under study.
      Figure thumbnail gr1
      Figure 1Number of records per temperature-humidity index (THI) class in early (5–105 DIM, stage 1), middle (106–205 DIM, stage 2), and late (206–305 DIM, stage 3) stages of lactation.

      Statistical Analyses

      At the phenotypic level, the effect of THI on milk production and milk FA traits was evaluated using the following statistical model 1:
      yijklmno=μ+HTDi+MFj+ASk+THIl(Sm)+Cn+eijklmno,
      [1]


      where yijklmno was the test-day record; µ was the overall mean; HTDi was the effect of ith herd-test-date; MFj was the jth milking frequency (2 or 3 times per day); ASk was the kth age-season of calving class; THIl(Sm) was the lth THI class nested within the mth stage of lactation (3 stages were defined: 5–105 DIM, 106–205 DIM, and 206–305 DIM); Cn was the random effect of the nth cow, and eijklmno was the random residual effect. From model 1, we computed and studied least squares means for milk production and milk FA traits in dependency of THI classes.
      Effects of THI on genetic parameters and breeding values were inferred via genomic RNM. In this regard, in a multiple-trait modeling approach, same traits from different lactation stages were considered as different traits. The RNM (model 2) was defined as follows:
      yijklmns=HTDis+MFjs+o=02αkoszo(t)+o=02βloszo(t)+o=02γmoszo(t)+eijklmns,
      [2]


      where yijklmns was the nth test-day record of the mth cow from the sth lactation stage; HTDis was the fixed effect of the ith herd-test-date for the sth lactation stage; MFjs was the jth milking frequency class; αkos was the oth fixed regression coefficient specific to the kth age-season of calving class; βlos was the oth random regression coefficient for the additive genetic effect of the lth animal at the sth stage of lactation; γmos was the oth random regression coefficient for the permanent environmental effect of the mth cow, zo(t) was a vector of covariates of size o of Legendre polynomials for fixed and random regressions evaluated at the tth THI class (10 classes); and eijklmns was the random residual effect. Fixed and random regressions were modeled using the second order of Legendre polynomials. To illustrate the 3-trait modeling approach for the 3 lactation stages, the model notation in matrix version was
      [y1y2y3]=[X1000X2000X3][b1b2b3]+[Z1000Z2000Z3][a1a2a3]+[W1000W2000W3][p1p2p3][e1e2e3],


      where ys was the vector of observations for the same trait at the sth stage of lactation; b was the vector for fixed effects as described above; a and p were the vectors of random regression coefficients for animal additive genetic and permanent environmental effects, respectively; Xs was the incidence matrix for fixed effects; Zs and Ws were the incidence matrices for additive genetic and permanent environmental effects, respectively, and es was the vector for random residual effects. The following (co)variance structure was assumed for random effects:
      Var[βγe]=[HW000ImW000InR],


      where β, γ, and e were additive genetic, permanent environmental, and residual effects, respectively; W and P were 9 × 9 covariance matrices of random regression coefficients for direct genetic and permanent environmental effects, respectively; R was a covariance matrix of random residual effects; Im was an identity matrix of size m × m for the permanent environmental effect (m is the number of cows with records); In was an identity matrix of size n × n for the residual (n is the number of test-day records); and ⊗ denotes the Kronecker product. Residual variances were assumed to be constant across THI classes because stable residual standard deviation across THI were estimated in preliminary phenotypic analyses for all traits. In a single-step genomic BLUP approach as developed by
      • Aguilar I.
      • Misztal I.
      • Johnson D.
      • Legarra A.
      • Tsuruta S.
      • Lawlor T.
      Hot topic: A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score.
      , the matrix H combines the pedigree-based numerator relationship matrix (A) with the genomic relationship matrix (G) matrix. The inverse of H was defined as H-1=A-1+[000Gw-1-A22-1], where Gw = 0.95G + 0.05A22 and A22 was the submatrix of A for the genotyped animals. The Gw matrix was constructed as proposed by
      • VanRaden P.M.
      Efficient methods to compute genomic predictions.
      to render it invertible.
      For the sth stage of lactation, additive and permanent environmental variances as a function of THI classes (i.e., σa(st)2 and σpe(st)2, respectively) were the diagonal elements of QWsQ′ and QPsQ′, respectively. In these equations, Ws and Ps were 3 × 3 covariance matrices of random regression coefficients, and Q was the matrix of Legendre polynomials. The heritability at the sth stage of lactation for the tth THI class ( h(st)2) with the assumption of homogeneous residual variance across THI classes ( σe(s)2) was h(st)2=σa(st)2σa(st)2+σpe(st)2+σe(s)2. The genetic correlation between same traits across THI classes at different stages of lactation was rg(i,j)=qtW(i,j)qtqtW(i)qt×qtW(j)qt, where W(i,j) was the genetic covariance matrix between the ith and the jth stages of lactation; W(i) and W(j) were the genetic covariance matrices for ith and jth stages of lactation, respectively, and qt was the vector of Legendre polynomials for the tth THI class. For the estimation of the average genetic correlation across stages of lactation, qt was a vector including the sum of Legendre polynomials for THI classes 1 to 10.
      Genetic parameters were estimated using the GIBBS2F90 program (

      Misztal, I., S. Tsuruta, T. Strabel, B. Auvray, T. Druet, and D. Lee. 2002. BLUPF90 and related programs. Communication no. 28–07 in Proc. of the 7th World Congress for the Genetic Applied Livestock Production, Montpellier, France.

      ). For each trait analysis, a single chain of 100,000 samples was run. The first 20,000 samples were considered as burn-in period and discarded. Posterior means of parameter estimates were calculated from every 50th sample. Convergence of the Gibbs chain was checked by visual inspection of the trace plots for genetic (co)variance components.
      Possible limitations implying potential biases on estimates for genetic (co)variance components address the modeling strategy for the environmental descriptor THI. Instead of modeling a continuous THI covariate, we created 10 different THI classes. In previous random regression or reaction norm model applications, we observed increased genetic variances at the extreme ends of the continuous descriptor scale, possibly due to limited number of observations for specific descriptor levels (e.g.,
      • Al-Kanaan A.
      • König S.
      • Brügemann K.
      Effects of heat stress on semen characteristics of Holstein bulls estimated on a continuous phenotypic and genetic scale.
      ). In consequence, in the present study, we created THI classes. However, such approach of data categorization does not fully reflect the trait-specific physiological background. For example, values during values during the thermoneutral zone may be influenced by the random regressions, causing a smoothed curve for the random regression coefficients.

      Assessment of Model Stability

      We focused on 2 strategies to create training sets [i.e., considering (1) a time and (2) a time plus an environmental perspective to assess model stability]. In the time approach, sires were divided into 2 groups according to their birth year (<2010 and ≥2010). In a next step, test-day records from 1,562 daughters of 96 genotyped sires born after 2010 (New_sire) were masked, comprising 14,571 records for test-day production traits and 8,328 records for FA.
      In the environmental approach, we additionally aimed on sire breeding value estimations under “new” extreme heat stress conditions. In this regard, test-day records from daughters of sires born after 2010 in THI classes 8 to 10 (New_env) were masked (comprising 4,067 and 1,819 records for test-day production traits and FA, respectively).
      Estimated breeding values from the model including all phenotypic records (full data set) were considered as the best estimates for the true breeding values of individuals. To assess the model stability of the RNM, we correlated the GEBV of new sires in the reduced data set with their GEBV in the full data set. Such validation was applied for both strategies (New_sire and New_env) in different stages of lactation. With regard to the 96 genotyped sires born after 2010, 18 sires had one genotyped daughter, 19 sires had 2 to 5 genotyped daughters, 25 sires had 6 to 20 genotyped daughters, and 34 sires had 21 to 122 genotyped daughters. The accuracy of GEBV for sires in the full data set was calculated as rit=1-PEVitσat2, where PEVit was the prediction error variance (PEV) for the ith sire in the tth THI class and σat2 was additive genetic variance in the tth THI class. According to
      • Tier B.
      • Meyer K.
      Approximating prediction error covariances among additive genetic effects within animals in multiple-trait and random regression models.
      , PEVit was the diagonal of QTiQ′, where Q was a 10 × 3 matrix with the values of 3 coefficients of the second-order Legendre polynomials for THI classes 1 to 10, and Ti was a 3 × 3 matrix of random regression coefficients of prediction error (co)variances for the ith sire. Accuracies of GEBV for sires from THI classes 8 to 10 were averaged (Supplemental Table S1, https://jlupub.ub.uni-giessen.de/handle/jlupub/51). Sires with averaged accuracies lower than 0.60 were excluded from the validation set. Accordingly, among production traits, the smallest and the largest validation sets were allocated for PY and MY with 72 and 93 sires, respectively. Among FA traits, the smallest (65 sires) and the largest (93 sires) validation sets were created for MUFA at the first stage of lactation and for SFA at the second stage of lactation, respectively.

      RESULTS AND DISCUSSION

      Determination of Climate Stress

      The number of test-day records across THI for different lactation stages is presented in Figure 1. Temperature-humidity index values lower than 51.4 (i.e., assigned to THI classes 1 and 2) were observed during winter months, and THI values higher than 60.9 mostly represent the summer period (i.e., the records assigned to THI classes 8 to 10).
      Least squares means for all traits in dependency of THI within lactation stage (results from model 1) were averaged across lactation stages (to reduce the number of figures) and are depicted in Figure 2. In THI classes 1 and 2, estimated daily average MY was lower than in THI classes representing the thermoneutral zone. The decline in MY and FY was also noticeable in THI class 10. A decline in PY was observed beyond THI class 7. The increase of THI was associated with a decrease in C16:0 and SFA contents. Fluctuating yield curves were observed for C18:0 and UFA (MUFA, PUFA, and UFA), with high contents at the extreme ends of the THI scale. The increase of average daily yield was substantial for PUFA in THI class 10.
      Figure thumbnail gr2
      Figure 2Least squares means averaged across stages of lactation with corresponding averaged SE for test-day milk yield (MY, in kg × 10−1), fat yield (FY, in kg), protein yield (PY, in kg), C16:0, MUFA, SFA, UFA, C18:0, and PUFA in dependency of the maximum hourly temperature-humidity index (THI). Units for fatty acids are g/100 g of milk.
      Changes in milk production, as well as in FA profiles, were identified under heat stress conditions in middle Europe (i.e., for THI larger than 62;
      • Hammami H.
      • Vandenplas J.
      • Vanrobays M.-L.
      • Rekik B.
      • Bastin C.
      • Gengler N.
      Genetic analysis of heat stress effects on yield traits, udder health, and fatty acids of Walloon Holstein cows.
      ;
      • Bohlouli M.
      • Alijani S.
      • Naderi S.
      • Yin T.
      • König S.
      Prediction accuracies and genetic parameters for test-day traits from genomic and pedigree-based random regression models with or without heat stress interactions.
      ).
      • Soyeurt H.
      • Dardenne P.
      • Dehareng F.
      • Bastin C.
      • Gengler N.
      Genetic parameters of saturated and monounsaturated fatty acid content and the ratio of saturated to unsaturated fatty acids in bovine milk.
      identified the lowest concentrations of SFA and the largest contents of UFA during summer.
      • Chilliard Y.
      • Ferlay A.
      • Doreau M.
      Dietary control of milk fat nutritional quality in the dairy cow: Trans and polyunsaturated fatty acids, and conjugated linoleic acid.
      reported seasonal effects on FA concentrations due to feeding particularities. The hot and humid climate caused a negative energy balance in early lactation, with ongoing effects on MY and milk compositions (
      • Penasa M.
      • Tiezzi F.
      • Gottardo P.
      • Cassandro M.
      • De Marchi M.
      Genetics of milk fatty acid groups predicted during routine data recording in Holstein dairy cattle.
      ). The mobilization of body fat depots in heat-stressed cows increased the content of UFA and, conversely, decreased SFA (
      • Gross J.
      • van Dorland H.A.
      • Bruckmaier R.M.
      • Schwarz F.J.
      Milk fatty acid profile related to energy balance in dairy cows.
      ;
      • Penasa M.
      • Tiezzi F.
      • Gottardo P.
      • Cassandro M.
      • De Marchi M.
      Genetics of milk fatty acid groups predicted during routine data recording in Holstein dairy cattle.
      ). Correspondingly, the significant increase of PUFA at high THI as identified in the present study, might be due to the negative energy balance of cows under heat stress conditions.
      About 9% of the phenotypic records belonged to test dates with air temperatures lower than 0°C.
      • Mader T.
      • Johnson L.
      • Gaughan J.
      A comprehensive index for assessing environmental stress in animals.
      associated very low air temperatures with cold stress and identified a pronounced MY decline. Similar observations for low temperatures were made in German Holstein cows from different production systems (
      • Brügemann K.
      • Gernand E.
      • König von Borstel U.
      • König S.
      Defining and evaluating heat stress thresholds in different dairy cow production systems.
      ). Low temperatures imply increased energy expenditure and feed intake in lactating cows to maintain homeothermy (
      • Collier R.J.
      • Renquist B.J.
      • Xiao Y.
      A 100-Year Review: Stress physiology including heat stress.
      ), with primary energy utilization for maintenance and restricted blood flow to the udder (
      • Johnson H.
      Bioclimate effects on growth, reproduction and milk production.
      ;
      • Brouček J.
      • Letkovičová M.
      • Kovalčuj K.
      Estimation of cold stress effect on dairy cows.
      ).

      Heritabilities and Genetic Variances

      For milk production and milk FA traits, estimates of additive genetic variances and heritabilities are presented in Supplemental Figure S1 (https://jlupub.ub.uni-giessen.de/handle/jlupub/51) and Figure 3, respectively, for all THI by parity effects. Heritabilities for MY ranged from 0.17 to 0.35 and were larger throughout than corresponding estimates for FY and PY. Heritabilities for C16:0 and SFA were larger than the estimates for UFA. Generally, heritabilities for FA are in line with estimates from test-day random regression models (
      • Soyeurt H.
      • Dardenne P.
      • Dehareng F.
      • Bastin C.
      • Gengler N.
      Genetic parameters of saturated and monounsaturated fatty acid content and the ratio of saturated to unsaturated fatty acids in bovine milk.
      ;
      • Bastin C.
      • Soyeurt H.
      • Gengler N.
      Genetic parameters of milk production traits and fatty acid contents in milk for Holstein cows in parity 1–3.
      ).
      • Bastin C.
      • Soyeurt H.
      • Gengler N.
      Genetic parameters of milk production traits and fatty acid contents in milk for Holstein cows in parity 1–3.
      reported average daily heritability estimates of 0.43 for C16:0, 0.25 for C18:0, 0.26 for MUFA, 0.31 for PUFA, 0.46 for SFA, and 0.27 for UFA.
      • Palmquist D.
      • Beaulieu A.D.
      • Barbano D.
      Feed and animal factors influencing milk fat composition.
      and
      • Dewhurst R.J.
      • Moorby J.M.
      • Vlaeminck B.
      • Fievez V.
      Apparent recovery of duodenal odd-and branched-chain fatty acids in milk of dairy cows.
      identified strong influence of genetic factors on variations for contents of SFA, implying moderate to large heritabilities for C16:0 and SFA. In contrast, diet compositions and environmental factors have a stronger effect on MUFA, PUFA, and UFA, explaining the moderate heritabilities. Generally, estimates for genetic variances and heritabilities indicate genetic variations and different genetic mechanisms for FA with alterations of lactation stages and of THI.
      Figure thumbnail gr3
      Figure 3Posterior estimates of heritabilities for test-day milk production and milk fatty acid traits in the first (solid line), second (dashed line), and third (dotted line) stage of lactation for distinct temperature-humidity index (THI) classes. Error bars represent SD of posterior estimates. MY = milk yield; PY = protein yield; FY = fat yield.

      Effect of Lactation Stage.

      For MY, PY, C16:0 and SFA, low genetic variances were estimated in the first stage of lactation (Supplemental Figure S1). Largest estimates were identified in the last stage of lactation. Very similar genetic variances across stages of lactation were estimated for C18:0, MUFA, PUFA, and UFA. For all traits, heritabilities in the first stage of lactation were lower than corresponding estimates for the second and the last lactation stage, except for C18:0 and MUFA, PUFA, and UFA in THI class 10. For MY, FA, and PY as well as for UFA in the first stage of lactation, heritabilities (Figure 3) were lower than 0.20. Heritabilities in the second and the last lactation stage were moderate to large for C16:0 and SFA, in the range from 0.42 to 0.61, and from 0.43 to 0.64, respectively. Regarding SFA, moderate to large heritabilities (0.26 in early lactation to 0.50 in mid-lactation) were reported by
      • Soyeurt H.
      • Dardenne P.
      • Dehareng F.
      • Bastin C.
      • Gengler N.
      Genetic parameters of saturated and monounsaturated fatty acid content and the ratio of saturated to unsaturated fatty acids in bovine milk.
      , considering Holstein dairy cattle from Belgium. For all traits, the larger heritability estimates in the second and the last stages of lactation might be due to the genetic component of lactation persistency (
      • Soyeurt H.
      • Dardenne P.
      • Dehareng F.
      • Bastin C.
      • Gengler N.
      Genetic parameters of saturated and monounsaturated fatty acid content and the ratio of saturated to unsaturated fatty acids in bovine milk.
      ).

      Effect of THI.

      Additive genetic and permanent environmental variances were quite constant for MY and PY over the entire range of THI (Supplemental Figure S1 and S2, respectively; https://jlupub.ub.uni-giessen.de/handle/jlupub/51). For FY, C16:0, and SFA, estimates for the additive genetic variance decreased with increasing THI. For C18:0 and UFA (MUFA, PUFA, and UFA) in the first stage of lactation, additive genetic and permanent environmental variances in THI class 1 and under heat stress (THI class 10) were substantially larger than under temperate climatic conditions (the middle THI classes). Increased genetic variances at the extreme ends of THI scale indicate that genotypes react differently because of environmental sensitivity (
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effects of heat stress in Italian Holstein dairy cattle.
      ). However, increased variances at the extreme ends of the THI scale also might be an artifact of random regression modeling based on second or third order Legendre polynomials. Similarly,
      • Gernand E.
      • König S.
      Random regression test-day model for clinical mastitis: Genetic parameters, model comparison, and correlations with indicator traits.
      reported increased genetic variances for health traits very early and late in lactation when modeling the time axis for DIM with second and third order Legendre polynomials. Hence, we assume 2 factors contributing to the increased trait variances: first, the artifacts of the statistical modeling approach, and second, pronounced genetic differentiation with increasing stress conditions. In such context,
      • Schierenbeck S.
      • Reinhardt F.
      • Reents R.
      • Simianer H.
      • König S.
      Identification of informative cooperator herds for progeny testing based on yield deviations.
      found stronger dispersion of daughter yield deviations for SCS in herds representing harsh environments. Nevertheless, C18:0, MUFA, PUFA, and UFA were identified as climatically sensitive traits. With regard to the large genetic variances at the extreme ends of THI scale, selection on C18:0, MUFA, PUFA, and UFA will contribute to thermotolerance improvements genetically.
      For all traits, the heritability patterns across THI (Figure 3) reflect the genetic variance alterations. Heritabilities for milk production traits, C16:0, and SFA decreased with increasing THI. At the first stage of lactation, heritabilities for C18:0 and UFA increased at the extreme ends of THI scale. The C16:0 and SFA as well as MY, FY, and PY displayed quite stable heritabilities in dependency of THI.
      • Bohlouli M.
      • Alijani S.
      • Naderi S.
      • Yin T.
      • König S.
      Prediction accuracies and genetic parameters for test-day traits from genomic and pedigree-based random regression models with or without heat stress interactions.
      identified very similar heritabilities for MY by THI in early lactation, but decreasing MY heritabilities with increasing THI in late lactation. Heritability alterations for milk FA traits across THI are comparable with results reported by
      • Hammami H.
      • Vandenplas J.
      • Vanrobays M.-L.
      • Rekik B.
      • Bastin C.
      • Gengler N.
      Genetic analysis of heat stress effects on yield traits, udder health, and fatty acids of Walloon Holstein cows.
      , who applied a linear reaction norm model with first-order Legendre polynomials for THI, but ignored within lactation genetic components. In the second and the last stages of lactation, heritabilities for SFA, MUFA, and UFA decreased with increasing THI. Comparably,
      • Hammami H.
      • Vandenplas J.
      • Vanrobays M.-L.
      • Rekik B.
      • Bastin C.
      • Gengler N.
      Genetic analysis of heat stress effects on yield traits, udder health, and fatty acids of Walloon Holstein cows.
      estimated lower heritabilities for C16:0 and SFA under heat stress conditions. Furthermore, they reported a slight heritability increase for PUFA with increasing THI, but quite constant heritabilities for UFA across the THI range.

      Genetic Correlations Across Lactation Stages

      Estimates of genetic correlations averaged across THI between different stages of lactation are presented in Table 2. The lowest genetic correlations were estimated between traits from the first lactation stage with the corresponding trait from stages 2 and 3. In this regard, the lowest genetic correlation was 0.69, estimated for UFA. All genetic correlations between same traits from the second and the third stages of lactation were larger than 0.96.
      Table 2Posterior estimates for genetic correlations averaged across temperature-humidity index and corresponding SD (in parentheses) between test-day milk yield (MY), fat yield (FY), protein yield (PY), C16:0, C18:0, MUFA, PUFA, UFA, and SFA from different lactation stages
      TraitBetween lactation stages
      Stage 1: 5–105 DIM; stage 2: 106–205 DIM; stage 3: 206–305 DIM.
      1 and 21 and 32 and 3
      MY0.91 (0.014)0.82 (0.026)0.97 (0.006)
      FY0.92 (0.014)0.84 (0.028)0.96 (0.009)
      PY0.90 (0.019)0.82 (0.033)0.97 (0.007)
      C16:00.97 (0.007)0.94 (0.014)0.99 (0.004)
      C18:00.92 (0.025)0.88 (0.036)0.99 (0.006)
      MUFA0.79 (0.069)0.73 (0.088)0.98 (0.008)
      PUFA0.87 (0.034)0.78 (0.056)0.98 (0.008)
      SFA0.97 (0.007)0.95 (0.013)0.99 (0.003)
      UFA0.76 (0.084)0.69 (0.107)0.97 (0.012)
      1 Stage 1: 5–105 DIM; stage 2: 106–205 DIM; stage 3: 206–305 DIM.
      The between-stage genetic correlations considering same traits for distinct THI classes are presented in Figure 4. Estimated genetic correlations between the second and the third stages of lactation were quite large for all milk production and milk FA traits. Between the first and remaining lactation stages, genetic correlations larger than 0.75 were estimated for MY, FY, PY, C16:0, and SFA. For C18:0 and UFA in the middle THI classes, genetic correlations between traits from the first stage of lactation with the same traits from stages 2 and 3 were lower than 0.50, but increased at the extreme ends of the THI scale.
      Figure thumbnail gr4
      Figure 4Posterior estimates for genetic correlations between test-day milk production traits from different lactation stages and between milk fatty acids from different lactation stages for distinct temperature-humidity index (THI) classes. Posterior SD of genetic correlations ranged from 0.005 to 0.096. MY = milk yield; PY = protein yield; FY = fat yield.
      During early lactation, cows experience a negative energy balance due to the decline in DMI. In consequence, cows mobilize their lipid reserves from adipose tissues to provide energy requirements for maintenance and production (
      • Bell A.W.
      Regulation of organic nutrient metabolism during transition from late pregnancy to early lactation.
      ). Body fat mobilization increases the plasma nonesterified FA concentration (
      • Bell A.W.
      Regulation of organic nutrient metabolism during transition from late pregnancy to early lactation.
      ;
      • de Vries M.
      • Veerkamp R.
      Energy balance of dairy cattle in relation to milk production variables and fertility.
      ), which is used by the mammary glands for milk fat synthesis (
      • Bielak A.
      • Derno M.
      • Tuchscherer A.
      • Hammon H.
      • Susenbeth A.
      • Kuhla B.
      Body fat mobilization in early lactation influences methane production of dairy cows.
      ). Accordingly, the content of preformed C18 FA (i.e., mainly 18:1, 18:2, and 18:3) increased in milk fat in the early stage of lactation (
      • Palmquist D.
      • Beaulieu A.D.
      • Barbano D.
      Feed and animal factors influencing milk fat composition.
      ;
      • Chilliard Y.
      • Ferlay A.
      • Doreau M.
      Dietary control of milk fat nutritional quality in the dairy cow: Trans and polyunsaturated fatty acids, and conjugated linoleic acid.
      ;
      • Bastin C.
      • Soyeurt H.
      • Gengler N.
      Genetic parameters of milk production traits and fatty acid contents in milk for Holstein cows in parity 1–3.
      ). In addition, for energy supply, cows are generally fed diets enriched with concentrates during the first stage of lactation (
      • Humer E.
      • Petri R.
      • Aschenbach J.
      • Bradford B.
      • Penner G.
      • Tafaj M.
      • Südekum K.-H.
      • Zebeli Q.
      Invited review: Practical feeding management recommendations to mitigate the risk of subacute ruminal acidosis in dairy cattle.
      ). In this regard,
      • Palmquist D.
      • Beaulieu A.D.
      • Barbano D.
      Feed and animal factors influencing milk fat composition.
      and
      • Krag K.
      • Poulsen N.A.
      • Larsen M.K.
      • Larsen L.B.
      • Janss L.L.
      • Buitenhuis B.
      Genetic parameters for milk fatty acids in Danish Holstein cattle based on SNP markers using a Bayesian approach.
      identified strong feed composition influences on UFA than on SFA concentrations in milk fat. The described metabolic pathways and different feeding managements across stages of lactation contribute to increasing intraanimal variations (between stages of lactation) for C18:0 and UFA. The given arguments may explain differences in genetic regulation mechanisms in early lactation for C18:0, MUFA, PUFA, and UFA syntheses, and thus, the low genetic correlations between the same traits from the first and remaining lactation stages. The low genetic correlations between same traits across lactation stages especially were observed in the temperate THI classes. However, with regard to both extreme THI classes, large genetic correlations were estimated for C18:0, MUFA, PUFA, or UFA across lactation stages, indicating similar genetic mechanisms for these traits under stress conditions.

      Genetic Correlations Across THI

      Genetic correlations between milk traits at THI class 10 with same milk traits at remaining THI classes, separately for lactations stages 1, 2, and 3, are presented in Figure 5. For MY and FY, genetic correlations larger than 0 0.80 in all lactation stages disproved any genotype by environment interactions across THI classes. Among milk production traits, PY responded most sensitive to heat stress. In the first stage of lactation, genetic correlations between PY from THI class 10 with PY from THI classes 1 to 5 were lower than 0.80. Regarding milk proteins,
      • Bernabucci U.
      • Basiricò L.
      • Morera P.
      • Dipasquale D.
      • Vitali A.
      • Cappelli F.P.
      • Calamari L.
      Effect of summer season on milk protein fractions in Holstein cows.
      ,
      • Hu H.
      • Wang J.Q.
      • Gao H.N.
      • Li S.L.
      • Zhang Y.D.
      • Zheng N.
      Heat-induced apoptosis and gene expression in bovine mammary epithelial cells.
      , and
      • Ma L.
      • Yang Y.
      • Zhao X.
      • Wang F.
      • Gao S.
      • Bu D.
      Heat stress induces proteomic changes in the liver and mammary tissue of dairy cows independent of feed intake: An iTRAQ study.
      found a decline in the contents of α-caseins and β-caseins under heat stress conditions.
      Figure thumbnail gr5
      Figure 5Posterior estimates for genetic correlations between test-day milk production traits and milk fatty acids from temperature-humidity index (THI) class 10 with respective traits from the remaining THI classes for lactation stages 1, 2, and 3. Posterior SD of genetic correlations ranged from 0.001 to 0.054. MY = milk yield; PY = protein yield; FY = fat yield.
      For C18:0, MUFA, and UFA in the first lactation stage, genetic correlations between same traits from THI class 10 and from THI classes 1 to 7 were in the range from −0.54 to 0.25, indicating notable genotype by environment interactions. The low genetic correlations indicate climate sensitivity for these traits. Climate sensitivity already was identified above, when studying the pronounced genetic variance heterogeneity for C18:0 and UFA along the THI gradient for thermoneutral and heat stress conditions (Supplemental Figure S1).
      The bioenergetics status under heat stress is comparable with physiological mechanisms, which cows experience during early lactation (
      • Moore C.
      • Kay J.
      • Collier R.J.
      • VanBaale M.
      • Baumgard L.
      Effect of supplemental conjugated linoleic acids on heat-stressed Brown Swiss and Holstein cows.
      ;
      • Hammami H.
      • Vandenplas J.
      • Vanrobays M.-L.
      • Rekik B.
      • Bastin C.
      • Gengler N.
      Genetic analysis of heat stress effects on yield traits, udder health, and fatty acids of Walloon Holstein cows.
      ).
      • Moore C.
      • Kay J.
      • Collier R.J.
      • VanBaale M.
      • Baumgard L.
      Effect of supplemental conjugated linoleic acids on heat-stressed Brown Swiss and Holstein cows.
      reported that heat-stressed dairy cows lose significant amounts of BW, indicating a negative energy balance (
      • Bastin C.
      • Soyeurt H.
      • Gengler N.
      Genetic parameters of milk production traits and fatty acid contents in milk for Holstein cows in parity 1–3.
      ). Therefore, with focus on the first lactation stage, we hypothesize that heat stress intensifies an increase of UFA due to the mobilization of adipose FA. Quantitative genetically, such physiological mechanisms might explain the low genetic correlations between UFA in THI class 10 with the same traits from THI classes 2 to 6.
      In the second lactation stage, genetic correlations between MY, FY, PY, and FA from THI class 10 with the same trait from remaining THI classes were generally larger than 0.80. The smallest genetic correlation was 0.75, estimated for UFA at THI classes 4 with UFA at THI class 10. The genetic correlations between same traits from different THI classes in the third lactation stage reflect the genetic correlation pattern from lactation stage 2. Again, the lowest genetic correlation with 0.6 was estimated for UFA considering THI classes 5 and 10.

      Quantifying Environmental Sensitivity

      Changes in GEBV of the 8 top-ranked and 8 bottom-ranked genotyped sires along THI classes for the 3 lactation stages are presented in Supplemental Figure S3 (https://jlupub.ub.uni-giessen.de/handle/jlupub/51). Ranking of sires was done according to the sum of THI GEBV for the respective lactation stage. For MY, FY, C16:0, and SFA, the marginal re-rankings among top sires across THI do not justify separate genetic evaluations for different climate conditions. In contrast, regarding the top sires for PY in the second and the last stages of lactation, some re-rankings were observed. For example, the top ranked sire at low THI had a GEBV below the average at high THI. For C18:0 and UFA, re-rankings of sires across THI classes were stronger in the first stage of lactation than in stages 2 or 3, reflecting the genetic correlation estimates as presented in Figure 5. Regarding Figure 5, genetic correlations between these traits at THI 10 and at remaining THI classes were only moderate in the first lactation stage, but quite high for the remaining stages. For UFA in the first stage of lactation, breeding values of top sires increased at the extreme ends of THI scale. In contrast, for the bottom sires, breeding values were highest in the middle THI classes. Such breeding value differences indicate different responses to environmental stress for sires with high or low genetic merit for milk FA. Daughters of top sires reacted more sensitive on THI alterations than the daughters of the bottom sires, especially in the stress period after calving in the first stage of lactation.

      Model Stability

      Model stability was the correlation between GEBV of new sires in the reduced data set with their GEBV in the full data set. The model stabilities for both validation strategies New_sire and New_env in different lactation stage strategies for milk production and FA traits are given in Table 3. Correlations as indicators for the model stability for MY, FY, and PY were larger than corresponding values for FA. Larger correlations for milk production traits than for milk FA may be due to the larger number of phenotypic records in genetic evaluations for MY, FY, and PY (
      • de Roos A.
      • Hayes B.
      • Goddard M.
      Reliability of genomic predictions across multiple populations.
      ;
      • Hayes B.J.
      • Bowman P.J.
      • Chamberlain A.
      • Goddard M.
      Invited review: Genomic selection in dairy cattle: Progress and challenges.
      ). For all milk FA traits, the model stability was quite poor in the first stage of lactation. With regard to the New_sire validation, the lowest correlations were 0.23 (PUFA), 0.33 (MUFA), and 0.33 (UFA). The heritability as a major component describing the genetic trait architecture is a further parameter influencing the model stability (
      • Hayes B.J.
      • Bowman P.J.
      • Chamberlain A.
      • Goddard M.
      Invited review: Genomic selection in dairy cattle: Progress and challenges.
      ). The correlations increased with increasing trait heritabilities in the present as well as in previous studies (e.g.,
      • Yao C.
      • De Los Campos G.
      • VandeHaar M.
      • Spurlock D.
      • Armentano L.
      • Coffey M.
      • De Haas Y.
      • Veerkamp R.
      • Staples C.
      • Connor E.
      • Wang Z.
      • Hanigan M.D.
      • Tempelman R.J.
      • Weigel K.A.
      Use of genotype × environment interaction model to accommodate genetic heterogeneity for residual feed intake, dry matter intake, net energy in milk, and metabolic body weight in dairy cattle.
      ;
      • Mota L.F.M.
      • Fernandes Jr., G.A.
      • Herrera A.C.
      • Scalez D.C.B.
      • Espigolan R.
      • Magalhães A.F.B.
      • Carvalheiro R.
      • Baldi F.
      • Albuquerque L.G.
      Genomic reaction norm models exploiting genotype × environment interaction on sexual precocity indicator traits in Nellore cattle.
      ).
      Table 3Correlations between genomic estimated breeding values (GEBV) of new sires in the reduced data set and their GEBV in the full data set for test-day milk yield (MY), fat yield (FY), protein yield (PY), C16:0, C18:0, MUFA, PUFA, UFA, and SFA at different lactation stages
      Stage 1: 5–105 DIM; stage 2: 106–205 DIM; stage 3: 206–305 DIM.
      for 2 validation strategies
      ItemValidation strategy
      New_sire
      New_sire: daughter records from genotyped sires born after 2010 with at least 2 genotyped daughters were masked in genomic predictions.
      New_env
      New_env: daughter records in temperature-humidity index classes 8 to 10 from genotyped sires born after 2010 were masked in genomic predictions.
      Stage 1Stage 2Stage 3Stage 1Stage 2Stage 3
      MY0.800.810.800.920.920.91
      FY0.750.740.720.880.890.88
      PY0.770.790.790.870.900.91
      C16:00.400.420.430.470.530.53
      C18:00.350.390.390.460.510.51
      MUFA0.330.360.380.430.440.49
      PUFA0.230.300.320.450.500.46
      SFA0.390.440.430.460.540.52
      UFA0.330.340.350.430.460.47
      1 Stage 1: 5–105 DIM; stage 2: 106–205 DIM; stage 3: 206–305 DIM.
      2 New_sire: daughter records from genotyped sires born after 2010 with at least 2 genotyped daughters were masked in genomic predictions.
      3 New_env: daughter records in temperature-humidity index classes 8 to 10 from genotyped sires born after 2010 were masked in genomic predictions.
      Dissimilar validation approaches used in genomic analyses complicate comparisons across studies. Regarding RNM, correlations between GEBV of the youngest sires with the phenotype adjusted for fixed effects (
      • Mota L.F.M.
      • Fernandes Jr., G.A.
      • Herrera A.C.
      • Scalez D.C.B.
      • Espigolan R.
      • Magalhães A.F.B.
      • Carvalheiro R.
      • Baldi F.
      • Albuquerque L.G.
      Genomic reaction norm models exploiting genotype × environment interaction on sexual precocity indicator traits in Nellore cattle.
      ) or with their daughter yield deviations (
      • Zhang Z.
      • Kargo M.
      • Liu A.
      • Thomasen J.R.
      • Pan Y.
      • Su G.
      Genotype-by-environment interaction of fertility traits in Danish Holstein cattle using a single-step genomic reaction norm model.
      ) were used as indicators for the model stability. Small correlations were found by
      • Zhang Z.
      • Kargo M.
      • Liu A.
      • Thomasen J.R.
      • Pan Y.
      • Su G.
      Genotype-by-environment interaction of fertility traits in Danish Holstein cattle using a single-step genomic reaction norm model.
      for low heritability female fertility traits.
      • Mota L.F.M.
      • Fernandes Jr., G.A.
      • Herrera A.C.
      • Scalez D.C.B.
      • Espigolan R.
      • Magalhães A.F.B.
      • Carvalheiro R.
      • Baldi F.
      • Albuquerque L.G.
      Genomic reaction norm models exploiting genotype × environment interaction on sexual precocity indicator traits in Nellore cattle.
      reported correlations smaller than 0.50 for early pregnancy in heifers and scrotal circumference traits.
      • Tiezzi F.
      • de Los Campos G.
      • Gaddis K.P.
      • Maltecca C.
      Genotype by environment (climate) interaction improves genomic prediction for production traits in US Holstein cattle.
      focused on comprehensive evaluations of cross-validation approaches. For the characterization of cow reactions on environmental alterations, they proposed climate variables for cross-validations.
      In the current study, model stabilities for all traits substantially improved (~plus 30%) when comparing results from the New_env strategy with results from the New_sire strategy. Accordingly,
      • Guo G.
      • Zhao F.
      • Wang Y.
      • Zhang Y.
      • Du L.
      • Su G.
      Comparison of single-trait and multiple-trait genomic prediction models.
      identified an increase of prediction accuracies for the trait of interest when considering information from correlated traits from other environments. In simulations, the prediction accuracy was improved when progressively including phenotyped cows from additional environments (
      • Yin T.
      • Pimentel E.
      • König v. Borstel U.
      • König S.
      Strategy for the simulation and analysis of longitudinal phenotypic and genomic data in the context of a temperature× humidity-dependent covariate.
      ;
      • Bohlouli M.
      • Alijani S.
      • Javaremi A.N.
      • König S.
      • Yin T.
      Genomic prediction by considering genotype× environment interaction using different genomic architectures.
      ). Gain in model stabilities due to increased daughter records (i.e., the comparison New_env versus New_sire) was very obvious for UFA in the first stage of lactation, especially for PUFA with 0.23 for New_sire and 0.45 for New_env. For production traits, prediction accuracies were quite large for both strategies New_sire and New_env. However, regarding UFA in the early lactation stage, the effect of the negative energy balance might bias genomic predictions. Generally, for the new sires, we suggest consideration of daughter records from temperate climates in genomic predictions under heat stress conditions.

      CONCLUSIONS

      Additive genetic variances for C18:0, MUFA, PUFA, and UFA were larger at extreme THI than under temperate climate conditions. Such indications for environmental sensitivity suggest selective improvement of heat tolerance based on these FA in environments representing heat stress. High genetic correlations across THI classes in different stages of lactation disprove possible genotype by environment interactions for MY, FY, C16:0, and SFA. Genetic correlations were quite small between C18:0 and UFA under thermoneutral with corresponding traits under thermal stress conditions. Finally, our results supported previous studies, that in the presence of environmental sensitivity, UFA contents in early lactation can be very affordable traits to improve thermotolerance using genomic RNM. For all FA traits, we showed a moderate genomic model stability of RNM for sires in heat stress environments when considering daughter records from temperate climates. For UFA, genetic (co)variance components were very similar for low and for high THI, which may reflect the same genetic mechanisms under stress conditions. Therefore, further research is needed to identify a combination of environmental descriptors or alternative functions for modeling continuous THI that can describe both cold and heat stresses.

      ACKNOWLEDGMENTS

      The authors gratefully acknowledge the financial support provided by the German Research Foundation, DFG, through grant number KO 3520/8-1. We also acknowledge funding for M. Bohlouli from the Alexander von Humboldt Foundation (Bonn, Germany) and his access throughout his project to the resources of the University of Liège – Gembloux Agro-Bio Tech partly supported by the Fonds de la Recherche Scientifique – FNRS under Grants n° T.0095.19 (PDR “DEEPSELECT”) and J.0174.18 (CDR “PREDICT-2”). The authors have not stated any conflicts of interest.

      REFERENCES

        • Aguilar I.
        • Misztal I.
        • Johnson D.
        • Legarra A.
        • Tsuruta S.
        • Lawlor T.
        Hot topic: A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score.
        J. Dairy Sci. 2010; 93 (20105546): 743-752
        • Aguilar I.
        • Misztal I.
        • Tsuruta S.
        • Legarra A.
        • Wang H.
        PREGSF90 – POSTGSF90: Computational tools for the implementation of single-step genomic selection and genome-wide association with ungenotyped individuals in BLUPF90 program. Proc. 10th WCGALP, Vancouver, Canada.
        • Al-Kanaan A.
        • König S.
        • Brügemann K.
        Effects of heat stress on semen characteristics of Holstein bulls estimated on a continuous phenotypic and genetic scale.
        Livest. Sci. 2015; 177: 15-24
        • Bastin C.
        • Soyeurt H.
        • Gengler N.
        Genetic parameters of milk production traits and fatty acid contents in milk for Holstein cows in parity 1–3.
        J. Anim. Breed. Genet. 2013; 130 (23496012): 118-127
        • Bell A.W.
        Regulation of organic nutrient metabolism during transition from late pregnancy to early lactation.
        J. Anim. Sci. 1995; 73 (8582872): 2804-2819
        • Bernabucci U.
        • Basiricò L.
        • Morera P.
        • Dipasquale D.
        • Vitali A.
        • Cappelli F.P.
        • Calamari L.
        Effect of summer season on milk protein fractions in Holstein cows.
        J. Dairy Sci. 2015; 98 (25547301): 1815-1827
        • Bernabucci U.
        • Biffani S.
        • Buggiotti L.
        • Vitali A.
        • Lacetera N.
        • Nardone A.
        The effects of heat stress in Italian Holstein dairy cattle.
        J. Dairy Sci. 2014; 97 (24210494): 471-486
        • Bielak A.
        • Derno M.
        • Tuchscherer A.
        • Hammon H.
        • Susenbeth A.
        • Kuhla B.
        Body fat mobilization in early lactation influences methane production of dairy cows.
        Sci. Rep. 2016; 6 (27306038)28135
        • Bohlouli M.
        • Alijani S.
        • Javaremi A.N.
        • König S.
        • Yin T.
        Genomic prediction by considering genotype× environment interaction using different genomic architectures.
        Ann. Anim. Sci. 2017; 17: 683-701
        • Bohlouli M.
        • Alijani S.
        • Naderi S.
        • Yin T.
        • König S.
        Prediction accuracies and genetic parameters for test-day traits from genomic and pedigree-based random regression models with or without heat stress interactions.
        J. Dairy Sci. 2019; 102 (30343923): 488-502
        • Bohmanova J.
        • Misztal I.
        • Cole J.
        Temperature-humidity indices as indicators of milk production losses due to heat stress.
        J. Dairy Sci. 2007; 90 (17369235): 1947-1956
        • Bohmanova J.
        • Misztal I.
        • Tsuruta S.
        • Norman H.
        • Lawlor T.
        Genotype by environment interaction due to heat stress.
        J. Dairy Sci. 2008; 91 (18218772): 840-846
        • Brouček J.
        • Letkovičová M.
        • Kovalčuj K.
        Estimation of cold stress effect on dairy cows.
        Int. J. Biometeorol. 1991; 35 (1917124): 29-32
        • Brügemann K.
        • Gernand E.
        • König von Borstel U.
        • König S.
        Defining and evaluating heat stress thresholds in different dairy cow production systems.
        Arch. Tierzucht. 2012; 55: 13-24
        • Brügemann K.
        • Gernand E.
        • Von Borstel U.
        • König S.
        Genetic analyses of protein yield in dairy cows applying random regression models with time-dependent and temperature x humidity-dependent covariates.
        J. Dairy Sci. 2011; 94 (21787948): 4129-4139
        • Cheruiyot E.K.
        • Nguyen T.T.T.
        • Haile-Mariam M.
        • Cocks B.G.
        • Abdelsayed M.
        • Pryce J.E.
        Genotype-by-environment (temperature-humidity) interaction of milk production traits in Australian Holstein cattle.
        J. Dairy Sci. 2020; 103 (31864748): 2460-2476
        • Chilliard Y.
        • Ferlay A.
        • Doreau M.
        Dietary control of milk fat nutritional quality in the dairy cow: Trans and polyunsaturated fatty acids, and conjugated linoleic acid.
        Productions Animales (France). 2001; 14: 323-335
        • Collier R.J.
        • Renquist B.J.
        • Xiao Y.
        A 100-Year Review: Stress physiology including heat stress.
        J. Dairy Sci. 2017; 100 (29153170): 10367-10380
        • de Roos A.
        • Hayes B.
        • Goddard M.
        Reliability of genomic predictions across multiple populations.
        Genetics. 2009; 183 (19822733): 1545-1553
        • de Vries M.
        • Veerkamp R.
        Energy balance of dairy cattle in relation to milk production variables and fertility.
        J. Dairy Sci. 2000; 83 (10659965): 62-69
        • Dewhurst R.J.
        • Moorby J.M.
        • Vlaeminck B.
        • Fievez V.
        Apparent recovery of duodenal odd-and branched-chain fatty acids in milk of dairy cows.
        J. Dairy Sci. 2007; 90 (17369218): 1775-1780
        • Falconer D.S.
        Selection in different environments—Effects on environmental sensitivity (reaction norm) and on mean performance.
        Genet. Res. 1990; 56: 57-70
        • Falconer D.S.
        • Mackay T.F.C.
        Introduction to Quantitative Genetics.
        4th ed. Longman, 1996
        • Gernand E.
        • König S.
        Random regression test-day model for clinical mastitis: Genetic parameters, model comparison, and correlations with indicator traits.
        J. Dairy Sci. 2014; 97 (24731633): 3953-3963
        • Gross J.
        • van Dorland H.A.
        • Bruckmaier R.M.
        • Schwarz F.J.
        Milk fatty acid profile related to energy balance in dairy cows.
        J. Dairy Res. 2011; 78 (21843394): 479-488
        • Guo G.
        • Zhao F.
        • Wang Y.
        • Zhang Y.
        • Du L.
        • Su G.
        Comparison of single-trait and multiple-trait genomic prediction models.
        BMC Genet. 2014; 15 (24593261): 30
        • Haile-Mariam M.
        • Carrick M.
        • Goddard M.
        Genotype by environment interaction for fertility, survival, and milk production traits in Australian dairy cattle.
        J. Dairy Sci. 2008; 91 (19038960): 4840-4853
        • Hammami H.
        • Bormann J.
        • M'hamdi N.
        • Montaldo H.H.
        • Gengler N.
        Evaluation of heat stress effects on production traits and somatic cell score of Holsteins in a temperate environment.
        J. Dairy Sci. 2013; 96 (23313002): 1844-1855
        • Hammami H.
        • Rekik B.
        • Bastin C.
        • Soyeurt H.
        • Bormann J.
        • Stoll J.
        • Gengler N.
        Environmental sensitivity for milk yield in Luxembourg and Tunisian Holsteins by herd management level.
        J. Dairy Sci. 2009; 92 (19700723): 4604-4612
        • Hammami H.
        • Vandenplas J.
        • Vanrobays M.-L.
        • Rekik B.
        • Bastin C.
        • Gengler N.
        Genetic analysis of heat stress effects on yield traits, udder health, and fatty acids of Walloon Holstein cows.
        J. Dairy Sci. 2015; 98 (25958288): 4956-4968
        • Hayes B.J.
        • Bowman P.J.
        • Chamberlain A.
        • Goddard M.
        Invited review: Genomic selection in dairy cattle: Progress and challenges.
        J. Dairy Sci. 2009; 92 (19164653): 433-443
        • Herbut P.
        • Angrecka S.
        • Walczak J.
        Environmental parameters to assessing of heat stress in dairy cattle—A review.
        Int. J. Biometeorol. 2018; 62 (30368680): 2089-2097
        • Hijmans R.J.
        • Williams E.
        • Vennes C.
        • Hijmans M.R.J.
        Package 'geosphere’.
        https://cran.r-project.org/web/packages/geosphere/index.html
        Date: 2016
        Date accessed: November 19, 2017
        • Hill D.L.
        • Wall E.
        Dairy cattle in a temperate climate: The effects of weather on milk yield and composition depend on management.
        Animal. 2015; 9 (25315451): 138-149
        • Hu H.
        • Wang J.Q.
        • Gao H.N.
        • Li S.L.
        • Zhang Y.D.
        • Zheng N.
        Heat-induced apoptosis and gene expression in bovine mammary epithelial cells.
        Anim. Prod. Sci. 2016; 56: 918-926
        • Humer E.
        • Petri R.
        • Aschenbach J.
        • Bradford B.
        • Penner G.
        • Tafaj M.
        • Südekum K.-H.
        • Zebeli Q.
        Invited review: Practical feeding management recommendations to mitigate the risk of subacute ruminal acidosis in dairy cattle.
        J. Dairy Sci. 2018; 101 (29153519): 872-888
        • Johnson H.
        Bioclimate effects on growth, reproduction and milk production.
        in: Johnson H.D. Bioclimatology and the Adaptation of Livestock. Elsevier, 1987: 35-57
        • Krag K.
        • Poulsen N.A.
        • Larsen M.K.
        • Larsen L.B.
        • Janss L.L.
        • Buitenhuis B.
        Genetic parameters for milk fatty acids in Danish Holstein cattle based on SNP markers using a Bayesian approach.
        BMC Genet. 2013; 14 (24024882): 79
        • Lynch M.
        • Walsh B.
        Genetics and Analysis of Quantitative Traits.
        Sinauer Assoc, 1998
        • Ma L.
        • Yang Y.
        • Zhao X.
        • Wang F.
        • Gao S.
        • Bu D.
        Heat stress induces proteomic changes in the liver and mammary tissue of dairy cows independent of feed intake: An iTRAQ study.
        PLoS One. 2019; 14 (30625175)e0209182
        • Mader T.
        • Johnson L.
        • Gaughan J.
        A comprehensive index for assessing environmental stress in animals.
        J. Anim. Sci. 2010; 88 (20118427): 2153-2165
        • Misztal I.
        • Strabel T.
        • Jamrozik J.
        • Mäntysaari E.A.
        • Meuwissen T.H.E.
        Strategies for estimating the parameters needed for different test-day models.
        J. Dairy Sci. 2000; 83 (10821589): 1125-1134
      1. Misztal, I., S. Tsuruta, T. Strabel, B. Auvray, T. Druet, and D. Lee. 2002. BLUPF90 and related programs. Communication no. 28–07 in Proc. of the 7th World Congress for the Genetic Applied Livestock Production, Montpellier, France.

        • Moore C.
        • Kay J.
        • Collier R.J.
        • VanBaale M.
        • Baumgard L.
        Effect of supplemental conjugated linoleic acids on heat-stressed Brown Swiss and Holstein cows.
        J. Dairy Sci. 2005; 88 (15829665): 1732-1740
        • Mota L.F.M.
        • Fernandes Jr., G.A.
        • Herrera A.C.
        • Scalez D.C.B.
        • Espigolan R.
        • Magalhães A.F.B.
        • Carvalheiro R.
        • Baldi F.
        • Albuquerque L.G.
        Genomic reaction norm models exploiting genotype × environment interaction on sexual precocity indicator traits in Nellore cattle.
        Anim. Genet. 2020; 51 (31944356): 210-223
        • NRC
        A Guide to Environmental Research on Animals.
        National Academies, 1971
        • Palmquist D.
        • Beaulieu A.D.
        • Barbano D.
        Feed and animal factors influencing milk fat composition.
        J. Dairy Sci. 1993; 76 (8326036): 1753-1771
        • Penasa M.
        • Tiezzi F.
        • Gottardo P.
        • Cassandro M.
        • De Marchi M.
        Genetics of milk fatty acid groups predicted during routine data recording in Holstein dairy cattle.
        Livest. Sci. 2015; 173: 9-13
        • Ravagnolo O.
        • Misztal I.
        Genetic component of heat stress in dairy cattle, parameter estimation.
        J. Dairy Sci. 2000; 83 (11003247): 2126-2130
        • Ravagnolo O.
        • Misztal I.
        Effect of heat stress on nonreturn rate in Holstein cows: Genetic analyses.
        J. Dairy Sci. 2002; 85 (12487476): 3092-3100
        • Schierenbeck S.
        • Reinhardt F.
        • Reents R.
        • Simianer H.
        • König S.
        Identification of informative cooperator herds for progeny testing based on yield deviations.
        J. Dairy Sci. 2011; 94 (21426998): 2071-2082
        • Segelke D.
        • Chen J.
        • Liu Z.
        • Reinhardt F.
        • Thaller G.
        • Reents R.
        Reliability of genomic prediction for German Holsteins using imputed genotypes from low-density chips.
        J. Dairy Sci. 2012; 95 (22916947): 5403-5411
        • Soyeurt H.
        • Dardenne P.
        • Dehareng F.
        • Bastin C.
        • Gengler N.
        Genetic parameters of saturated and monounsaturated fatty acid content and the ratio of saturated to unsaturated fatty acids in bovine milk.
        J. Dairy Sci. 2008; 91 (18765620): 3611-3626
        • Soyeurt H.
        • Dehareng F.
        • Gengler N.
        • McParland S.
        • Wall E.
        • Berry D.P.
        • Coffey M.
        • Dardenne P.
        Mid-infrared prediction of bovine milk fatty acids across multiple breeds, production systems, and countries.
        J. Dairy Sci. 2011; 94 (21426953): 1657-1667
        • Tian H.
        • Zheng N.
        • Wang W.
        • Cheng J.B.
        • Li S.L.
        • Zhang Y.D.
        • Wang J.Q.
        Integrated metabolomics study of the milk of heat-stressed lactating dairy cows.
        Sci. Rep. 2016; 6 (27048914)24208
        • Tier B.
        • Meyer K.
        Approximating prediction error covariances among additive genetic effects within animals in multiple-trait and random regression models.
        J. Anim. Breed. Genet. 2004; 121: 77-89
        • Tiezzi F.
        • de Los Campos G.
        • Gaddis K.P.
        • Maltecca C.
        Genotype by environment (climate) interaction improves genomic prediction for production traits in US Holstein cattle.
        J. Dairy Sci. 2017; 100 (28109596): 2042-2056
        • Umphrey J.
        • Moss B.
        • Wilcox C.
        • Van Horn H.
        Interrelationships in lactating Holsteins of rectal and skin temperatures, milk yield and composition, dry matter intake, body weight, and feed efficiency in summer in Alabama.
        J. Dairy Sci. 2001; 84 (11814024): 2680-2685
        • VanRaden P.M.
        Efficient methods to compute genomic predictions.
        J. Dairy Sci. 2008; 91 (18946147): 4414-4423
        • Wiggans G.
        • Sonstegard T.
        • VanRaden P.
        • Matukumalli L.
        • Schnabel R.
        • Taylor J.
        • Schenkel F.
        • Van Tassell C.
        Selection of single-nucleotide polymorphisms and quality of genotypes used in genomic evaluation of dairy cattle in the United States and Canada.
        J. Dairy Sci. 2009; 92 (19528621): 3431-3436
        • Yao C.
        • De Los Campos G.
        • VandeHaar M.
        • Spurlock D.
        • Armentano L.
        • Coffey M.
        • De Haas Y.
        • Veerkamp R.
        • Staples C.
        • Connor E.
        • Wang Z.
        • Hanigan M.D.
        • Tempelman R.J.
        • Weigel K.A.
        Use of genotype × environment interaction model to accommodate genetic heterogeneity for residual feed intake, dry matter intake, net energy in milk, and metabolic body weight in dairy cattle.
        J. Dairy Sci. 2017; 100 (28109605): 2007-2016
        • Yin T.
        • König S.
        Heritabilities and genetic correlations in the same traits across different strata of herds created according to continuous genomic, genetic, and phenotypic descriptors.
        J. Dairy Sci. 2018; 101 (29248231): 2171-2186
        • Yin T.
        • Pimentel E.
        • König v. Borstel U.
        • König S.
        Strategy for the simulation and analysis of longitudinal phenotypic and genomic data in the context of a temperature× humidity-dependent covariate.
        J. Dairy Sci. 2014; 97 (24485687): 2444-2454
        • Zhang Z.
        • Kargo M.
        • Liu A.
        • Thomasen J.R.
        • Pan Y.
        • Su G.
        Genotype-by-environment interaction of fertility traits in Danish Holstein cattle using a single-step genomic reaction norm model.
        Heredity. 2019; 123 (30760882): 202-214