All experimental procedures were approved by the Ethics Committee on Animal Use of the Federal University of Minas Gerais (UFMG) (Protocol 255/2010).
Database
The database used in this study consisted of records of Nellore cattle born between 1990 and 2018, which belong to the Novo Mundo Farm located in the municipality of Uberaba, state of Minas Gerais, Brazil. The animals were reared on Urochloa pasture and had free access to mineral supplement throughout the year. The complete description of the farm and main managements used can be found in Passafaro et al. (2015). The sexual precocity traits evaluated included PFC and scrotal circumference measured at 12 and 18 months of age (SC12 and SC18). The following resistance-related traits were analyzed: tick count (TC), gastrointestinal nematode egg count (NEC), and Eimeria spp. sporulated oocyst count (EOC). Weaning weight was used as an anchor trait in the analyses to reduce bias in the parameter estimates for the traits of interest.
Weaning weight
Weaning weight is the body weight adjusted to 205 days of age (W205). Animals with valid W205 records were born between 1990 and 2018. The animals were weaned between 145 and 265 days of age and average daily gain (ADG) from birth to weaning was estimated. The W205 was then obtained by summing the birth weight and weight gain at 205 days. Values outside the interval obtained by the mean of the trait ± 3 standard deviations were excluded. Data of all animals belonging to management groups with fewer than four animals were also excluded. The management groups were formed by animals of the same sex that were weaned in the same month and year. The W205 results and the correlation of this trait with the other traits will not be discussed in this article because W205 was only used as an anchor trait in the analyses.
Reproductive data
The female reproductive data comprised records of heifers born between 1999 and 2016. After weaning, the females were divided into batches of approximately 30 animals and exposed to a teaser bull. The animals remain in the batches for biostimulation until they are over yearling, when new batches are formed for the breeding season. The breeding periods began when the heifers were approximately 21 months old. These were subjected to natural mating for periods until 210 days. During the breeding season, the heifers were submitted monthly to transrectal palpation for the diagnosis of pregnancy.
The sexual precocity trait evaluated in females was PFC, which is a binary trait. Phenotype 1 is assigned to heifers that calved and 0 to those that did not calve after the first breeding season. The month and year of birth of the heifer were considered for formation of the management groups and groups with fewer than four animals were excluded.
The SC12 and SC18 data were collected from males born between 1990 and 2017. Scrotal circumference was measured with a specific tape measure placed transversely around the largest diameter of the scrotum. Measurements obtained between 280 and 430 days of age were considered for SC12 and measurements obtained between 500 and 625 days for SC18. Values outside the intervals obtained by the means ± 3 standard deviations were excluded. The records of all animals in management groups with fewer than four animals were also discarded. In this case, the management group consisted of animals born in the same month and year whose measurements were made on the same date.
Parasitological data
The resistance data were obtained for males and females born between 2010 and 2016. The animals were dewormed at birth and weaning with 4% ivermectin (1ml of Ivermectin per 50Kg of live body weight - Master LP, Ouro Fino Saúde Animal, Cravinhos, SP). No acaricides or anticoccidial agents were used during the data collection period. Natural parasite infestations were evaluated and paddock rotation allowed the batches to pass through different levels of exposure and challenge.
The parasitological data were collected from the animals during weight recording. Ticks were counted on the animal and feces were collected for obtaining NEC and EOC. In females, data were collected during weight recordings at 12 and 18 months of age. In males, data collection was started when the animals were close to 330 days of age and was repeated every 56 days, totaling up to five samplings per animal.
Ticks were counted using the technique proposed by Wharthorn and Utech (1970). Only engorged Rhipicephalus (Boophilus) microplus females present on the right side of the animal and measuring > 4.5 mm in length were counted (Roberts 1968; Wharton and Utech 1970).
For NEC and EOC per gram feces, fecal samples were collected directly from the rectal ampoule of the animal after stimulation and homogenization of the fecal bolus. Plastic bags identified with a sample number for each animal were used. The counts were performed at the Laboratory of Parasite Diseases, Veterinary School of UFMG, by the modified MacMaster technique (Ueno and Gonçalves 1998). Details of the method can be found in Ribeiro et al. (2021). For the analysis, the values counted within the imprinted reticles of the McMaster chamber were considered.
For statistical analysis, the counts were transformed in order to meet the assumptions of normality. Natural logarithmic transformation was chosen because its results are less biased when compared to the use of generalized linear models (St-Pierre et al. 2018).
Outliers were not excluded during database editing in order to preserve the natural biological variability of resistance-related traits. However, data of animals belonging to management groups with fewer than four observations were excluded. In this case, the management group consisted of animals of the same sex and born in the same month and year whose counts were performed on the same date.
After editing the databases for the seven traits, the pedigree data of animals with at least one valid phenotype, as well as of their ancestors, were extracted from the farm’s zootechnical archive. The relationship matrix (represented by \(A\) in the statistical models) consisted of the records of 26,749 animals.
Statistical analysis
Bayesian inference was used to obtain samples of posterior distributions of the (co)variance components and genetic parameters in multitrait analysis using the THRGIBBS3F90 program (Misztal et al. 2022). The analyses were performed for all pairs of precocity and resistance traits. Weaning weight was included as anchor trait in all analyses.
The statistical model for W205 included direct and maternal additive genetic effects, maternal permanent environmental effects, and error as random effects. Management group, linear and quadratic effects of age of cow at calving, and the linear effect of age of animal at weight measurement were included as fixed effects.
The PFC data were analyzed using a threshold model with a probit link function (Gianola and Foulley 1983). In this model, the direct additive genetic effect and error were included as random effects and the management group as fixed effect. For SC12 and SC18, the statistical models included the direct additive genetic effect and error as random effects, and the linear effect of age of animal at measurement and management group as fixed effects.
Repeatability models were used for the resistance traits. The models included the direct additive genetic effect, individual permanent environmental effect and error as random effects, and the management group and age of animal at counting (linear effect) as fixed effects.
Generically, the model can be written in matrix notation as follows:
$$\left[\begin{array}{c}{y}_{1}\\ {y}_{2}\\ {y}_{3}\end{array}\right]=\left[\begin{array}{ccc}{X}_{1}& 0& 0\\ 0& {X}_{2}& 0\\ 0& 0& {X}_{3}\end{array}\right]\left[\begin{array}{c}{\beta }_{1}\\ {\beta }_{2}\\ {\beta }_{3}\end{array}\right]+\left[\begin{array}{ccc}{Z}_{1}& 0& 0\\ 0& {Z}_{2}& 0\\ 0& 0& {Z}_{3}\end{array}\right]\left[\begin{array}{c}{a}_{1}\\ {a}_{2}\\ {a}_{3}\end{array}\right]+\left[\begin{array}{ccc}{M}_{1}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{m}_{1}\\ 0\\ 0\end{array}\right]+\left[\begin{array}{ccc}{C}_{1}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{c}_{1}\\ 0\\ 0\end{array}\right]+\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& {P}_{3}\end{array}\right]\left[\begin{array}{c}0\\ 0\\ {p}_{3}\end{array}\right]+\left[\begin{array}{c}{e}_{1}\\ {e}_{2}\\ {e}_{3}\end{array}\right]$$
where \({y}_{n}\) is the vector of phenotypic observations of trait \(n\); \({X}_{n}\) is the incidence matrix relating the fixed effect in β to the observations; \({Z}_{n}\) is the incidence matrix that assigns the breeding values of \({a}_{n}\) to each individual in the relationship matrix, and \({e}_{n}\) is the residual. The matrix components with subscript 1 refer exclusively to W205; specifically for W205, \({M}_{1}\) is the incidence matrix relating the maternal effect in \({m}_{1}\) to observations \({y}_{1}\), and \({C}_{1}\) is the incidence matrix relating the maternal permanent environmental effect in \({c}_{1}\) to observations \({y}_{1}\). Components with subscripts 2 and 3 indicate a reproductive and a resistance trait, respectively; in this case, \({P}_{3}\) is the incidence matrix relating the permanent environmental effect of the animal in \({p}_{3}\) to the parasite counts.
To obtain samples of posterior distributions, uniform prior distributions were assumed for the vector with solutions for the fixed effects \({\left[{\beta }_{1} {\beta }_{2} {\beta }_{3} \right]}^{t}\); normal distributions for the vectors with solutions for the random genetic effects \({\left[{a}_{1} {a}_{2} {a}_{3} {m}_{1} \right]}^{t}|G\), \({c}_{1}|{\sigma }_{{c}_{1}}^{2}\), \({p}_{3}|{\sigma }_{{p}_{3}}^{2}\), and \({\left[{e}_{1} {e}_{2} {e}_{3} \right]}^{t}|R\), where \(G={G}_{0}⨂A\) and \(R={R}_{0}⨂I\); inverse Wishart distributions for the additive genetic (co)variance matrix \({G}_{0} \text{w}\text{i}\text{t}\text{h} \text{h}\text{y}\text{p}\text{e}\text{r}\text{p}\text{a}\text{r}\text{a}\text{m}\text{e}\text{t}\text{e}\text{r}\text{s} {V}_{a} \text{a}\text{n}\text{d} {S}_{a}\), where \({G}_{0}=\left[\begin{array}{cccc}{\sigma }_{{a}_{1}}^{2}& {\sigma }_{{a}_{1}{a}_{2}}& {\sigma }_{{a}_{1}{a}_{3}}& {\sigma }_{{a}_{1}{m}_{1}}\\ {\sigma }_{{a}_{1}{a}_{2}}& {\sigma }_{{a}_{2}}^{2}& {\sigma }_{{a}_{2}{a}_{3}}& {\sigma }_{{a}_{2}{m}_{1}}\\ {\sigma }_{{a}_{1}{a}_{3}}& {\sigma }_{{a}_{2}{a}_{3}}& {\sigma }_{{a}_{3}}^{2}& {\sigma }_{{a}_{3}{m}_{1}}\\ {\sigma }_{{a}_{1}{m}_{1}}& {\sigma }_{{a}_{2}{m}_{1}}& {\sigma }_{{a}_{3}{m}_{1}}& {\sigma }_{{m}_{1}}^{2}\end{array}\right]\); and for the residual (co)variance matrix \({R}_{0}\) with hyperparameters \({V}_{e}\) and \({S}_{e}\), where \({R}_{0}=\left[\begin{array}{ccc}{\sigma }_{{e}_{1}}^{2}& {\sigma }_{{e}_{1}{e}_{2}}& {\sigma }_{{e}_{1}{e}_{3}}\\ {\sigma }_{{e}_{1}{e}_{2}}& {\sigma }_{{e}_{2}}^{2}& {\sigma }_{{e}_{2}{e}_{3}}\\ {\sigma }_{{e}_{1}{e}_{3}}& {\sigma }_{{e}_{2}{e}_{3}}& {\sigma }_{{e}_{3}}^{2}\end{array}\right]\); and inverse chi-squared distributions for variances \({\sigma }_{{c}_{1}}^{2}\) and \({\sigma }_{{p}_{3}}^{2}\).
Samples of the joint a posteriori distributions of the (co)variance components were generated with the Gibbs sampler using the THRGIBBS3F90 program (Misztal et al. 2022). Chains with 1,100,000 cycles were generated, with a burn-in period of 100,000 cycles and a thinning interval of 10 cycles. The convergence of the Gibbs sampling algorithm was checked by visual inspection of the graph generated with the BOA package (Smith 2005) of the R software. The Geweke test (Geweke 1991) was also used for analysis of chain convergence. Statistical measures of a posteriori distributions were obtained using the POSTGIBBSF90 program (Misztal et al. 2022), with estimation of posterior means and 90% highest posterior density (HPD) intervals.
Genetic (\({r}_{{a}_{n,{n}^{{\prime }}}}\)), residual (\({r}_{{e}_{n,{n}^{{\prime }}}}\)) and phenotypic (\({r}_{{p}_{n,{n}^{{\prime }}}}\)) correlations were estimated for each pair traits. To obtain the other parameters, the samples of a posteriori distributions of each multitrait analysis were grouped in such a way that the estimate for each trait was made from 500,000 samples. The following genetic parameters were calculated from these samples for each trait: additive genetic variance \(\left({\sigma }_{a}^{2}\right)\), individual permanent environmental variance (\({\sigma }_{pe}^{2}\)), residual variance \(\left({\sigma }_{e}^{2}\right)\), phenotypic variance \(\left({\sigma }_{p}^{2}\right)\), and heritability \(\left({h}^{2}\right)\).