Two approaches to account for genotype-by-environment interactions for production traits and age at first calving in South African Holstein cattle

If not accounted for, genotype x environment (G×E) interactions can decrease the accuracy of genetic evaluations and the efficiency of breeding schemes. These interactions are reflected by genetic correlations between countries lower than 1. In countries that are characterized by a heterogeneity of production systems, they are also likely to exist within country, especially when production systems are diverse, as is the case in South Africa. We illustrate several alternative approaches to assess the existence of G×E interactions for production traits and age at first calving in Holsteins in South Africa. Data from 257,836 first lactation cows were used. First, phenotypes that were collected in different regions were considered as separate traits and various multivariate animal models were fitted to calculate the estimates of heritability for each region and the genetic correlations between them. Second, a random regression approach using long-term averages of climatic variables at the herd level in a reaction norm model, was used as an alternative way to account for G×E interactions. Genetic parameter estimates and goodness-of-fit measures were compared. Genetic correlations between regions as low as 0.80 or even lower were found for production traits, which reflect strong G×E interactions within South Africa that can be linked to the production systems (pasture vs total mixed ration). A random regression model including average rainfall during several decades in the herd surroundings gave the best goodness-of-fit for production traits. This can be related to a preference for total mixed ration on farms with limited rainfall. For age at first calving, the best model was based on a random regression on maximum relative humidity and maximum temperature in summer. Our results indicate that G×E interactions can be accounted for when genetic evaluations of production traits are performed in South Africa, by either considering production records in different regions as different correlated traits or using a reaction norm model based on herd management characteristics. From a statistical point of view, climatic variables such as average rainfall over a long period can be included in a random regression model as proxies of herd production systems and climate.

climate in the interior plateau, a subtropical climate in the east and a small area in the northwest with a desert climate, all with summer rainfall. With an average rainfall of only 464 mm, it is a relatively dry country. Most of the country has warm, sunny days and cool nights. Winters in South Africa occur between June and August [1,2].
The total number of dairy cattle in South Africa is estimated at about 1.37 million [3] with a majority of Holstein cows. These cows are kept under a total mixed ration (TMR) or under a pasture-based production system or a combination of both, depending on the geographical region [4]. Currently between 65 and 75% of the milk production is based on pasture. However, many of these pasture-based systems increasingly incorporate additional feeding such as concentrates or forage crops as hay or silage [5].
In South Africa, breeding values for production traits are computed three times a year [6]. Currently, no consideration is given to the different production environments or climatic conditions when bulls are evaluated. In most national genetic evaluations, local genotype x environment (G×E) interactions are considered as negligible for production traits in dairy cattle [7]. However, if these interactions are present but ignored [8], they can substantially decrease the accuracy of the estimated breeding values (EBV) of bulls, generate biases according to the production environment of their daughters, and ultimately decrease the efficiency of selection at the farmer or country level, i.e. the lower true reliability of the genetic evaluations will reflect negatively on the genetic gain because bull selection for a particular environment is based on inadequate information.
G×E interactions in dairy cattle are most often studied by defining phenotypic expressions in various environments as separate traits [9][10][11][12][13][14][15] (say, N traits for N environments) and by estimating the genetic correlations between these N separate traits. An estimate that is significantly different from 1 is interpreted as a signature of a G×E interaction. This is, in particular, the basic assumption underlying the multiple across-country evaluations (MACE) that are routinely calculated by Interbull [10] for national dairy cattle traits. A drawback of the MACE approaches is that the number of (co)variance genetic components (CGC) to be estimated increases quadratically with N. For MACE of milk production in the Holstein breed, this represents 435 CGC for 29 countries [10], which requires significant postprocessing in order to obtain results that are both realistic and consistent.
An intermediate option is to consider only n (< N) underlying traits, in which case the number of GCG increases quadratically only with n. Particular strategies have been proposed to find an optimum number n of underlying traits in international genetic evaluations, such as rank reduction [16,17] and principal components and factor analyses [18,19].
At the national level, G×E interactions also exist (e.g., [11][12][13]), and a particular approach to reveal them consists in adding an interaction term to the traditional quantitative genetic model. For example, a random regression mixed linear model can be built including a reaction norm function of a continuous variable [11,20,21], such as a climatic or herd management variable, to reflect the fact that the observed phenotype is the result of a genetic component which differs according to the local environmental conditions [22]. Other strategies consist in constructing structural models in which environmental and genetic variables are used to describe and measure similarities (e.g., between geographical regions). In fact, when modelling G×E interactions, there is a continuum between the two extreme alternatives to describe the link between the genetic component and the environment of the animal, i.e. "one trait stands for one environment" up to "the genetic component is a continuous function of environmental parameters".
The purpose of this study was to illustrate the wide range of potential genetic evaluation models along this continuum in the particular context of G×E interactions in South African Holstein cows for production traits and age at first calving. The same phenotypes will be described as separate traits in different regions, or by including in the model a reaction norm involving climatic characteristics of the herds averaged over a long period. For simplicity, only first lactation records are considered in this study, to focus more easily on the detection of potential G×E interactions. This is in contrast with the current genetic evaluation model in South Africa, which is based on a 3-lactation fixed regression test-day model [6,23].

Methods
Data from 378,782 first lactation Holstein cows recorded over a period of 30 years (1982-2012) were used in the analysis. The traits analyzed were first lactation milk, fat and protein yields (corrected to a 305-day equivalent) and age at first calving (AFC) in months. Only herds in the seven major administrative regions (provinces) with dairy production were considered. Cows were required to have a known sire and an age at first calving between 20 and 40 months. Contemporary groups were defined within herd and included at least 20 cows that were the progeny of sires with at least 20 recorded daughters each in the complete data set. After editing, records from 257,836 first lactation cows from 702 herds were first grouped into four regions according to global climatic characteristics: Western Cape (1), Eastern Cape (2), Free State, Gauteng, North West and Limpopo (3), Kwazulu-Natal (4). Table 1 shows the number of observations and average performances in each region. The data of 10 global climatic variables (yearly rainfall, maximum and minimum temperature, maximum and minimum relative humidity, solar radiation, evapotranspiration, maximum temperature in summer, maximum humidity in summer and maximum solar radiation in summer) from the closest weather station, based on the GPS position of the herd, were averaged over a period of 50 years and added to each herd. To this list, it would have been interesting to add a composite climatic variable such as the average or extreme daily temperature-humidity index (THI) [24] but this was not possible as only yearly averages of temperature and humidity were available. Overall, these variables are used to reflect the long-term climatic characteristics of the herds and not temporary events such as a heat wave or a period of drought. Table 2 presents the average and standard deviation of each climatic variable in the four groups of regions. For all the herds, the feeding system that is categorized in three possibilities (pasture-based, total mixed ration (TMR) or a mixture between the two) had to be known.
In a first set of analyses, three linear models (summarized in Table 3) were fitted and analyzed using restricted maximum likelihood, as implemented in the Wombat software [25] The fixed effects included were the same for all analyses: herd-year, calving season (quarter) and age at first calving in months (except for the analyses considering age at first calving as the trait under study).
The first model was a univariate model: the phenotype was considered for the whole population as a unique trait with homogenous variance parameters across environments. This is basically the same assumption as in the current fixed regression test-day model used to estimate breeding values of Holstein cows in South Africa [6,14].
The second model investigated G×E interactions between regions using a multivariate model: the same phenotypes in each of the four regions were considered as four distinct traits, without preconceived hypotheses about the source of potential G×E interactions. Genetic correlations between traits (here, between regions) were estimated as well as the specific genetic and residual  variances from which a heritability could be derived for each region.
The third model was also a multivariate model, but with a reduction of the rank of the correlation matrix that was imposed a priori. The Wombat software [25] reparametrizes the genetic (co)variance matrix in order to force a matrix of rank inferior to the number of traits (i.e., 3, 2 or 1, in our case). This results in a more parsimonious description of the correlation between traits and forces an increased resemblance between regions. The extreme situation is a final rank equal to 1, which implies that all the genetic correlations between regions are equal to 1 but the genetic and residual variances-and therefore heritability-may vary. Goodness-of-fit of the various models were compared using the Bayesian information criterion (BIC), which is recommended when comparing non-nested mixed models with a large number of observations [26]. A lower BIC represents a better model after accounting for the number of parameters to be estimated.
A second group of analyses was set up in order to assess whether the observed differences between regions in terms of production could be related to the climatic conditions that prevailed in each herd. Only milk yield was analyzed and a random regression analysis on climatic variables was used. First, as a base model, a unique (homogeneous) genetic variance and heterogeneous residual variances (one per region) were estimated. This is equivalent to the univariate analysis with a specific residual variance for each region. Then, the model was extended to a random regression model by including each climatic variable as a reaction norm variable: as a result, the genetic parameters of the random regression model included two genetic effects: a fixed term and a regression coefficient on the climatic variable (one at a time). The comparison of these models should reveal which climatic variables contribute to a G×E interaction. Then, two climatic variables were included simultaneously together with the fixed term and the most relevant combination of the two climatic variables was chosen. These random regression analyses were initially run assuming a full rank genetic variance-covariance matrix. This rank was then reduced to 1. In such a case, the first eigenvector of the resulting genetic covariance matrix can be interpreted as the optimal linear combination of the climatic variables describing the environmental component of the genetic variability.

Results
A summary of the statistics describing first lactation milk, fat and protein yield and age at first calving in the four regions is in Table 1. The mean value of the three production traits differed between the four regions, but for age at first calving, it only differed between Western Cape and the three other regions. The climatic characteristics of the four (groups of ) regions are in Table 2. Average rainfall is the variable that varies the most between regions, with the lowest average in the Western Cape.
Genetic parameter estimates for production traits and age at first calving measured in each of the four regions are in Table 4, under six models: two univariate genetic models assuming either a common (U1) or four different (U4) residual variances for the four different regions, and four multivariate models assuming different genetic and residual variances in each region and different assumptions about the rank of the genetic variance matrix (from full (MTR4) to rank 1 (MTR1)). The heritability estimates are generally in agreement with the range of values found in the literature [27][28][29][30]. They are lower for protein yield and even more so for fat yield than for milk yield (see Additional file 1: Table S1), possibly because of less accurate on-farm measurements of fat and protein contents compared to milk yield (records are collected by the farmers themselves). Regardless of the trait, rank reductions of the genetic variance matrix did not modify the estimated heritabilities significantly. Heritability estimates varied between regions or models used, especially for production traits. All the genetic correlations estimated using the multivariate models without rank reduction (MTR4) ranged from 0.619 to 0.868 and were usually lower than 0.8. They were all statistically different from 1. These results underline the existence of important G×E interactions, not only as scale effects but also inducing strong bull re-rankings. The effect of rank reduction on genetic correlations was susbtantial, Table 3 Overview of the different models tested y im performance of animal i in region m ; f sum of fixed effects (herd-year, calving season and age at first calving for production traits); a im additive genetic value for animal i in region m ; j reaction norm coefficient for climatic variable j ; b ij standardized value of climatic variable j for animal i ; ε i ( ε im ): residual for animal i (in region m)

Analysis by region
Univariate Analysis including climatic variables Univariate with heterogeneous residual Reaction norm on 1 variable Reaction norm on > 1 variables y ir = f + a i + j b ij + ε im especially for reduction to a rank of 2 (MTR2), or-obviously-1 (MTR1). The inclusion of these interactions in the genetic evaluation model could be more precise by using climatic variables: they appear to define more accurately the variety of environments than the geographical limits, e.g., regions. It appears that including region   as a factor of heterogeneity for residual variances already greatly improved the goodness-of-fit of the model. Table 5 shows the respective importance of the two eigenvectors of the genetic variance matrix when a reaction norm model was used with one climatic variable included at a time for each studied trait (C1-C10). Given the eigenvector coordinates, it follows that the first eigenvalue roughly represents the part of the total genetic variance that is explained by the average genetic effect regardless of the environment. The second eigenvalue represents the part explained by the reaction norm on the climatic variable.
The most important climatic variable for all production traits is the average rainfall, contributing to about 10% of the total genetic variance (C1). In fact, the other climatic variables did not lead to any improvement of the BIC. In contrast, age at first calving was more influenced by the average maximum temperature in summer (C8) and by the average minimum relative humidity (C5), representing about a quarter of the total genetic variance, i.e., the two components related to THI.
To highlight the impact of average rainfall, cows were partitioned into four groups depending on the average rainfall of their herd and a multivariate analysis was performed for milk production (RF4). The results are in Table 6. This analysis showed that a higher level of rainfall was often associated with a lower milk production but also a reduced residual variance and a higher proportion of pasture-based herds (among the herds for which this information was available). However, and quite surprisingly, the heritability was higher in herds with an above-average rainfall level, in spite of their lower average production.
All the models used were compared according to their goodness-of-fit [BIC, see Table 7 and (Additional file 1: Table S2)]. It can be first concluded that including a heterogeneous residual variance according to the region (U4) strongly increased the goodness-of-fit for the univariate model. However, for milk yield, a reaction norm model including average rainfall (C1) led to a better fit (lower BIC) than models based on regions. Including more climatic variables in the model (C2V, C5V) did not necessarily increase the goodness-of-fit, but greatly increase computation time. The best model was clearly the multivariate model that groups herds according to rainfall level (RF4), with distinct residual variances by region. For age at first calving, the best model was the model that included a reaction norm on maximum relative humidity in summer, followed by the model with a reaction norm on maximum temperature in summer, suggesting again that a model including a reaction norm on temperature-humidity index  C2V: Reaction norm on average rainfall and maximum temperature with genetic correlation = 0 and reduced rank 956 C5V: Reaction norm on 5 variables with genetic correlation between "average production" and a linear regression term on climatic variables 967 RF4: Multivariate, grouping herds according to level of rainfall 0 a (THI)-not considered here because only yearly averages of temperature and of humidity were availablecould be even better.

Discussion
After more than 50 years of uninterrupted development, genetic (and genomic) evaluations continue to play a fundamental role in the sustainable improvement of dairy cattle. To avoid computational complications, the first animal models were simple, based on strong assumptions, such as cow selection performed only through their sire-to-cow path or successive lactations considered as repetitions of a same trait. These unrealistic assumptions were progressively dropped in many countries with the development of, e.g., multiple-trait animal models and/or models based on test-day records (e.g. [31]). However, some of these initial simplifications are still in use, for example in international evaluations, where daughter yield deviations (DYD) on a standard lactation scale are first derived in order to evaluate bulls on different national scales; these approximations are considered as minor and indispensable in order to properly account for large within-breed G×E interactions between countries, i.e. in international genetic evaluations, a trade-off is necessary between the accurate modeling of national performances and the exploration of G×E interactions. Our study is another illustration that strong G×E interactions do exist also within country-in our case, South-Africa-with genetic correlations between South-African environments as low as 0.8. In spite of the use of a quite simplistic lactation model, we could demonstrate that genetic correlations for lactation milk yield between regions with different climates could be of the same order of magnitude as the correlation considered in Interbull evaluations between South Africa as a whole and European or North American countries [6]. A likely explanation is the large contrast in average rainfall between farms in different parts of the country, which conditions the management of the cows, i.e. a higher average rainfall favors pasture-based feeding systems. The effect of the feeding systems on genetic parameters has already been demonstrated in previous studies [12,13,32]. The use of climatic variables averaged over many years in the genetic part of the model is an indirect way of accounting for the most likely feeding system when this information is not available for all farms. This is very different from, e.g., the inclusion of THI in the model to reflect the short-term influence of extreme climatic conditions (e.g., [33,34]).
The model that compared herds according to average rainfall over 50 years was statistically better than the model that opposed regions. Interestingly, in contrast with what is classically found in North America or Europe [22] where genetic variance and heritability tend to increase with level of production, the highest heritability was found in the group of herds with substantially lower average production levels (Table 6), due to a much lower estimate of the residual variance. A potential explanation is that large herds based on TMR are grouped and fed according to level of production, with less opportunity to express genetic differences between cows competing for grazing. For clarity, we used a rather trivial lactation model in this study but our conclusions can be extended to more complex genetic or genomic evaluation models, for example based on multivariate random regression models [34].

Conclusions
This study shows that important G×E interactions exist in Holstein cattle in South Africa that could result in large re-rankings of bulls. In absence of a generalized recording of the farms's feeding systems, these interactions could be taken into consideration in national genetic evaluations by defining a heterogeneous residual variance by region and a genetic term involving a reaction norm on long-term averages of climatic variables for each herd. In local breeding programs, calculating a different breeding value depending on the environment (region or feeding system) should be considered.