A general breeding program for the reproductive traits of the Brahman breed is carried out in the experimental station “La Cumaca”, Facultad de Ciencias Veterinarias of the Universidad Central de Venezuela, this represents a valuable source of genes for this population¹.

The methodology used for genetic evaluations has different crucial elements. First, the WA548, used to estimate the genetic values (GVs), may be biased, since it assumes growth is linear. The results published ^2,3,4,5 for B. indicus show growth variations throughout the behavioral tests. Furthermore, the sex of the animal is generally considered a fixed effect in the model, which implicitly assumes that the (co)variance components are the same in both sexes. This approach can incorporate another biased source for the GVs estimation, decreasing their accuracy and thus, affecting the breeding program development^6,7,8. The exposed elements can affect the genetic correlations of the same trait between both sexes (r_FM), which could manifest genotype-sex interaction effects.

The LW has been generally expressed at a fixed point, although it seems quite reasonable to examine these relationships throughout the age. If the information is available, the (co)variance components can be estimated using multitrait (MT) models, or, preferably, random regression (RR) models. Previous studies have compared the MT and RR models for the LW in Bos indicus cattle^2,4; their results show the advantages of RR. However, this longitudinal approach has not examined the relationships between the sex of the animals. Therefore, more evidence is required, particularly when considering possible relationships between LW and female reproductive behavior (RB).

The importance of RB in the beef cattle economy is well known. However, the seemingly low h² of most of the reproductive traits^9,10,11 has been a limiting factor for its use as a direct selection criterion. As an alternative, previous studies have reported the scrotal circumference (SC) or the measured LW in young males, and its response correlated with the RB of females measured by services per pregnancy, days to calving, and the PR at the first service. These encouraging results^12,13 correspond to a fixed age, but the evolution of these trends throughout the age and until first calving remains unknown.

This study aimed to estimate the heritabilities and genetic correlations between the WA548 and the PR of heifers in their first breeding season using a MT model; as well as the genetic (co)variance components of the LW of both sexes regarding age using a RR model; and to compare the genetic values for each i^th age (GVi), based on the RR, with the MT-based GVs.

The experimental station “La Cumaca”, located at 472 m asl., near San Felipe City, Yaracuy State, Venezuela, has an extension of 433 ha, with 300 ha cultivated with Guinea, Star, Swazzi, Pará, and Wire grasses. The annual mean precipitation is 1,650 mm, with a mean temperature of 24 to 31.9 °C, and a mean relative humidity of 84 %¹⁴. It has a herd of pure, registered Brahman cattle, with approximately 180 cows in production.

The LW adjusted to 548 d of age (WA548) of 3,120 animals, born between February 2000 and June 2011, was modified, eliminating records with pedigree inconsistencies, absence, or problems in the date of birth. Finally, there were 2,777 animal records available (1,377 females and 1,400 males). These animals were born from 984 mothers (729 in the data vector) and 107 sires (48 in the data vector). The pedigree file included 3,977 animals. A total of 94,752 individual LW records from 1,776 females and 1,864 males, born between February 1978 and June 2011, were used. These animals were born from 1,291 mothers (929 in the data vector) and 128 sires (58 in the data vector), and the pedigree file included 4,070 animals. These data were edited, eliminating those records with pedigree inconsistencies, absence or problems in the date of birth, and data recorded less than 30 or more than 570 d of age. Data outside the range of ± 3.2 standard deviations within a 30-d range age classes were removed. Finally, a total of 53,258 individual data were available from 1,737 females and 1,803 males.

Several models were analyzed using the SAS GLM procedure¹⁵. Table 1 shows some indicators of the studied data.

There were completed three block analyses using the ASReml3 program¹⁶:

Block 0. Multivariate (MU) model for WA548 and PR.

Where:

y_iis a vector that corresponds to the WA₅₄₈ and the PR analyzed at the same time;

b_iis a fixed-effects vector of the jk ^th combination of sex-year-month (with 275 levels for WA₅₄₈ and 124 for PR);

a_iis a random correlated vector due to the genetic additive effect of the i ^th animal with data and its predecessors without records (4,070 levels) for WA₅₄₈ and PR;

e_ilis a random residual vector correlated between trait 1 and 2;

X and Z are incidence matrices that connect the fixed and random effects with the vector of observations.

This model assumes that:

Var [aiei]=[Gi00Ri]

In which G _i =[σa12σa12σa21σa22]⊗A, where σa12 σa22 represent the genetic variances for both traits, and σa12 their covariance. The residual (co)variance R _i =[σe12σe12σe21σe22] includes σe12 and σe22, which represent the variances for both traits and σe12 their covariance. A is the relationship matrix between all the animals and ⊗ is the product symbol. With these parameters, it was estimated the h² for each trait (hai2) and the genetic correlations (r _gi) between both traits, using linear functions of the corresponding components and classical equations¹⁷. The GVs for each trait were estimated as a solution of the described model, and the accuracy (Acc_ij) of such estimates according to:

Acc_ij = 1−Pevσi2*100

Where: Pev is the prediction error variance (individual value for each animal and study trait), and σi2: is the genetic variance of the trait in the studied population. It was applied this same model in its univariate form in order to present similar parameters to the original program in the experimental station.

To analyze the LW at different ages, it was applied RR, using different models and without considering the maternal effects, which variated in the fitting order of the polynomial for random effects, as well as the estimates of the total or intrasexual (co)variance components of the animal. In total, two model blocks were made:

Block 1 - Assumes that the (co)variance components are the same in both sexes.

Block 2 - Assumes that the (co)variance components are not the same in both sexes.

In both blocks, Y_ijkl represents the different LW estimates in the l ak^th animal, of the j^th sex. The fixed effects (fixed_i) were year-month of control with 674 levels: sex-age at calving with 18 levels, represented in all the models, so that the results can be compared by applying the LogL information criteria; BIC and AIC. The six models only differ in the fitting order of the Legendre polynomial (Φi) for random effects and the residual variance (R), considered homogeneous for block 1 and intra j^th sex in block 2. The strategy applied in block 2 consisted in estimating the inter- and intrasexual (co)variance components of the animal. The Z₁ incidence matrix contains the elements 1 or 0 to connect each observation with the random effects of maternal permanent environment (q_m) with 91 levels. For both blocks, the population growth curve was modeled by a regression coefficient (b₁) dividing the age by the live weight, using the Φi coefficients of third-order, the random genetic effects of r =1, 2, 3 orders, and the individual permanent effect (p_i) of first-order, due to the repetitions of the same trait in the animals throughout the age scale. The expected variance components in both blocks were:

Where A is the numerator of the relationship matrix between the animals with data and their ancestors without records (n= 4,070 total animals). I_p is the identity matrix for the random effects of the individual permanent environment (p=3540 dimension for block 1 with σi2 variance, p_h= 1,737, and p_m=1,803 levels for females and males, respectively in block 2, with variances included in the random regression matrix of individual permanent effect intra j^th sex of the animal (K _P:j). I_q is the incidence matrix for the maternal permanent environment (q dimension=1,291 mothers and σm2  variance). R is the residual error with I_n as incidence matrix (n=53,258 records in block 1 and σe2  variance, and n_h= 25,781 and n_m=27,477 for females and males in block 2 with σe:j2 variance, respectively). G₀ have a random regression matrix (K_G) of (r+1)*(r+1) dimension, in the most complex models of block 2 the elements will be:

ASReml automatically produces the principal component analysis of this matrix, which facilitates the interpretation of the estimated (co)variances trajectories. Herein, K_G is a symmetrical matrix that consists of four submatrices with the same (co)variance components for the genetic effects in females (K_h); males (K_m), and their covariances (K_hm), with their corresponding variances of the intercept (σho2  and σmo2); slope (σhs2  and σms2 ), and covariances (σhso2 , σmso2, σmhs2,  and σhmos2). In these cases, the subscripts o and s indicate intercept and slope, respectively. The described matrix applies to a fitting order of r=1. Therefore, each submatrix has a 2x2 dimension; for r = 2 it will be 3x3, and for r = 3 it will be 4x4, and the additional components will be the quadratic and cubic terms, respectively. For block 1, the K_G matrix does not represent the estimates for each sex. For both blocks, manipulating the elements of these matrices, as well as the r-order Legendre polynomial coefficients (Φi), it is possible to estimate the (co)variance components throughout the age and for each sex¹⁸:σhi2=ΦiKhΦi';σmi2=ΦiKmΦi', andσhmij= ΦiKhmΦj'.

Generally, the genetic parameters of h² and r_g can be determined using classical equations¹⁷. The GVs of LW are determined for each sex using the best model solution where, for the k^th animal it will have:

and where Φi are the corresponding Legendre polynomial coefficients for each i^th point on the age scale. In this model, each animal (total; female or male) will be assigned a genetic function (a_k) linked to the effects of the intercept, slope, and other terms according to the fitting order of the chosen polynomial.

Table 2 shows the genetic parameters obtained from block 0, where the correlations between the estimated GVs determined by MU and MT were incorporated.

For the WA548 there were no differences between MU and MT. For the PR, the MT increased the h², improving the accuracy of the GVs. The r_g between both traits was positive (r_g =0.309), which indicates the absence of antagonism in the improvement of both traits. The correlations between the GVs, based on the models, were higher than 0.897, from which it is inferred that there will be no changes in the order of merit for both procedures. The MT model has additional advantages, manifested in a higher correlation with the GV for PR, as well as greater accuracy in the GV estimates for this last trait, whose h² value was low.

The fitness of the six models in blocks 1 and 2 was determined using the LogL, AIC, and BIC criteria, all three agreed that the third-order polynomial for the genetic effect is the best fit to the data. Block 2 models present better results, which demonstrate that there is a significant variation between sexes for the (co)variance genetic components.

The LW h² throughout the age in both sexes, as well as the genetic correlation between them is showed in Figure 1. The h² trends show slight increases as age progresses, being higher for females. The r_g reflect an inverse pattern with values ranging from 0.25 to 0.35. In contrast, the frequency distribution of the GVs for WA548, estimated according to block 0 models and RR, in Figure 2 shows an overlap of the three GVs estimates. The principal component analysis of the KG matrix for the chosen model 6 demonstrates that the first (vp1) and second (vp2) eigenvectors explained the 57 and 31 % of the genetic variation, respectively. The GVs of the best 200 animals in the MU (official current method) and the RR throughout age and for each sex are shown in Figure 3. This figure shows that in males the trend is positive, while in females we can find animals with negative GVs.

The merit evolution based on the year of birth of the animals is shown in Figure 4. The annual genetic progress was of 0.933 ± 0.021 kg/yr for the WA548, the principal trait in the applied breeding program; for the PR it was of 0.354 ± 0.010 %/yr.

The h² values for the WA548 are similar to several references published in Venezuela for this type of animal^19,20, as well as with the results published for B. indicus in different Latin American countries²¹. Genetic progress for the WA548 was lower than that published by other authors^20,22.

The h² estimates with the MU for PR were low, which is similar to most of the publications about reproductive traits^10,23. However, the h² levels for PR with the MT increased, (h²=0.109 vs 0.087), which boosts the average accuracy of the GVs in 15.7 %. The r_g levels between both traits suggest the absence of antagonism, which means that a selection process for LW and PR is possible, this approach has already been suggested in other studies^11,12,24. The greatest genetic variability (Figure 2) and different h² and r_g levels for the WA548 throughout the age-sex scales (Figure 1) indicate that the expression of this trait should not be considered as an expression of the same trait in both sexes. The latter agrees with other published results^7,8.

The sexual dimorphism (SD), evident in these data, has been studied in detail in the evolutionary context of the populations, creating a debate about the importance of the heterogeneous variance between sexes and its effects in the specialization and adaptability of the populations²⁵, while there are previous statements about changes in SD as a correlated response to fertility selection. In relation to this last point of view, these results present a new approach, this study presents the GVs in each sex for the WA548M and WA548F (model 6 solution results, block 2), which makes possible to estimate a SD of genetic origin like SDg=GVWA548M - GVWPA548F, and these estimates of SDg can be related with the GVs of the same animals for the PR (MT model solution, block 0). The analysis results indicate that a quadratic equation (TG= 0.574 + 0.1045*SDg + 0.000765*SDg² - 0 .0000244*SDg³, and R²= 96.1 %) was the best fit for the data, with an order increase of +1.2 % in PR for each 10 kg of SDg, with a maximum point when 40>SDg<60 with PR= +4.4 %. However, when -10>SDg<0, the PR was -1.4 %. These results are encouraging, but more research is required on this topic, which may have an important practical application in beef cattle production systems.

This study detected a wide genetic variety in the LW and PR of Brahman animals. It is suggested to use MT models, which allow substantial increases in h² values and the accuracy of the estimated GVs, particularly in the PR. The RR analysis indicated that the h² and r_g levels between the LW of females and males vary throughout the age scale, which means that they should not be considered as expressions of the same trait. Finally, this study identified an important genetic variability in sexual dimorphism, which is related to the PR, although this suggestion requires further investigation with a larger number of animals.