Heterosis in the components of lactation curves of Girolando cows

Abstract

The objective of this study was to evaluate the effect of heterosis on the lactation curve components of Girolando cattle obtained by fitting different mathematical models. Data consisted of 258,891 test day milk yield records of the first lactation from 37,965 cows of Minas Gerais State (Brazil) between 1998 and 2014. Those cows were from the Holstein breed (H), Gyr breed (G) and six genetic cross-breedings of Holstein × Gyr, (1/4H, 3/4G (1/4H); 3/8H, 5/8G (3/8H); 1/2H, 1/2G (1/2H), 5/8H, 3/8G (5/8H); 3/4H, 1/4G (3/4H); 7/8H, 1/8G (7/8H)), which is officially named as Girolando breed in Brazil. The Wood’s linear model (WD_lin), Wood’s non-linear model (WD_nlin), Wilmink’s model (WL) and Ali and Schaeffer’s model (ASH) were used for estimating the peak milk yield (PY), time to peak yield (PT), 305-day milk yield (TMY) and four different persistency measures (P, P_2:1, P_3:1 and P_3:2). Regardless of the fitted model, the highest estimates of PY and TMY were for the H group. The heterosis effect was significant (p < .001) for TMY and all components of the lactation curve, except for P_2:1. Girolando cattle presented a heterosis effect of 12.30% and 13.03% for PY and TMY, respectively. The magnitude of heterosis effect was larger for PT (24.18%), whereas the different persistency measures presented the smallest magnitude of heterosis values. The producers may use the different genetic groups to benefit from the heterosis mainly for the time to peak, peak yield and 305-day milk yield.HIGHLIGHTS• Girolando cows in production systems in Brazil has shown positive results. • The greater productive efficiency may be because of heterosis on the production.• Thus effect of heterosis on lactation curve can contribute to the improvement of the production system.

Keywords:

• Cross-breeding
• dairy cows
• mathematical model
• persistency of lactation

Introduction

In tropical countries, such as Brazil, cattle breeding is based mainly on cross-bred animals from Bos taurus taurus x Bos taurus indicus (Silva et al. 2015). The objective of cross-breeding is to obtain more adapted animals to tropical conditions by means of heterosis and breed complementarity (Canaza-Cayo et al. 2014). Thus, around the 1940’s Girolando Cattle arose when Brazilian farmers began to obtain different genetic compositions between Holstein and Gyr breeds (Facó et al. 2002). The individual performance of animals with different genetic compositions will depend on the prediction of the shape of the lactation curve (Pereira et al. 2016). The pattern of the lactation curve is very important in the study of dairy cattle (Bangar and Verma 2017). In general, the common pattern of lactation curves is that milk yield increases in the early stage of lactation reaches a peak yield and then decreases slowly until the end of lactation (Chang et al. 2001). The lactation curve presents three important components, namely milk yield at peak, time to peak yield and persistence of lactation (Wasike et al. 2014).

The use of mathematical models to describe the shape of lactation curves in genetic programmes allows establishing strategies to optimise selection of more efficient genotypes for the farmer in several production systems (Oliveira et al. 2007; Hossein-Zadeh 2017). In the last decade, many mathematical models were developed to describe milk yield along lactation (Wood 1967; Ali and Schaeffer 1987; Wilmink 1987). Those models were classified into linear and non-linear models. In linear models, the parameters are linear functions of days in milk or of their transformation, and can thus be estimated by linear regression (Torshizi et al. 2011). Non-linear models are non-linear in the parameters and they are more statistically and computationally demanding, due to non-linearity than their linear equivalent models (Glória et al. 2012). The advantage of such models is that the parameters represent the ascending or descending parts of lactation curve (Varanis et al. 2016), which allow attributing a biological interpretation to each parameter (Glória et al. 2012).

Considering that the shape of lactation curves is different between animals and breeds, obtaining the parameters of the curve is a requisite for the study of heterosis on the components of the lactation curve. The term heterosis, which was firstly described in 1914, may be defined as the increased vigour of cross-bred animals in comparison to the average performance of parental breeds (Shull 1948). When heterosis is correctly used it may be a bonus in the genetic gain of traditional breeding schemes of livestock species (Sorensen et al. 2008). In general, heterosis for total milk yield may range from 17.3% (Rege 1998) to 28% (Cunningham and Syrstad 1987) in crosses between European breeds and zebu. However, there is a lack of studies about the effect of heterosis on the components of the lactation curve.

The objective of this study was to use different models to better describe the lactation curve to evaluate the effect of heterosis on the components of the lactation curve of Girolando cows.

Materials and methods

Data

Data were from the Association of Holstein Breeders of Minas Gerais (ACGHMG), the Brazilian Association of Dairy Gyr (ABCGIL) and the Brazilian Association of Girolando (ABCG).

The data set consisted of 258,891 test day milk yield records of the first lactation from 37,965 cows of Minas Gerais State (Brazil) between 1998 and 2014. Those cows were from the Holstein breed (H), Gyr breed (G) and six genetic cross-breedings of Holstein × Gyr (1/4H, 3/4G (1/4H); 3/8H, 5/8G (3/8H); 1/2H, 1/2G (1/2H); 5/8H, 3/8G (5/8H); 3/4H, 1/4G (3/4H); 7/8H, 1/8G (7/8H)), which is officially named as Girolando breed in Brazil. A minimum of 4 and maximum of 10 test days, obtained from 5 to 305 days in milk, were considered for estimating the lactation curves by the mathematical models. Abnormal yield values or outliers were checked by graphical techniques such as normal probability plots and boxplots, as well as by median, mean, mode, skewness and kurtosis. Thus, the test day milk yield and the 305-day milk yield records were removed if milk yield were out of the range from 3 to 45 kg and from 686.07 to 11026.40 kg, respectively. The descriptive statistics of the edited data are presented in Table 1.

Table 1. Description of the database.

Genetic group	Proportion of Holstein gene, %	Number of herd	Number of cows	Number of test-day records	TMY^a, kg
H	100.00	288	17,148	99,802	7537.98
1/4H	25.00	119	392	2820	4794.49
1/2H	50.00	304	3712	27,473	5439.11
3/4H	75.00	290	4667	37,549	5380.94
3/8H	37.50	88	379	2914	4610.39
5/8H	62.50	280	3830	29,397	4661.99
7/8H	87.50	144	1212	9868	5375.60
G	0	92	6625	49,068	3767.69

^a305-day milk yield from the database.

H: Holstein breed; G: Gyr breed; 1/4H, 3/8H, 1/2H, 5/8H, 3/4H, 7/8H: six genetic cross-breedings of Holstein X Gyr which is officially named as Girolando breed in Brazil.

Lactation curve models

The linear and non-linear models used to fit test day milk yield along lactation of Holstein, Gyr and Girolando breeds were:

• (1) Wood’s non-linear model (WD_nlin) (Wood 1967): $$Y_t={at}^b{e}^{-ct}.$$

• (2) Wilmink’s model (WL) (Wilmink 1987): $$Y_t=a+\mathrm{ }{be}^{-kt}+ct.$$

• (3) Wood’s linear model (WD_lin) (Wood 1967): $$Y_t=\log \mathrm{ a}+\mathrm{b }\mathrm{log }t-ct.$$ In the models, Y_t is the average daily yield in the tth test day of lactation; a is the initial milk yield just after calving; b is the ascending slope parameter up to the peak yield; c is the descending slope parameter and t is the length of time since calving. The constant value k for the WL model was determined in the previous analysis and fixed as 0.05.

• (4) Ali and Schaeffer’s model (ASH) (Ali and Schaeffer 1987): $$Y_t=a+b\left(\frac{t}{305}\right)+{c\left(\frac{t}{305}\right)}^2+d\mathrm{ ln}\left(\frac{305}{t}\right)+f{\text{ ln}}^2\left(\frac{305}{t}\right).$$

In this model, Y_t is the average daily yield in the tth test day of lactation, a is associated with peak yield, b and c are associated with the decreasing slope, d and e are associated with the increasing slope. After the estimation of the parameters of each mathematical model, the mean peak yield (PY), the mean peak time (PT) and the mean persistency (P) were calculated. For persistency, the other three different measures (P_2:1, P_3:1 and P_3:2) were also calculated as proposed by Johansoon and Hansson (1940): $$P_{2:1}=\frac{\text{Milk yield between }101\text{ and }200\text{ days after parturition}}{\text{Milk yield in the first }100\text{ days of the lactation}}\times 100$$ $$P_{3:1}=\frac{\text{Milk yield between }201\text{ and }300\text{ days after parturition}}{\text{Milk yield in the first }100\text{ days of the lactation}}\times 100$$ $$P_{3:2}=\frac{\text{Milk yield between }201\text{ and }300\text{ days after parturition}}{\text{Milk yield between from }101\text{ to }200\text{ days after parturition}}\times 100$$

The different measures of persistency (P_2:1, P_3:1 and P_3:2) use the ratios between the average milk yields obtained in different parts of lactation and all three measures were expressed as percentage.

The average predicted 305-day milk yield (TMY) was estimated for each genetic group in the respective model using the equation of Vargas et al. (2000): $$\text{TMY}=\sum _{t=5}^{305}y\left(t\right),$$ where TMY denotes the average predicted 305-day milk yield and y (t) is the milk yield at day t (5, 6, 7, 8, 9 …, 305) estimated by the corresponding mathematical model.

Breed and heterosis effects

The proportion of genes from the ith breed was calculated for each cow as: $${\alpha }_i^p=\mathrm{ }\frac{\left({\alpha }_i^s+\mathrm{ }{\alpha }_i^d\right)}{2},$$ where ${\alpha }_i^p$ is the proportion of genes from breed $i$ in the progeny, ${\alpha }_i^s$ is the proportion of breed $i$ in the sire, and ${\alpha }_i^d$ is the proportion of breed $i$ in the dam.

Coefficients of specific heterosis in a given cross were calculated between any pair of the dairy breeds using the following identity (Dickerson 1973): $${\delta }_{ij}^p=\mathrm{ }{\alpha }_i^s{\alpha }_j^d+{\alpha }_j^s{\alpha }_i^d,$$ where ${\delta }_{ij}^p$ is the coefficient of expected heterosis between fractions of breeds $i$ and $j$ in the progeny, ${\alpha }_i^s$ and ${\alpha }_j^s$ are proportions of breeds $i$ and $j$ in the sire, and ${\alpha }_i^d$ and ${\alpha }_j^s$ are proportions of breed $i$ and $j$ in the dam.

Those specific effects of heterosis were used for the six genetic groups of Girolando because the distribution of cows across classes of coefficients of expected heterosis was suitable for this purpose (Penasa et al. 2010). The coefficient of general heterosis for each cow was obtained by summing coefficients of specific heterosis previously calculated.

Statistical analyses

Each model was fitted to test day milk yield records using NLIN, REG and AUTOREG procedures in SAS (Statistical Analysis System, version 9.3, Cary North Caroline, USA, 2002). The non-linear models were adjusted to the milk yield records as the iteration method of Gauss-Newton.

The models were tested for goodness of fit using the root means square error (RMSE), Akaike’s information criterion (AIC) and Bayesian information criterion (BIC).

RMSE was calculated as follows: $$\text{RMSE}=\sqrt{\frac{\text{RSS}}{n-p-1}},$$ where RSS is the residual sum of squares, n is the number of observations (data points) and p is the number of parameters in the equation.

AIC was calculated using the following equation: $$\text{AIC}=n\times \ln \left(\text{RSS}\right)+2\mathrm{p}.$$

Also, BIC was calculated as follows: $$\text{BIC}=n\text{ ln}\left(\frac{\text{RSS}}{n}\right)+p\mathrm{ ln}\left(n\right).$$

Smaller numerical values of RMSE, AIC and BIC indicate the better fit when comparing the different models.

The heterosis effect of the components of the lactation curve in Girolando cattle was estimated by the MIXED procedure in SAS. The heterosis effects were obtained after fitting the following mixed linear model: $$Y_{jkl}\mathrm{ }=\mathrm{ }\mathrm{\mu }+H_j+C_k+\sum _{q=1}^2{{\phi }_{q}a}^q+\beta f+\lambda h+e_{jkl},$$ where

• $Y_{jkl}$ is the observation $l$ taken in cow $k$ and herd $j$.

• $\mu$ is an constant.

• $H_j$ is the random effect of herd $j$.

• $C_k$ is the random effect of cow $k$.

• ${\phi }_q$ are regression coefficients associated with the linear ($q=1)$ and quadratic ($q=2$) effects of age of cow.

• $\beta$ is the regression coefficient associated with the linear effect of proportion of Holstein (f).

• $\lambda$ is the regression coefficient associated with the linear effect of heterosis (h) between Holstein and Gyr.

• $e_{jkl}$ is residual random error associated to observation $Y_{jkl}$.

Results

The parameters estimated for the different genetic groups by the various non-linear and linear models, and the criteria for evaluating the quality of fit of these models, are in Tables 2 and 3. WD_nlin estimated positive values for the a, b and c parameters, with b and c close to zero. The WL model presented positive values for parameter a and negative values for parameters b and c in all genetic groups. In linear models, WD_lin estimated positive values for the parameters a and b, and negative but near to zero for parameter c. The ASH model estimated the highest values for the parameter a when compared to the other models, negative values for b and d and values close to zero for f parameter.

Table 2. Estimated parameters (mean ± SE) from different non-linear models, root mean square error (RMSE), Bayesian information criterion (BIC) and Akaike information criterion (AIC).

Breed (composition)	Model	Parameters			RMSE	AIC	BIC
		a	b	c
1/2H	WDnlin	12.989 (0.308)	0.128 (0.006)	0.001 (0.00007)	6,63	387173,90	106398,20
	WL	21.296 (0.112)	–6.850 (0.374)	–0.018 (0.0006)	6,64	387193,50	106417,80
1/4H	WDnlin	11.210 (0.917)	0.134 (0.024)	0.001 (0.0002)	6,65	33335,97	10950,33
	WL	18.897 (0.345)	–6.866 (1.128)	–0.017 (0.0020)	6,64	33332,58	10946,95
3/4H	WDnlin	11.678 (0.227)	0.156 (0.005)	0.002 (0.00005)	6,07	534781,30	139288,10
	WL	21.353 (0.088)	–8.228 (0.297)	–0.020 (0.0004)	6,07	534834,70	139341,60
3/8H	WDnlin	12.623 (1.001)	0.095 (0.023)	0.001 (0.0002)	6,30	34258,65	11030,79
	WL	18.182 (0.328)	–4.669 (1.103)	–0.016 (0.0019)	6,30	34257,28	11029,41
5/8H	WDnlin	12.185 (0.285)	0.107 (0.006)	0.001 (0.00007)	5,94	410374,20	107943,70
	WL	18.395 (0.097)	–5.157 (0.324)	–0.017 (0.0005)	5,94	410367,60	107937,10
7/8H	WDnlin	12.593 (0.462)	0.133 (0.010)	0.002 (0.0001)	5,98	127105,40	36370,50
	WL	21.113 (0.170)	7.253 (0.569)	0.018 (0.0009)	5,98	127107,20	36372,29
G	WDnlin	10.291 (0.402)	0.128 (0.010)	0.002 (0.00009)	5,78	707699,80	177744,60
	WL	17.039 (0.087)	–7.952 (0.700)	–0.018 (0.0004)	5,78	707686,80	177731,60
H	WDnlin	20.202 (0.148)	0.173 (0.002)	0.002 (0.00002)	6,81	1540066,00	391279,80
	WL	29.765 (0.060)	–12.007 (0.200)	–0.026 (0.003)	6,82	1540191,00	391404,80

a, b, c: parameters that define the scale and shape of the curve in the model; WD_nlin: Wood’s non-linear model; WL: Wilmink’s model; H: Holstein breed; G: Gyr breed; 1/4H, 3/8H, 1/2H, 5/8H, 3/4H, 7/8H: six genetic cross-breedings of Holstein X Gyr which is officially named as Girolando breed in Brazil.

Table 3. Estimated parameters (mean ± SE) from different linear models, root mean square error (RMSE), Bayesian information criterion (BIC) and Akaike information criterion (AIC).

Genetic composition	Model	Parameters
		a	b	c	d	f	RMSE	AIC	BIC
1/2H	WDlin	2.517 (0.024)	0.124 (0.007)	0.001 (0.00007)	–	–	0,49	−39632,36	−39630,36
	ASH	28.398 (1.448)	–16.068 (2.998)	–3.070 (1.710)	–3.557 (0.570)	0.013 (0.005)	6,48	387172,24	106412,94
1/4H	WDlin	2.379 (0.081)	0.122 (0.024)	–0.001 (0.0002)	–	–	0,52	−3711,54	−3709,54
	ASH	31.012 (4.405)	–28.171 (9.198)	11.668 (5.299)	–5.193 (1.724)	0.020 (0.0164)	6,49	33336,43	10962,69
3/4H	WDlin	2.380 (0.019)	0.122 (0.057)	–0.001 (0.00005)	–	–	0,45	−59359,81	−59357,81
	ASH	26.963 (1.130)	–12.630 (2.331)	0.560 (1.321)	–3.098 (10.445)	0.004 (0.004)	6,07	534783,46	139307,42
3/8H	WDlin	2.535 (0.082)	0.164 (0.024)	–0.002 (0.0002)	–	–	0,51	−3869,39	−3867,39
	ASH	24.581 (4.241)	−17.271 (8.805)	6.109 (5.046)	–2.768 (1.668)	0.005 (0.0162)	6,14	34261,42	11045,50
5/8H	WDlin	2.457 (0.024)	0.102 (0.007)	–0.001 (0.00007)	–	–	0,49	−42059,24	−42057,24
	ASH	24.315 (1.236)	–15.746 (2.559)	4.706 (1.456)	–2.686 (0.486)	0.004 (0.004)	6,00	383513,33	102754,02
7/8H	WDlin	2.494 (0.038)	0.130 (0.011)	–0.002 (0.0001)	–	–	0,46	−15304,19	−15302,19
	ASH	30.270 (2.175)	−20.982 (4.487)	6.197 (2.542)	–4.282 (0.857)	0.016 (0.008)	5,81	127105,75	36385,22
G	WDlin	2.338 (0.0435)	0.108 (0.011)	−0.002 (0.00009)	–	–	0,49	−69745,83	−69743,83
	ASH	20.586 (4.513)	–12.865 (7.459)	3.769 (3.058)	–1.221 (2.3597)	0.003 (0.0773)	5,60	707687,72	177750,10
H	WDlin	2.666 (0.009)	0.180 (0.002)	–0.002 (0.00002)	–	–	0,41	−176444,62	−176442,62
	ASH	36.537 (1.392)	–15.773 (0.011)	0.364 (0.00001)	−3.818 (0.989)	−0.003 (0.197)	6,66	1540051,44	391283,81

a, b, c, d, f: parameters that define the scale and shape of the curve in the model; WD_lin: Wood’s linear model; ASH: Ali and Schaeffer’s model; H: Holstein breed; G: Gyr breed; 1/4H, 3/8H, 1/2H, 5/8H, 3/4H, 7/8H: six genetic cross-breedings of Holstein X Gyr which is officially named as Girolando breed in Brazil.

Differences between non-linear and linear models were found based on the different fit quality tests (Tables 2 and 3). In the non-linear models (WD_nlin and WL) the values of RMSE, AIC and BIC were similar and ranged from 5.77 to 6.81 when applied to test day milk yield records of different genetic groups. The WD_nlin model presented the smallest values of RMSE, AIC and BIC for the 1/2H, 3/4H, 7/8H and H groups, whereas WL model presented the smallest values for 1/4H, 3/8H, 5/8H and G groups. However, in the linear models (WD_lin and ASH), RMSE differed between models and groups, ranging from 0.41 to 0.51, for WD_lin, and from 5.60 to 6.66, for ASH. Similarly, AIC and BIC values were different between genetic models and groups and were smaller when adjusted by the WD_lin model.

The shape of lactation curves was similar when the curves calculated by WD_nlin, WL, WD_lin and ASH models are compared with the mean observed test day yield records, except for the observed shapes of the curves of the 1/4H, 3/8H and 7/8H genetic groups (Figure 1). The lactation curve of the genetic group H presented the best fit when the shape of the observed curve is compared with the curves calculated by WD_nlin, WL and ASH models.

Figure 1. Trajectory of lactation curves estimated by 305-day milk yield from database (TMY), Wood’s linear model (WDlin), Ali and Schaeffer’s model (ASH), Wood’s non-linear model (WDnlin) and Wilmink’s (WL) model for 1/2H (A), 1/4H (B), 3/4H (C), 3/8H (D), 5/8H (E), 7/8 (F), G (G) and H (H) genetic groups.

The average values of PY, PT and TMY presented larger differences between models and genetic groups than those estimated for the different persistency measures (P, P_2:1 P_3:1 and P_3:2) (Table 4). The values of PY and TMY for the genetic groups ranged from 14.30 to 36.75 kg and from 3380.43 to 7539.26 kg, respectively, when estimated by the different models. The highest estimates of PY and TMY estimated by the WD_nlin, WL and WD_lin models were for the genetic group H, followed by the group 1/2H. Regardless of the model, ASH presented the higher value of PY for H group, followed by 3/4H. Similarly, the highest value of TMY estimated by ASH was also for H group, followed by 1/2H group. Independently of the model, the smaller PY and TMY were obtained for animals belonging to the genetic group G.

Table 4. Average estimates of peak yield (PY), peak time (PT), persistency (P, P_2:1; P_3:1, and P_3:2) and 305-day milk yield (TMY) for different genetic groups estimated by linear and non-linear models.

Model		Genetic composition
		1/2H	1/4H	3/4H	3/8H	5/8H	7/8H	G	H
WDnlin	PY	21.26	18.88	19.71	17.79	18.05	21.12	15.41	36.75
	PT	128.00	134.00	78.00	95.00	107.00	133.00	64.00	86.50
	P	7.79	7.83	7.18	7.56	7.64	7.82	7.01	7.28
	P2:1	101.75	102.11	102.74	99.15	99.62	101.96	95.11	103.46
	P3:1	90.41	90.67	90.36	87.88	87.92	90.45	114.64	91.40
	P3:2	88.84	88.79	87.95	88.63	88.26	88.79	120.53	88.36
	TMY	5438.78	4801.83	5382.73	4611.83	4644.91	5369.65	3762.42	7538.93
WL	PY	20.46	18.13	20.44	17.48	17.65	20.28	17.03	28.53
	PT	49.80	56.30	39.49	114.69	82.00	47.98	43.78	25.14
	P	6.85	6.86	8.22	4.66	5.15	7.25	7.95	12.00
	P2:1	100.48	101.09	101.07	98.60	90.83	100.69	94.49	101.73
	P3:1	90.56	90.95	90.44	88.04	88.03	90.53	82.29	91.41
	P3:2	90.12	89.96	89.44	89.29	89.07	89.91	87.08	89.84
	TMY	5443.5	4806.20	5387.18	4614.01	4647.01	5373.20	3765.42	7539.26
WDlin	PY	19.89	17.16	17.18	22.05	16.89	18.29	14.30	33.38
	PT	124.00	122.00	122.00	82.00	102.00	65.00	54.00	63.33
	P	7.760	7.75	7.75	7.23	7.61	7.02	6.88	6.91
	P2:1	100.26	100.34	101.76	96.47	97.78	101.51	93.14	100.57
	P3:1	87.98	88.23	88.13	84.31	85.89	88.14	79.29	88.14
	P3:2	87.74	87.93	86.60	87.39	87.83	86.82	85.13	87.63
	TMY	5016.91	4350.20	5019.70	4185.49	4265.06	5010.43	3380.43	5010.44
ASH	PY	28.39	31.01	26.96	24.51	24.31	30.27	20.58	36.53
	PT	5.23	2.41	22.55	2.82	3.34	3.38	3.41	43.33
	P	–	–	–	–	–	–	–	–
	P2:1	101.85	100.47	103.00	98.21	98.93	101.49	94.26	103.51
	P3:1	90.30	91.00	90.32	88.17	88.01	90.40	81.93	91.57
	P3:2	88.66	90.57	87.68	89.77	88.96	89.07	86.91	88.46
	TMY	5439.76	4811.29	5381.73	4617.19	4647.94	5372.84	3766.69	7536.41

P: persistency calculated by linear or non-linear parameters; P_2:1, P_3:1 and P_3:2: different milk yield persistency measures proposed by Johansson and Hansson (1940); WD_lin: Wood’s linear model; ASH: Ali and Schaeffer’s model; WD_nlin: Wood’s non-linear model; WL: Wilmink’s model; H: Holstein breed; G: Gyr breed; 1/4H, 3/8H, 1/2H, 5/8H, 3/4H, 7/8H: six genetic cross-breedings of Holstein X Gyr which is officially named as Girolando breed in Brazil.

The values of time to peak yield (PT) estimated by WD_nlin, WL, WD_lin and ASH for the cross-bred animals ranged from 78 to 134, from 39.5 to 114, from 65 to 124 and from 2 to 22 days in milk, respectively. The PT estimated by the various models were between 3 and 64 days in milk for Gyr and between 25 and 86 days for Holstein. .For the various combinations of genetic group and lactation curve model, the highest PT values were estimated for 1/4H (WD_nlin), 1/4H (WD_lin), 3/8H (WL) and H (ASH). For the same models, the smallest values for PT were obtained in the G, G, H and 5/8H genetic groups, respectively, with WD_nlin, WD_lin, WL and ASH models, respectively.

The persistency values estimated by WD_nlin, WL and WD_lin showed the lowest values of all models, ranging from 7.01 to 7.83, 5.15 to 12.00 and 6.91 to 7.76, respectively. For P_2:1, P_3:1 and P_3:2, the estimates ranged from 90.83 to 103.51%, 79.29 to 114.64% and 85.13 to 120.53%, respectively.

In regard to the non-linear models, the heterosis effect was evaluated only by WD_nlin, whose criteria of choice were based on the comparisons of the RMSE, AIC and BIC values between models (Table 2). In addition, the best quality of fit of curves was obtained for the different genetic groups by fitting the WD_nlin model (Figure 1). The heterosis effect estimated by WD_nlin was significantly different from zero (p < .001) for most traits, except for P_2:1 (Table 5). The heterosis effect of parameter a was positive and significantly different of zero (p < .001). For the parameters b and c, the estimates of heterosis effect were negative and significant (p < .001). The highest magnitudes of heterosis effect were observed for the PT (24.18%), followed by TMY (13.03%), parameter a (12.72%) and PY (12.30%). The lowest values of heterosis effect were observed for the different measurements of persistency, which were 0.20, 0.10, 0.55 and 0.42, for P, P_2:1, P_3:1 and P_3:2, respectively. The lactation curves of the different genetic groups of Girolando cattle were 2.79% more persistent than the average of persistency of parental breeds.

Table 5. Heterosis effect (mean ± SE) for Wood parameters (a, b and c), peak yield (PY), peak time (PT), persistency measures (P, P_2:1, P_3:1, and P_3:2) and 305-day milk yield (TMY) estimated by Wood’s non-linear model.

	Heterosis	Heterosis, %	SE	p Value
a	1.9400	12.72	0.020000	<.001
b	–0.0050	–3.32	0.000800	<.001
c	–0.0006	–30.00	0.000015	<.001
PY	3.2100	12.30	0.054000	<.0001
PT	18.200	24.18	1.009000	<.0001
P	0.2000	2.79	0.012000	<.0001
P2:1	0.1000	0.10	0.059000	.0684
P3:1	0.5500	0.53	0.052000	<.0001
P3:2	0.4200	0.40	0.015000	<.0001
TMY	736.36	13.03	14.728000	<.0001

Discussion

The application of mathematical models in genetic evaluations of lactation curve depends on the accuracy of the estimates, as well as the biological interpretation of the shape of the curve (Pollott and Gootwine 2000). After choosing the best model to fit the lactation curve, the emphasis of selection can be focused only on the special parts of the lactation curve (Hossein-Zadeh 2016). However, regarding the description of the lactation curves and its components, literature studies should be analysed with caution when interpreting the results they found (Pereira et al. 2016), as the mathematical models are tools that may affect the genetic evaluation and, consequently, the selection of best sires and dams in genetic selection schemes (Şeahin et al. 2015).

First, lactation curves were described by linear and non-linear models. Then the heterosis effect on the components of the lactation curve as average milk yield at peak, average peak time and average milk yield persistency were studied. The lactation curves estimated by linear and non-linear models in the different genetic groups presented a typical pattern. Such a pattern may be attributed to the values of the estimated parameters (regression coefficients) in the models. In the Wood model (WD_nlin and WD_lin), typical curves were estimated when positive values were found for a, b and c parameters, with b between 0 and 1 and always higher than 0 (Vadillo et al. 2012). The negative values of the parameters b and/or c in this model are considered a problematic or atypical result (Pollott and Gootwine 2000), which makes it impractical to calculate peak yield, peak time and total milk yield (Wood 1967). The estimated values of the parameters in the WL model showed that the shape of the lactation curve in the different genetic groups presented a typical pattern. According to Macciotta et al. (2005), when the values of the parameters b and c are negative, the lactation curve presents a typical pattern when fitted by WL model. In ASH model, a similar interpretation of the parameters is not feasible because the parameters do not present a biological interpretation (Macciotta et al. 2011). Graphical visualisation of the average lactation curves, fitted by ASH model in the different genetic groups, showed that the present estimated curves in this study presented a typical shape. Atypical shapes of the curve are characterised by the absence of a peak in the lactation curve (Macciotta et al. 2005).

As the estimation of atypical lactation curves is indicative of a decrease in the quality of fit, criteria were used to evaluate the quality of fit of linear (WD_lin and ASH) and non-linear (WD_nlin and WL) models. Comparisons between the values of RMSE, AIC and BIC showed that WL model presented better quality of fit for 1/4H, 3/8H, 5/8H and G groups than that of WD_nlin model (Table 2). However, WD_nlin indicated the best fit for 1/2H, 3/4H, 7/8H and H groups. Torshizi et al. (2011) and Bangar and Verma (2017), using similar criteria, reported the superiority of the non-linear Wood’s model compared to Wood’s linear model when fitting test day milk yield records of Holstein and Gyr breeds, respectively. In fact, Wood’s model in the non-linear form has been considered one of the best and most popular functions to describe the lactation curve in dairy cattle (Macciotta et al. 2011).

Among linear models, WD_lin was indicated the best fit of lactation curves of genetic groups because it presented lower values of RMSE, AIC and BIC than the values in ASH model (Table 3). However, despite the indication of the best fit by criteria, the shape of the curve estimated by WD_lin presented worse fit (lower precision) than the curves estimated by WD_nlin, WL and ASH when compared the observed curve (Figure 1).

In several studies, WD_nlin, ASH and WL models were successfully used in the adjustment of individual lactation curves of dairy cattle (Macciotta et al. 2005; Silvestre et al. 2006; Torshizi et al. 2011). The use of these models in the adjustment of the test day milk records is the most recommended option for primiparous Holstein cows, when there is not a lack of records in a population or subpopulation (Konchagul and Yazgan 2008) or abrupt reduction in the number of records in certain periods of lactation (Costa et al. 2005). That reduced number of animals and, consequently, the low number of test day records, may explain the fact that the models did not show the best quality of fit for the 1/4H, 3/8H and 7/8H groups.

WD_nlin and WL models underestimated the estimated milk yields in the first stage of the lactation curve, except for G group by WL model and 1/2H group by WD_nlin and WL models (Figure 1). That indicated that both models tend to underestimate test day milk yields before the peak of lactation. Torshizi et al. (2011), using the same models, found that the milk yields between 5 and 100 days in milk by Wood non-linear model and between 101 and 200 days by Wilmink model were underestimated, whereas the yields were overestimated at the final stage of lactation by both models. Similarly, the overestimation of milk yields at the final stage of lactation were found for the 1/4H, 3/4H and 3/8H calculated by WL and WD_nlin models which could, to some extent, be due to the small number of test day records at the end of lactation in some of the groups. Some studies also have reported the influence of the distribution of the available test day records along days in milk on the precision or quality of the adjustment of lactation curves when using Wood and Wilmink models (Silvestre et al. 2006; Konchagul and Yazgan 2008). Pollott and Gootwine (2000) and Torshizi et al. (2011) reported that those models overestimated the milk yields at final stage of lactation of dairy sheep and dairy cows, respectively.

The PY values estimated by ASH model in the genetic groups were higher than the values estimated by WD_nlin, WL and WD_lin (Table 4), whereas the values estimated for PY by the WD_nlin, WL and WD_lin models were much closer to each other on average compared with ASH. Although ASH model may estimate higher values than the estimated by Wood and Wilmink models, it also means that it may produce negative values or too high values at the beginning or end of lactation curve (Macciotta et al. 2005; Silvestre et al. 2006). In WD_nlin, WL and WD_lin models, the highest values of PY were for H group, followed by 1/2H. In addition, those genetic groups presented the highest estimates for TMY, regardless of the fitted model. According to Hossein-Zadeh (2017), there is a positive association between the milk yield at the peak and the cumulative milk yield at 305 days. Cows with high milk yield at the peak of lactation may present higher total milk yield at 305 days than cows with lower peak yield. It suggested that the selection of animals could be based on the milk yield at the peak of lactation (Hossein-Zadeh 2014).

The estimates of TMY were similar between WD_nlin, WL and ASH models and between these models and the TMY values estimated by the Brazilian official model because the differences between them were lower than 4 kg in average (Table 1). Among all models, WD_lin estimated the lowest values (–14%) compared with the other models, which means about –660 kg in average. According to Wood (1967), WD_lin model tends to overestimate the yields at the first stage of lactation, as well as to underestimate the yields at the final stage of lactation. The underestimation of milk yields at final days in milk of lactation curve, in this study, could be an explanation for the lowest values of TMY for the different genetic groups.

WD_nlin estimated the highest values of PT, followed by WD_lin and WL. The PT values estimated by ASH model for the genetic groups were between 67% and 87%, lower than those estimated by other models, which represents from 46 to 92 days of difference. It allows inferring that ASH was not a suitable model, because it underestimated the values for the PT trait in the different genetic groups. In general, the time to peak yield (PT) estimated by WD_nlin, WL and WD_lin for the cross-bred animals were from 30% to 159% and 49% to 94% higher than those estimates for Holstein and Gyr, respectively. It indicates that cross-bred animals take between 21 and 49 days longer to reach the peak than pure breeds. However, ASH model estimated values for cross-bred animals, which were 85% lower than the 43 days for Holstein but was close to the 3 days estimated for the Gyr breed.

There are not many studies comparing the time to peak of Holstein, Gyr and cross-bred Holstein × Gyr animals in the literature. Some studies reported the time to peak between 60 and 90 days for Holstein (Torshizi et al. 2011), around 60 days for Gyr (Herrera et al. 2008) and from 26 to 56.5 for the Holstein × Gyr cross-bred animals (Glória et al. 2010; Jacopini et al. 2016). In fact, the literature is not clear about whether the time to peak of cross-bred animals is shorter than that of pure breeds or not. The large variation in the values of time to peak estimated for the genetic groups compared with those reported in the literature could be attributed to many factors as the model or function of estimation or the lower number of test day records (or the lack of) in the first stage of lactation, which, in consequence, affects the shape of lactation curve (Oliveira et al. 2007; Glória et al. 2010; Torshizi et al. 2011). Besides, such metabolic and behaviour effects, as well as age and weight of animals may also affect the results (Borges et al. 2015).

The fact is that the shape of the lactation curve, as well as the day in milk, may affect the results of different measures of persistency. Therefore, the choice of the best persistency measure will depend on the model used for fit the lactation curves.

The differences of persistency values between the models within the same genetic groups were in the average of 3 kg for the majority of the estimated values, considering P, P_2:1, P_3:1 and P_3:2 measures. It means, in general, that persistency was the component of the lactation curve that showed the lower variation between models. The persistency of milk yield is one of the most important traits because cows with high persistent milk yield can maintain a high level of milk yield after the peak of lactation (Cobuci et al. 2003). That increase in the amount of milk yield after the peak yield could decrease the costs of production.

The values of P estimated by WD_nlin and WD_lin of the cross-bred genetic groups were around 1% larger than that in Holstein or Gyr group, whereas the values by WL model were in average between 1 and 6% lower than that in the pure breeds. For P_3:1, P_2:1, and P_3:2, the persistency values of cross-bred genetic groups estimated by almost all models were between 0.9 and 3% lower than that in the Holstein group. The cross-bred animals presented values between 2 and 8% higher than that in Gyr group, except for P_3:1 by WD_nlin model. For the P measure of persistency, it could indicate that the cross-bred animals showed more persistency of milk yields than the pure breeds. However, the other three persistency measures indicated that cross-bred genetic groups had less persistent milk yields than Holstein, but more persistent milk yields than Gyr breed. Thus, this could suggest that there is a positive association between TMY and persistency, because Holstein breed presented higher values of 305 day milk yield and more persistency of milk yield than the cross-bred groups for almost all models. According to Gengler (1996), persistency is influenced by the level of milk yield. Similarly, to this study, Hossein-Zadeh (2016) reported a positive association between the different measures of persistency (P_2:1, P_3:1 and P_3:2) and milk yield up to 305 days using different non-linear models.

Many authors have studied the heterosis effect on milk yield (Cunningham and Syrstad, 1987; Rege, 1998; Penasa et al. 2010). In Brazil, the heterosis effect has had an important contribution on the performance of the different genetic groups of Girolando cattle (Facó et al. 2002; Facó et al. 2008). WD_nlin model was chosen to estimate the heterosis effect of all traits in the genetic groups of cross-bred cows of Girolando cattle based on the quality of fit criteria and the estimated lactation curves. Such a choice was made because WD_nlin model presented a better fit quality for a larger number of animals within the genetic groups (1/2H, 3/4H, 7/8H and H), when compared with the WL model that presented better fit for the groups (1/4H, 3/8H, 5/8H and G). In addition, the WD_nlin model presented better accuracy in the adjustment of the lactation curve of the different genetic groups compared with the linear model (WD_lin) with better adjustment according to RMSE, AIC and BIC.

The heterosis effect was significant for TMY and the most components of lactation curves. In general, those results were expected, because the best results from the cross-bred animals are obtained when the parental pure breeds are more genetically distant. If parental breeds have different alleles or frequency of different alleles, offspring will show more heterogeneity and heterosis compared with crosses between parental breeds with similar allele frequencies (Sorensen et al. 2008).

There was heterosis effect of parameter a estimated by WD_nlin, which is associated to the average yield along the initial stage of lactation (Oliveira et al. 2007). The cows of the different genetic groups of the Girolando breed presented 12.72% more milk at the beginning of lactation compared with the average of their parental breeds. In contrast, the b and c parameters of that model had a negative heterosis effect (–0.005 and –0.0006, respectively). According to Wood (1967), the parameter b represents the average rate of increase in the milk yield after the peak yield and c the average rate of decline of milk yield after the peak. Although the parameter c had a negative value for the heterosis effect, it could be favourable for the farmers, milk industry and health cows. It means that a smoother decline in the milk yield after the peak yield is associated with the flatter shape of the lactation curve, that is, to the milk yield persistency. Consequently, more persistent lactating cows present a more balanced distribution of milk yield along lactation (Gengler 1996).

TMY and PY had expressive values of heterosis effect of 13.03% and 12.30%, respectively. The heterosis effect of TMY and PY followed the same positive association with the average yields estimated in Table 4. It means that the higher the TMY, the higher PY. Although the increase in the milk yield of cross-bred cows compared with the average of their parental pure breeds was an aspect of economic interest for farmers, higher peak yields are associated to metabolic disorders (Remppis et al. 2011). Cows with smaller PY but higher persistency are more desirable for dairy cattle (Atashi et al. 2013). In Brazil, it is a fact that the selection of top bulls from pure breeds has been focussed on the 305-day milk yield. Consequently, it might have influenced the heterosis effect of PY. It also could suggest that previous selection of parental breeds with higher TMY and smaller PY would bring more benefits to Girolando cattle.

Among all components of the lactation curve, the lower effect of heterosis was found for the different measurements of persistency (P, P_3:1 and P_3:2). That result may have occurred because cross-bred animals with a lower percentage of Holstein breed genes present shorter lactation duration when compared to the Holstein. Thus, it may have affected the lower heterosis effect of persistency.

The highest estimate of heterosis effect among all traits was for PT (24.18%). It means that cross-bred cows from the genetic groups of Girolando cattle may present a later lactation peak compared with the Holsteins. Several studies have reported that the increase of time to peak yield was followed by the increase in the 305-day milk yield of Holstein cows (Muir et al. 2004; Němečková et al. 2015). The main benefit of a longer time to reach the peak yield compared with pure breeds could be the decrease in the probability of metabolic problems caused by a high peak yield and/or an earlier time to peak. It is known that feed intake increases in the initial phase of the lactation, but peak milk yield usually occurs before the maximum feed intake is reached. Therefore, in this early stage cows must use the energy of their body storage to produce milk, and as a consequence enters in a negative energy balance (Remppis et al. 2011). Therefore, a later time to peak yield could decrease the magnitude of this negative energy balance, which could bring not only economic benefits for farmers but also benefits for the health and welfare of cows.

Conclusions

On the basis of criteria of goodness of fit, the results of this study showed that Wood linear model provided the best fit of lactation curve for different genetic groups. However, Wood’s non-linear model described most appropriately the lactation curve compared with the linear model.

The heterosis of 305-day milk yield is more associated to the heterosis of components of the initial stage of the lactation curve (initial milk yield, peak yield and time to peak) than that in the intermediate and final stage of lactation (persistency). It could be suggested for further studies in the genetic programmes the introduction of persistency as a criterion of selection of progenitors from pure breeds to increase persistency and 305-day milk yield in Girolando.

Ethical approval

All authors declare that this study follows the principles of the Declaration of Helsinki.

Disclosure statement

No potential conflict of interest was reported by the authors.