Advertisement

Genotype by environment interaction due to heat stress in Brown Swiss cattle

Open AccessPublished:December 29, 2022DOI:https://doi.org/10.3168/jds.2021-21551

      ABSTRACT

      Due to its geographical position and a highly variable orography, Italy is characterized by several climatic areas and thus, by many different dairy cow farming systems. Brown Swiss cattle, in this context, are a very appreciated genetic resource for their adaptability and low metabolic requirement. The significant heterogeneity in farming systems may consist of genotype by environment (G × E) interactions with neglected changes in animals' rank position. The objective of this study was to investigate G × E for heat tolerance in Brown Swiss cattle for several production traits (milk, fat, and protein yield in kilograms; fat, protein, and cheese yield in percentage) and 2 derivate traits (fat-corrected milk and energy-corrected milk). We used the daily maximum temperature-humidity index (THI) range, calculated according to weather stations' data from 2008 to 2018 in Italy, and 202,776 test-day records from 23,396 Brown Swiss cows from 639 herds. Two different methodologies were applied to estimate the effect of the environmental variable (THI) on genetic parameters: (1) the reaction norm model, which uses a continuous random covariate to estimate the animal additive effect, and (2) the multitrait model, which splits each production pattern as a distinct and correlated trait according to the first (a thermal comfort condition), third (a moderate heat stress condition), and fifth (a severe heat stress condition) mean THI value quintile. The results from the reaction norm model showed a descending trend of the additive genetic effect until THI reached the value of 80. Then we recorded an increase with high extreme THI values (THI 90). Permanent environmental variance at increasing THI values revealed an opposite trend: The plot of heritability and the ratio of animal permanent environmental variance to phenotypic variance showed that when the environmental condition worsens, the additive genetic and permanent environmental component for production traits play a growing role. The negative additive genetic correlation between slope and linear random coefficient indicates no linear relationship between the production traits or under heat stress conditions, except for milk yield and protein yield. In tridimensional wireframe plots, the extreme margin decreases until a minimum of ~0.90 of genetic correlation in the ECM trait, showing that the magnitude of G × E interaction is greater than the other traits. Genetic correlation values in Brown Swiss suggest the possibility of moderate changes in animals' estimated breeding value in heat stress conditions. Results indicated a moderate G × E interaction but significant variability in sire response related to their production level.

      Key words

      INTRODUCTION

      Modern breeding programs in animals selected for milk production have focused on quantitative and qualitative production, indirectly leading to an increase in the animals' sensitivity to the environment (
      • Bryant J.R.
      • Lopez-Villalobos N.
      • Pryce J.E.
      • Holmes C.W.
      • Johnson D.L.
      • Garrick D.J.
      Environmental sensitivity in New Zealand dairy cattle.
      ;
      • Herbut P.
      • Angrecka S.
      • Walczak J.
      Environmental parameters to assessing of heat stress in dairy cattle—A review.
      ;
      • Martinez-Castillero M.
      • Toledo-Alvarado H.
      • Pegolo S.
      • Vazquez A.I.
      • de Los Campos G.
      • Varona L.
      • Finocchiaro R.
      • Bittante G.
      • Cecchinato A.
      Genetic parameters for fertility traits assessed in herds divergent in milk energy output in Holstein-Friesian, Brown Swiss, and Simmental cattle.
      ). Although nowadays, all selection schemes include traits related to the resilience of cows and their longevity (
      • Stocco G.
      • Cipolat-Gotet C.
      • Gasparotto V.
      • Cecchinato A.
      • Bittante G.
      Breed of cow and herd productivity affect milk nutrient recovery in curd, and cheese yield, efficiency and daily production.
      ;
      • De Vries A.
      • Marcondes M.I.
      Review: Overview of factors affecting productive lifespan of dairy cows.
      ), little has been done regarding animals' tolerance to environmental temperature changes (
      • Rovelli G.
      • Ceccobelli S.
      • Perini F.
      • Demir E.
      • Mastrangelo S.
      • Conte G.
      • Abeni F.
      • Marletta D.
      • Ciampolini R.
      • Cassandro M.
      • Bernabucci U.
      • Lasagna E.
      The genetics of phenotypic plasticity in livestock in the era of climate change: A review.
      ). However, the wide use of cooling systems in intensive milking production systems to mitigate the negative effects of heat stress has considerably increased. Cooling systems represent an additional cost and increase water and energy demand for milk production in warm climate conditions. Therefore, mainly from a progressive global warning perspective, the mid- and long-term overall aim of livestock production systems is to select more tolerant and effectively acclimating animals (
      • Dillon P.
      • Buckley F.
      • O'Connor P.
      • Hegarty D.
      • Rath M.
      A comparison of different dairy cow breeds on a seasonal grass-based system of milk production: 1. Milk production, live weight, body condition score and DM intake.
      ).
      The temperature-humidity index (THI) allows for the estimation of when dairy cows become stressed due to heat and humidity levels, mainly in intensive farming systems. It can be calculated with different formulas, even if the most suitable ones give greater weight to humidity (
      • Bohmanova J.
      • Misztal I.
      • Cole J.B.
      Temperature-humidity indices as indicators of milk production losses due to heat stress.
      ). In the first papers dealing with this topic (
      • Ravagnolo O.
      • Misztal I.
      • Hoogenboom G.
      Genetic component of heat stress in dairy cattle, development of heat index function.
      ;
      • Aguilar I.
      • Misztal I.
      • Tsuruta S.
      Genetic components of heat stress for dairy cattle with multiple lactations.
      ), the average daily maximum THI values of the 3 d before the test day were suitable for explaining the highest degrees of the observed variance due to thermal stress. Moreover, it is interesting that
      • Negri R.
      • Aguilar I.
      • Feltes G.L.
      • Machado J.D.
      • Braccini Neto J.
      • Costa-Maia F.M.
      • Cobuci J.A.
      Inclusion of bioclimatic variables in genetic evaluations of dairy cattle.
      observed that the diurnal temperature variation was also a relevant factor in assessing animal welfare.
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effects of heat stress in Italian Holstein dairy cattle.
      used daily mean THI in the 15 d before the test days. Because it is challenging to directly measure stress on individuals, it must be quantified as the relationship between the genetic component of the animal and the environment. Different approaches have been proposed to estimate the genetic parameters for different traits related to heat stress.
      The reaction norm (RN) models fitted through random regression models are a powerful tool to identify and quantify the genotype × environment (G × E) interaction. An RN model is a mixed-effect model in which additive genetic (AG) effects are modeled using random regression coefficients (
      • Martin J.G.A.
      • Nussey D.H.
      • Wilson A.J.
      • Réale D.
      Measuring individual differences in reaction norms in field and experimental studies: A power analysis of random regression models.
      ;
      • Windig J.J.
      • Mulder H.A.
      • Bohthe-Wilhelmus D.I.
      • Veerkamp R.F.
      Simultaneous estimation of genotype by environment interaction accounting for discrete and continuous environmental descriptors in Irish dairy cattle.
      ). A second approach considers the phenotypes expressed in different environments as different but correlated traits (
      • 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.
      ;
      • Martinez-Castillero M.
      • Toledo-Alvarado H.
      • Pegolo S.
      • Vazquez A.I.
      • de Los Campos G.
      • Varona L.
      • Finocchiaro R.
      • Bittante G.
      • Cecchinato A.
      Genetic parameters for fertility traits assessed in herds divergent in milk energy output in Holstein-Friesian, Brown Swiss, and Simmental cattle.
      ).
      The choice depends on using the environmental gradient as a continuous or discrete factor, even if both situations may coincide (
      • Windig J.J.
      • Mulder H.A.
      • Bohthe-Wilhelmus D.I.
      • Veerkamp R.F.
      Simultaneous estimation of genotype by environment interaction accounting for discrete and continuous environmental descriptors in Irish dairy cattle.
      ).
      Previous studies of genetic parameters have been carried out in several European countries (
      • Brügemann K.
      • Gernand E.
      • von Borstel U.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.
      ;
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effects of heat stress in Italian Holstein dairy cattle.
      ;
      • Carabaño M.J.
      • Logar B.
      • Bormann J.
      • Minet J.
      • Vanrobays M.L.
      • Díaz C.
      • Tychon B.
      • Gengler N.
      • Hammami H.
      Modeling heat stress under different environmental conditions.
      ), Australia (Nguyen et al.; 2016;
      • 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.
      ), New Zealand (
      • Bryant J.R.
      • Lopez-Villalobos N.
      • Pryce J.E.
      • Holmes C.W.
      • Johnson D.L.
      • Garrick D.J.
      Environmental sensitivity in New Zealand dairy cattle.
      ), China (
      • Shi R.
      • Brito L.F.
      • Liu A.
      • Luo H.
      • Chen Z.
      • Liu L.
      • Guo G.
      • Mulder H.
      • Ducro B.
      • van der Linden A.
      • Wang Y.
      Genotype-by-environment interaction in Holstein heifer fertility traits using single-step genomic reaction norm models.
      ), and South America (
      • Santana Jr., M.L.
      • Bignardi A.B.
      • Pereira R.J.
      • Menéndez-Buxadera A.
      • El Faro L.
      Random regression models to account for the effect of genotype by environment interaction due to heat stress on the milk yield of Holstein cows under tropical conditions.
      ), though they focused on the Holstein Friesian breeds and even Bos indicus breeds in tropical climate conditions (
      • Santana Jr., M.L.
      • Pereira R.J.
      • Bignardi A.B.
      • Filho A.E.
      • Menendez-Buxadera A.
      • El Faro L.
      Detrimental effect of selection for milk yield on genetic tolerance to heat stress in purebred Zebu cattle: Genetic parameters and trends.
      ). A few studies tested the association between THI values and production traits in other breeds, such as Jersey (
      • Smith D.L.
      • Smith T.
      • Rude B.J.
      • Ward S.H.
      Short communication: Comparison of the effects of heat stress on milk and component yields and somatic cell score in Holstein and Jersey cows.
      ) and Brown Swiss (
      • Maggiolino A.
      • Dahl G.E.
      • Bartolomeo N.
      • Bernabucci U.
      • Vitali A.
      • Serio G.
      • Cassandro M.
      • Centoducati G.
      • Santus E.
      • De Palo P.
      Estimation of maximum thermo-hygrometric index thresholds affecting milk production in Italian Brown Swiss cattle.
      ). However, the estimation of genetic parameters needs to be deepened, as it plays a pivotal role in understanding the differences among dairy cattle breeds. The present study aimed to investigate the G × E interaction in Italian Brown Swiss cattle using an RN analysis model and comparing the results with a classical approach of a multitrait model in different productive traits.

      MATERIALS AND METHODS

      Ethics Statement

      Animal welfare and use committee approval was not needed for this study because data sets were obtained from pre-existing databases based on routine animal production recording procedures.

      Data Set

      Data were provided by the Italian Brown Swiss Breeders Association (Verona, Italy). They comprised 202,776 test-day records (from 2008 to 2018) of 23,396 Brown Swiss cows from 639 herds, including the following traits: milk yield (MY) (kg/d), protein percentage (PP), daily protein yield (PY), fat percentage (FP), and daily fat yield (FY). Additionally, the 4% FCM yield was calculated for each test-day record according to the formula reported by
      • Gaines W.L.
      • Davidson F.A.
      Relation between Percentage Fat Content and Yield of Milk: Correction of Milk Yield for Fat Content.
      :
      4%FCM=0.4×MY+15×FY.


      The ECM yield was calculated for each test-day record according to the formula reported by
      • Yan M.-J.
      • Humphreys J.
      • Holden N.M.
      An evaluation of life cycle assessment of European milk production.
      :
      ECMkg/d=MY×[0.25+(0.122×FP)+(0.077×PP)].


      Cheese yield percentage (CHY) was calculated according to the formula reported by
      • Van Slyke L.L.
      • Publow C.A.
      The Science and Practice of Cheese-Making: A Treatise on the Manufacture of American Cheddar Cheese and Other Varieties.
      , modified for Parmigiano-Reggiano cheese (
      • Formaggioni P.
      • Summer A.
      • Malacarne M.
      • Franceschi P.
      • Mucchetti G.
      Italian and Italian-style hard cooked cheeses: Predictive formulas for Parmigiano-Reggiano 24-h cheese yield.
      ):
      CHY%=[(0.93×FP+(casein%0.1))×1.09]/(100M),


      in which the fat recovery factor is 0.93, the factor accounting for casein loss is 0.1; 1.09 represents the constituents of cheese solids (noncasein and nonfat), corresponding to 9% of the sum; and M is the moisture of fresh cheese, assumed to be 39.5%, as for 24 h Parmigiano-Reggiano cheese (
      • Formaggioni P.
      • Summer A.
      • Malacarne M.
      • Franceschi P.
      • Mucchetti G.
      Italian and Italian-style hard cooked cheeses: Predictive formulas for Parmigiano-Reggiano 24-h cheese yield.
      ).
      The data relating to the structure of the population are presented in Table 1 and traits' descriptive statistics are reported in Table 2. Full pedigree observation was used in the estimation.
      Table 1Descriptive statistics of data used in the study
      ItemValue
      Number of herds583
      Number of records174,875
      Herd test date/number of herds12,068
      Number of cows21,242
      Number of sires1,507
      Average number of cows per herd test date14.48
      Average number of daughters per sire14.09
      Table 2Overall description of the different production traits used in the present paper, clusterized according to each of the 3 quintiles used in the multitrait model
      PP = protein percentage; PY = protein yield; FP = fat percentage; FY = fat yield; MY = milk yield; CHY = cheese percentage; min = minimum; max = maximum; N = number of test-day records.
      SetItemPP (%)PY (kg/d)FP (%)FY (kg/d)MY (kg/d)4% FCM (kg/d)ECM (kg/d)CHY (%)
      OverallMin2.300.051.190.0021.21.391.450.062
      Max5.031.906.712.3455.255.7055.410.17
      Mean3.670.944.031.0325.8925.8626.570.12
      SD0.380.270.720.337.657.757.750.015
      N172,835173,111171,817172,401174,040172,769172,87182,772
      First quintileMin2.3500.0821.4600.0352.0002.0722.1800.064
      Max5.0301.8856.6002.29853.80053.72855.4170.172
      Mean3.7320.9524.1381.05325.74726.12426.8870.120
      SD0.3940.2700.7160.3417.7787.9337.9330.015
      N30,25030,35430,11630,18230,49530,26530,29216,550
      Third quintileMin2.3200.0851.4300.0021.9001.7581.9570.062
      Max5.0001.8516.6302.28054.80054.34254.9870.172
      Mean3.6880.9444.0501.03425.80325.86926.6000.118
      SD0.3760.2610.7050.3297.5447.6437.6350.015
      N32,02632,01331,76531,89732,10131,93931,97416,550
      Fifth quintileMin2.3800.0521.1900.0011.2001.3891.4500.062
      Max5.0201.8826.5702.27352.80054.58553.2790.170
      Mean3.5760.9253.8851.00026.03025.45726.0950.115
      SD0.3690.2570.7070.3207.5067.4857.4690.014
      N35,31935,38335,18835,33335,69335,39035,35417,155
      1 PP = protein percentage; PY = protein yield; FP = fat percentage; FY = fat yield; MY = milk yield; CHY = cheese percentage; min = minimum; max = maximum; N = number of test-day records.
      A first editing was performed fitting a linear model for each trait separately:
      Yijkl=DIMi+PARj+YMCk+rijkl,


      where Y is the single-trait records l for each test day at i class of DIM defined as follows: from 0 to 30, from 31 to 90, from 91 to 180, from 181 to 270, and from 271 to 365; PAR is the parity class j (1, 2, 3, and >4); YMC the k year-month of calving class; and r the residual. Records of test dates with residuals falling outside 3 standard deviations from the mean were removed from the data set. Model and data manipulation were performed in R software (
      • R Core Team
      R: A Language and Environment for Statistical Computing.
      ) using lm function and the tidyverse package (
      • Wickham H.
      • Averick M.
      • Bryan J.
      • Chang W.
      • D'Agostino McGowan L.
      • François R.
      • Grolemund G.
      • Hayes A.
      • Henry L.
      • Hester J.
      • Kuhn M.
      • Pedersen T.L.
      • Miller E.
      • Bache S.M.
      • Müller K.
      • Ooms J.
      • Robinson D.
      • Seidel D.P.
      • Spinu V.
      • Takahashi K.
      • Vaughan D.
      • Wilke C.
      • Woo K.
      • Yutani H.
      Welcome to the tidyverse.
      ).

      Statistical Methods

      All statistical analyses were carried out separately for each trait. Data preparation, editing, and graphics were performed using the R programming environment v.4.0.1 (
      • R Core Team
      R: A Language and Environment for Statistical Computing.
      ).
      Age at calving (AC) was included as 2 mo class from birth date to each parity (59 classes). The contemporary group effect (HTD) was constructed by concatenating the herd to the test day date (year, month, and day), retaining classes with at least 5 records.
      The THI was calculated with the commonly used formula (
      • NRC
      A Guide to Environmental Research on Animals.
      ) using the records from 76 weather stations located throughout the Italian territory, based on maximum daily temperature and relative humidity. As described in a previous work (
      • Maggiolino A.
      • Dahl G.E.
      • Bartolomeo N.
      • Bernabucci U.
      • Vitali A.
      • Serio G.
      • Cassandro M.
      • Centoducati G.
      • Santus E.
      • De Palo P.
      Estimation of maximum thermo-hygrometric index thresholds affecting milk production in Italian Brown Swiss cattle.
      ), the farms that had no weather stations within 10 km, were located higher than 700 m above sea level, were linked to a weather station at a different altitude (>50 m), or had significant orography or hydrographic bodies separating them from the weather station were discarded. The THI included in the model was the average maximum daily THI value of the 5 d before test day (THIm). It was considered as were used as environmental variable, using second-order Legendre polynomial (LP) coefficients. We used the 5 d value to test the effect of duration of heat stress rather than its intensity only, as applied by other authors (
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effects of heat stress in Italian Holstein dairy cattle.
      ). The effects of heat stress and THI thresholds have been previously investigated in Brown Swiss, and MY traits begin to decline at THI >70 (
      • Maggiolino A.
      • Dahl G.E.
      • Bartolomeo N.
      • Bernabucci U.
      • Vitali A.
      • Serio G.
      • Cassandro M.
      • Centoducati G.
      • Santus E.
      • De Palo P.
      Estimation of maximum thermo-hygrometric index thresholds affecting milk production in Italian Brown Swiss cattle.
      ,
      • Maggiolino A.
      • Landi V.
      • Bartolomeo N.
      • Bernabucci U.
      • Santus E.
      • Bragaglio A.
      • De Palo P.
      Effect of heat waves on some Italian Brown Swiss dairy cows' production patterns.
      ). Based on these results, the THI threshold was set at 70 (i.e., if THI <70, then THI = 70) for all traits. Tests with THI values >90 were arbitrarily given a value of 90. This was done to avoid possible artifacts of variance estimation using RN models, which might lead to unexpected trajectories at extreme THIm values as environmental descriptor due to few data points (
      • Misztal I.
      • Strabel T.
      • Jamrozik J.
      • Mantysaari E.A.
      • Meuwissen T.H.
      Strategies for estimating the parameters needed for different test-day models.
      ).
      To estimate (co)variance components both for random regression and for multitrait analysis, the Gibbs sampling procedure was carried out by the THRGIBBS1F90 program (
      • Tsuruta S.
      • Misztal I.
      THRGIBBS1F90 for estimation of variance components with threshold-linear models.
      ). A flat prior was used as fixed effects, and an inverted Wishart distribution was used as prior for the random effects. For each analysis, 500,000 samples were run and every 100th sample was saved after discarding a burn-in of 100,000 iterations, for a total of 4,000 samples. The outputs were obtained using the Postgibbsf90 software (
      • Tsuruta S.
      • Misztal I.
      THRGIBBS1F90 for estimation of variance components with threshold-linear models.
      ). Convergence was calculated from a visual inspection of trace plots. The analyses were run on the ReCaS server at the University of Bari (https://www.recas-bari.it/) using a 10 core Intel Xeon CPU E5–2650L v.2 (Intel Corporation) at 1.70 gigahertz with 37 gigabytes of random-access memory.

      Reaction Norm Analysis Model

      Reaction norm is a random regression mixed model more suitable to a continuous descriptor coefficient. A single-trait animal model was fitted with the following RN model:
      Yijklmo=μ+HTDi+PARj+DIMk+ACl+YMCm+n=0lpazonaon+n=0lppzonpeon+eijklmo,


      where Yijklmo is the measurement of each trait (separately) of each cow o at herd-test-date class i; µ is overall mean; HTD is the fixed effect of herd-test-date class i;PAR is the fixed effect of j parity class; DIM the fixed effect day in milk class k; AC is the fixed effect of age at calving class l; YMC is the fixed of class m effect of year of calving nested with month of calving; aon and peon are the nth random regression of AG and permanent environmental (PE) effects for the oth cow, respectively; zon is the nth-order of LP for cow o; lpa and lpp denote the order of LP (intercept, first and second LP coefficient) for AG and PE effects, respectively; and eijklmo is the random residual. The variance or covariance structure for additive animal effects were the following, given a and p are the vectors of random regression additive and PE coefficient of cow o at parity j andr the residual variance vector, then
      Var[apr]=[AG0000IP0000IR0],


      where A is the numerator relationship matrix; I is the identity matrix; G0 and P0 are 3 × 3 matrices of (co)variances for additive and permanent environmental effects, respectively; R0 is residual variances corresponding to each trait; and ⊗ denotes the Kronecker product of matrices.
      Heritability (ha2) for each trait in the RN model was calculated as follows:
      ha2=σa2σa2+σpe2+σr2,


      where σa2, σpe2, σr2 are the additive genetic, permanent environmental, and residual variance.
      The estimate of genetic correlation, RGs,i, and PE correlation, RPEs,i, between intercept and linear regression coefficient (i and s) were calculated as
      RGs,i=σaa1σa×σa1andRPEs,i=σpepe1σpe×σpe1.


      Finally, the estimate of random regression coefficients had been used for the calculation of response of animals across the THI variation scale. The variance structure for the AG components for all animals (a) was assumed to be
      var(a)=GoA,


      with Go representing (co)variances between random regression coefficients for the AG component and A is the AG relationships matrix.
      Genetic (co)variances for performance under thermal loads T and T′ were calculated as follows, where Z and Z′ are the vector of Legendre polynomial covariables associated with each regression coefficient evaluated at temperature humidity T:
      GT,T=Z(T)GoZ(T).


      Similarly, PE (co)variances under thermal loads T and T′ have been calculated as
      PT,T=Z(T)PoZ(T).


      Overall phenotypic variance was obtained as
      VT=GT,T+PT,T+RV,


      where RV is the estimated residual variance (assumed homogeneous along the THI values).
      Then, h2 and ratio of animal permanent environmental variance on the phenotypic variance (PEE) of the trait at THI value were obtained as
      hT2=GT,T'VTandPEET,T'=PT,T'VT


      Genetic and PE correlations (RGTT' and RPETT') between performance under THI values T and T′ were obtained as follows, where sqrt is the square root:
      RGT,T=GT,Tsqrt(GT×GT)andRPET,T=PT,Tsqrt(PT×PT).


      The EBV deviation for animal i at a given T (THI) were computed from the estimated solutions for genetic regression coefficients and the estimated (co)variance matrices:
      gkT=Z(T)ak, 


      where gT represents the AG value for performance of animal k at the thermal load T, Z′(T) is the vector of LP covariables associated with each regression coefficient evaluated at T. Then we selected bulls with at least 90 daughters, resulting in a set of 88 subjects for all traits (68 for cheese yield production in percentage) for rank calculations and plots. Three ranks were compared based on traditional breeding (e.g., intercept of RN model) values at THI 76 and THI 90, respectively.

      Multitrait Analysis

      In this section, we investigated the efficiency and applicability of multitrait multienvironment (MT), considering combined different lactations. In Australian Holsteins, according to
      • Hayes B.J.
      • Carrick M.
      • Bowman P.
      • Goddard M.E.
      Genotype × environment interaction for milk production of daughters of Australian dairy sires from test-day records.
      , milk production traits are affected at THI values >60. Mild signs of heat stress are observed at a THI of 68 to 74, and a THI ≥75 caused drastic decreases in production efficiency (
      • Aguilar I.
      • Misztal I.
      • Tsuruta S.
      Short communication: Genetic trends of milk yield under heat stress for US Holsteins.
      ;
      • De Rensis F.
      • Garcia-Ispierto I.
      • López-Gatius F.
      Seasonal heat stress: Clinical implications and hormone treatments for the fertility of dairy cows.
      ). Stress condition limits are variable among traits, breeds, and countries, so it is difficult to set a general threshold limit. Therefore, in our work, we first divided the observation according to the THIm quintiles and considered the first within 34≤ THIm ≤55, where cows are generally considered to be in thermal comfort conditions; the third within 63≤ THIm ≤69, where cows are considered to be in low or mild heat stress; and the fifth within 78≤ THIm ≤95, in which cows are considered to be under severe heat stress conditions, following the general THI thresholds described in
      • Renaudeau D.
      • Collin A.
      • Yahav S.
      • de Basilio V.
      • Gourdine J.L.
      • Collier R.J.
      Adaptation to hot climate and strategies to alleviate heat stress in livestock production.
      . The chosen THI range is larger when compared with other studies (
      • 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.
      ), due, in our case, to the smallest data set, which did not allow the selection of a narrow class. The considered quintiles, containing 33,306, 35,311, and 39,469 test-day records, respectively (Figure 1 and Table 2), showed a mean number of records per animal of 8 ±5 (minimum 1 and maximum 42). Each of the 3 selected quintiles was considered a different but correlated trait. Then, the following MT animal model (MT model) was fitted:
      Figure thumbnail gr1
      Figure 1Distribution of temperature-humidity index (THI) values in the data used in this study divided in quintiles.
      Yijklmot=μ+HTDi+PARj+DIMk+ACl+YMCm+cot+pot+eijklmot,


      where Yijklmot is a measurement of the 3 traits consisting in the first, third, and fifth quintile of each original trait separately, referring to each cow o;µ is the overall mean of the trait; fixed effects were the same above reported in the RN model; cot was the random additive effect of the cow o for each trait t, and pot was the random PE effect for cow o for each trait t; eijklmot is the residual effect for each trait t.
      In the covariance structure, cp, pp, and rp, were the additive, PE, and residual effect of each cow for the trait t, respectively:
      Var[cppprp]=[AG0000IP0000IR0], 


      where A is the numerator relationship matrix; I is the identity matrix; and G0 and P0 are 3 × 3 matrices of (co)variances for additive and PE effects, respectively. R0 is a diagonal matrix of residual variances corresponding to each trait (first, third, and fifth quintile), and ⊗ denotes the Kronecker product of matrices. Residual covariance was fixed at zero because each animal is performing in only one environment at the same time.
      Heritability and ratio of animal permanent environmental variance to phenotypic variance (PEE) for each trait (quintile 1, 3, or 5) were respectively calculated as follows:
      hiq2=σaiq2σaiq2+σpeiq2+σriq2andPEEiq=σpeiq2σaiq2+σpeiq2+σriq2,


      where σ2aiq, σ2peiq, and σ2riq are the AG, PE, and residual variance, respectively, of trait i at quintile q (with q = 1, 3, or 5, alternatively).
      Estimates of genetic correlation and PE RGiq and RPEiq were calculated as
      RGiq=σaiqsqrt(σaiq2×σaiq2)andRPEiq=σpeiqsqrt(σpeiq2×σpeiq2),


      where σaiq and σpeiq are the additive and permanent (co)variance of trait i with trait i′ at quintile q (q = 1, 3, or 5).

      RESULTS

      Additive Genetic and PE Variances

      The estimates for the MT model shown in Table 3, Table 4 are comparable with PP values and slightly higher for PY for RN (Table 5).
      Table 3Additive genetic (AG), permanent environmental (PE), and residual (RE) variances; h2 estimates, and lower and upper bounds of the 95% highest posterior density (HPD95) region from the multitrait model for all investigated traits
      Trait
      CHY= cheese percentage.
      First quintileHPD95Third quintileHPD95Fifth quintileHPD95
      AGPEREh2AGPEREh2AGPEREh2
      Protein (%)0.0250.0150.0430.300.276–0.3380.0250.0100.0400.330.304–0.3570.0230.0110.0360.330.291–0.359
      Protein yield (kg/d)0.0030.0110.0200.080.061–0.1050.0030.0100.0200.080.064–0.1070.0030.0100.0190.080.062–0.107
      Fat (%)0.0590.0420.2760.170.132–0.1770.0610.0260.2870.160.139–0.1840.0510.0330.2880.140.117–0.159
      Fat yield (kg/d)0.0050.0140.0430.070.055–0.0950.0040.0120.0460.070.051–0.0880.0030.0120.0430.060.042–0.074
      Milk (kg/d)2.7909.09414.3520.110.079–0.1302.3007.96815.7160.090.069–0.1102.0998.96014.3880.080.061–0.103
      4% FCM (kg/d)2.3648.85919.0150.090.059–0.1032.0267.06720.4430.070.051–0.0881.5237.84518.8690.050.039–0.068
      ECM (kg/d)2.3959.15318.2210.080.056–0.1022.0167.37419.5120.070.051–0.0861.7307.98317.7140.060.047–0.082
      CHY (% × 100)35.24322.09693.1250.230.193–0.28031.62317.09197.7090.220.184–0.25028.81216.39694.2650.210.170–0.242
      1 CHY= cheese percentage.
      Table 4Genetic and permanent environmental correlation in MT model for first, third, and fifth quintiles
      Trait
      CHY = cheese percentage.
      Additive genetic correlationPermanent environmental correlation
      First vs. third (HPD95)
      HPD95 = lower and upper bounds of the 95% highest posterior density region.
      Third vs. fifth (HPD95)First vs. fifth (HPD95)First vs. third (HPD95)Third vs. fifth (HPD95)First vs. fifth (HPD95)
      Protein (%)0.996 (0.987–0.999)0.985 (0.972–0.995)0.973 (0.950–0.988)0.973 (0.966–0.979)0.785 (0.740–0.827)0.621 (0.548–0.686)
      Protein yield (kg/d)0.995 (0,990–0.999)0.983 (0.967–0.997)0.965 (0.933–0.993)0.961 (0.954–0.968)0.874 (0.857–0.889)0.707 (0.675–0.742)
      Fat (%)0.997 (0.995–0.999)0.966 (0.947–0.988)0.961 (0.938–0.989)0.898 (0.852–0.939)0.911 (0.868–0.957)0.65 (0.518–0.759)
      Fat yield (kg/d)0.986 (0.979–0.997)0.948 (0.904–0.989)0.951 (0.905–0.982)0.982 (0.972–0.989)0.839 (0.808–0.873)0.726 (0.680–0.770)
      Milk (kg/d)0.995 (0.981–0.999)0.977 (0.957–0.994)0.953 (0.896–0.990)0.969 (0.958–0.979)0.862 (0.838–0.886)0.714 (0.682–0.743)
      4% FCM (kg/d)0.993 (0.988–0.998)0.965 (0.932–0.989)0.977 (0.958–0.993)0.968 (0.956–0.975)0.865 (0.843–0.884)0.713 (0.671–0.747)
      ECM (kg/d)0.989 (0.972–0.999)0.975 (0.954–0.996)0.944 (0.886–0.994)0.966 (0.954–0.974)0.873 (0.854–0.891)0.719 (0.681–0.759)
      CHY (% × 100)0.995 (0.992–0.998)0.977 (0.955–0.996)0.954 (0.913–0.988)0.954 (0.923–0.987)0.844 (0.776–0.914)0.661(0.551–0.784)
      1 CHY = cheese percentage.
      2 HPD95 = lower and upper bounds of the 95% highest posterior density region.
      Table 5Additive genetic, permanent environmental, and residual variance for random regression coefficients and h2 for the investigated traits obtained using reaction norm model
      σa02 = intercept additive genetic variance; σa12 and σa22 = first- and second-order slope additive genetic variance, respectively; σpe02 = intercept permanent environmental variance; σpe12 and σpe22 = first- and second-order slopes permanent environmental variance; σe2 = residual variance; ha2 = heritability; HPD95 = lower and upper bounds of the 95% highest posterior density region.
      Trait
      CHY = cheese percentage.
      σa02σa1σa22σpe02σpe1σpe22σe2ha2 (HPD95)
      Protein (%)0.0473.66E-044.89E-050.0203.33E-043.67E-040.0400.44 (0.404–0.475)
      Protein yield (kg/d)0.0057.80E-053.56E-052.19E-040.0200.12 (0.097–0.147)
      Fat (%)0.1040.0012.76E-040.0671.29E-032.05E-030.2840.23 (0.199–0.257)
      Fat yield (kg/d)0.0067.13E-053.90E-050.0244.01E-042.93E-040.408 (0.065–0.106)
      Milk (kg/d)4.3800.0131.65E-0216.7150.3880.16314.9350.12 (0.093–0.150)
      4% FCM (kg/d)3.0240.0252.54E-0215.0590.2230.17319.5070.08 (0.062–0.100)
      ECM (kg/d)3.0990.0562.19E-0215.5540.2120.20518.5360.08 (0.066–0.103)
      CHY (% × 100)0.5780.0023.43E-030.3270.0110.0050.9720.31 (0.266–0.349)
      1 σa02 = intercept additive genetic variance; σa12 and σa22 = first- and second-order slope additive genetic variance, respectively; σpe02 = intercept permanent environmental variance; σpe12 and σpe22 = first- and second-order slopes permanent environmental variance; σe2 = residual variance; ha2 = heritability; HPD95 = lower and upper bounds of the 95% highest posterior density region.
      2 CHY = cheese percentage.
      Reaction norm model estimates of the AG and PE variance, with the increasing of THI, are shown in Figure 2, Figure 3, respectively.
      Figure thumbnail gr2
      Figure 2Plot of additive genetics (top) and permanent environmental variance (bottom) posterior mean estimation along temperature-humidity index (THI) values, for protein percentage (A), protein yield in kilograms per day (B), fat yield percentage (C), and fat yield in kilograms per day (D).
      Figure thumbnail gr3
      Figure 3Plot of additive genetics (top) and permanent environmental variance (bottom) along temperature-humidity index (THI) posterior mean estimation, for milk yield in kilograms per day (A), 4% FCM in kilograms per day (B), ECM in kilograms per day (C), and cheese yield percentage (D).
      The PP curve showed a decreasing trend from THI 70 toward extreme values (THI 90), similar to that recorded in FCM and CHY, although in the latter 2 traits, the decreasing trend was more evident. In FP, FY, MY, and ECM, the curve's trend drew a U shape, with minimum variance values reached at THI 80. Unlike in the PY, the curve had a slightly increasing trend up to THI 85 and then grew more rapidly up to THI 90.
      The PE graph showed the same bell trend in all the investigated traits, with the maximum variance value at THI 80 and slightly higher values at THI ranging from 70 to 75, compared with the extreme opposite, between THI 85 and 90. In PP, a different trend was observed: The curve rapidly dropped toward minimum values at THI 75 and then rose toward THI 80 at approximately 50% of the variance observed at THI 70, remaining almost linear up to THI 90.
      The residual variance in the 2 models was always higher than that of AG and PE, except for PP in the RN model.

      Heritability and Ratio of Animal PE Variance to Phenotypic Variance

      The trend of h2 and the PEE estimate for the RN model are shown in Figure 4, Figure 5, respectively. The PEE curve followed a similar trend to that of PE in Figure 2, Figure 3, and the curves described the trend of h2 in all traits. Generally, a change of the curve shape was observed at THI 80, when the h2 value increases toward extreme THI values. Only PP, FCM, and CHY decreased at THI values of 80 up to the highest ones.
      Figure thumbnail gr5
      Figure 5Plot of h2 (top) and permanent environment effect (bottom) along temperature-humidity index (THI) posterior mean estimation, for milk yield in kilograms per day (A), 4% FCM in kilograms per day (B), ECM in kilograms per day (C), and cheese yield percentage (D).
      In Table 3, the h2 trend in the MT model is overall superimposable to that observed in the RN model, with slight differences in the last quantile, at THI values ranging from 80 to 90. Notably, for PY, FP, and MY, higher h2 values were recorded in the RN model in the range of THI 88–90, compared with the fifth quintile of the MT model (THI 78–95). All the other traits showed h2 values similar to or slightly lower in the MT model than the RN model.

      Genetic and PE Correlations

      Supplemental Table S1 (https://figshare.com/articles/journal_contribution/Supplementary_Table_1_docx/21728531;
      • Landi V.
      • Maggiolino A.
      • Cecchinato A.
      • Mota L.
      • Bernabucci U.
      • Rossoni A.
      • De Palo P.
      Supplementary Table 1.docx. figshare. Journal contribution.
      ) shows the RG and RPE between linear and quadratic RR coefficients and the intercept of the RN model. The values were always negative, except for RG12, whereas RG13 was positive for PY, FP, MY, and EMC. RG23 was negative only in PP, MY, and CHY. RPE, on the contrary, showed negative values except for PP in RPE12. Table 4 shows the correlations RGiq and RPEiq between the 3 correlated traits in the MT model. In the case of RGiq, the values were all higher than 0.94. The values of the RGiq of the first versus the third quintile were, on average, 0.992 against a lower value observed between the first and fifth quintiles (0.959).
      On the contrary, intermediate values were observed between the second and fifth (0.971) quintile, and CHY traits recorded the lowest values (0.944 and 0.945, respectively), followed by FY and MY (0.951 and 0.953, respectively). For PY and FP, values of 0.965 and 0.961 were observed, respectively, and finally, PP and FCM were located toward the highest values (0.973 and 0.977, respectively). The RPEiq, indicating the relationship between the different quintiles of the PE, recorded lower values (average 0.959, 0.855, and 0.691 from quintiles 1 to 5, respectively), and the difference between RGiq and RPEiq in the same quintile increased considerably from the first to the fifth (0.033, 0.116, and 0.286, respectively). The RPEiq between the first and fifth quintiles was the lowest for FP (0.898) and the highest for PP (0.973). On the contrary, between the third and fifth quintile, PP was the lowest value (0.785) and FP was the highest (0.911). Similarly, between the second and third quintile, PP showed the lowest value (0.621) and FY the highest (0.726).
      The graphic representation of the nTHI × nTHI order matrix of the RGTT' and RPETT' correlations along the THI scale is shown in Figures 6, 7, 8, and 9. Results confirmed the outcomes of the MT model. The graphs showed extreme values overall lower than the MT model. Similar 3-dimensional surfaces are observed in all the investigated traits for the RPETT, which specifically showed a rapid drop of the correlation values starting from THI 70 up to the highest values. The concavity in the center of the surface suggests a decrease in the correlation values before the extreme values of THI (around THI 80–86). The lower values ranged from approximately ~83 (FCM, MY, PY, FY, and ECM) to ~80 for PP, and lower values, approximately 77–78, for CHY and FP. The plot of RGTT showed a continuous and linear decrease for all the traits toward extreme values of THI. PP and PY showed similar trends with a high correlation (~0.99), gradually decreasing to correlation values of ~0.94 and ~0.92, respectively. The ECM trait showed a lower correlation value with 0.87 at a THI ranging from 76 to 90.
      Figure thumbnail gr6
      Figure 6Tridimensional wireframe plot of estimated posterior means of the additive genetic and permanent environment correlation across different temperature-humidity index (THI) values for fat percentage and fat yield in kilograms per day.
      Figure thumbnail gr7
      Figure 7Tridimensional wireframe plot of estimated posterior means of the additive genetic and permanent environment correlation across different temperature-humidity index (THI) values for protein percentage and protein yield in kilograms per day.
      Figure thumbnail gr8
      Figure 8Tridimensional wireframe plot of estimated posterior means of the additive genetic and permanent environment correlation across different temperature-humidity index (THI) values for milk yield in kilograms per day and 4% FCM yield in kilograms per day.
      Figure thumbnail gr9
      Figure 9Tridimensional wireframe plot of estimated posterior means of the additive genetic and permanent environment correlation across different temperature-humidity index (THI) values for ECM yield in kilograms per day and cheese yield in percentage.

      Estimation of the EBVs Including the THI Effect

      Table 6 shows the ranking experienced by the first 50 sires; for brevity, only the PP and MY traits based on RN model results are shown. We chose 2 target points within the THI scale (THI = 76 and THI = 90) to represent moderate and severe heat stress, respectively. In PP 26 and 37, sires experienced a change in rank position at THI 76 and 90, respectively, while in MY 7 and 12 bulls, respectively. As expected, at THI 90, more sires experienced a change in rank position. At THI 90, the highest change observed in a bull was 5 positions in PP and 1 in MY. However, the PY trait showed that a more significant number of sires changed position in rank (43), and the widest position change for a single sire (21 positions).
      Table 6Bull ranking comparison based on reaction norm model results using protein percentage and milk yield
      Rank 1 = traditional breeding (e.g., intercept of reaction norm model) values rank; Rank 2 = ranking at temperature-humidity index (THI) 76; Rank 3 = ranking at THI 90; P2 = number of the position change experienced by the bull at THI 76; P3 = number of the position change experienced by the bull at THI 90.
      BullProtein (%)Milk yield (kg/d)
      Rank 1Rank 2P2Rank 3P3BullRank 1Rank 2P2Rank 3P3
      103,96111170,343111
      70,6072231104,054222
      87,353332176,750333
      44,1364451109,840444
      88,3315541104,686555
      48,794671665,470666
      65,387761771,711777
      87,84988871,712888
      83,05999988,348999
      53,884101011170,864101010
      29,4731112110159,264111111
      48,8051211113170,616121212
      110,227131312171,7181313141
      71,71414141459,4621414131
      87,176151516148,794151515
      54,093161615194,388161616
      35,2471718120359,271171717
      83,0791817121354,121181818
      44,317191918164,879191919
      54,034202019171,6532020211
      70,7762124317470,8522121201
      48,8002221125348,8062223122
      29,49023221221104,68323221241
      53,8792423123154,05024251231
      70,4832526124170,2662524125
      48,793262512688,3462626271
      39,475272728171,0342727261
      59,462282827165,114282828
      64,66329292983,0682929301
      77,914303113065,4733030291
      44,1433130133288,331313131
      87,005323231148,804323232
      98,637333334183,064333333
      48,792343438459,672343434
      70,25835353570,258353535
      64,6493638232444,000363636
      71,0343736136148,795373737
      48,9593837140276,756383838
      43,99939393970,483393939
      87,488404042254,105404040
      70,8524145437471,709414141
      65,473424241177,9234243142
      110,2214341246387,8494342143
      60,6764443143170,7764444451
      35,4514546144139,2954545441
      35,4844648245144,009464646
      39,369474924754,110474747
      32,234484714854,093484848
      54,1184944554565,387494949
      48,8085051149159,6735051150
      1 Rank 1 = traditional breeding (e.g., intercept of reaction norm model) values rank; Rank 2 = ranking at temperature-humidity index (THI) 76; Rank 3 = ranking at THI 90; P2 = number of the position change experienced by the bull at THI 76; P3 = number of the position change experienced by the bull at THI 90.
      In Figure 10, Figure 11, bulls' EBVs belonging to the first (bottom-ranked sires), third (medium-ranked sires), and fifth quintile (top-ranked sires) of the initial selection is displayed, resulting in ~18 bulls for class. A similar pattern is repeated for all traits where the bulls with an average EBV (third quintile) are clustered in a more homogeneous group. On the contrary, the top sires (fifth quintile) are more heterogeneous and dispersed over a greater range of values.
      Figure thumbnail gr10
      Figure 10EBVs for protein yield in percentage (A) and kilograms per day (B), fat yield in percentage (C) and kilograms per day (D), along with the value of temperature-humidity index (THI) of first (bottom), third (medium), and fifth (top) standard (intercept of reaction norm model) EBV quintile of sires with at least 90 daughters. Black dotted lines are the within-quintile average trend.
      Figure thumbnail gr11
      Figure 11EBVs for milk yield in kilograms per day (A), 4% FCM yield in kilograms per day (B), ECM yield in kilograms per day (C) and cheese yield in percentage (D) along with the value of temperature-humidity index (THI) of first (bottom), third (medium), and fifth (top) standard (intercept of reaction norm model) EBV quintile of sires with at least 90 daughters. Black dotted lines are the within-quintile average trend.
      The EBVs for production traits and the studied temperature range generally declined at high temperatures for the top-ranked animals but with slight differences. The top-rank sires are most heterogeneous in the PY trait, with several bulls increasing their EBVs toward extreme THI values. All top-rank bulls tended to decline for FY and MY traits until THI 85, but the trend flexed toward higher values. For FCM and ECM, the most evident decreasing trend of top-rank bulls has been recorded.

      DISCUSSION

      Climatic Variables and Modeling

      To the best of our knowledge, the study of genetic components for resistance to heat stress in Brown Swiss cattle has been estimated for the first time. Previous studies showed that the maximum THI levels and the number of continuous days over the threshold before the test day negatively affect production (
      • Maggiolino A.
      • Dahl G.E.
      • Bartolomeo N.
      • Bernabucci U.
      • Vitali A.
      • Serio G.
      • Cassandro M.
      • Centoducati G.
      • Santus E.
      • De Palo P.
      Estimation of maximum thermo-hygrometric index thresholds affecting milk production in Italian Brown Swiss cattle.
      ,
      • Maggiolino A.
      • Landi V.
      • Bartolomeo N.
      • Bernabucci U.
      • Santus E.
      • Bragaglio A.
      • De Palo P.
      Effect of heat waves on some Italian Brown Swiss dairy cows' production patterns.
      ). So, upper critical THI thresholds are higher than those reported in Italian Holstein Friesians (
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effects of heat stress in Italian Holstein dairy cattle.
      ). The limited literature on this specific breed reported that the Brown Swiss breed is more rustic and resistant to environmental conditions (
      • El-Tarabany M.S.
      • Roushdy E.M.
      • El-Tarabany A.A.
      Production and health performance of Holstein, Brown Swiss and their crosses under subtropical environmental conditions.
      ).
      Many studies identified and reported THI thresholds beyond which production traits worsen, and these values have been used to characterize and model the THI effect. The heat stress threshold has been established worldwide at 72 units of THI (
      • NRC
      A Guide to Environmental Research on Animals.
      ,
      • Collier R.
      • Hall L.
      • Rungruang S.
      • Zimbelman R.
      Quantifying heat stress and its impact on metabolism and performance.
      ). This threshold is consistent with the values reported by
      • Ravagnolo O.
      • Misztal I.
      Genetic component of heat stress in dairy cattle, parameter estimation.
      in US Holsteins.
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effects of heat stress in Italian Holstein dairy cattle.
      , using test-day information from Italian Holsteins, calculated a different threshold by parity class and production trait, with a mean value of THI 71. In Italian Brown Swiss,
      • Maggiolino A.
      • Dahl G.E.
      • Bartolomeo N.
      • Bernabucci U.
      • Vitali A.
      • Serio G.
      • Cassandro M.
      • Centoducati G.
      • Santus E.
      • De Palo P.
      Estimation of maximum thermo-hygrometric index thresholds affecting milk production in Italian Brown Swiss cattle.
      found a mean threshold of THI 75, with a minimum of THI 64 (fat kg/d) and a maximum of THI 75 (FCM at 4%). Based on this previous information, we considered a unique upper-thermal comfort threshold for all the production traits of THI 70 to evaluate the variability of this parameter according to each production trait and cow factors.
      In G × E interaction assessment, a classical approach splits a trait into several correlated traits depending on the physical environment or farm characteristics (
      • Windig J.J.
      • Calus M.P.L.
      • Beerda B.
      • Veerkamp R.F.
      Genetic correlations between milk production and health and fertility depending on herd environment.
      ;
      • Martinez-Castillero M.
      • Toledo-Alvarado H.
      • Pegolo S.
      • Vazquez A.I.
      • de Los Campos G.
      • Varona L.
      • Finocchiaro R.
      • Bittante G.
      • Cecchinato A.
      Genetic parameters for fertility traits assessed in herds divergent in milk energy output in Holstein-Friesian, Brown Swiss, and Simmental cattle.
      ). The MT model used in this study, similar to what was reported by
      • 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.
      , exploits this strategy. In this approach lies the disadvantage that it is challenging to obtain enough observations in small temperature ranges in real conditions.
      Previous studies that used random regression models for similar investigations included the environmental variable as a function of the degree of stress of the animal, starting from a threshold of thermoneutrality, including it as a fixed coefficient or as random regression covariate (
      • Aguilar I.
      • Misztal I.
      • Tsuruta S.
      Genetic components of heat stress for dairy cattle with multiple lactations.
      ;
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effects of heat stress in Italian Holstein dairy cattle.
      ;
      • Negri R.
      • Aguilar I.
      • Feltes G.L.
      • Cobuci J.A.
      Selection for test-day milk yield and thermotolerance in Brazilian Holstein cattle.
      ). Some authors introduced more than a unique threshold for traits and calving class (
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effects of heat stress in Italian Holstein dairy cattle.
      ). Other studies highlight significant differences between Italian Brown Swiss and Holstein Friesians (
      • Maggiolino A.
      • Dahl G.E.
      • Bartolomeo N.
      • Bernabucci U.
      • Vitali A.
      • Serio G.
      • Cassandro M.
      • Centoducati G.
      • Santus E.
      • De Palo P.
      Estimation of maximum thermo-hygrometric index thresholds affecting milk production in Italian Brown Swiss cattle.
      ).
      • Carabaño M.J.
      • Logar B.
      • Bormann J.
      • Minet J.
      • Vanrobays M.L.
      • Díaz C.
      • Tychon B.
      • Gengler N.
      • Hammami H.
      Modeling heat stress under different environmental conditions.
      reported the difficulty of using the thresholds, given the impossibility of accurately describing the performance of each animal in interaction with the environment. Moreover, inconsistent results in some traits, such as FP and MY in kilograms, have been previously described (
      • 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.
      ;
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effects of heat stress in Italian Holstein dairy cattle.
      ;
      • Maggiolino A.
      • Dahl G.E.
      • Bartolomeo N.
      • Bernabucci U.
      • Vitali A.
      • Serio G.
      • Cassandro M.
      • Centoducati G.
      • Santus E.
      • De Palo P.
      Estimation of maximum thermo-hygrometric index thresholds affecting milk production in Italian Brown Swiss cattle.
      ). Polynomials covariate as LP allow more flexible modeling of the phenotypes trend according to THI changes, where multiple breaking points should be assumed, and the passage from the comfort zone to the stress zone can be gradually highlighted (
      • Carabaño M.J.
      • Logar B.
      • Bormann J.
      • Minet J.
      • Vanrobays M.L.
      • Díaz C.
      • Tychon B.
      • Gengler N.
      • Hammami H.
      Modeling heat stress under different environmental conditions.
      ).
      The RN model used the second-order LP of the mean values of THI (from the test day to 5 d before it) as a random regression coefficient on the animal additive effect. In many articles, an LP of order 3 (
      • 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.
      ;
      • Carabaño M.J.
      • Logar B.
      • Bormann J.
      • Minet J.
      • Vanrobays M.L.
      • Díaz C.
      • Tychon B.
      • Gengler N.
      • Hammami H.
      Modeling heat stress under different environmental conditions.
      ) or of order 4 (
      • Gernand E.
      • König S.
      • Kipp C.
      Influence of on-farm measurements for heat stress indicators on dairy cow productivity, female fertility, and health.
      ) has been used, whereas
      • Negri R.
      • Aguilar I.
      • Feltes G.L.
      • Cobuci J.A.
      Selection for test-day milk yield and thermotolerance in Brazilian Holstein cattle.
      , in a model comparison study, reported how the model with a coefficient of order 2 is the one with the best fitting parameters. In Polish Holstein Friesians, higher order coefficient models allowed a better fitting (
      • Strabel T.
      • Szyda J.
      • Ptak E.
      • Jamrozik J.
      Comparison of random regression test-day models for Polish Black and White cattle.
      ). Some authors used the first-order coefficient with a thermoneutrality zone established to THI 60 (
      • Nguyen T.T.T.
      • Bowman P.J.
      • Haile-Mariam M.
      • Pryce J.E.
      • Hayes B.J.
      Genomic selection for tolerance to heat stress in Australian dairy cattle.
      ;
      • 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.
      ). Negative implications have also been reported in the use of higher order mathematical models of LP because the AG and PE estimates can take on a more erratic and oscillatory trend and lead to an overestimation of the value, especially at extreme THI values, which, in addition, are even those with the fewest records and therefore more susceptible to erroneous estimates (
      • Pool M.H.
      • Meuwissen T.H.E.
      Reduction of the number of parameters needed for a polynomial random regression test day model.
      ;
      • López-Romero P.
      • Carabaño M.J.
      Comparing alternative random regression models to analyse first lactation daily milk yield data in Holstein–Friesian cattle.
      ;
      • Strabel T.
      • Szyda J.
      • Ptak E.
      • Jamrozik J.J.I.B.
      Comparison of random regression test-day models for production traits of dairy cattle in Poland.
      ).

      Variance Components, h2, and Ratio of Permanent Environmental Variance to Phenotypic Variance

      Except for the intercept, which indicates the production level when the covariate is null, the biological interpretation of each regression coefficient in polynomials of greater than linear order is not easy. The linear and quadratic coefficients show how the form of the individual response pattern changes as the thermal load varies. In our case, linear and quadratic coefficients show a decrease in variance value.
      Many authors reported a gradual and slight increase of AG for production traits at THI values higher than 80; we showed similar results for the RN model, only for PY, FP, FY, MY, and ECM, with a U-shaped AG variance pattern for production traits due to heat stress (
      • Ravagnolo O.
      • Misztal I.
      Genetic component of heat stress in dairy cattle, parameter estimation.
      ;
      • Brügemann K.
      • Gernand E.
      • von Borstel U.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.
      ;
      • 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.
      ). The AG variance for all the milk traits from the RN model was not stable at the different THI and traits (Figure 2, Figure 3), suggesting that a different response to selection is expected regardless of the THI environment.
      In the RN model, results of PE variance for MY, PY, and FY in Brown Swiss were more significant than compared with AG variance; other studies on Holsteins confirmed this evidence (
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effects of heat stress in Italian Holstein dairy cattle.
      ;
      • 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.
      ). In the range of THI 60–75, the shape of the PE curve in the Brown Swiss breed was similar to that observed by other authors (
      • 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.
      ). A decreasing trend was observed for PE variances at increasing values of THI or increasing heat stress, which agrees with other authors (
      • Brügemann K.
      • Gernand E.
      • von Borstel U.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.
      ).
      Higher PEE values were found at THI ranging from 75 to 85, suggesting that the animals' performances are strongly influenced by nongenetic factors (
      • Martinez-Castillero M.
      • Toledo-Alvarado H.
      • Pegolo S.
      • Vazquez A.I.
      • de Los Campos G.
      • Varona L.
      • Finocchiaro R.
      • Bittante G.
      • Cecchinato A.
      Genetic parameters for fertility traits assessed in herds divergent in milk energy output in Holstein-Friesian, Brown Swiss, and Simmental cattle.
      ). Larger PEE compared with h2 at THI 75–85 suggests that heat stress sensitivity is cow-specific and mainly acquired rather than genetic (
      • Boonkum W.
      • Misztal I.
      • Duangjinda M.
      • Pattarajinda V.
      • Tumwasorn S.
      • Sanpote J.
      Genetic effects of heat stress on milk yield of Thai Holstein crossbreds.
      ).
      The increasing values of h2 at the warmest environmental variable values were also observed by
      • Carabaño M.J.
      • Bachagha K.
      • Ramón M.
      • Díaz C.
      Modeling heat stress effect on Holstein cows under hot and dry conditions: Selection tools.
      ,
      • Carabaño M.J.
      • Pineda-Quiroga C.
      • Ugarte E.
      • Díaz C.
      • Ramón M.
      Genetic basis of thermotolerance in 2 local dairy sheep populations in the Iberian Peninsula.
      ), studying Spanish Holstein Friesians and some Spanish sheep breeds. PY, FP, FY, MY, FCM, and ECM quadratic polynomials tended to show a steep increase in h2 estimates at the extreme values of the considered independent variable (THI) due to artifacts of polynomial functions.
      • Brügemann K.
      • Gernand E.
      • von Borstel U.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.
      detected a trend of higher h2 at lower THI values, as recorded in the present study, except for PY, where the trend was the opposite.
      Unlike other studies (
      • Carabaño M.J.
      • Bachagha K.
      • Ramón M.
      • Díaz C.
      Modeling heat stress effect on Holstein cows under hot and dry conditions: Selection tools.
      ,
      • Carabaño M.J.
      • Pineda-Quiroga C.
      • Ugarte E.
      • Díaz C.
      • Ramón M.
      Genetic basis of thermotolerance in 2 local dairy sheep populations in the Iberian Peninsula.
      ;
      • 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.
      ), PY showed an almost sinusoidal trend (Figure 4B), where the estimates increased slightly until THI 85 and then more quickly until THI 90. The shape of the h2 curve found in many studies is variable, and it could depend on the different regression coefficient structures used and the strategies researchers adopted in setting values considered as the thermoneutrality zone. In Australian Holstein cattle,
      • 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.
      , using a first-order LP coefficient, a U-shaped curve was obtained, setting the neutrality limit at THI 60 and forcing all extreme (until THI 80) environmental data to 74.
      Figure thumbnail gr4
      Figure 4Plot of h2 (top) and permanent environment effect (bottom) along temperature-humidity index (THI) posterior mean estimation, for protein percentage (A), protein yield in kilograms per day (B), fat percentage (C), and fat yield in kilograms per day (D).

      Genetic and PE Correlation

      Our findings revealed that thermal tolerance varies significantly at the THI gradient's extremities (THI 85–90). Furthermore, we differentiated the results of production trait variation slopes along the THI gradient under heat stress (THI >70) based on the level of traditional EBVs value quintiles. Several studies on cattle (
      • Ravagnolo O.
      • Misztal I.
      Genetic component of heat stress in dairy cattle, parameter estimation.
      ;
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effects of heat stress in Italian Holstein dairy cattle.
      ;
      • Carabaño M.J.
      • Bachagha K.
      • Ramón M.
      • Díaz C.
      Modeling heat stress effect on Holstein cows under hot and dry conditions: Selection tools.
      ) and sheep (
      • Finocchiaro R.
      • van Kaam J.B.
      • Portolano B.
      • Misztal I.
      Effect of heat stress on production of Mediterranean dairy sheep.
      ;
      • Carabaño M.J.
      • Pineda-Quiroga C.
      • Ugarte E.
      • Díaz C.
      • Ramón M.
      Genetic basis of thermotolerance in 2 local dairy sheep populations in the Iberian Peninsula.
      ) showed an antagonistic relationship between high production level and heat tolerance. Antagonism between production and fitness attributes regularly makes long-term productivity improvement difficult (
      • Rauw W.M.
      • Kanis E.
      • Noordhuizen-Stassen E.N.
      • Grommers F.J.
      Undesirable side effects of selection for high production efficiency in farm animals: A review.
      ). Selection toward higher production will lead to a deterioration of the adaptation capacity of the animals to heat stress events. As the breeding scheme aimed mainly to improve productive traits, we had remarkable progress in productivity of the populations, and several undesirable correlated genetic responses for thermotolerance. The linear regression coefficient and intercept correlations were always negative, except for the MY and PY.
      • Maggiolino A.
      • Dahl G.E.
      • Bartolomeo N.
      • Bernabucci U.
      • Vitali A.
      • Serio G.
      • Cassandro M.
      • Centoducati G.
      • Santus E.
      • De Palo P.
      Estimation of maximum thermo-hygrometric index thresholds affecting milk production in Italian Brown Swiss cattle.
      recorded that MY is not affected by heat stress in the Italian Brown Swiss population.
      G × E interaction is directly indicated by a nonunit genetic correlation between extreme THI class values (
      • Carabaño M.J.
      • Pineda-Quiroga C.
      • Ugarte E.
      • Díaz C.
      • Ramón M.
      Genetic basis of thermotolerance in 2 local dairy sheep populations in the Iberian Peninsula.
      ).
      • Kolmodin R.
      • Bijma P.
      Response to mass selection when the genotype by environment interaction is modelled as a linear reaction norm.
      reported that a negative genetic correlation between intercept and random slope of an RN model results from selection for high yield in a nonlimiting environment, which leads to increased environmental sensitivity. Relevant G × E interaction in terms of genetic correlation is intended when the value is lower than 0.8 (
      • Robertson A.
      The sampling variance of the genetic correlation coefficient.
      ). Genetic productivity correlations under different thermal load regions were always above the 0.9 thresholds in our study.
      The evolution of Brown Swiss to different environmental conditions might explain the apparent better adaptation to extreme weather. The ECM and FY traits showed a lower value in agreement with results reported by
      • 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.
      .
      In the MT model, the correlations between different quintiles were consistently high (>0.90), which is in line with what was observed in studies that used a similar procedure (
      • Hayes B.J.
      • Carrick M.
      • Bowman P.
      • Goddard M.E.
      Genotype × environment interaction for milk production of daughters of Australian dairy sires from test-day records.
      ;
      • 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.
      ). However,
      • 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.
      recorded lower values between the 5th and 95th percentile, probably also due to the larger data set size, which allowed the authors to narrow the range of THI used to build the separate environments.
      Although the genetic correlation is strong, implying modest G × E interactions, the breeding value index may be used to select milk traits and heat tolerance (
      • Ravagnolo O.
      • Misztal I.
      • Hoogenboom G.
      Genetic component of heat stress in dairy cattle, development of heat index function.
      ). Furthermore, the PE association revealed that controlling the environment for high milk production cows will enhance the cows' unfavorable long-term environmental consequences of higher heat stress (
      • Boonkum W.
      • Duangjinda M.
      Estimation of genetic parameters for heat stress, including dominance gene effects, on milk yield in Thai Holstein dairy cattle.
      ). According to
      • Aguilar I.
      • Misztal I.
      • Tsuruta S.
      Genetic components of heat stress for dairy cattle with multiple lactations.
      , cows' responses to heat stress are mainly influenced by the environment rather than genetics. Moreover, evidence of the carryover effect of heat stress during a cow's life, including the in utero stage, has been demonstrated both in experimental and field conditions (
      • Laporta J.
      • Ferreira F.C.
      • Ouellet V.
      • Dado-Senn B.
      • Almeida A.K.
      • De Vries A.
      • Dahl G.E.
      Late-gestation heat stress impairs daughter and granddaughter lifetime performance.
      ;
      • Kipp C.
      • Brugemann K.
      • Yin T.
      • Halli K.
      • Konig S.
      Genotype by heat stress interactions for production and functional traits in dairy cows from an across-generation perspective.
      ).

      EBVs Change Implications

      Resilient animals are minimally affected by environmental disturbances (
      • Colditz I.G.
      • Hine B.C.
      Resilience in farm animals: Biology, management, breeding and implications for animal welfare.
      ;
      • Berghof T.V.L.
      • Poppe M.
      • Mulder H.A.
      Opportunities to improve resilience in animal breeding programs.
      ), and their performances are expected to be comparatively consistent under different environmental conditions. Our results suggested a proportion of resilient sires among the first and third quantile of traditional EBVs. However, sires with decreasing performances are appreciable in the fifth quintile. Top-ranked sires, if sensitive, will generate daughters with the worst performance under heat stress conditions. This behavior was described previously in several works, e.g.,
      • Carabaño M.J.
      • Bachagha K.
      • Ramón M.
      • Díaz C.
      Modeling heat stress effect on Holstein cows under hot and dry conditions: Selection tools.
      found a higher decrease in performance in top-ranked sires in Spanish Holsteins. In the Holstein Friesian breed, some authors reported that the antagonism of high productive performance under heat stress might be due to high metabolic activity implied by high production levels (
      • Bernabucci U.
      • Biffani S.
      • Buggiotti L.
      • Vitali A.
      • Lacetera N.
      • Nardone A.
      The effects of heat stress in Italian Holstein dairy cattle.
      ). Some authors suggested that heat stress in the Holstein breed could increase metabolic heat (
      • Carabaño M.J.
      • Bachagha K.
      • Ramón M.
      • Díaz C.
      Modeling heat stress effect on Holstein cows under hot and dry conditions: Selection tools.
      ). Evidence that the importance of the ability to dissipate heat is also indirectly a function of body size would partially justify the greater resilience of Brown Swiss (
      • Maggiolino A.
      • Dahl G.E.
      • Bartolomeo N.
      • Bernabucci U.
      • Vitali A.
      • Serio G.
      • Cassandro M.
      • Centoducati G.
      • Santus E.
      • De Palo P.
      Estimation of maximum thermo-hygrometric index thresholds affecting milk production in Italian Brown Swiss cattle.
      ). These results could be useful for breeding decisions related to bull choice, even if the extent of re-ranking for most of the sires and G × E magnitude in the Brown Swiss population may not justify separate genetic evaluations for high heat stress environments.
      On the one hand, in sire daughters in which heat stress sensitivity is present, a more controlled environment is necessary, such as one with cooling systems and suitable feeding strategies for lowering the core body temperature. On the other hand, an alternative solution may be the unique genetic evaluation and ranking of sires using environmental variables, suggesting to breeders that sires perform better in regions where the climate conditions can often exceed the thermal comfort thresholds (
      • Kolmodin R.
      • Strandberg E.
      • Madsen P.
      • Jensen J.
      • Jorjani H.
      genotype by environment interaction in Nordic dairy cattle studied using reaction norms.
      ;
      • Bryant J.R.
      • Lopez-Villalobos N.
      • Pryce J.E.
      • Holmes C.W.
      • Johnson D.L.
      • Garrick D.J.
      Environmental sensitivity in New Zealand dairy cattle.
      ;
      • Nguyen T.T.T.
      • Bowman P.J.
      • Haile-Mariam M.
      • Pryce J.E.
      • Hayes B.J.
      Genomic selection for tolerance to heat stress in Australian dairy cattle.
      ).
      • 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.
      suggested using the presence of heat stress tolerance to make the breeding program more flexible. However, a more extensive use should consider other phenotypes related to fertility and health.

      CONCLUSIONS

      The results shown in the present paper revealed a moderate G × E effect in Italian Brown Swiss breed cows. We observed changes in ranking position rather than a general trend in bulls' EBVs. However, several sires in the top position for productive performance of the daughters could experience a change in rank position in heat stress environments. Moreover, re-ranking many sires may not justify separate genetic evaluations for high heat stress environments. On the contrary, including environmental data in genetic evaluations could allow breeders to choose resilient sires in specific environmental areas where it is more likely to exceed the thermal comfort thresholds during the year, or where the conditions of high temperature and humidity have a longer duration.

      ACKNOWLEDGMENTS

      This study received no external funding. The authors are grateful to Stefano Biffani (CNR-IBBA, Milano) and Mayra Gomez Carpio (Italian Buffalo Breeders Association, Caserta, Italy) for their technical support. Moreover, the authors are grateful to Marco Giazzi and Marco Tadini, the president and vice president, respectively, of the Meteonetwork association (www.meteonetwork.it), and to the Italian Air Force weather service (www.aeronautica.difesa.it) for giving us access to all weather station data. The authors have not stated any conflicts of interest.

      REFERENCES

        • Aguilar I.
        • Misztal I.
        • Tsuruta S.
        Genetic components of heat stress for dairy cattle with multiple lactations.
        J. Dairy Sci. 2009; 92 (19841230): 5702-5711
        • Aguilar I.
        • Misztal I.
        • Tsuruta S.
        Short communication: Genetic trends of milk yield under heat stress for US Holsteins.
        J. Dairy Sci. 2010; 93 (20338455): 1754-1758
        • 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
        • Berghof T.V.L.
        • Poppe M.
        • Mulder H.A.
        Opportunities to improve resilience in animal breeding programs.
        Front Genet. 2019; 9: 692
        • 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
        • Bohmanova J.
        • Misztal I.
        • Cole J.B.
        Temperature-humidity indices as indicators of milk production losses due to heat stress.
        J. Dairy Sci. 2007; 90 (17369235): 1947-1956
        • Boonkum W.
        • Duangjinda M.
        Estimation of genetic parameters for heat stress, including dominance gene effects, on milk yield in Thai Holstein dairy cattle.
        Anim. Sci. J. 2015; 86 (25226870): 245-250
        • Boonkum W.
        • Misztal I.
        • Duangjinda M.
        • Pattarajinda V.
        • Tumwasorn S.
        • Sanpote J.
        Genetic effects of heat stress on milk yield of Thai Holstein crossbreds.
        J. Dairy Sci. 2011; 94 (21183060): 487-492
        • Brügemann K.
        • Gernand E.
        • von Borstel U.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
        • Bryant J.R.
        • Lopez-Villalobos N.
        • Pryce J.E.
        • Holmes C.W.
        • Johnson D.L.
        • Garrick D.J.
        Environmental sensitivity in New Zealand dairy cattle.
        J. Dairy Sci. 2007; 90 (17297127): 1538-1547
        • Carabaño M.J.
        • Bachagha K.
        • Ramón M.
        • Díaz C.
        Modeling heat stress effect on Holstein cows under hot and dry conditions: Selection tools.
        J. Dairy Sci. 2014; 97 (25262182): 7889-7904
        • Carabaño M.J.
        • Logar B.
        • Bormann J.
        • Minet J.
        • Vanrobays M.L.
        • Díaz C.
        • Tychon B.
        • Gengler N.
        • Hammami H.
        Modeling heat stress under different environmental conditions.
        J. Dairy Sci. 2016; 99 (26923054): 3798-3814
        • Carabaño M.J.
        • Pineda-Quiroga C.
        • Ugarte E.
        • Díaz C.
        • Ramón M.
        Genetic basis of thermotolerance in 2 local dairy sheep populations in the Iberian Peninsula.
        J. Dairy Sci. 2021; 104 (33612212): 5755-5767
        • 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
        • Colditz I.G.
        • Hine B.C.
        Resilience in farm animals: Biology, management, breeding and implications for animal welfare.
        Anim. Prod. Sci. 2016; 56: 1961-1983
        • Collier R.
        • Hall L.
        • Rungruang S.
        • Zimbelman R.
        Quantifying heat stress and its impact on metabolism and performance.
        in: Proc. Florida Ruminant Nutrition Symp. Univ. of Florida, Gainesville, FL2012: 74-83
        • De Rensis F.
        • Garcia-Ispierto I.
        • López-Gatius F.
        Seasonal heat stress: Clinical implications and hormone treatments for the fertility of dairy cows.
        Theriogenology. 2015; 84 (26025242): 659-666
        • De Vries A.
        • Marcondes M.I.
        Review: Overview of factors affecting productive lifespan of dairy cows.
        Animal. 2020; 14 (32024570): s155-s164
        • Dillon P.
        • Buckley F.
        • O'Connor P.
        • Hegarty D.
        • Rath M.
        A comparison of different dairy cow breeds on a seasonal grass-based system of milk production: 1. Milk production, live weight, body condition score and DM intake.
        Livest. Prod. Sci. 2003; 83: 21-33
        • El-Tarabany M.S.
        • Roushdy E.M.
        • El-Tarabany A.A.
        Production and health performance of Holstein, Brown Swiss and their crosses under subtropical environmental conditions.
        Anim. Prod. Sci. 2017; 57: 1137-1143
        • Finocchiaro R.
        • van Kaam J.B.
        • Portolano B.
        • Misztal I.
        Effect of heat stress on production of Mediterranean dairy sheep.
        J. Dairy Sci. 2005; 88 (15829679): 1855-1864
        • Formaggioni P.
        • Summer A.
        • Malacarne M.
        • Franceschi P.
        • Mucchetti G.
        Italian and Italian-style hard cooked cheeses: Predictive formulas for Parmigiano-Reggiano 24-h cheese yield.
        Int. Dairy J. 2015; 51: 52-58
        • Gaines W.L.
        • Davidson F.A.
        Relation between Percentage Fat Content and Yield of Milk: Correction of Milk Yield for Fat Content.
        University of Illinois Agricultural Experiment Station, 1923
        • Gernand E.
        • König S.
        • Kipp C.
        Influence of on-farm measurements for heat stress indicators on dairy cow productivity, female fertility, and health.
        J. Dairy Sci. 2019; 102 (31128870): 6660-6671
        • 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.
        • 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.
        • Carrick M.
        • Bowman P.
        • Goddard M.E.
        Genotype × environment interaction for milk production of daughters of Australian dairy sires from test-day records.
        J. Dairy Sci. 2003; 86 (14672205): 3736-3744
        • 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
        • Kipp C.
        • Brugemann K.
        • Yin T.
        • Halli K.
        • Konig S.
        Genotype by heat stress interactions for production and functional traits in dairy cows from an across-generation perspective.
        J. Dairy Sci. 2021; 104 (34099290): 10029-10039
        • Kolmodin R.
        • Bijma P.
        Response to mass selection when the genotype by environment interaction is modelled as a linear reaction norm.
        Genet. Sel. Evol. 2004; 36 (15231233): 435-454
        • Kolmodin R.
        • Strandberg E.
        • Madsen P.
        • Jensen J.
        • Jorjani H.
        genotype by environment interaction in Nordic dairy cattle studied using reaction norms.
        Acta Agric. Scand. A Anim. Sci. 2002; 52: 11-24
        • Landi V.
        • Maggiolino A.
        • Cecchinato A.
        • Mota L.
        • Bernabucci U.
        • Rossoni A.
        • De Palo P.
        Supplementary Table 1.docx. figshare. Journal contribution.
        • Laporta J.
        • Ferreira F.C.
        • Ouellet V.
        • Dado-Senn B.
        • Almeida A.K.
        • De Vries A.
        • Dahl G.E.
        Late-gestation heat stress impairs daughter and granddaughter lifetime performance.
        J. Dairy Sci. 2020; 103 (32534930): 7555-7568
        • López-Romero P.
        • Carabaño M.J.
        Comparing alternative random regression models to analyse first lactation daily milk yield data in Holstein–Friesian cattle.
        Livest. Prod. Sci. 2003; 82: 81-96
        • Maggiolino A.
        • Dahl G.E.
        • Bartolomeo N.
        • Bernabucci U.
        • Vitali A.
        • Serio G.
        • Cassandro M.
        • Centoducati G.
        • Santus E.
        • De Palo P.
        Estimation of maximum thermo-hygrometric index thresholds affecting milk production in Italian Brown Swiss cattle.
        J. Dairy Sci. 2020; 103 (32684476): 8541-8553
        • Maggiolino A.
        • Landi V.
        • Bartolomeo N.
        • Bernabucci U.
        • Santus E.
        • Bragaglio A.
        • De Palo P.
        Effect of heat waves on some Italian Brown Swiss dairy cows' production patterns.
        Front Anim Sci. 2022; 2800680
        • Martin J.G.A.
        • Nussey D.H.
        • Wilson A.J.
        • Réale D.
        Measuring individual differences in reaction norms in field and experimental studies: A power analysis of random regression models.
        Methods Ecol Evol. 2011; 2: 362-374
        • Martinez-Castillero M.
        • Toledo-Alvarado H.
        • Pegolo S.
        • Vazquez A.I.
        • de Los Campos G.
        • Varona L.
        • Finocchiaro R.
        • Bittante G.
        • Cecchinato A.
        Genetic parameters for fertility traits assessed in herds divergent in milk energy output in Holstein-Friesian, Brown Swiss, and Simmental cattle.
        J. Dairy Sci. 2020; 103 (33222858): 11545-11558
        • Misztal I.
        • Strabel T.
        • Jamrozik J.
        • Mantysaari E.A.
        • Meuwissen T.H.
        Strategies for estimating the parameters needed for different test-day models.
        J. Dairy Sci. 2000; 83 (10821589): 1125-1134
        • Negri R.
        • Aguilar I.
        • Feltes G.L.
        • Machado J.D.
        • Braccini Neto J.
        • Costa-Maia F.M.
        • Cobuci J.A.
        Inclusion of bioclimatic variables in genetic evaluations of dairy cattle.
        Anim. Biosci. 2021; 34 (32777914): 163-171
        • Negri R.
        • Aguilar I.
        • Feltes G.L.
        • Cobuci J.A.
        Selection for test-day milk yield and thermotolerance in Brazilian Holstein cattle.
        Animals (Basel). 2021; 11 (33430092): 128
        • Nguyen T.T.T.
        • Bowman P.J.
        • Haile-Mariam M.
        • Pryce J.E.
        • Hayes B.J.
        Genomic selection for tolerance to heat stress in Australian dairy cattle.
        J. Dairy Sci. 2016; 99 (27037467): 2849-2862
        • NRC
        A Guide to Environmental Research on Animals.
        Vol. 12. National Academy of Sciences, 1971
        • Pool M.H.
        • Meuwissen T.H.E.
        Reduction of the number of parameters needed for a polynomial random regression test day model.
        Livest. Prod. Sci. 2000; 64: 133-145
        • R Core Team
        R: A Language and Environment for Statistical Computing.
        R Foundation for Statistical Computing, 2019
        • Rauw W.M.
        • Kanis E.
        • Noordhuizen-Stassen E.N.
        • Grommers F.J.
        Undesirable side effects of selection for high production efficiency in farm animals: A review.
        Livest. Prod. Sci. 1998; 56: 15-33
        • 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.
        • Hoogenboom G.
        Genetic component of heat stress in dairy cattle, development of heat index function.
        J. Dairy Sci. 2000; 83 (11003246): 2120-2125
        • Renaudeau D.
        • Collin A.
        • Yahav S.
        • de Basilio V.
        • Gourdine J.L.
        • Collier R.J.
        Adaptation to hot climate and strategies to alleviate heat stress in livestock production.
        Animal. 2012; 6 (22558920): 707-728
        • Robertson A.
        The sampling variance of the genetic correlation coefficient.
        Biometrics. 1959; 15: 469-485
        • Rovelli G.
        • Ceccobelli S.
        • Perini F.
        • Demir E.
        • Mastrangelo S.
        • Conte G.
        • Abeni F.
        • Marletta D.
        • Ciampolini R.
        • Cassandro M.
        • Bernabucci U.
        • Lasagna E.
        The genetics of phenotypic plasticity in livestock in the era of climate change: A review.
        Ital. J. Anim. Sci. 2020; 19: 997-1014
        • Santana Jr., M.L.
        • Bignardi A.B.
        • Pereira R.J.
        • Menéndez-Buxadera A.
        • El Faro L.
        Random regression models to account for the effect of genotype by environment interaction due to heat stress on the milk yield of Holstein cows under tropical conditions.
        J. Appl. Genet. 2016; 57 (26155774): 119-127
        • Santana Jr., M.L.
        • Pereira R.J.
        • Bignardi A.B.
        • Filho A.E.
        • Menendez-Buxadera A.
        • El Faro L.
        Detrimental effect of selection for milk yield on genetic tolerance to heat stress in purebred Zebu cattle: Genetic parameters and trends.
        J. Dairy Sci. 2015; 98 (26476953): 9035-9043
        • Shi R.
        • Brito L.F.
        • Liu A.
        • Luo H.
        • Chen Z.
        • Liu L.
        • Guo G.
        • Mulder H.
        • Ducro B.
        • van der Linden A.
        • Wang Y.
        Genotype-by-environment interaction in Holstein heifer fertility traits using single-step genomic reaction norm models.
        BMC Genomics. 2021; 22 (33731012): 193
        • Smith D.L.
        • Smith T.
        • Rude B.J.
        • Ward S.H.
        Short communication: Comparison of the effects of heat stress on milk and component yields and somatic cell score in Holstein and Jersey cows.
        J. Dairy Sci. 2013; 96 (23498016): 3028-3033
        • Stocco G.
        • Cipolat-Gotet C.
        • Gasparotto V.
        • Cecchinato A.
        • Bittante G.
        Breed of cow and herd productivity affect milk nutrient recovery in curd, and cheese yield, efficiency and daily production.
        Animal. 2018; 12 (28712377): 434-444
        • Strabel T.
        • Szyda J.
        • Ptak E.
        • Jamrozik J.
        Comparison of random regression test-day models for Polish Black and White cattle.
        J. Dairy Sci. 2005; 88 (16162544): 3688-3699
        • Strabel T.
        • Szyda J.
        • Ptak E.
        • Jamrozik J.J.I.B.
        Comparison of random regression test-day models for production traits of dairy cattle in Poland.
        in: Proceedings of the 2003 Interbull Meeting. Interbull Centre, Uppsala, Sweden2003: 197-201
        • Tsuruta S.
        • Misztal I.
        THRGIBBS1F90 for estimation of variance components with threshold-linear models.
        in: Proc. 8th World Congress Gen. Appl. Livest. Prod. World Congress on Genetics Applied to Livestock Production, Belo Horizonte, Brazil. 2006: 27-31
        • Van Slyke L.L.
        • Publow C.A.
        The Science and Practice of Cheese-Making: A Treatise on the Manufacture of American Cheddar Cheese and Other Varieties.
        O. Judd Co., New York, NY1921
        • Wickham H.
        • Averick M.
        • Bryan J.
        • Chang W.
        • D'Agostino McGowan L.
        • François R.
        • Grolemund G.
        • Hayes A.
        • Henry L.
        • Hester J.
        • Kuhn M.
        • Pedersen T.L.
        • Miller E.
        • Bache S.M.
        • Müller K.
        • Ooms J.
        • Robinson D.
        • Seidel D.P.
        • Spinu V.
        • Takahashi K.
        • Vaughan D.
        • Wilke C.
        • Woo K.
        • Yutani H.
        Welcome to the tidyverse.
        J. Open Source Softw. 2019; 41686
        • Windig J.J.
        • Calus M.P.L.
        • Beerda B.
        • Veerkamp R.F.
        Genetic correlations between milk production and health and fertility depending on herd environment.
        J. Dairy Sci. 2006; 89 (16606748): 1765-1775
        • Windig J.J.
        • Mulder H.A.
        • Bohthe-Wilhelmus D.I.
        • Veerkamp R.F.
        Simultaneous estimation of genotype by environment interaction accounting for discrete and continuous environmental descriptors in Irish dairy cattle.
        J. Dairy Sci. 2011; 94 (21605783): 3137-3147
        • Yan M.-J.
        • Humphreys J.
        • Holden N.M.
        An evaluation of life cycle assessment of European milk production.
        J. Environ. Manage. 2011; 92 (21055870): 372-379