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

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.


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 et al., 2007;Herbut et al., 2018;Martinez-Castillero et al., 2020).Although nowadays, all selection schemes include traits related to the resilience of cows and their longevity (Stocco et al., 2018;De Vries and Marcondes, 2020), little has been done regarding animals' tolerance to environmental temperature changes (Rovelli et al., 2020).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 et al., 2003).
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 et al., 2007).In the first papers dealing with this topic (Ravagnolo et al., 2000;Aguilar et al., 2009), 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 et al. (2021a) observed that the diurnal temperature variation was also a relevant factor in assessing animal welfare.Bernabucci et al. (2014) 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 et al., 2011;Windig et al., 2011).A second approach considers the phenotypes expressed in different environments as different but correlated traits (Cheruiyot et al., 2020;Martinez-Castillero et al., 2020).
The choice depends on using the environmental gradient as a continuous or discrete factor, even if both situations may coincide (Windig et al., 2011).
Previous studies of genetic parameters have been carried out in several European countries (Brügemann et al., 2011;Bernabucci et al., 2014;Carabaño et al., 2016), Australia (Nguyen et al.;2016;Cheruiyot et al., 2020), New Zealand (Bryant et al., 2007), China (Shi et al., 2021), and South America (Santana et al., 2016), though they focused on the Holstein Friesian breeds and even Bos indicus breeds in tropical climate conditions (Santana et al., 2015).A few studies tested the association between THI values and production traits in other breeds, such as Jersey (Smith et al., 2013) and Brown Swiss (Maggiolino et al., 2020).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.

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 and Davidson (1923): The ECM yield was calculated for each test-day record according to the formula reported by Yan et al. (2011): Cheese yield percentage (CHY) was calculated according to the formula reported by Van Slyke and Publow (1921), modified for Parmigiano-Reggiano cheese (Formaggioni et al., 2015): 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 et al., 2015).
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.
A first editing was performed fitting a linear model for each trait separately: 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, 2019) using lm function and the tidyverse package (Wickham et al., 2019).

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, 2019).
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, 1971) 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 et al., 2020), 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 et al., 2014).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 et al., 2020(Maggiolino et al., , 2022)).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 et al., 2000).
To estimate (co)variance components both for random regression and for multitrait analysis, the Gibbs sampling procedure was carried out by the THRG-IBBS1F90 program (Tsuruta and Misztal, 2006).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 and Misztal, 2006).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
where A is the numerator relationship matrix; I is the identity matrix; G 0 and P 0 are 3 × 3 matrices of (co)variances for additive and permanent environmental effects, respectively; R 0 is residual variances corresponding to each trait; and ⊗ denotes the Kronecker product of matrices.Heritability h a 2 ( ) for each trait in the RN model was calculated as follows: where σ a 2 , σ pe 2 , σ r 2 are the additive genetic, permanent environmental, and residual variance.
The estimate of genetic correlation, RG s,i , and PE correlation, RPE s,i , between intercept and linear regression coefficient (i and s) were calculated as 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 with G o 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: Similarly, PE (co)variances under thermal loads T and T′ have been calculated as Overall phenotypic variance was obtained as

Landi et al.: GENETICS OF HEAT STRESS RESILIENCE IN BROWN SWISS CATTLE
where RV is the estimated residual variance (assumed homogeneous along the THI values).Then, h 2 and ratio of animal permanent environmental variance on the phenotypic variance (PEE) of the trait at THI value were obtained as , .
′ ′ ′ and Genetic and PE correlations (RG TT' and RPE TT' ) between performance under THI values T and T′ were obtained as follows, where sqrt is the square root: 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: where g T 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 et al. (2003), 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 et al., 2010;De Rensis et al., 2015).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 et al. (2012).The chosen THI range is larger when compared with other studies (Cheruiyot et al., 2020), 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: where Y ijklmot 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; c ot was the random additive effect of the cow o for each trait t, and p ot was the random PE effect for cow o for each trait t; e ijklmot is the residual effect for each trait t.
In the covariance structure, c p , p p , and r p , were the additive, PE, and residual effect of each cow for the trait t, respectively: where A is the numerator relationship matrix; I is the identity matrix; and G 0 and P 0 are 3 × 3 matrices of (co)variances for additive and PE effects, respectively.R 0 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: ) 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).

Additive Genetic and PE Variances
The estimates for the MT model shown in Tables 3 and 4 are comparable with PP values and slightly higher for PY for RN (Table 5).
Reaction norm model estimates of the AG and PE variance, with the increasing of THI, are shown in 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,  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 h 2 and the PEE estimate for the RN model are shown in Figures 4 and 5, respectively.The PEE curve followed a similar trend to that of PE in Figures 2 and 3, and the curves described the trend of h 2 in all traits.Generally, a change of the curve shape was observed at THI 80, when the h 2 value increases toward extreme THI values.Only PP, FCM, and CHY decreased at THI values of 80 up to the highest ones.
In Table 3, the h 2 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 h 2 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 h 2 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 et al., 2022) 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 RG 12 , whereas RG 13 was positive for PY, FP, MY, and EMC.RG 23 was negative only in PP, MY, and CHY.RPE, on the contrary, showed negative values except for PP in RPE 12 .Table 4 shows the correlations RG iq and RPE iq between the 3 correlated traits in the MT model.In the case of RG iq , the values were all higher than 0.94.The values of the RG iq 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 RPE iq , 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 RG iq and RPE iq in the same quintile increased considerably from the first to the fifth (0.033, 0.116, and 0.286, respectively).The RPE iq 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 RPE TT' correlations along the THI scale is shown in Figures 6-9  plot of RG TT 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.

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).
In Figures 10 and 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.
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.

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 et al., 2020(Maggiolino et al., , 2022)).So, upper critical THI thresholds are higher than those reported in Italian Holstein Friesians (Bernabucci et al., 2014).The limited literature on this specific breed reported that the Brown Swiss breed is more rustic and resistant to environmental conditions (El-Tarabany et al., 2017).
Many studies identified and reported THI thresholds beyond which production traits worsen, and these val- ues have been used to characterize and model the THI effect.The heat stress threshold has been established worldwide at 72 units of THI (NRC, 1971, Collier et al., 2012).This threshold is consistent with the values reported by Ravagnolo and Misztal (2000) in US Holsteins.Bernabucci et al. (2014), 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 et al. ( 2020) 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 et al., 2006;Martinez-Castillero et al., 2020).The MT model used in this study, similar to what was reported by Cheruiyot et al. (2020), 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 et al., 2009;Bernabucci et al., 2014;Negri et al., 2021b).Some authors introduced more than a unique threshold for traits and calving class (Bernabucci et al., 2014).Other studies highlight significant differences between Italian Brown Swiss and Holstein Friesians (Maggiolino et al., 2020).Carabaño et al. (2016) 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 et al., 2013; Berna-  , 2014;Maggiolino et al., 2020).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 et al., 2016).
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 et al., 2015;Carabaño et al., 2016) or of order 4 (Gernand et al., 2019) has been used, whereas Negri et al. (2021b), 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 et al., 2005).Some authors used the first-order coefficient with a thermoneutrality zone established to THI 60 (Nguyen et al., 2016;Cheruiyot et al., 2020).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 and Meuwissen, 2000;López-Romero and Carabaño, 2003;Strabel et al., 2003).

Variance Components, h 2 , 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 in- dividual 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 and Misztal, 2000;Brügemann et al., 2011;Cheruiyot et al., 2020).The AG variance for all the milk traits from the RN model was not stable at the different THI and traits (Figures 2 and 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 et al., 2014;Cheruiyot et al., 2020).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 et al., 2020).A decreasing trend was observed for PE variances at increasing values of THI or increasing heat stress, which agrees with other authors (Brügemann et al., 2011).
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 et al., 2020).Larger PEE compared with h 2 at THI 75-85 suggests that heat stress sensitivity is cow-specific and mainly acquired rather than genetic (Boonkum et al., 2011).
The increasing values of h 2 at the warmest environmental variable values were also observed by Carabaño et al. (2014Carabaño et al. ( , 2021)), 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 h 2 estimates at the extreme values of the considered independent variable (THI) due to artifacts of polynomial functions.Brügemann et al. (2011) detected a trend of higher h 2 at lower THI values, as recorded in the present study, except for PY, where the trend was the opposite.
Unlike other studies (Carabaño et al., 2014(Carabaño et al., , 2021;;Cheruiyot et al., 2020), 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 h 2 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 et al. (2020), 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.

Genetic and PE Correlation
Our findings revealed that thermal tolerance varies significantly at the THI gradient's extremities .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 and Misztal, 2000;Bernabucci et al., 2014;Carabaño et al., 2014) and sheep (Finocchiaro et al., 2005;Carabaño et al., 2021) 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 et al., 1998).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 et al. ( 2020) 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 et al., 2021).Kolmodin and Bijma (2004) 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, 1959).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 et al. (2020).
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 et al., 2003;Hammami et al., 2015).However, Cheruiyot et al. (2020) 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 et al., 2000).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 and Duangjinda, 2015).According to Aguilar et al. (2009), 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 et al., 2020;Kipp et al., 2021).

EBVs Change Implications
Resilient animals are minimally affected by environmental disturbances (Colditz and Hine, 2016;Berghof et al., 2019), 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 et al. (2014) 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 et al., 2014).Some authors suggested that heat stress in the Holstein breed could increase metabolic heat (Carabaño et al., 2014).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 et al., 2020).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 et al., 2002;Bryant et al., 2007;Nguyen et al., 2016).Cheruiyot et al. (2020) 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.
peiq 2 , and σ riq 2 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 RG iq and RPE iq were calculated as Figures 2 and 3, respectively.
Figure 1.Distribution of temperature-humidity index (THI) values in the data used in this study divided in quintiles.

Figure 2 .
Figure 2. Plot of additive genetics (top) and permanent environmental variance (bottom) posterior mean estimation along temperaturehumidity 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 3 .
Figure 3. Plot 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).
. 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 RPE TT , 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

Figure 4 .
Figure 4. Plot of h 2 (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).

Figure 5 .
Figure 5. Plot of h 2 (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).

Figure 6 .
Figure 6.Tridimensional 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 7 .
Figure 7. Tridimensional 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 8 .
Figure 8. Tridimensional 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 9 .
Figure 9. Tridimensional 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.

Figure 10 .
Figure 10.EBVs 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 11 .
Figure 11.EBVs 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.

Table 1 .
Landi et al.: GENETICS OF HEAT STRESS RESILIENCE IN BROWN SWISS CATTLE Descriptive statistics of data used in the study

Table 3 .
Additive genetic (AG), permanent environmental (PE), and residual (RE) variances; h 2 estimates, and lower and upper bounds of the 95% highest posterior density (HPD95) region from the multitrait model for all investigated traits

Table 4 .
Genetic and permanent environmental correlation in MT model for first, third, and fifth quintiles Trait 1 1 CHY = cheese percentage. 2 HPD95 = lower and upper bounds of the 95% highest posterior density region.

Table 5 .
Landi et al.:GENETICS OF HEAT STRESS RESILIENCE IN BROWN SWISS CATTLE Additive genetic, permanent environmental, and residual variance for random regression coefficients and h 2 for the investigated traits obtained using reaction norm model 1 2CHY = cheese percentage.

Table 6 .
Landi et al.: GENETICS OF HEAT STRESS RESILIENCE IN BROWN SWISS CATTLE Bull ranking comparison based on reaction norm model results using protein percentage and milk yield 1