Comparison of non-linear models to describe the growth in the Andalusian turkey breed

Abstract

Seven non-linear growth models were compared in the Andalusian turkey, an endangered native breed. To this aim, turkeys were weekly weighted until they reached 35 weeks. The goodness-of-fit and flexibility criteria of Brody, Von Bertalanffy, Verhulst, Logistic, Gompertz, Richards, and Sinusoidal growth models were evaluated to quantify their ability to describe the biological growth curve. Goodness-of-fit criteria were assessed comparing the mean squared error (MSE) and adjusted determination coefficient (Pseudo-R²), while the flexibility criteria of Akaike (AIC) and Bayesian information criteria (BIC) were evaluated to quantify the explanatory and predictive ability of the models tested. Afterward, all criteria were considered in a combined index to determine the most efficient model to describe and predict growth patterns. The best-fitting model for males growth was Logistic (MSE: 250,349.87; Pseudo-R²: 0.97) which also reported the best explanatory and predictive properties (AIC: 18,949.25; BIC: 18,963.24), while best goodness-of-fit criteria, explanatory and predictive capacity in females were reported for the Richards model (MSE: 144,432.45; Pseudo-R²: 0.95; AIC: 17,529.83; BIC: 17,549.02, respectively). Von Bertalanffy and Richards models underestimated the weight at early age stages, contrary to Logistic and Verhulst models. The asymptotic weight was higher in males than in females at all evaluated models, being 11,085.37 g for Logistic and 5,706.38 g for Richards, respectively. In conclusion, a marked sexual dimorphism is evident, with females reaching maturity earlier than males. The higher inflection point in males may enable their relatively easier commercial standardisation than in turkey hens.

Highlights

• Logistic was the best fitting model for males' growth and Richards for females.

• Females reach maturity earlier than males with higher rates of maturity.

• Andalusian turkey breed shows an intense sexual dimorphism.

Keywords:

• Autochthonous genotypes
• mathematical functions
• poultry
• sex
• zootechny

Introduction

The Andalusian turkey is an endangered Spanish native turkey breed distributed across the Southeast of Iberian Peninsula. As described by Canales Vergara et al. (2019), Andalusian turkeys probably descend from the first imported turkeys from Mexico in the early sixteenth century, during the colonisation period, when conquerors disembarked them in the port of Seville (Spain) upon their arrival in the Old World to then distribute these animals across the villages flanking the Guadalquivir river (Fernández et al. 2009). Nowadays its breeding and husbandry practices related are generally limited to rural areas of Andalusia and Extremadura, where Andalusian turkeys are mainly raised under semi-grazing conditions by backyard producers, taking advantage of their high rusticity, paired with their lightness and agility which confers them a high capacity to feed themselves by foraging the brush for insects and vegetation, diminishing their supplementary food intake requirements. The Andalusian turkey is mainly characterised by its black and roan feather patterns. Traditionally, a meat aptitude has been ascribed to the breed. However, by-product production of eggs weighing between 81 and 86 g (Arando et al. 2019) is regarded as an extra contribution to the familiar economy, out of the breeding season.

However, during the twentieth century, like many other native breeds, their populations significantly decreased as a consequence of the introduction of commercial hybrid lines, which displaced them, and in many cases, derived in the hybridisation seeking improved production parameters (González-Ariza et al. 2019).

One of the main differences between native turkey breeds and commercial meat-type turkey lines is the bodyweight that individuals reach, which is approximately about two or three times lower in native breeds (10 kg in Andalusian male turkey; 5–6.8 kg in Mexican backyard male turkeys and 20 kg in commercial line males) (López-Zavala et al., 2008; Zampiga et al. 2019). These bodyweight variations are usually associated with low growth rates and worst feed conversion in native breeds (Mata-Estrada et al. 2020). However, native turkey breeds have developed a genetic resistance to some diseases as a consequence of greater diversity which provides them with greater resistance than commercial lines (Kamara et al. 2007; Camacho‐Escobar et al. 2008), since heterogeneous populations are less susceptible to viral infections (Mpenda et al. 2019). In addition, in the last years, there has been an increasing interest in food produced under environmentally sustainable conditions, for which native breeds have a great potential (Nogales et al. 2017).

Animal growth is defined as the changes in volume, size, or shape over time until individuals reached maturity (Ulutas et al. 2010). Growth can be represented using models which provide the opportunity to mathematically describe growth patterns, with biologically interpretable parameters which have become a useful tool to estimate growth evolution at different ages (Lupi et al. 2015). Besides, growth parameters are highly heritable (Narinc et al. 2014; Kaplan et al. 2016). Growth describes a sigmoid-type curve with three differentiated stages until mature weight is reached: initial slow growth, steady acceleration, and an afterward period of deceleration (Lawrence et al. 2012).

In addition, fitting growth curves to growth functions not only allows to determine different management factors, such as the most suitable slaughtering age, sexual maturity age, health status, dairy dietary requirements but also it helps to measure the effect of selection on growth traits (Sariyel et al. 2017; Mata-Estrada et al. 2020).

As described by Narinc et al. (2017), the most widely used functions to model the poultry growth performance are Gompertz, Logistic, Richards, and Von Bertalanffy models and have been applied to fit growth in different species, such as guinea fowl (Nahashon et al. 2006), partridge (Sariyel et al. 2017), chicken (Mata-Estrada et al. 2020), ostrich (Ramos et al. 2013), quail (Narinc et al. 2010), and turkey (Porter et al. 2010).

Considering interspecies variability together with the characteristics and limitations of each particular non-linear growth model, it is necessary to carefully assess goodness-of-fit criteria to determine which best model to define the growth performance is. Furthermore, since growth has been reported to be affected by sex, breed, diet, lighting program, or rearing period (Farghly et al. 2018; Ciurescu et al. 2020; Van Staaveren et al. 2020), additional factors must be considered.

Therefore, the present study aimed to evaluate goodness-of-fit and flexibility criteria in seven growth models to determine the non-linear function that best describes and predicts the growth patterns of the Andalusian turkey.

Materials and methods

Animal data and management

This study was conducted in the Agropecuary Provincial Centre of Diputacion of Cordoba (Spain, 37°54′50.9″N–4°42′40.4″W). The eggs used in the study were obtained from the Andalusian Turkey base population (Figure 1) comprising in the conservation program of the breed (18 female and three male, reared in three different groups and whose ages ranged from 12 to 16 months). Base population animals were randomly collected from the backyards of rural areas around Andalusia to gather and maintain the highest levels of diversity considering that the severe endangerment status of the Andalusian turkey population.

Figure 1. Andalusian turkey breed males and females according de feather colour. (A) Black female; (B) Black male; (C) Roan female; (D) Roan male.

The eggs were daily collected during the breeding season (February to April 2019) and were stored at 17 °C and 85% of humidity until they were placed in the incubator (7 days maximum since collection). The eggs were incubated for 26 days in an incubator (Masalles, M240-I, Barcelona, Spain) at 37.2 °C, 55% humidity, and followed an automatic turning programme. On the 26th day of incubation, all eggs were placed in a tray at the hatchery cabinet (Masalles, 25-N HLC, Barcelona, Spain) maintained at 36.7 °C and 60% of humidity until hatching (two more days).

Considering the low egg production rate of turkey hens, four different incubations were performed to ensure that a sufficient number of eggs was obtained for the present study. A total of 140 eggs were incubated out of which 116 finally hatched (58 for each sex) and were considered in the study. Animals were wing-banded and individually weighed for the whole period during which the study took place. As a result, 3,000 observations were obtained once the database was filtered. The data file was purged and values that were twice the standard deviation above and below the average weight were removed to ensure normality of residuals assumption was met (Lupi et al. 2015). As a result of this purging process, data of 98% of the total animals was considered for the statistical analyses performed in the study. As suggested by Pizarro et al. (2019), normality was tested on all the observations, but also across the individuals of the same life status (age group or productive moment). Each bird was weighed at hatching and weekly until it reached 35 weeks of life using an electronic balance (measurement precision = 0.01 g, Cobos, CSB-600C, Barcelona, Spain) for birds weighing up to 600 g. For animals exceeding this weight (>600 g), a suspended electronic balance was used (measurement precision = 5 g, Kern CH50K100, Kern & Sohn, Balingen, Germany). The birds were subjected to the hatchery prophylaxis procedure (Newcastle, fowl pox, chronic respiratory disease, and fowl cholera vaccines). Electric heaters (Copele LGA, Copele, Murcia, Spain) were placed in the room where the animals were hosted for their first six weeks of life to help them to regulate body temperature. Afterward, turkeys were free-ranged (temperature, humidity, and photoperiod) until the end of the experimental period, which lasted for 35 weeks on the whole.

Food and water supplies were provided ad libitum. However, animals were supplemented with two different commercial diets depending on their physiologic status as reported in Table 1. A starter diet was supplied for the first four weeks, while a grower diet was provided to the animals from the end of the fourth week until the experimental period finished. Animal management was conducted by European Union Direction regulations (2010/63/EU) as transposed to Spanish Royal Decree Law 53/2013. The study protocol was submitted to a legally constituted ethics committee and deemed exempt from review (The Ethics Committee of Animal Experimentation of the University of Cordoba, Spain).

Table 1. Recipe and chemical composition of experimental diets.

	Starter (0–4 weeks)	Grower (>4 weeks)
Recipe	Shelled toasted soybeans flour	Wheat
	Flaked wheat	Shelled toasted soybeans flour
	Corn	Corn
	Monocalcium phosphate	Monocalcium phosphate
	Calcium carbonate	Calcium carbonate
	Sodium bicarbonate	Sodium bicarbonate
	Soybean oil	Soybean oil
Analysed chemical composition (%)
Crude protein	20.50	18.00
Crude lipids	2.10	3.00
Crude fibre	2.80	3.00
Crude ash	6.30	5.70
Calcium	0.90	1.00
Phosphorus	0.69	0.50
Sodium	0.15	0.18
Methionine	0.48	0.32
Lysine	1.13	0.90
Analysed gross energy (kcal/kg)	4303.69	4530.20

Model development

Table 2 presents the functions for the seven non-linear models selected in the present study which were namely; Richards, Brody, Logistic, Verhulst (a variant of Logistic model), Gompertz, Sinusoidal, and Von Bertalanffy (Ersoy et al. 2006; Nogales et al. 2017; Kuhi et al. 2019). The fitting procedure was developed using multiple observations of the dependent variable (weight) for each point of the independent variable (day of life). The assumption of residuals normality was tested using Shapiro–Wilk Francia’s tests (recommended for samples from 50 < n < 2,500) with StataCorp Stata version 14.2 (StataCorp, College Station, TX). Autocorrelation analysis was performed by testing the Durbin-Watson method. Residual values derive from the difference between observed and predicted values. Well-fitted models do not describe variable patterns when the residuals are plotted against the fitted values. If the variance of the residuals is non-constant then the residual variance is said to be heteroscedastic. Durbin Watson and Scatter plots were performed using SPSS Statistics for Windows statistical program, Version 25.0.

Table 2. Mathematical description of growth models, biological parameters and growth evaluators.

	Mathematical expression	Inflection weight	Inflection age	Growth rate
Brody	y = a*(1−b*exp(−k*t))	–	–	$V_c=ka\left(1-\frac{y}{a}\right)$
Von Bertalanffy	y = a*(1−b*exp(−k*t))**3	$y_i=\frac{8a}{27}$	$t_i=\frac{ln(3b)}{k}$	$V_c=3ky\left[{\left(\frac{a}{y}\right)}^{1/3}-1\right]$
Verhulst	y = a/(1+b*exp(−k*t))	$y_i=\frac{a}{2}$	$t_i=\frac{ln(b)}{k}$	$V_c=ky\left(1-\frac{y}{a}\right)$
Logistic	y = a*(1 + exp(−k*t))**(−m)	$y_i=\frac{a}{2}$	$t_i=\frac{-ln(2^{\frac{1}{m}}-1)}{k}$	$V_c=mk\frac{y_i}{2}\left(\frac{e^{-kt}}{1+e^{-kt}}\right)$
Gompertz	y = a*exp (−b*exp(−k*t))	$y_i=\frac{a}{e}$	$t_i=\frac{ln(b)}{k}$	$V_c=\mathit{kyln}\left(\frac{a}{y}\right)$
Richards	y = a*(1−b*exp(−k*t))**m	$y_i=a{\left(\frac{m-1}{m}\right)}^m$	$t_i=k^{-1}\mathit{\text{ln mb}}$	$V_c=\frac{0.5amk}{(2m-1)}$
Sinusoidal	y = a1+a2*((sin (ct−θ))^n	a1	$t_i=\frac{\theta }{c}$

y = weight, in g, at age t; t = age in days; a, a₁, a₂, c, n, k, and m = parameters.

Weight data were used to fit growth models. Due to the detection of high numbers of outliers, Pseudo-R² was used to correct for the potential inflation of variance explanatory ability to avoid the use of data that may distort the results (Lupi et al. 2015). The curve shape parameters considered (three for all models except from Richards model with four) were as follows: ‘a’, mature or asymptotic weight reached when time tends to infinity; ‘b’, integration constant that allows to calculate the inflexion age, that allows to calculate the inflexion age; ‘k’, which represents the exponential growth rate and ‘m’ which denotes the shape of the growth curve and determines the inflection point, when the animal starts the auto-deceleration stage until it reaches adult (Darmani Kuhi et al. 2010). For the sinusoidal model, parameters were: ‘a₁’, weight at inflexion; ‘a₂’, the overall rate parameter of the total growth development; ‘θ’, the phase shift (Kuhi et al. 2019). The most appropriate model was chosen as suggested by Lupi et al. (2015) and Pizarro et al. (2020a ) considering the following best fitting and flexibility criteria:

1. The lowest value of mean square error (MSE) for each corresponding function, as a measure of the variability of other potentially non-considered factors.

2. The highest value of the adjusted determination coefficient (Pseudo-R²), taking into account that the determination coefficient (R²) in non-linear models could overestimate results (Lupi et al. 2015). Pseudo-R² can be computed as: $$\mathrm{Pseudo}-{\mathrm{R}}^2=1-\frac{\mathit{SSResidual}}{{\mathit{SSTotal}}_{\mathit{Corrected}}}$$

3. The lowest value of the Akaike information criterion (AIC). This coefficient was estimated as AIC = n * ln (SSE/n) +2p, where n is the number of observations; SSE is the sum of squares of the error and p is the number of parameters + 1. AIC coefficient estimates variations in the quality of fit for present data and the differences in the number of parameters between both models.

4. The lowest value of Bayesian information criterion (BIC, Leonard and Hsu 2001), which is a model-order selection criterion: $\mathit{BIC}=n\mathit{ln}\left(\frac{\mathit{SSResidual}}{N}\right)+k*\mathit{ln}\left(n\right),$ where n is the number of observations; SSResidual is the sum of squares of the residuals and k is the number of parameters.

The goodness-of-fit criteria of Pseudo-R², MSE, and the flexibility criteria, of AIC and BIC were individually considered and used to rank the seven fitting models depending on their best performance. Those models presenting a higher ability to capture data variability (Pseudo-R²) and better explanatory (AIC) and predictive properties (BIC and MSE) scored higher in the rank. Computationally, those models presenting the lowest scores for AIC, BIC, and MSE, and the highest scores for Pseudo-R² presented higher scores in the rank. When the same value for a specific criterion was reported for two or more models, these were ranked equally.

Afterward, as goodness-of-fit and flexibility criteria may differ in terms of which their most desirable values are and what their magnitude is, a combined selection index (ICO) was developed following the premises in Van Vleck (1993) to summarise the position in the rank for each of the goodness-of-fit and flexibility criteria determined for each model. The combined index used (ICO) was as follows: $$\mathit{ICO}=\frac{\mathit{Pseudo}-R^2\mathit{Rank}\mathit{Position}*W1+\mathit{MSE}\mathit{Rank}\mathit{Position}*W2+\mathit{AIC}\mathit{Rank}\mathit{Position}*W3+\mathit{BIC}\mathit{Rank}\mathit{Position}*W4}{4},$$ where, W1: the weight for Pseudo-R² rank position; W2: the weight for MSE rank position; W3: the weight for AIC rank position;W4: the weight for BIC rank position. All criteria were given the same relevance in the ICO, hence, no coefficient was used, that is the proportion of 1: 1: 1: 1 was followed. The model presenting superior ICO values was selected as a model with the best fitting, explanatory and predictive properties.

Statistical analysis

NLIN procedure of the SAS statistical package (v. 9.0 software, SAS Institute Inc., Cary, NC, USA) was used for data processing, using individual observations to fit the population and considering the Marquardt iterative method as the estimation procedure to reduce the residuals and improve the goodness-of-fit.

Results

Residuals for growth observations of all animals and animals of the same age were normally distributed (p > .01) and heteroscedastic p < .01). Consequently, fitting non-linear regression models were suggested as the most appropriate approach. A summary of the results derived from biological curve data model fitting and flexibility criteria evaluation is reported in Table 3. As it can be observed, the adjusted determination coefficient (Pseudo-R²) was slightly higher in males than in females across all studied models, with similar values being reported between models within each sex, except for the Brody model, for which the results are lower. After the goodness-of-fit and flexibility criteria were considered and the ICO index was used to rank the models, the best goodness-of-fit criteria were reported for Richards model in the case of females, while for males, the best goodness-of-fit criteria were reported by the Logistic model. For males, Gompertz model offered similar Pseudo-R² values to those reported by Logistic model, however, MSE, AIC, and BIC values in the later were slightly higher, while the Richards model did not converge with data. Similarly, the Sinusoidal model did not converge for both sexes and was discarded.

Table 3. Biological growth curve estimation parameters for each model by sex in Andalusian turkey.

SEX	Model	a (s.e)	b (s.e)	m (s.e)	k (s.e)	Pseudo-R^2	MSE	AIC	BIC	ICO
F	Brody	16,447.92 (1,702.55)	1.024 (0.004)	–	0.0018 (0.0002)	0.93	182,918.84	17,877.27	17,891.16	1.00
	Von Bertalanffy	5,706.38 (81.07)	0.943 (0.020)	–	0.0130 (0.0004)	0.94	144,569.52	17,530.23	17,544.12	4.75
	Verhulst	4,841.85 (36.57)	22.596 (1.068)	–	0.0296 (0.0006)	0.94	160,290.50	17,682.49	17,696.38	2.00
	Logistic	5,192.19 (50.00)	–	6.003 (0.140)	0.0193 (0.0004)	0.94	146,626.63	17,551.07	17,564.96	3.00
	Gompertz	5,328.39 (58.58)	4.625 (0.121)	–	0.0171 (0.0004)	0.94	145,118.49	17,535.82	17,549.71	4.00
	Richards	5,553.76 (119.71)	0.727 (0.173)	4.597 (1.604)	0.0144 (0.0010)	0.94	144,432.45	17,529.83	17,549.02	5.00
M	Brody	1,618,960.42 (10,126,234.52)	1.001 (0.004)	–	0.0000 (0.0002)	0.95	443,993.79	19,822.43	19,836.42	1.00
	Von Bertalanffy	13,304.69 (255.51)	0.989 (0.015)	–	0.0094 (0.0002)	0.97	258,779.78	18,999.73	19,013.71	2.75
	Verhulst	9,716.63 (75.16)	34.124 (1.162)	–	0.0265 (0.0004)	0.97	268,406.80	19,055.39	19,069.38	2.00
	Logistic	11,085.37 (124.53)	–	6.847 (0.110)	0.0155 (0.0002)	0.97	250,349.87	18,949.25	18,963.24	4.50
	Gompertz	11,559.01 (152.15)	5.270 (0.097)	–	0.0137 (0.0003)	0.97	250,764.45	18,951.78	18,965.76	3.75

F: female; M: male; s.e.: standard error; pseudo-R²: non-linear determinative coefficient; MSE: mean square error; AIC: Akaike information coefficient; BIC: Bayesian information criteria; ICO: combined selection index. Richards model for males and Sinusoidal model for both sexes did not converged.

In regards to curve shape parameters present in the mathematical expressions, mature or asymptotic weight ‘a’ was higher in males than in females in all evaluated models, being 11,085.37 and 5,553.76 g, respectively, when best-fitting models were considered. However, the ‘a’ values obtained when Brody model was fitted did not report any coherence with the observed weights, hence this model was discarded.

With respect to the ‘b’ parameter, in terms considerably higher values were generally observed in Verhulst model. By contrast, ‘k’ parameter values were higher for females in all studied models.

The observed and predicted weights for the different models fitted are presented in Figure 2. As it can be appreciated, observed weights of males were higher than those of females, except for hatching (58.7 and 59.7 g, respectively). Predicted weights by Richards model (Figure 2), underestimated females weight at early stages (14.1 g at hatching and 40.3 g at day 7), same as Von Bertalanffy model (0 g at hatching and 5.5 g at day 7 for females and 1 g at hatching and 15.3 g at day 7 for males), By contrast, the Logistic and Verhulst models overestimated weight at early stages (Figure 2). Contrarily, when mature predicted and observed weights were compare, all models overestimated these weights, except Verhulst model for males.

Figure 2. Growth curves for both sexes of Andalusian turkey breed, observed and predicted with each model (F: female; M: male).

Finally, information about the inflection weight and age, and the growth rate at this specific moment for the studied models are presented in Table 4. Considering the most appropriate models, the inflection point was achieved at 12 weeks (83.68 days) and 1,798.26 g in females and 20 weeks (144.08 days) and 5,542.69 g in males, observing a difference in the time to reach an inflection weight of 61 days between both sexes.

Table 4. Inflection weight, inflection age and growth rate at the inflection point of fitting models for the biological curve of Andalusian turkey breed.

Model		Inflection weight (g)	Inflection age (days)	Growth rate (g/day)
Brody	F	–	–	–
Brody	M	–	–	–
Von Bertalanffy	F	1,690.78	79.96	32.99
Von Bertalanffy	M	3,942.13	115.39	55.73
Verhulst	F	2,420.93	105.3	35.84
Verhulst	M	4,858.32	133.40	64.28
Logistic	F	2,596.10	108.80	16.41
Logistic	M	5,542.69	144.08	28.39
Gompertz	F	1,960.20	89.34	33.60
Gompertz	M	4,252.32	121.65	58.09
Richards	F	1,798.26	83.68	22.47
Richards	M	–	–	–
Sinusoidal	F	–	–	–
Sinusoidal	M	–	–	–

F: female; M: male.

Discussion

Model choice decision rationale

Non-linear models have been widely contrasted as suitable methods to fit for growth in turkey hens and sires (Porter et al. 2010; Pérez-Lara et al. 2013; Segura-Correa et al. 2017). However, other alternatives, such as mixed models with random and fixed effects, have been formerly suggested in the literature to fit for the same aim (Testik et al. 1999). Mixed models (random and fixed) have several desirable properties, but their use needs to stick to certain assumptions (El Halimi 2005).

As suggested by Harrison et al. (2018), mixed models, either it is using fixed or random effects, need at least five ‘levels’ (groups) for a random intercept term to achieve robust estimates of variance. This was also supported by Gelman and Hill (2007) and Harrison (2015), who suggested the fact that fixed or random effects, with lower than five levels may not accurately estimate between population variance, due to variance estimate collapsing to zero, which in turn makes the model equivalent to non-linear modelling (Gelman and Hill 2007) or be non-zero, but incorrect when the small number of levels (groups) from which samples were taken is not representative of the true distribution of means. This was also addressed to be a potential source for variance and covariance distortion as a direct consequence of too few levels of a fixed or random effect as suggested by (Silk et al., 2020). Furthermore, fixed effects are generally inconsistent in the non-linear model as n grows (Newey 2007). This rationale may make it compulsory to apply integration methods, such as the adaptive Gaussian Quadrature method to prevent distortion from occurring (Afrouziyeh et al. 2021), which in turn increases the complexity of the methods used.

In our study, sex only accounts for two levels, strongly preventing against the use of non-linear mixed (random) models due to the aforementioned reasons. Sex must be considered a fixed effect rather than a random effect. In this context, randomising sex may lead to model degeneracy, biased estimation of random effect variance, inaccurate estimation of random effect variance, and high error potential for questions related to random effects. The random effects of other factors, such as the animal, could have been considered as well. However, as the same animals were repeatedly measured along the course of the study, those random effects corresponding to the same animal are correlated, hence, this acts against the fulfillment of the independence of observations assumption (Chan and Kuk 1997). As a result, and given the only additional factor considered was sex and animals were measured repeatedly, we decided to run non-linear procedures for both levels (male and females) separately.

Second, the present study aims to serve as a model for the evaluation of growth in native turkey populations. Randomly mixed models can be unstable if sample sizes across groups are highly unbalanced, which is likely to occur in native poultry populations for which sex censuses may be highly skewed, with male groups containing very little data or for which female/male ratio favours one of the sexes as it has been reported in the literature (Mata-Estrada et al. 2020; González-Ariza et al. 2021), for which separate sex modellisation has been suggested the most appropriate alternative.

Third, although fitting a factor as a fixed effect may provide a statement of the significance variation (differences across groups) in a relatively easy manner, testing differences among levels of a random effect is much more difficult for frequentist analyses, which may call for the application for Bayesian inference approaches to be considered (Kéry 2010). Last but not least, incorrect parameterisation of the random effects in the model could yield unreliable model estimates. Among others, failing to identify dependency structures, that is meet the assumption of non-independence (multiple clutch measures from a single bird), failing to specify random slopes to prevent constraining slopes of predictors to be identical across clusters in the data (Barr et al. 2013); and testing the significance of fixed effects at the wrong level may eventually lead to pseudoreplication and inflated Type I error rates which could be prevented by using residual degrees of freedom for fixed effects.

As above explained, non-linear models were chosen to fit for the time evolution of growth in male and female turkeys. No additional term was considered in the model, such as sex as a fixed effect, but sexes were modelled separately. The different models used comprised a variable number of curve shape parameters, thus they differed in terms of model complexity. The number of parameters considered in a model and the nature of its mathematical function may indeed be responsible for the better fitting properties of certain models against others as suggested by Pizarro et al. (2020a, 2020b). Consequently, and given ICO was used to sort models depending on their better fitting and flexibility properties, the discussion was divided into the following two sections.

Fitting and flexibility criteria model comparison

The present results agree with those in Ersoy et al. (2006), who reported Richards model to be the best-fitting model for turkey growth. Ersoy et al. (2006) suggested male turkeys matured more slowly and reached mature body weight later than female turkeys. Determination coefficient values were high (>0.98) but slightly differed across sexes. The same authors reported residual mean square error (RMSE) to be higher for male turkeys, for which a lower record number was used. By contrast, Richards model fails to converge for males in the present study, which resembles which was previously reported by other authors, such as Segura-Correa et al. (2017), for male and female hybrid turkeys and Porter et al. (2010) for turkey hens who considered larger samples during longer study periods. This suggests sample limitations in the aforementioned studies may have been a source for erroneous convergence in Richards model.

The Sinusoidal model did not converge for both sexes in the present study, even if contextually, recent studies conducted in turkeys have reported the Sinusoidal model to be the best-fitting model for growth (Kuhi et al. 2019). These differences may be a consequence of the properties of the data that was used, since Kuhi et al. (2019) used average values, while in the current study, the results were obtained from individual observations taken from the studied population.

Recent studies in turkeys reported better fitting properties, such as for Morgan Mercer Floding, polynomial regression, or Von Bertalanffy (Porter et al. 2010; Pérez-Lara et al. 2013; Segura-Correa et al. 2017). Such differences may derive from the specific growth patterns described by the species used in each particular case, breeds considered, or the rearing period during which the study took place (Segura-Correa et al. 2017).

Observed and predicted curve shape parameter comparison

Asymptotic weight (‘a’) in Andalusian turkeys was around 50% lower in females at 35 weeks of age, which suggested an intense sexual dimorphism across sexes, with females maturing faster than males, respectively. Similarly, recent studies reported the same trend in native Mexican bronze turkeys at 22 week of age. Specifically, a female’s body weight was reported to be 75% lower than that of male turkeys at the same age. (Pérez-Lara et al. 2013). The same trend was reported for mature weight, with it being 50–85% lower in females than in males.

In addition, native breeds’ weights at 35 weeks (5,553.76 and 11,085.37 g, for females and males, respectively) were lower than those of hybrid genotypes (14.63 and 25.88 g, at 22 weeks for female and males, respectively) (Segura-Correa et al. 2017). Higher mature weights in hybrid genotypes may be a consequence of intensive genetic selection.

As it can be inferred from the results, ‘a’ values when Brody model was fitted did not agree with observed values, hence the model was discarded. This could have been expected as Brody model has been addressed to be effective when describing growth patterns when growth rate starts at a maximum point to then decrease.

In this regard, recent studies have reported Brody model not to converge (Segura-Correa et al. 2017), which suggested this model not to be appropriate to estimate growth performance in turkeys. Similarly, Brody model has been deemed as inappropriate to fit for growth performance in other species (Lupi et al. 2015). Conversely, Brody model may be rather efficient to accurately predict growth in other species at older ages (Gómez et al. 2007), such as cattle (Nogales et al. 2017).

Most of the models used to fit for growth in poultry, are asymptotic and have a sigmoid structure (Yang et al. 2006). Brody model was typically exponential, while Gompertz, Logistic, Von Bertalanffy and Richards models are classical sigmoid functions (Gbangboche et al. 2008; Ware and Power 2017). Growth accelerates during the early stages of life, thus a greater increase in weight occurs which progressively decreases as the animal approaches adulthood. As a result, when the growth curve is graphically described, an ascending sigmoid line can be observed.

In this regard, the growth rate becomes altered as the animal reaches its adult size, which translates into an appreciable change in curve shape where the highest growth rate occurs (inflection point). Afterward, the growth rate decreases and stabilises describing a continuous trend, which mathematically coincides with a horizontal asymptote (Narinc et al. 2017; Nogales et al. 2017).

Von Bertalanffy presented the highest ‘a’ values for both sexes, while Verhulst reported the lowest ones. Similar results were reported by Segura-Correa et al. (2017), who observed higher mature weights when Von Bertalanffy model was fitted. By contrast, Şengül and Kiraz (2005) described higher ‘a’ values when Verhulst model was used. Verhulst model presented higher values of ‘b’ parameter as previously reported by Segura-Correa et al. (2017).

Another evidence for sexual dimorphism in turkeys was the maturity rate ‘k’ found in our study. Higher values of ‘k’ are indicative of turkeys reaching maturity earlier. As described by several authors and agreeing with the results derived from the present study, the highest values for ‘k’ parameter were reported in females (0.0296) and males (0.0265) in this study when Verhulst model was fitted, while the lowest values were reported for Von Bertalanffy model (0.0130 and 0.0094, for females and males, respectively) (Şengül and Kiraz 2005; Segura-Correa et al. 2017).

When Andalusian turkey’s observed weights at 35 weeks of age were compared to those in other native breeds, such as Mexican bronze turkeys (hatching weight 46.7 g males/43 g females; mature weight 5,827 g males/3,487 g females), hatching and mature weights were higher (58.7 g males/59.7 g females; 9,465 g males/4,762 g females) (Pérez-Lara et al. 2013).

Richards and Von Bertalanffy models tended to underestimate female and conjoined (female and males) predicted hatching body weight, respectively. This supports the results reported by Segura-Correa et al. (2017).

As described by Segura-Correa et al. (2017) while Von Bertalanffy model tends to underestimate hatching body weight in hybrid turkeys, the Logistic and Verhulst models tend to overestimate hatching body weight. This different behaviour may derive from the similarities and differences in their mathematical functions. For instance, Richards, Von Bertalanffy, Logistic, and Verhulst fix the origin through which the curve originates. In this context, the Gompertz model reported the best-predicted values for weight at hatching across sexes, and does not fix the curve to the origin, but uses the ‘a’ parameter instead.

By contrast, some studies conducted in turkeys and non-selected native chicken breeds suggest Gompertz model overestimate hatching body weight (Magothe et al. 2010; Segura-Correa et al. 2017; Mata-Estrada et al. 2020). Additionally, for females, a slight decrease in observed weights at 210 days of age (30 weeks) was observed. Such a slight decrease probably coincided with the beginning of the laying period of turkey hens, when the energetic requirements of females likely increased. For instance, the beginning of the laying period in native turkeys has been reported to occur later than in commercial turkey hen lines, for which the earliest ages for egg production that have been reported have been dated around 22–25 weeks (Siopes 2010), indicative of the lower precocity found in native breeds.

The results showed sharp inflection weight, ages, and growth rate fluctuations across sexes and models, with sharper inflection weights being reported in males for the Logistic model. Exceptionally, Verhulst model reported superior growth rates than the rest of the models fitted. This may evidence females may reach puberty before males, even if it is at lower weights. These results agree with those reported in recent studies conducted in hybrid and native turkey breeds (Ersoy et al. 2006; Pérez-Lara et al. 2013; Segura-Correa et al. 2017), and other species (Fernandes et al. 2020).

The economic importance of the inflection point bases on the fact that it represents the point of maximum growth rate and marks the moment when growth acceleration ceases. This point frequently coincides with puberty and its control in commercial turkey hens is crucial. The excessive energy production requirements derived as a consequence of fat tissue deposition collaterally occur in cases of premature puberty in turkey hens, which is associated with fertility disorders and deregulation in the follicle hierarchy (Dewez et al. 2018).

Therefore, it is essential to define this inflection point to design adequate feeding programmes or to accurately set slaughtering ages. The definition of slaughtering age is closely linked to turkey meat quality traits, such as pH, colour, and water holding capacity. As described by some authors, animals slaughtered at a lower age presented higher breast pH values, less dark, and lower drip and cooking losses (Sarica et al. 2011; Mauri et al. 2017). These variations in meat quality are associated with myoglobin concentration, collagen cross-links, or muscle fibre diameter (Lepetit 2008; Listrat et al. 2016).

In Andalusian turkeys, the defined slaughtering age for males at 20 weeks and 5,542.69 g could be an interesting alternative, which may recommend the commercialisation of turkey carcasses as whole individual pieces, provided their size may be suitable to fit in a home oven, hence making the product rather accessible.

However, lower weights reached by females at the inflection point (12 weeks and 1,798.26 g) may reduce their market opportunities, since their lower carcase weights may be insufficient to seek turkey hens standardised commercialisation. Consequently, the production of females may better focus on these females producing the individuals that supply for multiplication farms, as it occurs in other poultry meat lines.

Conclusions

In conclusion, limited samples may be a source for erroneous convergence in Richards model, which especially affects growth in turkey males. Differences in Sinusoidal model fitting may also be ascribed to how data is collected or considered. Contextually, when individual weight observations are used against when population data is averaged. A marked sexual dimorphism is evident, as denoted by higher values of maturity rate ‘k’, with females reaching maturity earlier but with lower weights than males. Lower mature weights may be a consequence of the lack of intensive genetic selection in local breeds. Models presenting a Sigmoid profile better graphically described growth in turkeys than exponential ones. These models represent growth acceleration in early stages as individuals approach adulthood up to the points of maximum weight. Afterward, the growth rate asymptotically decreases and stabilises. Underestimation or overestimation of hatching body weight may derive from origin being fixed or ‘a’ parametrised in functions. The higher inflection point in Andalusian male turkeys may enable their relatively easier commercial standardisation than in Andalusian turkey hens.

Ethical approval

The study follows the national guidelines and premises described in the Declaration of Helsinki. The Spanish Ministry of Economy and Competitivity through the Royal Decree Law 53/2013 and its credited institution the Ethics Committee of Animal Experimentation from the University of Córdoba permitted the application of the protocols. This Decree Law follows the European Union Directive 2010/63/UE, from the 22nd of September of 2010.

Acknowledgements

The authors would like to acknowledge the Association of Andalusian Turkey Breeders, Agropecuary Provincial Centre of Diputacion of Cordoba, and PAIDI AGR 218 research group for providing the animals and their support.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

The authors confirm that the data supporting the findings of this study are available within the article.