Introduction
The development of a breeding programme in Simmental breed is a decisive step in the evolution of a crossbreeding system in tropical regions, because of the need to use animals of superior genetic merit. However, Brazil is a country characterized by a large territory and bioclimatic diversity that culminate in the development of numerous production systems and increase the importance of including the effect of genotype–environment interaction (GEI) in genetic evaluations. For this purpose, reaction norm models have been widely used in animal breeding studies of cattle (Chiaia et al., Reference Chiaia, Lemos, Venturini, Aboujaoude, Berton, Feitosa, Carvalheiro, Albuquerque, Oliveira and Baldi2015; Ribeiro et al., Reference Ribeiro, Eler, Pedrosa, Rosa, Ferraz and Balieiro2015; Araujo Neto et al., Reference Araujo Neto, Pegolo, Aspilcueta-Borquis, Pessoa, Bonifácio, Lobo and Oliveira2018), poultry (Veloso et al., Reference Veloso, Abreu, Mota, Castro, Silva, Pires, Lima and Boari2015; Caetano et al., Reference Caetano, Silva, Calderano, Silva, Ribeiro, Oliveira and Mota2017), pigs (Knap and Su, Reference Knap and Su2008) and buffaloes (Rodrigues et al., Reference Rodrigues, Carneiro, Ramos, Ambrosini and Malhado2015; Freitas et al., Reference Freitas, Hurtado-Lugo, Santos, Aspilcueta-Borquis, Pegolo, Tonhati and Araujo Neto2019).
One of the main selection criteria used in beef cattle breeding programmes are weights measured at standard ages because of their easy measurement and moderate heritabilities (Araujo Neto et al., Reference Araujo Neto, Lobo, Mota and Oliveira2011). A matter of discussion in the analysis of these traits is the existence of genetic sexual dimorphism (Garrick et al., Reference Garrick, Pollak, Quaas and Van Vleck1989; Van Vleck and Cundiff, Reference Van Vleck and Cundiff1998; Näsholm, Reference Näsholm2004). Oikawa et al. (Reference Oikawa, Hammond and Tier1999), Diaz et al. (Reference Diaz, Araujo Neto, Marques and Oliveira2011) and Pegolo et al. (Reference Pegolo, Albuquerque, Lôbo and Oliveira2011) discuss the existence of sexual dimorphism in the phenotypic responses of animals to environmental changes, considering either discrete or continuous environments. Within this scenario, the aim of this study was to evaluate the effect of GEI on the yearling weight of Simmental cattle raised in Brazil, including models with sex dimorphism.
Materials and methods
For this study, records of weight at 365 days (W365) measured in 35 000 Simmental animals were provided by the Brazilian Association of Simmental and Simbrasil Cattle Breeders (ABCRSS). First, the contemporary groups were formed by the concatenation of farm, year and season of birth, sex, and management group at 365 days. The following consistency criteria were applied to the database: (a) sires with at least ten offspring were maintained; (b) contemporary groups with at least three animals were considered and (c) using the median and interquartile range, outliers within each contemporary group were eliminated. After consistency analysis, the database provided data for 17 464 animals divided into 2889 contemporary groups and a relationship matrix containing 26 258 animals (Fig. 1).
Structure of the database used in the reaction norms analysis: (a) number of male records; (b) number of female records and (c) average weights at 365 days of age for both sexes.
A key step in the analysis of reaction norms is the creation of an environmental gradient (EG) on which the trait of interest is regressed. For this purpose, the average W365 of the contemporary groups were first calculated. These values were standardized to the mean and standard deviation of the trait, multiplied by 10 and truncated. This procedure was repeated separately for each sex as proposed by Pegolo et al. (Reference Pegolo, Albuquerque, Lôbo and Oliveira2011). Figure 1 shows the distribution of males and females by EG, as well as the mean W365 by sex and these gradients.
Two approaches were adopted in this study to evaluate reaction norms for W365: a single-trait model that did not consider the heterogeneity of variances by sex for regression on the EGs, and a multitrait model in which the data for males and females were separated and considered different traits for the analysis of sexual dimorphism. Random regression models adjusted by first-order Legendre polynomials as a basis function were used. The model can be written in the matrix form (generalization for the two approaches) as follows:
where y, b, a and $\varepsilon$
are the vector of observations, systematic effects (contemporary groups), solutions for the additive genetic random regression coefficients and residuals, respectively; X and Z are incidence matrices that associate vectors b and a to the vector of observations y. In this analysis, the heterogeneity of the residual variance was modelled in four equidistant classes (determined in previous analyses, data not shown). Bayesian methods and Gibbs sampling were used for the estimation of variance components. In the Bayesian context, uniform and multivariate normal distributions were used for the systematic effects and genetic regression coefficients, respectively. For the covariance matrices (matrix of covariance functions and matrix of residual covariance for multitrait analysis), a scaled inverse Wishart distribution was assumed as prior.
The size of the Gibbs chain was determined considering the number of samples necessary for convergence and a minimum effective sample size of 200 samples (for each parameter). Based on these criteria, in both analyses, it was necessary to generate 1 700 000 samples and to consider a convergence period of 200 000 samples. The genetic parameters were estimated using the Gibbs3f90 program (Misztal et al., Reference Misztal, Tsuruta, Strabel, Auvray, Druet and Lee2002). The CODA package (Plummer et al., Reference Plummer, Best, Cowles and Vines2006) available in the R program (R Core Team, 2013) was used for analysis of convergence of the chain and estimation of the effective sample size.
The heritability for weight along the gradient and genetic correlations were calculated in the two analyses. Complementing the previous results, for the evaluation of sexual dimorphism, a statistical test was first applied to determine the difference between the variance and covariance parameters estimated for males and females with the multitrait model as follows: (a) a posterior distribution of the difference was obtained; (b) the 90 and 95% high-density intervals (HDI) of probability were estimated and (c) significance was verified when the HDI did not contain zero. Next, other genetic parameters that describe sexual dimorphism along the EG were calculated. The heritability for sexual dimorphism (delta – Δ) was calculated from the covariance of the difference as reported by Chapuis et al. (Reference Chapuis, Tixier-Boichard, Delabrosse and Ducrocq1996) and Mignon-Grasteau et al. (Reference Mignon-Grasteau, Beaumont, Poivey and Rochambeau1998), which can be expressed as:
where $\sigma _{m( g\_i) }^2$
and $\sigma _{f( g\_i) }^2$
are the genetic variances for males and females in the ith environment, respectively; $\sigma _{m( e\_i) }^2$
and $\sigma _{f( e\_i) }^2$
are residual variances for males and females in the ith environment, respectively and $\sigma _{mf( g\_i) }$
is the genetic covariance between sexes in the ith environment (it is important to note that the model assumes the absence of residual covariance between sexes).
Results
In single-trait analysis, the posterior means and HDI (between parentheses) for the variances of the intercept and slope coefficient were 995.82 (776.80–1227.00) and 214.00 (102.70–335.70) kg2, respectively. The posterior mean for the correlation between coefficients was estimated at 0.34 (0.05–0.63). For the residual variance classes, values of 1117.06 (979.60–1256.00), 1423.68 (1326.00–1509.00), 1684.26 (1578.00–1800.00) and 1027.27 (840.50–1208.00) kg2 were estimated.
In multitrait analysis, values of 1193.70 and 872.20 for the intercept and of 292.06 and 245.54 for the slope coefficient were obtained for males and females, respectively. For the residual variance classes, the estimated values for males and females were 1141.70 and 969.47 for class 1, 1619.60 and 1184.20 for class 2, 1833.30 and 1513.00 for class 3 and 1114.70 and 958.85 for class 4. Despite point differences, similarities in the posterior distributions of the parameters were observed between sexes (Figs 2(a)–(g)).
Posterior distribution of genetic parameters of males and females, using the multitrait model: (a) genetic variance of the intercept coefficient; (b) genetic variance of the slope coefficient; (c) genetic covariance between intercept and slope coefficient; (d) residual variance (class 1); (e) residual variance (class 2); (f) residual variance (class 3) and (g) residual variance (class 4).
In the analysis of the difference between sexes with the significance test based on the HDI (Table 1), only the residual variance in classes 2 and 3 exhibited a significant difference, suggesting the existence of residual heteroscedasticity due to sex. For the genetic variance components, a difference was only observed for the intercept when the 90% HDI was considered.
Mean and 95 and 90% HDI for the posterior distribution of parameter differences between males and females (male–female), for the multitrait model
HDI, high-density intervals; LL, low limit; HL, high limit.
Analysis of the trend of the heritability estimates obtained with the single-trait model (Fig. 3(c)) along the EG revealed a value of about 0.33 (EG: −21) in the worst environments, which decreased in the intermediate environments and reached a value of 0.24 in EG: −8, with a subsequent increase of the estimates up to 0.513 in EG: +23. Using the multitrait model, similar trends were observed for the heritability estimates, which ranged from 0.25 to 0.54 for males and from 0.23 to 0.50 for females (Fig. 4(c)). Regarding the genetic correlations between extremes of the EG obtained with the single-trait model, an estimate of 0.23 was obtained (Fig. 3(d)). The genetic correlations obtained with the multitrait model (Figs 4(d)–(f)) were also low between extreme environments (−0.05 and 0.18 for males and females, respectively). The genetic correlations between sexes (Fig. 4(d)) within each environment tended to increase, ranging from 0.64 (EG: −21) to 0.98 (EG: +23). The heritability estimates for delta (Fig. 5(a)) exhibited a decreasing trend, with most values being of low magnitude ranging from 0.01 (EG: +12) to 0.20 (EG: −21).
Estimates of genetic variance (a), residual variance (b), heritability (c) and genetic correlation between gradients (d) obtained in unitrait reaction norms analysis.
Estimates of genetic variance (a), residual variance (b), heritability (c), genetic correlation between sexes in environments (d) and genetic correlation between environment in sexes (e: males; f: females) obtained in multitrait reaction norms analysis.
Estimation of heritability of sexual dimorphism in Simmental breed cattle, using results from multitrait analysis.
Discussion
GEI for weight at 365 days (single-trait and multitrait analysis)
A general trend in recent years for the evaluation of GEI has been the consideration of quantitative EGs using reaction norm models. In these models, each individual is evaluated using two regression coefficients: one related to the genetic level of the animals and the other to their plasticity in response to environmental changes (Shariati et al., Reference Shariati, Su, Madsen and Sorensen2007). Thus, the first evidence of the importance of GEI for the weight of Simmental cattle in this study is the significance of the effect of the slope (based on the HDI). In addition, the correlation of moderate magnitude between the two coefficients, which was also reported in other studies on cattle using reaction norms (Pegolo et al., Reference Pegolo, Oliveira, Albuquerque, Bezerra and Lôbo2009; Araujo Neto et al., Reference Araujo Neto, Pegolo, Aspilcueta-Borquis, Pessoa, Bonifácio, Lobo and Oliveira2018), indicates changes in the ranking of animals for the different EGs (Santana et al., Reference Santana, Eler, Bignardi, Menéndez-Buxadera, Cardoso and Ferraz2015).
In all analyses performed using either the single- or multitrait model, there was a decreasing trend in the estimates until the intermediate environments, followed by another increase in better environments, demonstrating changes in the genetic parameters of the population according to the environment to which it is exposed, which characterizes the presence of GEI (Corrêa et al., Reference Corrêa, Dionello and Cardoso2009). Many studies analysing weights have reported an increasing trend of heritability estimates along the EG (Mattar et al., Reference Mattar, Silva, Alencar and Cardoso2011; Santana et al., Reference Santana, Bignardi, Eler, Cardoso and Ferraz2013; Chiaia et al., Reference Chiaia, Lemos, Venturini, Aboujaoude, Berton, Feitosa, Carvalheiro, Albuquerque, Oliveira and Baldi2015), in agreement with the premise established by Hammond (Reference Hammond1947) in which better environmental conditions allow expression of the genetic potential of the animal and a consequent increase in genetic variability. However, since estimates of high magnitude were also observed for extreme unfavourable environments (with very negative values), the greater variability in these gradients may be related to the adaptation potential of the animals to tropical environments and/or unfavourable production conditions.
The magnitude of the heritability parameters estimated with the two models ranged from medium to high. These values agree with those reported in the literature in which many studies found heritability estimates for W365 close to those obtained here: 0.17–0.42 for Nellore cattle (Pegolo et al., Reference Pegolo, Albuquerque, Lôbo and Oliveira2011), 0.21–0.58 for Tabapuã cattle (Souza et al., Reference Souza, Malhado, Neto, Martins Filho and Carneiro2016) and 0.19–0.78 for Mediterranean buffaloes (Rodrigues et al., Reference Rodrigues, Carneiro, Ramos, Ambrosini and Malhado2015). Furthermore, the genetic correlations between extreme environmental groups are within the range of −0.29 to 0.38 reported for weight in cattle (Chiaia et al., Reference Chiaia, Lemos, Venturini, Aboujaoude, Berton, Feitosa, Carvalheiro, Albuquerque, Oliveira and Baldi2015; Ribeiro et al., Reference Ribeiro, Eler, Pedrosa, Rosa, Ferraz and Balieiro2015; Araujo Neto et al., Reference Araujo Neto, Pegolo, Aspilcueta-Borquis, Pessoa, Bonifácio, Lobo and Oliveira2018). The agreement between the present results and the literature suggests the feasibility of the analyses performed here for the evaluation of GEI in the Simmental population in Brazil.
The present findings highlight the importance of considering the results of reaction norms in a breeding programme of Simmental cattle in Brazil. In the case of this breed, one factor that increases the complexity of adopting an animal breeding programme for robustness/plasticity is the frequent use of imported semen from sires not yet tested in Brazil. In this scenario, a complex effect is added to the relationship between environments and genotypes, the genotype–country of origin interaction, which could increase the bias in parameter estimates since this effect would not be considered in the evaluations.
Sexual dimorphism in reaction norm models (multitrait analysis of differences between sexes)
Using multitrait models that included each sex as a different trait, the existence of heteroscedasticity due to sex was only observed in two residual classes (Table 1). These differences might be attributed to specific unregistered or difficult to categorize managements between sexes. In addition, the response of individuals to the environment differs between sexes, as demonstrated by Iung et al. (Reference Iung, Neves, Mulder and Carvalheiro2017) who studied the genetic influence on residual variation. These authors described greater differences between sexes in modelling residual variances than in the average between the two categories.
Evidence of genetic heteroscedasticity is provided by the intercept coefficient, which exhibited a significant 90% HDI. Thus, the genetic average of the animals differed between sexes, which was expected since males are heavier than females. Additionally, the genetic correlation between sexes in lower-level environments and the magnitude of heritability for delta can be considered complementary evidence of heteroscedasticity. For the Simmental breed, several studies have investigated the possible effect of sexual dimorphism on weights at different ages but no consensus exists regarding the significance of this effect (Garrick et al., Reference Garrick, Pollak, Quaas and Van Vleck1989; Van Vleck and Cundiff, Reference Van Vleck and Cundiff1998). In this study, we were thus able to assess the significance of this effect and to quantify it in different environments with a multitrait model and estimation using a Bayesian approach.
Genetic correlations between sexes less than one (as observed in this study) can be attributed to genes located on the sex chromosomes (Ghafouri-Kesbi et al., Reference Ghafouri-Kesbi, Mianji, Pirsaraei, Hafezian, Baneh and Soleimani2014), although this type of analysis was not performed in this study. In a study on swine, Wittenburg et al. (Reference Wittenburg, Teuscher and Reinsch2011) partitioned the additive effects of autosomes and sex chromosomes and reported that the Y chromosome contributed significantly to the genetic variation in birth weight. Genetic correlation estimates less than unity between traits measured in males and females (as described in this study) may be due to limitations of the experimental design, differences in the environments of the performance tests, or differences in the genetic bases and/or selection objectives (Raidan et al., Reference Raidan, Porto-Neto and Reverter2019).
The estimates of genetic correlation between the sexes (Fig. 4(d)) demonstrate greater differences between males and females in lower production level. In a genetic evaluation of the performance of male and female Simmental cattle in two birth seasons, Diaz et al. (Reference Diaz, Araujo Neto, Marques and Oliveira2011) observed lower genetic correlations between sexes in the dry season, providing evidence of an antagonism between genetic sexual dimorphism and production level, similar to the present results. Oikawa et al. (Reference Oikawa, Hammond and Tier1999), in an Angus population, described positive correlations of high magnitude in environments of high technological level and of low magnitude in poor environments, and suggested the existence of a genotype–sex interaction at low production levels. The heritability estimates for the difference between sexes, corroborated by the genetic correlations between males and females, indicate that the lower the production level, the more marked is the sexual dimorphism at the genetic level. Since most beef cattle herds in Brazil use weights at standard ages as the main selection criterion, the lower selection intensity in females may also contribute to this difference. However, Oikawa et al. (Reference Oikawa, Hammond and Tier1999) argued that attributing the cause of sexual dimorphism to this difference in domestic populations is a simplistic explanation. What may explain the phenomenon observed in this study is the fact that in environments with less resources, smaller (and therefore less demanding) females are more likely to become fecundated and to produce offspring. Thus, a lower genetic gain is found in females when compared to males and an antagonistic evolution between sexes is established.
In general, there was little evidence of sexual dimorphism across environments (high correlations). According to Van Der Heide et al. (Reference Van Der Heide, Lourenço, Chen, Herring, Sapp, Moser, Tsuruta, Masuda, Ducro and Misztal2016), this finding can be explained by the phenomenon of dimorphism at advanced ages during the life of domestic animals. In addition, we used linear reaction norms to model GEI. Although this model is frequently used to analyse traits of economic importance in cattle, recent data show that higher-order models have a greater capacity to model the effect of GEI (Carvalheiro et al., Reference Carvalheiro, Costilla, Neves, Albuquerque, Moore and Hayes2019). Thus, a question to be raised is that linear reaction norm models may be unable to compute the different responses of sexes to environmental changes. Pereira et al. (Reference Pereira, Serodio, Ventrua, Araujo Neto and Pegolo2018) proposed the use of cubic reaction norms for the analysis of weight at 450 days in Nellore animals, with the results showing greater genetic stability in females. However, despite providing greater flexibility in the modelling of GEI curves, a larger number of parameters can pose problems related to convergence in the estimation process (Barr et al., Reference Barr, Levy, Scheepers and Tily2013; Bates et al., Reference Bates, Kliegl, Vasishth and Baayeb2015; Matuschek et al., Reference Matuschek, Kliegl, Vasishth, Baayen and Bates2017), overfitting (Hawkins, Reference Hawkins2004; Lever et al., Reference Lever, Krzywinski and Altman2016) and/or collinearity (Budescu, Reference Budescu1980; Brauner and Shacham, Reference Brauner and Shacham1998).
An important aspect that needs to be considered is how to include the concept of sexual dimorphism in reaction norms for animal breeding programmes. Sexual dimorphism could be used to maximize the productive efficiency of breeding females without compromising the weight gain of males by selecting females whose size is consistent with the technological level of the property. Thus, in stress-free environments with abundant food, females with a larger biotype are desired, while in situations of stress or resource scarcity smaller females are preferable (Ritchie, Reference Ritchie1995; Rosa et al., Reference Rosa, Lôbo, Oliveira and Borjas2000). The genetic gains of this selection strategy would be greater in environments with lower production levels not only because of the heritability for delta but also because of the positive correlation obtained with the selection of females.
Other considerations
A first point to be commented on is the contemporary groups size and the distribution of records number along the EG, used in the analyses of reaction norms. Considering that the mean contemporary group for W365 was used as basis for formation of the EG, the CG size may interfere in the results estimation, increasing the amplitude of the HDIs of the parameters and masking the differences between the sexes. Thus, the results presented should be evaluated with criteria, although the point results are interesting.
Another point to be commented on is the genomic tools application, as already used in Simmental breed animals (Wang et al., Reference Wang, Miao, Chang, Xia, An, Li, Lingyang, Lupei, Xue, Junya and Huijiang2019; Cesarani et al., Reference Cesarani, Garcia, Hidalgo, Degano, Vicario, Macciotta and Lourenco2021). In studies of reaction norms, either to investigate GEI or resilience, the single nucleotide polymorphism markers application using the single-step GBLUP method has enabled higher prediction accuracy compared to traditional methods (Oliveira et al., Reference Oliveira, Lourenço, Tsuruta, Misztal, Santos, Araujo Neto, Aspilcueta-Borquis, Baldi, Carvalheiro, Camargo, Albuquerque and Tonhati2018; Araujo Neto et al., Reference Araujo Neto, Santos, Arce, Aspilcueta Borquis, Santos, Guimarães, Nascimento, Oliveira and Tonhati2022). In addition, through genomic analysis it is possible to detect candidate genes that participate in biological processes which turn animals more tolerant to challenging environmental conditions and that can be considered in statistical models to assist in the inclusion of GEI in animal breeding programmes (Braz et al., Reference Braz, Rowan, Schnabel and Decker2021).
Conclusion
In view of the results presented, it can be concluded that there is a need to consider the effects of GEI in the genetic evaluation of 365-day weight in Simmental cattle, but without the need to apply sexual dimorphism models (in view of the reduced evidence presented in the study).
Acknowledgements
We acknowledge the Associação Brasileira dos Criadores das Raças Simental e Simbrasil (ABCRSS) which made the data available for the study.
Author contributions
FRAN, NTP and HNO contributed to the study conception and design. Material preparation and data collection were performed by LFAM and HNO. The statistical analysis was performed by FRAN, RRAB and CDSA. The revision of references and the first draft of the manuscript was written by GFM, APCG, DJAS and JCGS. All authors read and approved the final manuscript.
Financial support
This research received no specific grant from any funding agency, commercial or not-for-profit sectors.
Conflict of interest
The authors declare there are no conflicts of interest.
Ethical standards
Animal care and use committee approval was not necessary for this study because the data were obtained from an existing database.
