Evaluation of genotype by environment interactions on milk production traits of Holstein cows in southern Brazil

Objective This study assessed the possible existence of genotype by environment interactions for milk, fat and protein yields in Holstein cattle raised in one of the most important milk production basins in Brazil. Methods Changes in the genetic parameters and breeding values were evaluated for 57,967 animals from three distinct regions of southern Brazil, divided according to differences in climate. The genotype by environment interaction was determined by genetic correlations between regions, estimated by the restricted maximum likelihood, considering the animal model. Bull rankings were investigated to verify the ratio of coincident selected animals between regions for each trait. Results The estimates of heritability coefficients were similar between two regions, but were lower in the third evaluated area, for all traits. Genetic correlations between regions were high, ranging from 0.91 to 0.99 for milk, fat and protein yields, representing the absence of a genotype by environment interaction for productive traits. The percentage of selection error between regions for the top 10% of animals ranged from 0.88% to 2.07% for milk yield, 0.99% to 2.46% for fat yield and 0.59% to 3.15% for protein yield. Conclusion A slight change in genotype between areas was expected since no significant genotype by environment interactions were identified, facilitating the process of selecting Holstein cattle in southern Brazil.


INTRODUCTION
Reproductive biotechnology has resulted in the extensive and global distribution of genetic material from animals with high productive potential. For this reason, genotypes that are chosen from areas dissimilar to the region in which they are to be introduced may not satisfy breeder expectations. This aspect may be related to the genotype lacking an adaptation to the territory in which it is introduced; this in turn will affect the genotype by reducing the expression of its maximum productive potential [1]. The genotype×environment interaction (GEI) has several implications that may significantly affect dairy production, including factors responsible for the production of milk, fat and protein. Apart from this, the individual breeding values can be affected, finally culminating in a reclassification of the breeding rankings [2].
In Brazil, a major proportion of milk production comes from the Holstein breed [3]. However, very few studies have focused on productive traits in response to the variety of climates prevailing in the different regions of the country. Most temperate countries, particularly the United States, Canada and Netherlands, which have preferred the Holstein breed for a long time, are well-known as big exporters of semen from the breed [4]. They are now expanding their markets, making successful entries into the tropical and subtropical zones of the world, including Brazil. It therefore becomes crucial to assess whether the imported genetic material will meet the performance expected by breeders in various environmental contexts and emphasises the significance of GEI studies in dairy cattle.
The Paraná Basin in southern Brazil is well-known for having the highest dairy productivity per animal in the country. However, it is characterised by sharp climatic diversity, with some regions experiencing temperate climate conditions, while others have subtropical climate [5]. Therefore, such climatic differentiations occurring within the state do exert some influence that can alter the traits of interest to some degree. Furthermore, in Brazil, findings regarding the presence of GEI in the Holstein breed are meagre, with reports being available but limited only in terms of milk production [6] and not fat and protein production. Therefore, the present study focused on assessing the influence of the GEI on the production of milk, fat and protein in Holstein herds in Paraná, Brazil.

MATERIALS AND METHODS
Data on 305 d milk yield (MY), fat yield (FY), and protein yields (PY) recorded from 57,967 primiparous cows in 375 herds were taken from the database of the Holstein Association in Paraná, Brazil, considering the years between 1990 and 2015. A pedigree file containing the animal's identification, sire and dam was used, totalling 79,387 animals in the relationship matrix.
The effects included in the model were the fixed effects of contemporary group and the cow age at calving as a covariate (linear and quadratic), beyond the random additive genetic effect. The contemporary groups were created considering the interactions of herd-year-season, with four seasons of the partum being considered, i.e., December to February, March to May, June to August, and September to November. The data were checked and records that included errors, insufficient information, animals of unknown parentage, progenies of bulls which only pertain to one herd, and contemporary groups containing fewer than three animals were removed.
The state of Paraná in southern Brazil is situated between 22° 30′58″ and 26° 43′00″ south latitude and 48° 05′37″ and 54° 37′08″ west longitude and extends across 199,307 km 2 . Based on the climatic classification of the Secretariat of Agriculture of the State of Paraná, Brazil [7], the following different regions were easily distinguished: region 1 (R1): mesothermal climate humid and super humid (lacking a dry season, but experiencing severe winter, severe and frequent frost, rainy and mild summers, at altitudes higher than 850 to 900 m); region 2 (R2): mesothermal climate with no dry season (harsh winter with average rainfall and some frost, rainy summers and high temperature, altitudes normally below 850 to 900 m); region 3 (R3): mesothermal climate with a distinct dry season (hot summers and low occurrence of frost, mostly lower than 850 to 900 m).
Multi-trait analyses were performed using the GEI for a group that included the three regions (R1+R2+R3), for each individual trait, considering the same trait in each region as a distinct characteristic. The connectivity between the herds was assured by connecting sires, maintaining only those that produced at least one cow simultaneously, in at least two of the three regions in the study. A total of 1,016 sires were assessed in this study, genetically connected by region as follows: R1 and R2 (951 sires); R1 and R3 (528 sires); R2 and R3 (480 sires); R1, R2, and R3 (471 sires).
In the general matrix format, the model was represented as given: y = Xb+Za+e, where y is the vector of the analysed trait, b refers to the vector of solutions for the fixed effects containing the contemporary group and covariate of the age at birth, a is the vector of the solutions for the additive genetic random effect, X and Z are the incidence matrices for the fixed and additive genetic effects, respectively, and e is the vector of the random residuals. For the multi-trait analysis, regarding the joint analysis of the different regions, the matrix model was explained as given:  where the G and R matrices are defined as G = G0⊗A and R = R0⊗I, respecti 120 genetic covariance matrix of the three regions and R0 is the residual matrix fo 121 will be observed.

122
The presence or absence of the GEI was determined by the genetic c 123 explained by Falconer [8]. The variance components and genetic paramete 124 maximum likelihood method using VCE6.0 software according to Groenevel

125
The breeding values were obtained using PEST software [10]. Additionally, th 126 of the bulls were submitted to the Pearson and Spearman correlations, respect 127 where y 1 is the vector of observations of MY, FY, and PY in R1, y 2 is the vector of the observations of MY, FY, and PY in R2 and y 3 is the vector of the observations of MY, FY, and PY in R3, respectively; X 1 , X 2 , and X 3 , are the fixed effects incidence matrices, contained in vectors β 1 , β 2 , and β 3 , respectively; Z 1 , Z 2 , and Z 3 are the incidence matrices of the additive genetic effects contained in vectors a 1 , a 2 , and a 3 , respectively; e 1 , e 2 , and e 3 are the random error vectors associated with vectors y 1 , y 2 , and y 3 ; 1, 2 and 3 are the R1, R2, and R3 regions, respectively. where the G and R matrices are defined as G = G0⊗A and R = R0⊗I, respectively 120 genetic covariance matrix of the three regions and R0 is the residual matrix for th 121 will be observed.

122
The presence or absence of the GEI was determined by the genetic corre 123 explained by Falconer [8]. The variance components and genetic parameters w 124 maximum likelihood method using VCE6.0 software according to Groeneveld [9

125
The breeding values were obtained using PEST software [10]. Additionally, the b 126 of the bulls were submitted to the Pearson and Spearman correlations, respectivel 127 [11]. Regarding the bulls, the coincident animals between R1, R2, and R3 were es where the G and R matrices are defined as G = G0⊗A and R = R0⊗I, resp 120 genetic covariance matrix of the three regions and R0 is the residual matri 121 will be observed.

122
The presence or absence of the GEI was determined by the geneti 123 explained by Falconer [8]. The variance components and genetic param 124 maximum likelihood method using VCE6.0 software according to Groene

125
The breeding values were obtained using PEST software [10]. Additionall

126
of the bulls were submitted to the Pearson and Spearman correlations, resp 127 [11]. Regarding the bulls, the coincident animals between R1, R2, and R3 128 with the highest predicted transmission ability (PTA) were chosen for e 129 method employed by Pedrosa et al [12].

130
A and R = R0

112
are the fixed effects incidence matrices, contained in vectors β1, β2, and β3, respectively; Z1, Z2, and Z3 are the 113 incidence matrices of the additive genetic effects contained in vectors a1, a2, and a3, respectively; e1, e2, and e3 are 114 the random error vectors associated with vectors y1, y2, and y3; 1, 2 and 3 are the R1, R2, and R3 regions, where the G and R matrices are defined as G = G0⊗A and R = R0⊗I, respectively; where G0 refers to the additive 120 genetic covariance matrix of the three regions and R0 is the residual matrix for the three regions where the animal 121 will be observed.

122
The presence or absence of the GEI was determined by the genetic correlations between the regions, as 123 explained by Falconer [8]. The variance components and genetic parameters were estimated by the restricted 124 maximum likelihood method using VCE6.0 software according to Groeneveld [9], considering the animal model.

125
The breeding values were obtained using PEST software [10]. Additionally, the breeding values and classification 126 of the bulls were submitted to the Pearson and Spearman correlations, respectively, adopting the CORR procedure 127 [11]. Regarding the bulls, the coincident animals between R1, R2, and R3 were established, when 10% of the bulls 128 with the highest predicted transmission ability (PTA) were chosen for each trait in each location, based on the 129 method employed by Pedrosa et al [12].
I, respectively; where G0 refers to the additive genetic covariance matrix of the three regions and R0 is the residual matrix for the three regions where the animal will be observed.
The presence or absence of the GEI was determined by the genetic correlations between the regions, as explained by Falconer [8]. The variance components and genetic parameters were estimated by the restricted maximum likelihood method using VCE6.0 software according to Groeneveld [9], considering the animal model. The breeding values were obtained using PEST software [10]. Additionally, the breeding values and classification of the bulls were submitted to the Pearson and Spearman correlations, respectively, adopting the CORR procedure [11]. Regarding the bulls, the coincident animals between R1, R2, and R3 were established, when 10% of the bulls with the highest predicted transmission ability (PTA) were chosen for each trait in each location, based on the method employed by Pedrosa et al [12].

RESULTS
The additive, residual and phenotypic variances are shown in Table 1. Table 2 provides the heritabilities for each region and genetic correlations between each region for the three traits. R2 revealed equal heritability for MY in relation to R1. For R1 and R2, the value was 0.21, and for R3 it was 0. 16

DISCUSSION
From the variances gathered for the three traits in the regions under investigation, it was demonstrated that the environmental variance accounted for a major portion of the total variance  [13]. In this instance, selection increased the additive genetic variance in the herds, and thus raised the estimate of heritability of the traits over the long term.
The heritability values for MY, FY in R1 and R2 were found to be moderate in magnitude and were similar to those reported by Campos et al [14] in their evaluation of data from Holstein cows in Brazilian herds. Huquet et al [15] also found values of moderate to high heritability (from 0.39 to 0.47) for MYs and FYs in their evaluation of Holstein cattle in France. Such high values for heritability recorded in some studies may occur because these herds undergo significant genetic selection in the programs widely adopted in their countries and form adapted genotypes. Additionally, the environmental control practiced in these countries is more intense, utilising more homogeneous production systems and thus minimising the environmental variance. The R3 was the region that presented the lower heritabilities. Apart from this, the heritability results obtained in all three regions for PY concur with the work of Bernabucci et al [16] and Campos et al [14], i.e. 0.15 and 0.17, respectively. The results of 0.17, 0.16, and 0.10 reported for R1, R2, and R3, respectively, reiterate that PY has low heritability. Considering this, it can be concluded that selection for FY has a tendency to reveal greater genetic gains in smaller generations, as compared to PY.
From the genetic correlations applied in the estimation of the presence of the GEI for MY, the changes in the expression of the genotype between R1 and R3 (0.99) because of environmental differences, are less in relation to the comparative   [17] stated that the presence of the GEI is evident when the coefficients of genetic correlations are below 0.80. Therefore, no interaction for MY was evident, because all the regions reported values higher than the one mentioned above. Zwald et al [2] highlighted the role of environment as the possible predominant factor in the effect on production, provided that the manifestations of the genetic components are influenced by the variables such as temperature differences between the regions, for example, characterising the presence of environmental genotype interactions between the different regions. However, the absence of any interaction among the regions reveals that factors including temperature, herd size and influence of genetic material from outside were inadequate to bring about significant alterations in the phenotype expression among the regions being investigated.
The results of the genetic correlations for FY revealed the absence of any significant influence of a GEI as its value fell within the range of 0.93 to 0.99 between the regions. It has already been established that fat is one of the constituents that induces greater instability in milk; variability in FY occurs due to several factors. However, the results reiterate that the changes caused by distinctive genetics among the animals was of no significance. According to the investigation of Montaldo et al [13], the interactions between the FY in Canada, the United States and Chile is influenced by climatic and regional differences, which tend to induce a notable GEI for FY, confirming the reported results.
For PY, variations in the genetic correlations among the evaluated regions were observed from 0.91 between R1 and R2 to 0.99 between R1 and R3. Therefore, it was accepted that R1 appeared to show more similarity in the effects to R3 in terms of the expression of this characteristic, relative to R2. However, from the magnitude of the results shown, the environmental influence was insufficient to induce a significant interaction effect. This occurs because, according to Kolmodin  et al [18], alterations in the average rainfall exert little or no effect on the PY, reducing the effects of environmental variance on the total component. However, latitudinal differences have been cited in the literature as the principal reason for the alterations in temperature and time of day [19]. Carabaño et al [20] pointed out that temperature can be a determinant factor that influences the interaction on the production of protein in milk. However, in this study, all the regions of the Paraná were situated within the same latitude range (-30° to -20°). These factors thus reveal that latitudinal differences between the three regions did not exert a noteworthy effect on the response of this characteristic. Finally, this fact can also justify the lack of any remarkable effect on MY and FY, which showed correlations close to 1, as mentioned earlier.
The present study relied on Pearson's correlations to confirm the veracity of the genetic correlations and further the understanding of the absence of GEI in the Paraná Basin. The results recorded for MY, FY, and PY revealed only very slight distinctions between the genetic values of the animals from one region to another, because the Pearson correlation ranged from 0.98 to 1.00 for the characteristics, between the three regions. Further, the Spearman correlations for the three traits reiterated that the alterations among the bull ranks was low in magnitude between the regions. The values followed Pearson' s standard, ranging from 0.98 to 1.00. Calus and Veerkamp [21] suggested that the GEI might either induce a change in the classifications of the animals, termed reclassification, or reveal only the presence of differences in the breeding values of the animals, without necessitating any classification changes, otherwise termed the scale effect. This emphasises the fact that the animals assessed here showed only slight changes in breeding values between the different regions, with negligible changes in the classification rankings. Figure 1 shows that the selection error for MY between the regions was in the range of 0.88% to 2.07%, and thus of no significance. This means that the number of bulls selected in  one region, when they would be evaluated in another, should not be insignificant. In practice, this confirms the absence of a GEI for MY between the regions, based on the genetic correlations given earlier in Table 2. The best method would be to choose a specific selection program that could include the progenies evaluated in R1 and then the results of PTAs could be extrapolated for the selection of the same animals in the other two regions; then, the percentage of certainty would be around 99.12% for R3 and 97.93% for R2. This claim was also confirmed for fat and protein production. When considering FY, as revealed in Figure 2, the selection error ranged from 0.99% to 2.46%, and for PY ( Figure 3) from 0.59% to 3.15%. Therefore, for these two characteristics, our confidence in a selection program done in R1 and adopted in the other regions would vary from 96.85% to 99.41%, showing up as highly reliable. It is also noteworthy that the regions that exhibited genetic correlations of 0.99, like R1 and R3 for MY, R2 and R3 for FY and R1 and R3 for PY, showed greater uniformity in the bull PTA dispersions.
According to Mulder et al [22], when high genetic correlations are evident between different regions, the highest average genetic gain occurs with only a single selection program, rather than several different programs. The present study supports this assertion because, as previously mentioned, high genetic correlations demonstrate differences of low magnitude in the animal PTAs. Furthermore, as stated above, most of the animals in the Paraná herds were imported from regions experiencing different climatic conditions, which frequently decreased the selection efficiency. Considering these aspects, a cooperative selection program by breeders specifically tailored for such traits would be very acceptable, as it would create animals adapted to these climatic conditions, which could be used by producers in Paraná. It is also the authors' view that such a methodology, besides increasing genetic progress, could also induce very positive results due to the drop in semen price, thus proving to be advantageous to small producers as well, who form a big part of the primary milk chain in the Paraná Basin.
Based on the findings of this study, it is evident that the three climatic regions assessed in the Paraná Basin in Brazil showed no significant GEI for the traits assessed in the Holstein breed. The breeding values of the animals and their classifications revealed no alterations, irrespective of the regions where they were produced. This makes it clear that the genetic predictions concluded in these regions assessed can, without significant bias, be utilised between regions. Additionally, regarding the absence of any influence exerted by the GEI, a genetic selection program alone will be required for the Holstein dairy herd in Paraná, indicating heightened efficiency and lowered additional expenditure for the implementation of development techniques.

CONFLICT OF INTEREST
We certify that there is no conflict of interest with any financial organization regarding the material discussed in the manuscript.