INTRODUCTION
In developing countries, rabbit meat is used as an alternative healthy source of protein in the human diet, as this meat is low in calories and has a low fat and sodium content. Furthermore, the consumers’ awareness about the health benefits of consuming rabbit meat has increased significantly [
1]. In Brazil, the consumption of this product is still negligible, when compared to other meat categories (e.g., beef, pork, and chicken). Also, the Brazilian consumption of rabbit meat is low in the international context, for example, in comparison to Europe and China [
2].
Growth is one of the most important traits used for evaluation of animal production, as a low growth rate may result in low weight at a certain age or delay in slaughter, which would result in animals with low or no market value. For these reasons, studies focusing on growth curves of animals and the development of robust methods of predicting them have increased significantly in animal sciences [
3].
Carcass is the main component of live weight and has the highest market value in animals used for meat production. Therefore, muscles positively correlated with the animal finishing (e.g.,
Longissimus dorsi) can be good predictors of carcass composition. In this sense, the technique of real-time ultrasonography has been widely used for the evaluation of carcass and body composition, especially in beef cattle [
4], small ruminants [
5], and pigs [
6]. This technique has been relatively less used for the description of growth patterns in rabbits and, thus, there is still much to be done.
For a description of the development and growth of carcass and its components, nonlinear regression models can be used to evaluate strategies that allow the improvement of animal performance, especially regarding the weight gain and feed efficiency. Parameters of nonlinear models estimated for growth of rabbits can provide information for selection for carcass traits associated with the best slaughter age and prediction of productive efficiency indexes.
The lack of scientific studies regarding the use of mixed models and nonlinear regression applied to the description of ponderal and carcass growth trajectories in New Zealand rabbits (Oryctolagus cuniculus) makes it impossible to obtain more reliable definitions of the right age for the slaughter of these animals. Therefore, this study aimed to identify nonlinear regression models, under the context of mixed models, to describe the growth trajectory of New Zealand White rabbits based on weight records and carcass measurements obtained using ultrasonography.
MATERIALS AND METHODS
Animal care
All the experimental procedures carried out in this study were approved by the Committee on Ethics in the Use of Animals (CEUA) of the Federal University of Piauí, Brazil (protocol number 328/17).
Experiment animals and sample collection
The data used in this study were collected from 66 male New Zealand White rabbits (Oryctolagus cuniculus), between 2018 and 2019. During the experimental period, the animals were raised in a didactic-productive module of cuniculture located at the Technical College of Bom Jesus of the Federal University of Piauí, in the municipality of Bom Jesus, Piauí, Brazil (latitude 9°04′57.8″S, longitude 44°19′36.8″W at an altitude of 277 m a.s.l.).
After weaning (30-d-old), the animals were maintained in individual galvanized steel cages (80 cm long, 75 cm wide, and 67 cm high) each containing a plastic feeder and a plastic drinker. The cages were kept in a place with temperature regulated by fan ventilation. Water was provided ad libitum and 150 g/shift/rabbit of pellet commercial ration was offered to the animals (minimum levels ensured: 88% of dry matter; 12% of moisture; 17% of crude protein; 3.37% of ether extract; 15% of crude fiber; 12% of mineral matter; 2% of calcium; 0.75% of total phosphorus; 0.94% of lysine; 0.63% of methionine + cystine; and 2,300 kcal/kg of digestible energy).
Phenotypic records of live body weight (BW) and loin eye area (LEA) collected every seven days from weaning until 150 days of age were used to determine the model that best describes the growth trajectory of the rabbits. All animals used in the study were weighed during the morning shift using a digital scale.
The measurements of LEA were performed using an ultra-sound machine equipped with a linear transducer (probe) with a frequency of 3.5 MHz. Pictures were collected between the 12th and 13th thoracic vertebrae, at the left side of the animal’s body, perpendicular to the long axis of the Longissimus dorsi muscle (loin eye). This anatomical region is widely used for the evaluation of finishing and muscling traits in different species.
Statistical analyses
The traits BW and LEA were calculated as a function of the animal age (days) using each of the following nonlinear models: Brody [
7]; von Bertalanffy [
8]; Richards [
9]; Logistic [
10]; Gompertz [
11]; Meloun 1 [
12]; modified Michaelis-Menten [
13]; and Santana [
14]. Four possibilities were compared for the evaluated models (
Table 1): a model without random effects named as nonlinear fixed effects model; two models with inclusion of random effect; and a model with inclusion of random effect on the parameters
A and
k simultaneously.
Firstly, all the models mentioned above were tested considering the parameters as fixed. In this step, we used the following goodness of fit indicators for the tested models: coefficient of determination (R
2), calculated as the square of the correlation between the predicted and observed values; mean squared error (MSE); the percentage of convergence (%C) of each model; Akaike information criterion (AIC), which is given by
AIC = −
log(
L) + 2
p, where
log(
L) is the logarithm of the likelihood function of the probability density function; Bayesian information criterion (BIC), expressed as
BIC = −2log(
L) +
plog(
n); and the mean absolute deviation of residuals (MAD), as proposed by Sarmento et al [
15]. MAD is given by the following equation:
where: Yi is the observed value; Ŷi is the estimated value; and n is the sample size.
It is important to mention that higher values of R
2 or %C, as well as lower values of MSE, AIC, BIC, or MAD indicate a better model fit. After choosing the model that best described the growth trajectories for LEA and BW of the studied animals, this model was used under the context of mixed models. For this, we considered two parameters that admit biological interpretation (
A and
k), individually or combined. In this step, we used the following criteria to determine the best model: AIC; BIC; and the residual variance (
σe2). AIC and BIC were used because these criteria use the analysis of residual independence and the degree of parameterization to improve the precision of fitting of the nonlinear models evaluated. All analyses were performed using the SAS software [
16].
Under the context of mixed models, considering yij as the measurement j (BW or LEA) of the individual i and tij as the age of this animal (days), the regression model has residuals that follow a normal distribution, with mean zero and constant variance
σe2. The parameters A and k were considered as random, with normal distribution, whereas B was considered as a fixed parameter.
The observations (yij) are independent regarding the index i, but not with respect to j, because at fixing i, the measurements (yij) are taken longitudinally for a single animal. Therefore, it is necessary to include intra-individual variance components in the model.
Considering
Ai~N(A;σa2) and
Ki~N(K;σk2), where
σa2 and
σk2 are variance components of the parameters A and k; if M = (Ai; Ki) is a vector of random effects and
N′=(A,B,K;σe2) is a vector of fixed effects, we have f(yi|ti, N, Mi) ω (Mi;Σ) is the joint probability density function, where
yi′=(yi1; yi2;: : : ; yiJi),ti′=(ti1; ti2;: : : ; tiJi), Σ′=(σa2; σk2), and ω is the joint density of Ai and Ki. The marginal likelihood function is given by
L(N;Σ)=∏I-1n∫f(yi∣ti,N,Mi)ω(M;Σ)dMi.
The estimators for parameters in
N and Σ are obtained maximizing L(
N; Σ) in relation to these numbers. The PROC NLMIXED procedure of the SAS software [
16] minimizes numerically −L(
N; Σ) regarding the parameters
N and Σ, so that the variance-covariance matrix is approximated to the estimators obtained by the inverse of the Hessian matrix.
For the joint analysis, considering simultaneously A and k as random, the assumptions for these parameters are described as
[AiKi]~NM([AiKi],[σa2σa,kσk,aσk2]), where σa,k = σk,a is the covariance between A and k.
For the biological interpretability of the estimates of parameters, we calculated the absolute growth rate (AGR), inflection point (IP), and relative growth rate (RGR) for each selected model. AGR was obtained from the first derivative of the model with respect to time (∂y/∂t). AGR estimates the increase in the animal weight for each time unit t (days or months) in the growth trajectory (i.e., the average growth rate of animals in the population). Regarding IP, this indicates the point where the body growth rate of the animal is the maximum (i.e., the point in which the growth switches from a fast phase to a slower phase). RGR, in turn, is the ratio of AGR to the ŷi (predicted BW or LEA) of the selected model.
RESULTS AND DISCUSSION
The variability of BW and LEA tended to decrease with age (
Table 2), probably because the management provided higher uniformity as the animals aged. Several reports have shown that the BW of rabbits tends to increase with their age and older rabbits have higher meat incorporation [
17,
18]. Therefore, it is expected that male rabbits of the same breed submitted to the same nutritional management, as in the present study, have similar weight gain and carcass development over time.
All the studied models converged; however, the Richards model did not estimate the AIC and BIC values for BW, probably due to the parameterization of this model (
Table 3). In general, the comparison of the models used in this study is based on the quality of fit and computational difficulties. We used the percentage of convergence considering the maximum number of iterations of 2,000. Thus, the best model was the one that best fit the data and had the most appropriate results for the studied traits. It is important to mention that the higher the number of criteria considered, the most reliable is the indication of which is the best model in each situation [
19].
For BW, we observed that the von Bertalanffy model was the most appropriate (
Table 3). On the other hand, the Logistic model was the most appropriate for LEA. The results of this type of comparison rely on the dataset used, which allows one to infer that there is not a single model that is the best in all situations, but there is indeed a model appropriate for each case.
Observing the R
2, it is possible to note that this criterion was not so informative for LEA, as R
2 values are considered low for adjustment in this case (
Table 3). Regarding BW, R
2 can be considered only to exclude some models (Richards, Meloun 1, modified Michaelis-Menten, and Santana). Therefore, the coefficient of determination is not a good indicator for model selection, because the other models (Brody, Gompertz, Logistic, and von Bertalanffy) also showed high R
2 values (0.98). Similar values of R
2 in different models have also been reported by other authors that used this criterion for the selection of models to describe the growth curve in chickens [
20,
21]. In these cases, it is necessary to use different criteria of adjustment, due to the low power of decision of R
2 [
21].
Our results showed expressive differences in the variability of the estimates of MSE in functions of the models used. As mentioned before, lower MSE values indicate better adjustments of the model in function of the trait evaluated. MAD is another parameter that indicates the quality of fit of the growth curves; however, this parameter must not be considered exclusively [
15].
For the von Bertalanffy model, the inclusion of the random effect in
A (σa2) resulted in a reduction of approximately 11% of the AIC and BIC values (
Table 4), in comparison to estimates obtained without the inclusion of random effect.
In the graphical analysis, we can observe that the best adjustment was the one that included the random effect in
A (
Figure 1). On the other hand, the inclusion of two effects in
A and
k in the von Bertalanffy model did not show a satisfactory result, probably because the higher number of parameters makes more difficult the convergence of the model, which complicates the estimation of the error.
The parameter
A represents the estimate of the asymptotic weight, which is interpreted as the adult weight; however, there are controversies about the optimal adult weight, which relies on the species, breed, previous selection, management system, and weather conditions. The parameter
A indicates the weight that the animal can reach at maturity and is useful for prediction of results and planning of the whole activity regarding the raising and breeding [
22,
23].
In a study in New Zealand rabbits, Santos et al [
24] reported that the Gompertz model had the best fit for growth from birth until 150 days of age. However, the mean adult BW (2,026.84 g) reported by these authors was lower than our findings.
The model selection in the present study did not differ from the findings of Ferreira et al [
25], in a study in New Zealand rabbits. In both studies, the von Bertalanffy model had the highest estimate for the parameter
A.
The parameter
k represents the growth speed to reach the asymptotic weight (at maturity). In our study, the estimate of
k obtained using the von Bertalanffy model was 0.01904. Animals with a high maturity rate are considered more precocious than those that have a lower maturity rate. More precocious, prolific, productive, and resistant breeds have a higher market value, as these animals can be slaughtered earlier and with higher carcass yield. This optimizes the costs with feeding (which represents 70% of the total cost of the production) without affecting the environmental and animal welfare rules [
26].
With the inclusion of two random effects, the AIC and BIC values also decreased, in comparison to the results obtained without the inclusion of random effect for BW (
Table 4). In the dataset used in the present study, the most recommended would be to include only one random effect. This can be justified by the difficulty of convergence of the model, as this model became more parameterized when the two random effects were included. Furthermore, it is important to mention that the lowest residual variance for the von Bertalanffy model was observed when the random effect was included in
A.
It is important to mention that the parameter
B is a constant of integration that is related to the initial weight of the animal and represents its degree of maturity at birth. Higher
B values are associated with lower birth weights in the von Bertalanffy model. In the Logistic growth model,
B has a fixed value of one; therefore, in this case, there is no biological interpretation [
27].
When the random effect (
σa2) was included in
A in the Logistic model used for LEA, the AIC, BIC, and residual variance values decreased by 13.0, 13.5, and approximately 50%, respectively, in comparison to the estimates obtained without the inclusion of random effect (
Table 5).
Longitudinal data derived from studies on growth may present different variances during an animal’s life. Furthermore, repeated measurements of the same individual generate correlated residuals, which compromise the efficiency of fixed models [
28,
29]. In the analyses that include random effects in the model, we assumed that the variation in response (BW or LEA) for all individuals follow the same trend (
Figure 2); however, the animals may have different behaviors. Random effects represent individual effects as random deviations from the fixed treatment effects and allow each animal to have its own growth trajectory [
30].
Regarding the residual variance observed in all models, the use of models containing high residual variance can result in a higher distribution of errors at describing the growth of rabbits.
The IP for BW was at 40 days of age, when the average BW of the animals was 752.98 g (
Figure 3). Therefore, the maximum body growth rate occurred at 40-d-old (21.51 g/d) (i.e., at this time, the growth switched from a fast phase to a slower phase). It is important to mention that the individuals did not stop growing after 40 days of age.
Regarding LEA, the IP was observed at 16-d-old, when the animals had 3.16 cm
2 of LEA (
Figure 3). This indicates that the muscle growth of the studied animals became slower at 16 days of age (0.04 cm
2/d). This can be justified by the fact that growth has allometric characteristics (i.e., the tissues have different growth rates that change in different phases of the animal life) [
17].
The IP was fixed in approximately 30% of the asymptotic weight in the von Bertalanffy model, and approximately 35% of the asymptotic growth of LEA in the Logistic model (
Figure 3). Note that the IP of the trajectories coincides with the point of maximum growth rate, which can be useful to define the best moment for slaughter, from the economic point of view. It is important to mention that the age of animals at slaughter is also strongly influenced by the production costs and consumer preferences.
The precocity of ponderal and muscle growths observed in the pre- and post-weaning phases, respectively, indicates that there is a need for adjustments in the feed management of the rabbits after weaning (
Figure 3). For the muscle growth at maturity (90-d-old), we observed almost null values (0.003 mm
2/d). It is important to note that, at this age, the muscle mass reaches its maximum point, thus, the weight gain comes only from gaining fat. Regarding the maturity weight, the animal growth was only 9.09 g/d. From 90 days of age onwards, the weight gain decreased and was lower than that observed at 2-d-old (9.67 g/d). Thus, under conditions similar to those evaluated in this study, it would not be economically advantageous to maintain New Zealand White rabbits in the herd after 90-d-old and animals at this age should be slaughtered for meat production.
The RGR values decreased with age, with estimates that reached 0.8% and 0.1% for BW and LEA using the von Bertalanffy and Logistic models, respectively, at 90 days of age. RGR is expressed as the proportion of the increase in the animal weight or muscle growth for each day with respect to the predicted weight (i.e., the amount of BW that the animal gained in a certain age in relation to its BW at that age). The RGR values in the IP for each trait were 0.028 and 0.014 for BW (von Bertalanffy model) and LEA (Logistic model), respectively (
Figure 3). Thus, the gain was proportional to 2.8% and 1.4% for BW and LEA, respectively, in relation to body mass and muscling of the animals.