If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
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.
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 (
). 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 (
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 (
), 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
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 (
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.
). 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
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 (
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
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.
Set
Item
PP (%)
PY (kg/d)
FP (%)
FY (kg/d)
MY (kg/d)
4% FCM (kg/d)
ECM (kg/d)
CHY (%)
Overall
Min
2.30
0.05
1.19
0.002
1.2
1.39
1.45
0.062
Max
5.03
1.90
6.71
2.34
55.2
55.70
55.41
0.17
Mean
3.67
0.94
4.03
1.03
25.89
25.86
26.57
0.12
SD
0.38
0.27
0.72
0.33
7.65
7.75
7.75
0.015
N
172,835
173,111
171,817
172,401
174,040
172,769
172,871
82,772
First quintile
Min
2.350
0.082
1.460
0.035
2.000
2.072
2.180
0.064
Max
5.030
1.885
6.600
2.298
53.800
53.728
55.417
0.172
Mean
3.732
0.952
4.138
1.053
25.747
26.124
26.887
0.120
SD
0.394
0.270
0.716
0.341
7.778
7.933
7.933
0.015
N
30,250
30,354
30,116
30,182
30,495
30,265
30,292
16,550
Third quintile
Min
2.320
0.085
1.430
0.002
1.900
1.758
1.957
0.062
Max
5.000
1.851
6.630
2.280
54.800
54.342
54.987
0.172
Mean
3.688
0.944
4.050
1.034
25.803
25.869
26.600
0.118
SD
0.376
0.261
0.705
0.329
7.544
7.643
7.635
0.015
N
32,026
32,013
31,765
31,897
32,101
31,939
31,974
16,550
Fifth quintile
Min
2.380
0.052
1.190
0.001
1.200
1.389
1.450
0.062
Max
5.020
1.882
6.570
2.273
52.800
54.585
53.279
0.170
Mean
3.576
0.925
3.885
1.000
26.030
25.457
26.095
0.115
SD
0.369
0.257
0.707
0.320
7.506
7.485
7.469
0.014
N
35,319
35,383
35,188
35,333
35,693
35,390
35,354
17,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:
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 (
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 (
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 (
) 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 (
), 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 (
). 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 (
To estimate (co)variance components both for random regression and for multitrait analysis, the Gibbs sampling procedure was carried out by the THRGIBBS1F90 program (
). 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 (
). 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:
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
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
for each trait in the RN model was calculated as follows:
where
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
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 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:
Similarly, PE (co)variances under thermal loads T and T′ have been calculated as
Overall phenotypic variance was obtained as
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
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:
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 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
, 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 (
). 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
), 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 1Distribution of temperature-humidity index (THI) values in the data used in this study divided in quintiles.
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:
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:
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
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
Table 5Additive genetic, permanent environmental, and residual variance for random regression coefficients and h2 for the investigated traits obtained using reaction norm model
Reaction norm model estimates of the AG and PE variance, with the increasing of THI, are shown in Figure 2, Figure 3, respectively.
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 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 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.
) 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 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 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 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 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.
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.
Bull
Protein (%)
Milk yield (kg/d)
Rank 1
Rank 2
P2
Rank 3
P3
Bull
Rank 1
Rank 2
P2
Rank 3
P3
103,961
1
1
1
70,343
1
1
1
70,607
2
2
3
1
104,054
2
2
2
87,353
3
3
2
1
76,750
3
3
3
44,136
4
4
5
1
109,840
4
4
4
88,331
5
5
4
1
104,686
5
5
5
48,794
6
7
1
6
65,470
6
6
6
65,387
7
6
1
7
71,711
7
7
7
87,849
8
8
8
71,712
8
8
8
83,059
9
9
9
88,348
9
9
9
53,884
10
10
11
1
70,864
10
10
10
29,473
11
12
1
10
1
59,264
11
11
11
48,805
12
11
1
13
1
70,616
12
12
12
110,227
13
13
12
1
71,718
13
13
14
1
71,714
14
14
14
59,462
14
14
13
1
87,176
15
15
16
1
48,794
15
15
15
54,093
16
16
15
1
94,388
16
16
16
35,247
17
18
1
20
3
59,271
17
17
17
83,079
18
17
1
21
3
54,121
18
18
18
44,317
19
19
18
1
64,879
19
19
19
54,034
20
20
19
1
71,653
20
20
21
1
70,776
21
24
3
17
4
70,852
21
21
20
1
48,800
22
21
1
25
3
48,806
22
23
1
22
29,490
23
22
1
22
1
104,683
23
22
1
24
1
53,879
24
23
1
23
1
54,050
24
25
1
23
1
70,483
25
26
1
24
1
70,266
25
24
1
25
48,793
26
25
1
26
88,346
26
26
27
1
39,475
27
27
28
1
71,034
27
27
26
1
59,462
28
28
27
1
65,114
28
28
28
64,663
29
29
29
83,068
29
29
30
1
77,914
30
31
1
30
65,473
30
30
29
1
44,143
31
30
1
33
2
88,331
31
31
31
87,005
32
32
31
1
48,804
32
32
32
98,637
33
33
34
1
83,064
33
33
33
48,792
34
34
38
4
59,672
34
34
34
70,258
35
35
35
70,258
35
35
35
64,649
36
38
2
32
4
44,000
36
36
36
71,034
37
36
1
36
1
48,795
37
37
37
48,959
38
37
1
40
2
76,756
38
38
38
43,999
39
39
39
70,483
39
39
39
87,488
40
40
42
2
54,105
40
40
40
70,852
41
45
4
37
4
71,709
41
41
41
65,473
42
42
41
1
77,923
42
43
1
42
110,221
43
41
2
46
3
87,849
43
42
1
43
60,676
44
43
1
43
1
70,776
44
44
45
1
35,451
45
46
1
44
1
39,295
45
45
44
1
35,484
46
48
2
45
1
44,009
46
46
46
39,369
47
49
2
47
54,110
47
47
47
32,234
48
47
1
48
54,093
48
48
48
54,118
49
44
5
54
5
65,387
49
49
49
48,808
50
51
1
49
1
59,673
50
51
1
50
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 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 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 (
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 (
, 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,
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 (
, 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 (
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 (
). 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 (
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 (
, 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 (
). 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 (
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 (
). 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 (
), 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.
), 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,
, 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 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).
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 (
) showed an antagonistic relationship between high production level and heat tolerance. Antagonism between production and fitness attributes regularly makes long-term productivity improvement difficult (
). 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.
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 (
). 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
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 (
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 (
). 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 (
, 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 (
), 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.,
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 (
). 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 (
). 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 (
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.
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.
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.
00753-6&id=gr1.jpg){kind=link}
00753-6&id=gr2.jpg){kind=link}
00753-6&id=gr3.jpg){kind=link}
00753-6&id=gr5.jpg){kind=link}
00753-6&id=gr6.jpg){kind=link}
00753-6&id=gr7.jpg){kind=link}
00753-6&id=gr8.jpg){kind=link}
00753-6&id=gr9.jpg){kind=link}
00753-6&id=gr10.jpg){kind=link}
00753-6&id=gr11.jpg){kind=link}
00753-6&id=gr4.jpg){kind=link}