Heteroscedastic Reaction Norm Models Improve the Assessment of Genotype by Environment Interaction for Growth, Reproductive, and Visual Score Traits in Nellore Cattle

Simple Summary In animal breeding, the different climatic conditions and production systems in beef cattle can result in the genotype by environment interaction. These interactions also need to be considered in genetic evaluations. In this study, five reaction norm models were tested on eleven traits related to growth, reproduction, and visual score measured in Nellore cattle breeding programs to check changes in the estimated genetic values of sires according to the selection environment. Using the best fitted statistical model, the presence of genotype by environment interaction was observed for some traits such as age at first calving, scrotal circumference, weaning to yearling weight gain, and yearling weight. For these traits, it is possible to select the best sires to increase productivity and reduce environmental sensitivity. Overall, the reaction norms trajectories for these traits seem to be affected by a non-linear component, and selecting robust animals for these traits is an alternative to increase production and reduce environmental sensitivity. Abstract The assessment of the presence of genotype by environment interaction (GxE) in beef cattle is very important in tropical countries with diverse climatic conditions and production systems. The present study aimed to assess the presence of GxE by using different reaction norm models for eleven traits related to growth, reproduction, and visual score in Nellore cattle. We studied five reaction norm models (RNM), fitting a linear model considering homoscedastic residual variance (RNM_homo), and four models considering heteroskedasticity, being linear (RNM_hete), quadratic (RNM_quad), linear spline (RNM_l-l), and quadratic spline (RNM_q-q). There was the presence of GxE for age at first calving (AFC), scrotal circumference (SC), weaning to yearling weight gain (WYG), and yearling weight (YW). The best models were RNM_l-l for YW and RNM_q-q for AFC, SC, and WYG. The heritability estimates for RNM_l-l ranged from 0.07 to 0.20, 0.42 to 0.61, 0.24 to 0.42, and 0.47 to 0.63 for AFC, SC, WYG, and YW, respectively. The heteroskedasticity in reaction norm models improves the assessment of the presence of GxE for YW, WYG, AFC, and SC. Additionally, the trajectories of reaction norms for these traits seem to be affected by a non-linear component, and selecting robust animals for these traits is an alternative to increase production and reduce environmental sensitivity.


Introduction
The main concern of modern breeding programs in the tropics is the selection and mating of animals aiming to obtain progeny able to produce more efficiently under different production systems. Improvements in animal performance are directly related to genetic aspects and environmental conditions in which the animals are raised. Under tropical production systems, Nellore cattle are raised mainly on pasture under a wide range of climatic regions with seasonal variations in forage production and quality with differences in supplementation strategy [1] and temperature variation [2]. Such heterogeneous environmental factors have been recognized to directly affect animal performance and can provide a specific response from animals to changes in production systems [3][4][5].
The effect of genotype by environment interaction (GxE) has received special attention from several studies regarding its effect on traits related to growth, reproduction, and health [6][7][8][9][10]. These authors have pointed out the re-ranking of animals as the main interaction effect of GxE in tropical production systems for traits with economic importance in breeding programs. In Nellore cattle, dissecting the GxE interaction is commonly assessed by reaction norm models [9,11,12]. The reaction norm model allows fitting the phenotypic information as a function of a continuous environmental gradient by including covariates, regulated by environmental and genetic factors, which are shared with the target trait [13]. The reaction norm approach traditionally regresses the phenotypic value on the environmental gradients, assuming linear trajectories partitioning the estimated breeding value of animals into an independent component (intercept) and an environmental gradient dependent component (slope) associated with environmental sensitivity [14][15][16][17][18]. However, there is evidence that environmental sensitivity exhibits a non-linear trajectory across the environmental gradient. The GxE interaction evaluation using reaction norm models that contemplate the linear or non-linear trajectory is of utmost importance in finding the best model that fits the data based on environmental gradients to provide a more accurate genetic evaluation prediction [9].
The knowledge of the effects caused by GxE for traits of growth, reproduction, and visual scores, in addition, to the identification of linearity or not of the reaction norms trajectories, can be relevant to adjust strategies in breeding programs, considering the differences between environmental conditions of production systems. Therefore, more robust animals, i.e., slopes closer to zero and intercepts with high estimated breeding value (EBV) can be selected for herds with variation in the production levels [19]. Hence, the objective was to assess the presence of GxE by using different reaction norm models, contemplating linear and non-linear reaction norm trajectories, for growth, reproduction, and visual score traits evaluated in Nellore cattle taking into account goodness of fitting measures.

Dataset Description
Approximately one million phenotypic records, collected from 1984 to 2016, from animals raised in different commercial farms located in four regions of Brazil (Southeast, Midwest, North, and Northeast) and Paraguay, were used in the present study. This data comprises a subset of the Alliance Nellore dataset (www.gensys.com.br (accessed on 24 September 2022)). The traits considered were birth to weaning weight gain (BWG), conformation (WC), finishing precocity (WP), muscling (WM) at weaning, yearling weight (YW), weaning to yearling weight gain (WYG), conformation (YC), finishing precocity (YP), muscling (YM) at yearling, scrotal circumference (SC), and age at first calving (AFC), which are part of the genetic evaluation of breeding programs.
Visual scores (WC, WP, WM, YC, YP, YM) were obtained from the evaluation of animals belonging to the same contemporary group (CG), i.e., animals that have had an equal opportunity to perform, by trained technicians. The technicians assess the scores of the animals within the CG applying relative scores from 1 to 5 for each trait (conformation, finishing precocity and muscling), where 1 means that the animal has low expression of the trait and 5 means the animal has maximum expression of the trait. Initially, the technicians determine the extremes 1 and 5 and the average animal (3) within each CG. Next, each animal is evaluated individually and scores are assigned.
When measuring the conformation, the technician observes the meat production capacity in the carcass, imagining the slaughter at the time of evaluation, being observed mainly through the animal length, depth, and thoracic arch. Precocity is observed as the degree of finishing of the carcass, in which animals with good finishing precocity are those with good chest opening, good rib depth, heavy groin, and are at the beginning of subcutaneous fat deposition, especially at the base of the tail, combined with good body development. In muscularity, the development of muscular mass is evaluated as a whole.
Summary statistics of the dataset used for each trait after data editing are shown in Table 1. The data editing (removing the missing data and outliers) was carried out using R software [20]. The CG's composition for each trait is described in Table 2. Only data from CGs with at least 20 animals were considered for the analyses. Birth year refers to the agricultural year (from July to June) in which the animal was born. The birth season was separated into two seasons, dry (autumn and winter), in which there is low rainfall, and rainy season (spring and summer), in which there is greater precipitation and, consequently, better pasture conditions. Besides considering CG as fixed effect, the model for weaning traits included the covariate age of animal at recording (as linear effect) and age of dam at calving (as linear and quadratic effects). For SC, the model included the covariate age of animal at recording as linear effect and for YW, YC, YP, and YM as linear and quadratic effects. For WYG, the model included the covariate period in days between weaning and yearling as linear and quadratic effects.

Environment Descriptor
When studying the GxE through reaction norm models, it is necessary to consider a continuous environmental descriptor in order to verify the animal response in each gradient. One way to define the environment descriptor is to use CG solutions previously estimated, once CG estimates reflect the environmental condition in which the contemporaries were raised [21]. Therefore, the environmental gradients (EGs) which were used to describe the level of production were based on CG solutions (obtained from model 1 or 2 described below) from birth to weaning weight gain for the weaning traits, on CG solutions from weaning to yearling weight gain for the visual score at yearling traits and WYG, and on CG solutions from yearling weight for YW, AFC, and SC. To estimate EGs, animal model 1 was used for weaning traits and animal model 2 for yearling traits.
where: y w and y y are the vectors of observations for BWG and WYG (or YW), respectively; b w and b y are the vector of fixed effects (contemporary group and covariates) at weaning and at yearling, respectively; a w and a y are the vectors of genetic additive effects; m w is a vector of random maternal genetic effects; c w is a vector of random maternal permanent environmental effects, e w and e y are the vectors of random residual effects; and X w , X y , Z w , Z y , M w and W w are incidence matrices related to b w , b y , a w , a y , m w and c w , respectively. The following assumptions were made for animal model 1: where: A is the numerator relationship matrix between animals; I is the identity matrix; σ 2 a is the additive genetic variance, σ 2 m is the maternal genetic variance; σ 2 c is the maternal permanent environmental variance, σ 2 e is the residual variance and Aσ am is the genetic covariance between the additive genetic effect and the maternal genetic effect. For animal model 2, the following assumptions were made: a ∼ N(0, Aσ 2 a ); e ∼ N(0, Iσ 2 e ). The parameters of interest were estimated using the restricted maximum likelihood method implemented in AIREMLF90 program [22].
The CG solutions were standardized for a mean 0 and variance 1. Only records belonging to CG with standardized EG solutions between −3 and 3 sd were kept, in which after filtering, the minimum, average, and maximum standardized EG values for BWG and YW were equal to −3.0, 0.0, and 3.0 sd, and for WYG was equal to −2.77, 0.0, and 3.0 sd, corresponding to the respective effects of the CG solutions equal to 98.00, 152.75, and 207. 36

Reaction Norm Models (RNM)
A total of five reaction norm models were tested for the traits evaluated, as in [9]. The first model (RNM_homo) assumes homogeneity of residual variance for all EGs, as below: where: y ij is the phenotypic data of animal j in given environment i, β is the vector of fixed effects and covariates; X j is the incidence matrix related to β; ∅ is the overall linear fixed regression coefficient of y ij inÊG i ;ÊG i is the standardized EG solution of environment i, estimated previously; b 0 j is the random direct additive genetic effect or the intercept of animal j for an average EG, b 1 j is the coefficient of random regression or slope of the animal j in the environment represented byÊG i ; m j is the maternal genetic effect (considered only for the weaning traits), and e ij is residual effects. An attempt was made to also consider the maternal permanent environmental effect in the RNM for weaning traits, but the analyses did not converge. However, it is important to mention that the additive genetic and residual variances were little influenced by using a reduced model, without modeling the maternal permanent environmental effect (Table S1). The following assumptions were assumed for RNM_homo: where: A is the relationship matrix based on pedigree information (⊗ is the Kronecker product); and σ b 0 b 1 are the variances of the intercept, the slope and their covariance, respectively; σ 2 m is the variance of the maternal genetic effect for weaning traits; I is an identity matrix and σ 2 e is the residual variance. The second model (RNM_hete) differed in relation to RNM_homo by assuming that the residual variance is heterogeneous across EGs, using a linear regression onÊG i , with the intercept and slope coefficients being modeled using the log-residual function [23].
The third model (RNM_quad) also assumed heteroscedasticity, and considered a polynomial quadratic regression to model the fixed curve and the reaction norm of the random effects (additive genetic and residual) instead of linear regression, according to the equation: where: ∅ 2 is the quadratic fixed regression coefficient of y ij onÊG i ;ÊG 2 i isÊG i squared; b 2 j is the quadratic effect of the additive genetic effect of the reaction norm of animal j on EG i expressed as a deviation from ∅ 2 ; the other terms were as described for RNM_homo.
The fourth model (RNM_l-l) also assumed heteroscedasticity, and used a linear-linear spline function to model the genetic merit of animals across the EG, using the same knot in the average EG (knot = 0) for all animals, according to the following equation: where: ∅ * 1 is the difference of the regression coefficient between the first and the second segments of the linear-linear spline function of y ij onÊG * i , whereÊG * i is equal to zero if EG i is less than or equal to zero, ifÊG i is greater than zero,ÊG * i is equal toÊG i ; b * 1 j is the difference between the first and the second segments of the linear-linear spline function of the random additive genetic effect of animal j onÊG * i expressed as a deviation from ∅ * 1 ; the other terms were described previously.
The fifth model (RNM_q-q) is similar to RNM_l-l, but considering a quadratic-quadratic spline function, according to the equation: where: ∅ * 2 is the difference of the regression coefficient between the first and the second segments of the quadratic-quadratic spline function of y ij on ∅ * 2Ê G 2 * i , whereÊG 2 * i is equal to zero ifÊG i is less than or equal to zero, ifÊG i is greater than zero,ÊG 2 * i is equal toÊG 2 i ; b * 2 j is the difference between the first and the second segments of the quadratic- quadratic spline function of the random additive genetic effect of animal j onÊG * i expressed as a deviation from ∅ * 2 ; the other terms were described previously. The models RNM_quad, RNM_l-l, and RNM_q-q were also tested considering homogeneity of residual variance, but the heteroscedastic models outperformed them (results not shown). Because of this, was decided to keep only one homoscedastic model (RNM_homo).
In addition, as mentioned before, some analyses were performed considering different maternal random effects to understand the influence of including or omitting these effects in the animal model (Table S1).

Model Comparison
The models were compared based on the transformed values of AIC and Akaike weight (AICw), in which the interpretation of the values becomes straightforward, as it uses what would be the probability of choice of a model in relation to the other proposed models [24]. In addition, the plausibility of the biological interpretation of the results was verified, observing the trajectories of heritability estimates along the environmental gradient and the correlations between intercept and slopes. The heritability estimates were determined as the additive variance in each EG divided by the phenotypic variance (additive plus residual variances) in each EG.
All sires with reaction norms further inspected have at least fifty offspring and among these at least ten offspring in the low EG (−3 to −1) and ten in the high EG (1 to 3) or, for AFC, five offspring in the low and five in the high EG.
Ten sires (top10) with the highest EBV in each environmental gradient (low, medium, and high) were selected to have their reaction norm described along the environmental gradient for YW, WYG, AFC, and SC. Additionally, forty sires with the highest EBV in the medium environmental gradient were selected. Of these, the ten most robust sires (top10 robust), i.e., with a slope closer to zero, were shown ( Figures 5D-8D).
To assess the magnitude of GxE in different environments, Spearman correlations were calculated to compare the sire's classification in the model with the best fit.

Environmental Sensitivity
A plasticity scale was assumed based on the absolute individual value of the slope (b 1 j ) and standard deviation of the population slope (σ b 1 ). The animals were classified as:

Reaction Norm Models
The slope/intercept ratios for slope and intercept variance estimates were low for all weaning traits (Table 3) and for yearling visual scores (Table S2). In addition, it is possible to note that for all weaning traits there is a low correlation between intercept and slope, especially for heteroscedastic models, with the exception of RNM_l-l and RNM_qq, in which a high and negative correlation is observed between the first and second slope, suggesting a difference in the pattern of sensitivity from the first to the second segment (Table 3). According to model comparison criteria (AIC and AICw), for weaning traits (Table 3) and yearling visual scores (Table S2), the linear models (RNM_homo and RNM_hete) were the models that best fit the data.
For YW, WYG, SC, and AFC the slope/intercept ratio and the correlation between slope and intercept were higher than for the other traits (Table 4). For these traits, it can also be noted that there is a high and negative correlation between intercept and slope for RNM_l-l and RNM_q-q, which suggests a difference in the pattern of sensitivity between the first and second segments. The models that best fit the data according to AIC and AICw were RNM_l-l for YW and RNM_q-q for WYG, AFC, and SC (Table 4), which suggests that there is a non-linear component in the reaction norm for these traits. Table 3. Estimates of variance (diagonals), covariance (upper triangular; in italic), and correlation (lower triangular; in bold) among coefficients of reaction norm models (RNM) for the additive genetic effect of weaning traits in Nellore cattle, with respective residual variance estimates, Akaike information criterion (AIC), and AIC weight (AICw).

Heritability Estimates
Weaning traits (BWG, WC, WP, and WM) showed similar heritability estimates (h 2 ) for the different models between the intermediate gradients (−2 to 2 sd). However, a tendency for the overestimation of heritability estimates in the extreme gradients was observed for quadratic models. Furthermore, homoscedastic and linear heteroscedastic models presented similar heritability estimates for weaning traits (Figure 1; Table S3). als 2022, 12, x FOR PEER REVIEW 12 of

Heritability Estimates
Weaning traits (BWG, WC, WP, and WM) showed similar heritability estimates ( for the different models between the intermediate gradients (−2 to 2 sd). However, a te dency for the overestimation of heritability estimates in the extreme gradients was o served for quadratic models. Furthermore, homoscedastic and linear heteroscedas models presented similar heritability estimates for weaning traits (Figure 1; Table S3). For WYG, there are differences between models in the heritability estimates, main at the extreme gradients ( Figure 2). For RNM_homo, the heritability estimates increas with the increase in EGs and tended to be underestimated in poor gradients and over timated in better gradients when compared to models that assume heterogeneity of res ual variance. Quadratic models (RNM_quad and RMN_q-q) resulted in an overestimati of heritability in extreme gradients for WYG. For YC, YP, and YM, similar heritabil estimates were observed for all models (Figure 2), and the same occurred for SC (Figu 3). For YW, heritability estimates of different models followed a similar pattern as WYG, where quadratic models overestimated heritability in extreme environments a RNM_homo underestimated heritability in poor gradients and overestimated it in bet gradients (Figure 3). The greatest discrepancies in heritability estimates between mod were observed for AFC (Figure 3). For this trait, RNM_homo overestimated heritabil for better environments and quadratic models (RNM_quad and RMN_q-q) overestimat heritability for poor environments. For WYG, there are differences between models in the heritability estimates, mainly at the extreme gradients ( Figure 2). For RNM_homo, the heritability estimates increased with the increase in EGs and tended to be underestimated in poor gradients and overestimated in better gradients when compared to models that assume heterogeneity of residual variance. Quadratic models (RNM_quad and RMN_q-q) resulted in an overestimation of heritability in extreme gradients for WYG. For YC, YP, and YM, similar heritability estimates were observed for all models (Figure 2), and the same occurred for SC (Figure 3). For YW, heritability estimates of different models followed a similar pattern as for WYG, where quadratic models overestimated heritability in extreme environments and RNM_homo underestimated heritability in poor gradients and overestimated it in better gradients ( Figure 3). The greatest discrepancies in heritability estimates between models were observed for AFC (Figure 3). For this trait, RNM_homo overestimated heritability for better environments and quadratic models (RNM_quad and RMN_q-q) overestimated heritability for poor environments.

Environmental Sensitivity
Sires' EBVs for medium and high EGs were highly correlated, even for traits with more pronounced GxE effect (WYG, YW, AFC, SC), but low correlations were observ in contrasting EBVs for low with EBVs for medium and high EGs, especially for the tra WYG and AFC (Figure 4). High correlations among EBVs for different EGs were observ

Environmental Sensitivity
Sires' EBVs for medium and high EGs were highly correlated, even for traits with more pronounced GxE effect (WYG, YW, AFC, SC), but low correlations were observ in contrasting EBVs for low with EBVs for medium and high EGs, especially for the tra WYG and AFC (Figure 4). High correlations among EBVs for different EGs were observ

Environmental Sensitivity
Sires' EBVs for medium and high EGs were highly correlated, even for traits with a more pronounced GxE effect (WYG, YW, AFC, SC), but low correlations were observed in contrasting EBVs for low with EBVs for medium and high EGs, especially for the traits WYG and AFC (Figure 4). High correlations among EBVs for different EGs were observed for the weaning and visual score traits ( Figure S2). als 2022, 12, x FOR PEER REVIEW 14 of Figure 4. Estimates of Spearman correlation between sire's EBVs for low, medium, and high gra ents along the environmental gradient for weaning to yearling weight gain (A), yearling weight ( age at first calving (C), and scrotal circumference (D).
Because of the evidence of a non-linear component affecting the reaction norms a the overestimation of heritability for quadratic models (RNM_quad and RNM_q-q) in t extreme environments, it was decided to use results from RNM_l-l to further inspect t reaction norms of the top10 sires for WYG, YW, AFC, and SC ( Figures 5-8, respectivel It is possible to observe that there is a re-ranking of the top10 sires along the environme tal gradients when selecting the top10 sires in the low (−3), medium (0), or high (+3) en ronmental gradients. The strategy of selecting sires tacking into account their EBVs for t average environment and also the variation of their EBVs between EGs (top10 robu sires) resulted in lower reranking, in comparison with selection for target environme (Figures 5-8). As no evidence of GxE was found for weaning traits and yearling visu scores, the reaction norms for all sires in all models are presented as supplementary m terial ( Figure S3-S13), as well top 10 corresponding figures for these traits (Figure S1 S20). Because of the evidence of a non-linear component affecting the reaction norms and the overestimation of heritability for quadratic models (RNM_quad and RNM_qq) in the extreme environments, it was decided to use results from RNM_l-l to further inspect the reaction norms of the top10 sires for WYG, YW, AFC, and SC ( Figures 5-8, respectively). It is possible to observe that there is a re-ranking of the top10 sires along the environmental gradients when selecting the top10 sires in the low (−3), medium (0), or high (+3) environmental gradients. The strategy of selecting sires tacking into account their EBVs for the average environment and also the variation of their EBVs between EGs (top10 robust sires) resulted in lower reranking, in comparison with selection for target environments (Figures 5-8). As no evidence of GxE was found for weaning traits and yearling visual scores, the reaction norms for all sires in all models are presented as Supplementary Material ( Figure S3-S13), as well top 10 corresponding figures for these traits ( Figure S14-S20).      The percentage of robust, plastic, and extremely plastic animals for YW, WYG, SC, and AFC are shown in Figure 9. For all traits, there was a predominance of robust animals according to the adopted criterium. AFC presented the highest proportion of plastic and extremely plastic animals among the studied traits.  The percentage of robust, plastic, and extremely plastic animals for YW, WYG, SC, and AFC are shown in Figure 9. For all traits, there was a predominance of robust animals according to the adopted criterium. AFC presented the highest proportion of plastic and extremely plastic animals among the studied traits. The percentage of robust, plastic, and extremely plastic animals for YW, WYG, SC, and AFC are shown in Figure 9. For all traits, there was a predominance of robust animals according to the adopted criterium. AFC presented the highest proportion of plastic and extremely plastic animals among the studied traits. .

Discussion
The identification of the presence of GxE can be performed by using multi-trait or reaction norms models [26]. In this study, five reaction norm models were evaluated for eleven traits (BWG, WC, WP, WM, YW, WYG, YC, YP, YM, SC, and AFC) in Nellore cattle raised under pasture conditions. The models were compared by the goodness of fitting considering the AIC and residual variance. Heritability estimates, reaction norms of EBVs, and environmental sensitivity were also analyzed.

Reaction Norm Models
For weaning traits and visual scores at yearling, the linear reaction norm models presented the best fit, suggesting that the animals present a single pattern of sensitivity along the environmental gradient. Ambrosini et al. [27], also observed that for weaning weight in Nellore cattle, the homoscedastic model outperformed the heteroscedastic model. However, as can be seen in the Table S1, the animal model surpassed the reaction norm models for WC, WP, and WM according to the AIC criteria, suggesting that there is no presence of GxE for these traits.
For YW, WYG, AFC, and SC, all heteroscedastic models outperformed the homoscedastic models, showing that heterogeneity of variances in different environments can be indicative of the presence of GxE [28]. Residual heterogeneity in the evaluation of GxE should be considered in the genetic evaluation of these traits in Nellore cattle to improve the partition of phenotypic variance in genetic and environmental components, enabling an increase in the response to selection [29]. For YW, AFC, and SC, Chiaia et al. [3] also found the superiority of the heteroscedastic model when compared with the homoscedastic.
For YW, the RNM_l-l model had the best fit and the RNM_q-q model had the best fit for WYG, SC, and AFC, suggesting that animal sensitivity will depend on the level of the environmental gradient, since some animals may be less sensitive to variation in more severe environmental gradients and more sensitive to variation in less challenging environmental gradients and vice versa [9].
According to Streit et al. [30], including residual variance heterogeneity in the model is important for obtaining unbiased genetic parameters. For YW, WYG, AFC, and SC, regardless of the model order in this study, heteroskedastic models were superior to homoscedastic models, mainly for AFC, where the genetic variance was highly overestimated as the environmental gradient improved in the homoscedastic model. According to Meyer [31], choosing the best model is a compromise between what we might want to capture and the amount of information available in the data. Although quadratic models (RNM_quad e RNM_q-q) presented a better fit for the majority of the traits most affected by GxE (WYG, AFC, and SC), they tended to overestimate heritability for extreme environmental gradients. For this reason, RNM_l-l seems to be a good alternative model for a compromise between the goodness of fit and the plausibility of heritability estimates.

Heritability Estimates
Results obtained in this study illustrated that the heritability estimates varied along the environmental gradient (Figures 1-3), due to variation in additive genetic and residual variance (for heteroscedastic models) as a function of the environmental gradient. Heritability estimates varying along EGs were also observed by Ambrosini et al. [32], Chiaia et al. [3], and Ribeiro et al. [33] for YW in Nellore cattle; by Chiaia et al. [3], Lemos et al. [34], and Mota et al. [12] for SC in Nellore cattle; by Chiaia et al. [3], Lemos et al. [34], and Araujo Neto et al. [4] for AFC in Nellore cattle; by Carvalheiro et al. [9] for adjusted weight gain from weaning to yearling in Nellore cattle; by Cardoso and Tempelman [35] for postweaning weight gain in Hereford cattle; by Macneil, Cardoso, and Hay [36] for preweaning gain in Hereford cattle; and by Ambrosini et al. [27] for weaning weight in Nellore cattle.
In general, in low environmental conditions, there was an overestimation of the genetic variance for RNM_quad (Figures 1-3). Trying to overcome this, an RNM_q-q was modeled as spline functions tend to be more robust against problems related to "end-of-range" estimates [37]. Nevertheless, the quadratic spline also overestimated heritability in extreme environments, which may be related to the low number of phenotypes observed in these environments ( Figure S1), impairing the estimates at the extremes [37,38].
Heritability estimates for weaning traits and yearling visual scores, considering reaction norm models in the intermediate gradients (−2 to 2), were very similar to heritability estimates using the animal model ( Figure 2 and Table S3), being another evidence of the lack of GxE for these traits. In addition, visual scores are attributed to CG, which may have impaired the estimation of GxE for those traits. Considering only the additive genetic effect in the reaction norm for weaning traits and yearling visual scores, heritability estimates were close to those found by Vargas et al. [39] using an animal model and a subset of the data used in the present study.
For weaning traits, it should be noted that the direct heritability estimates do not differ when using the complete model or the model with the exclusion of maternal permanent environment effect (Table S1). The maternal genetic variance absorbs maternal permanent environment variance, making this parameter overestimated, but the direct heritability estimate is not affected by the exclusion of this parameter.
Heritability estimates for SC had very similar trajectories for all models along the environmental gradient, indicating a higher opportunity for response to selection in better environments. This, however, would not necessarily result in a better-correlated response to sexual precocity, as there is evidence in the literature that SC is not a good indicator of sexual precocity in better environments [3].
In general, the environmental variation affects heritability estimates as a consequence of the GxE, and the understanding of the behavior of heritability estimates along the EGs is of great importance, because there may be a greater genetic gain for each environment.

Environmental Sensitivity
No difference was found between the environmental gradients when observing the genetic correlation between them (all traits >0.8), suggesting that there is no difference in the phenotypes of this population due to the environmental gradient. However, this can happen because it is a parameter of the entire population studied, in which many animals do not have enough information to estimate the EBV in several environmental gradients. However, when the Spearman correlation between a sire's EBV is observed for YW, WYG, AFC, and SC (Figure 4), we can identify GxE in the studied gradients [40,41]. For visual score traits, this high correlation can be explained by the fact that the animals are evaluated within each CG, reducing the difference between each group.
Correlation estimates between the intercept and the slope(s) for the traits most affected by GxE were of medium to low magnitude, suggesting that animals are reclassified in different environments and that it is possible to carry out a combined selection in order to increase performance and reduce sensitivity [42,43]. However, even if the correlation between the intercept and the slope was high, the re-ranking of animals can still occur, as long as the EG interval is large enough and the slope variance is significantly different from zero because of the different trajectories of animal performance along the EGs are not parallel [44]. For YW (Table 3), WYG (Table 2), AFC (Table 3), and SC (Table 3), the correlations between intercept and slope were higher for the homoscedastic model and lower for the heteroscedastic models. Furthermore, it can be noted that the correlation between intercept and slope of the second segment for YW, WYG, AFC, and SC (Table 4) decreased compared with the correlation between intercept and slope of the first segment.
When a plastic animal is selected in the most favorable environment, its performance changes along the environmental gradient, as seen in Figures 5-8, for YW, WYG, AFC, and SC. Consequently, an animal with lower EBV will be used in a medium or low environment, which would result in lower genetic gains. When trying to get around this, the selection of sires that will be parents in the next generation should consider a sire's EBV for the environment in which its offspring will be raised.
Even though the re-ranking of sires among environmental gradients was expressive for some traits, most sires were classified as robust for YW, WYG, AFC, and SC (Figure 9), following the criteria adopted by Santana Jr et al. [25]. This highlights that observing the absolute value of the slope in relation to its standard deviation may not be an efficient criterion for identifying robust animals.
Mota et al. [5] suggested that selection is most effective when selecting the animal based on its EBV in the environment in which its progeny will be raised. Nevertheless, environmental conditions can vary substantially within a farm in different years, due to the seasonality of rainfall and, consequently, pasture conditions and differences in management. To overcome this, selection based on reaction norms becomes an important alternative for beef cattle breeding programs in the tropics, allowing joint selection for increased productivity and reduced sensitivity.

Conclusions
In Nellore cattle, it is not necessary to consider the presence of GxE in genetic evaluations for BWG and visual score traits (at weaning and yearling). On the other hand, for YW, WYG, SC, and AFC, accounting for GxE is important. Considering heteroskedasticity in reaction norm models improves the assessment of GxE for YW, WYG, AFC, and SC. In addition, the reaction norms for these traits seem to be affected by a non-linear component. The RMN_q-q model (best model according to the AIC criteria for WYG, AFC, and SC) overestimated the variances in the low environmental gradient. Therefore, among the tested models, RM_l-l seems the best option for genetic evaluation of YW, WYG, AFC, and SC in Nellore cattle. Furthermore, selecting robust animals for these traits is an alternative for increasing production and reducing environmental sensitivity.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/ 10.3390/ani12192613/s1, Table S1: variance estimates and heritability for weaning and yearling traits in Nellore cattle according to the animal model used; Table S2: estimates of variance (diagonals), covariance (upper triangular; in italic), and correlation (lower triangular; in bold) among coefficients of reaction norm models (RNM) for the additive genetic effect of yearling visual score traits in Nellore cattle, with respective residual variance estimates, Akaike information criterion (AIC) and AIC weight (AICw); Table S3: average (minimum and maximum) heritability estimates for weaning, yearling, and reproductive traits in Nellore cattle for different reaction norm models; Figure S1: number of phenotypes by environmental gradients for birth to weaning weight gain (A), weaning to yearling weight gain (B), weaning to yearling weight gain used for yearling visual scores (C), yearling weight (D), yearling weight used for age at first calving (E), and yearling weight used for scrotal circumference (F); Figure S2: estimates of Spearman correlation between sire's EBVs for low, medium, and high gradients along the environmental gradient for birth to weaning weight gain (A), conformation at weaning (B), finishing precocity at weaning (C), muscling ate weaning (D), conformation at yearling (E), finishing precocity at yearling (F), and muscling at yearling (G); Figure S3