Heterosis in the lactation curves of Girolando cows with emphasis on variations of the individual curves

ABSTRACT

Linear and nonlinear models were used to evaluate the effect of heterosis on the components of the lactation curves of different crossbred groups of cows from Girolando cattle. Data consisted of 233,587 test-day milk yield records from 33,995 primiparous cows from 1998 and 2014. The Wood’s linear model (WD_lin), Wood’s nonlinear model (WD_nlin), Wilmink’s model (WL) and Ali and Schaeffer’s model (ASH) were used to estimate individual peak milk yield (PY), day of peak milk yield (PT), 305-day milk yield (TMY), different persistency measures and parameters of lactation curves. The quality of fit of models was different for the genetic groups. The WL model was used to the estimation of heterosis because of the best fit of lactation curve. The heterosis effect was significant (P < .001) for TMY, PY and for the ‘a’ parameter of WL model (initial production). For TMY and PY, crossbred cows presented 14.64% and 20.60% larger yields than the average of the parental breeds, respectively. The heterosis effect from crossbreeding presents more benefits for the components of initial stage of milk yield and milk yield at the peak and with a lesser extent with the two final stages (persistency).

KEYWORDS:

• Crossbreeding
• dairy cattle
• daily milk yield
• lactation model

1. Introduction

Milk production in Brazil has made considerable increased considerably in the last decades, from 11.1 million litres in 1980 to more than 35 billion litres of milk in 2014 (Prata et al. 2015). It was estimated that about 80% of the milk produced in Brazil comes from cows that have Holstein and Gyr genes in their genetic composition (Silva et al. 2015). The crossbreeding between those two breeds (Holstein and Gyr) has been an important tool in order to increase milk yield and reproductive efficiency and improve the adaptation of animals in the tropical conditions, which uses the benefits from the heterosis expression and complementarity between breeds (Canaza-Cayo et al. 2014; Prata et al. 2015). Such benefits have contributed to the beginning of Girolando cattle formation in 1940 (Silva et al. 2015).

Animals with different genetic compositions from different crossbreedings present large genetic variation. In such a heterogenous group, individual lactation curves are useful for predicting the performance of milk yield (Pereira et al. 2016). Many mathematical models have been developed in order to describe the shape of the lactation curve and its graphical trajectory along days in milk (Ali and Schaeffer 1987; Wood 1967, Wilmink 1987). The modelling of lactation curve may predict the level of production with high accuracy as well as permit understanding the pattern of milk yield under different environmental conditions (Şeahin et al. 2015). Thus, functions that describe milk production in time can be very applicable in genetic breeding programmes, since through the mathematical model the milk production of cows can be predicted and the results can be applied in animal breeding programmes (Hossein-Zedeh 2017). Most studies on the lactation curve evaluate the average pattern between homogeneous groups of animals, even when the most important are the individual curves (Şeahin et al. 2015). On the other hand, considering that the shape of lactation curve may differ between breeds and animals within the breed, the study of different genetic groups is fundamental to verify the heterosis on the traits associated to the lactation curve traits. The heterosis effect from European-Zebu crossbred populations may range between 17% (Rege 1998) and 28% (Cunningham and Syrstad 1987).

Studies have verified the expression of heterosis on the milk yield traits and on the lactation length (Facó et al. 2005; Bryant et al. 2007; Lembeye et al. 2015). However, there are few studies about the heterosis effect on the shape of the lactation curve of Girolando cows, mainly for the different genetic groups of this breed.

The aim of this study was to verify the heterosis effect on the lactation curve characteristics of different compositions of genetic groups from Girolando cattle, using linear and nonlinear models.

2. Materials and methods

2.1. Data

Data from test-day records, which were collected between 1998 and 2014, consisted of 233,587 test day milk yield records (MY) of the first lactation from 33,995 cows of Minas Gerais State (Brazil). Those cows were from Holstein breed (H), Gyr breed (G), and six genetic crossbreedings of Holstein x 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 4–10 test day records per cow, 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 as normal probability plots and boxplots, as well as by median, mean, mode, skewness and kurtosis values. Then 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 to 11026 kg, respectively. The descriptive analysis of the edited data is presented in Table 1.

Table 1. Number of cows, herds and test day milk yield records as well as 305-day milk yield records (TMY) of genetic groups with different proportion of Holstein genes obtained from the Brazilian databases.

Genetic group	Proportion of Holstein-gene, %	Number of herds	Number of cows	Number of test-day records	TMY^a, kg
H	100	278	13,990	79,578	7537.98
1/4H	25.00	116	374	2733	4794.49
1/2H	50.00	299	3568	26,656	5439.11
3/4H	75.00	288	4526	36,691	5380.94
3/8H	37.50	87	368	2861	4610.39
5/8H	62.50	273	3696	28,574	4661.99
7/8H	87.50	144	1177	9636	5375.60
G	0	92	6296	46,858	3767.69

^a305-day milk yield from 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.

2.2. Lactation curve models

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

1. Wood’s nonlinear model (WD_nlin) (Wood 1967):$Y_t=a{t}^b{e}^{-ct}$

(2) Wilmink’s model (WL) (Wilmink 1987):$Y_t=a+b{e}^{-0.05t}+ct$(3) Wood’s linear model (WD_lin) (Wood 1967):$\log Y_t=\log a+b\log t-ct$

In the models, Y_t is the average daily yield in the t ^th 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. Ali and Schaeffer’s model (ASH) (Ali and Schaeffer 1987):$Y_t=a+b\left({\displaystyle \frac{t}{305}}\right)+c{\left({\displaystyle \frac{t}{305}}\right)}^2+d\ln \left({\displaystyle \frac{305}{t}}\right)+el{n}^2\left({\displaystyle \frac{305}{t}}\right)$

In this model, Y_t is the average daily yield in t ^th test day of lactation, and 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, peak milk yield (PY) and the day of maximum milk yield (PT) were estimated for each model using the mathematical functions as referred in the original papers. The persistency measures (P_2:1, P_3:1 e P_3:2 and P_weller) used were ratios between average milk yields obtained in different parts of the lactation and all four measures were expressed as percentage. The P_2:1, P_3:1 and P_3:2 were calculated as proposed by Johansoon and Hansson (1940):${\mathrm{P}}_{2:1}=\left({\displaystyle \frac{Milkyieldbetween101and200daysafterparturition}{Milkyieldinthefirst100daysofthelactation}}\right)\times 100\mathrm{\%}$ ${\mathrm{P}}_{3:1}=\left({\displaystyle \frac{Milkyieldbetween201and300daysafterparturition}{Milkyieldinthefirst100daysofthelactation}}\right)\times 100\mathrm{\%}$ ${\mathrm{P}}_{3:2}=\left({\displaystyle \frac{Milkyieldbetween201and300daysafterparturition}{Milkyieldbetweenfrom101to200daysafterparturition}}\right)\times 100\mathrm{\%}$

The Pw_eller were defined by Weller et al. (2006) as estimated milk production at 180 d after peak divided by estimated peak production in percent was described as:$P_{Weller}={\displaystyle \frac{Milkyield(270)}{Milkyield(90)}\times 100\mathrm{\%}}$where $\mathrm{Milk}\mathrm{yield}$ (270) and $\mathrm{Milk}\mathrm{yield}$ (90) are milk production at 270 and 90 days in milk, respectively.

Predicted 305-d MY (TMY) were obtained for each model using the following equation described by Vargas et al. (2000):$TMY=\sum _{t=5}^{305}y(t)$where TMY denotes predicted 305-d MY and y(t) represents MY at day t estimated by corresponding lactation models.

2.3 Breed and heterosis effects

The three breeds were considered with enough records to estimate the breed effect for all traits. The proportion of genes was calculated for each cow as:${\alpha }_i^p={\displaystyle \frac{({\alpha }_i^s+{\alpha }_i^d)}{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 were calculated between any pair of the dairy breeds using the following identify (Dickerson 1973):${\delta }_{ij}^p={\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.

The coefficients of specific heterosis effect were used for the six genetic groups in Girolando cattle, because the distribution of cows across classes of coefficients of expected heterosis was suitable for this purpose. The coefficient of general heterosis for each cow was obtained by summing coefficients of specific heterosis previously calculated.

2.4 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). The nonlinear models were fitted to the milk yield records as the iteration method of Gauss–Newton.

The models were tested for goodness of fit using the adjusted coefficient of determination (R ² _adj), root means square error (RMSE) and Akaike’s information criterion (AIC).

R ² _adj was calculated using the following formula:$R_{a}d{j}^2=1-[((n-1))/((n-p))](1-R^2),$where R ² is the coefficient of determination$R^2=1-{\displaystyle \frac{RSS}{TSS}}$

TSS is total sum square, RSS is residual sum of square, n is the number of observations (data points) and p is the number of parameters in the equation.

RMSE is kind of generalized standard deviation and was calculated as follows:$RMSE=\sqrt{{\displaystyle \frac{RSS}{n-p-1}}}$where RSS is the square root of residual sum of squares, n is the number of observations (test day records) and p is the number of parameters in the equation. The best model is the one with the lowest RMSE.

AIC was calculated using the following equation (Burnham and Anderson 2002):$AIC=n\times \ln (\mathrm{RSS})+2\mathrm{p}$

A smaller numerical value of AIC indicates a better fit when comparing models.

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

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

$\mu$ is na constant.

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

$C_k$ 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$ the regression coefficient associated with the linear effect of proportion of Holstein (f).

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

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

3. Results

The values of the parameters estimated by the non-linear (WD_nlin and WL) and linear models (WD_lin and ASH) and the criteria of the quality of fit of these models for test day milk yield records of cows from different genetic groups are associated to the pattern of the typical shape of the lactation curves calculated in this study (Table 2 and 3, Figure 1). The values of the parameters estimated by WD_nlin, considering the different genetic groups, ranged from 13.377–16.923 for parameter ‘a’, from 0.210–0.376 for parameter ‘b’ and from 0.0033–0.0052 for ‘c’. The WL model estimated higher values for parameter ‘a’ (16.831–22.312) than that in the WD_nlin, and negative values for parameters ‘b’ (−5.544–12.740) and ‘c’ (−0.0207 to −0.0003) in most genetic groups. In the linear models, the values of the lactation curve parameters estimated by WD_lin model for the different genetic groups ranged from 2.254–2.797 for the parameter ‘a’, from 0.040–0.210 for the ‘b’ and −0.001 to −0.003 for parameter ‘c’. The ASH model estimated positive values for the parameter ‘a’ (3.522–14.729) and negative for the other parameters ‘b’ (−0.006 to −0.025), ‘c’ (−0.0001 to −0.00009), ‘d’ (2.973 to −10.842) and ‘f’ (−0.049 to −1.236).

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 nonlinear model (WDnlin) and Wilmink model (WL) for 1/2H (A), 1/4H (B), 3/4H (C), 3/8H (D), 5/8H (E), 7/8H (F), G (G) and H (H) genetic groups.

Table 2. Estimated parameters (mean ± SE) of the lactation curve obtained from different nonlinear models, adjusted coefficient of determination (R ² _Adj), root mean square error (RMSE) and Akaike information criterion (AIC).

		Parameters			Statistics
Model	Breed (Composition)	a	b	c	R ^2 Adj	RMSE	AIC
WDnlin	1/2H	14.363 ± 3.838	0.258 ± 0.618	0.0036 ± 0.072	0.606	0.242	101,496.20
	1/4H	13.377 ± 3.904	0.261 ± 0.732	0.0036 ± 0.372	0.615	0.667	10,066.64
	3/4H	14.038 ± 3.841	0.254 ± 0.603	0.0033 ± 0.069	0.578	0.231	142,867
	3/8H	13.928 ± 3.837	0.238 ± 0.614	0.0041 ± 0.005	0.643	0.591	10,097.12
	5/8H	14.103 ± 3.823	0.210 ± 0.594	0.0035 ± 0.071	0.612	0.189	102,308.50
	7/8H	14.161 ± 3.834	0.249 ± 0.604	0.0036 ± 0.070	0.589	0.457	37,657.30
	G	12.125 ± 4.115	0.376 ± 0.736	0.0052 ± 0.082	0.699	0.094	148,354.70
	H	16.923 ± 3.728	0.224 ± 0.516	0.0028 ± 0.064	0.626	0.417	58,534.35
WL	1/2H	21.739 ± 2.893	−7.577 ± 4.441	−0.0223 ± 0.212	0.615	0.237	100,883.60
	1/4H	19.407 ± 2.292	−6.990 ± 4.277	−0.0234 ± 0.230	0.607	0.681	10,126.21
	3/4H	21.426 ± 2.759	−7.986 ± 4.351	−0.0207 ± 0.216	0.572	0.234	143,375.20
	3/8H	19.257 ± 2.773	−6.107 ± 4.213	−0.0293 ± 0.219	0.607	0.651	10,376.66
	5/8H	18.833 ± 2.694	−5.544 ± 4.113	−0.0235 ± 0.189	0.626	0.182	101,261.10
	7/8H	21.643 ± 2.733	−7.927 ± 4.236	−0.024 ± 0.192	0.597	0.448	37,478.76
	G	16.831 ± 2.404	−6.561 ± 4.888	−0.0241 ± 0.167	0.723	0.086	144,316.30
	H	22.312 ± 2.803	−12.740 ± 4.053	−0.0003 ± 0.223	0.652	0.389	57,475.23

Notes: a, b, c: parameters that define the scale and shape of the curve in the model; WD_nlin: Wood’s nonlinear model and 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 quality of fit criteria indicated that there were differences between non-linear and linear models as well as between genetic groups (Table 2 and 3). WD_nlin model presented the highest mean values of R² _adj for the 1/4H (0.615), 3/4H (0.578), 3/8H (0.643) groups, while WL presented the highest mean values for the 1/2H (0.615), 5/8H (0.626), 7/8H (0.597), G (0.723) and H (0.652) groups (Table 2). Similarly, the lowest values of RMSE and AIC fitted by WD_nlin were for the 1/4H, 3/4H, 3/8H groups, where as the WL presented the lowest mean values for the 1/2H, 5/8H, 7/8H, G and H groups.

The mean values of R ² _adj differed between linear models and different genetic groups. In the linear models, the ASH model presented the highest mean values of R ² _adj (0.500–0.629) considering all the genetic groups while for the WD_lin model the values ranged from 0.413–0.539 (Table 3). The mean RMSE and AIC values differed substantially between the linear models. In addition, the lowest values were observed for WD_lin regardless of the genetic group. In regard to the mean values of R ² _adj, the ASH model presented the best fit of the lactation curves of all genetic groups. The WD_lin was considered the most suitable model to describe the lactation curves according to RMSE and AIC criteria.

Table 3. Estimated parameters (mean ± SE) of the lactation curve obtained from different linear models, adjusted coefficient of determination (R² _Adj), root mean square error (RMSE) and Akaike information criterion (AIC).

		Parameters					Statistics
Model	Breed (Composition)	a	b	c	d	f	R ^2 Adj	RMSE	AIC
WDlin	1/2H	2.787 ± 1.019	0.076 0.601	−0.002 ± 0.067	––	–	0.430	0.166	−26.99
	1/4H	2.669 ± 1.052	0.078 0.626	−0.002 ± 0.070	–	–	0.465	0.175	−26.05
	3/4H	2.754 ± 1.005	0.082 0.583	−0.001 ± 0.063	–	–	0.413	0.170	−28.76
	3/8H	2.754 ± 1.039	0.040 0.622	−0.002 ± 0.070	–	–	0.505	0.169	−27.58
	5/8H	2.254 ± 1.144	0.177 0.629	−0.003 ± 0.070	–	–	0.472	0.164	−28.03
	7/8H	2.797 ± 1.001	0.070 0.589	−0.002 ± 0.064	–	–	0.423	0.172	−28.95
	G	2.625 ± 1.288	0.068 0.753	−0.003 ± 0.077	–	–	0.539	0.163	−27.91
	H	2.544 ± 1.036	0.210 0.558	−0.002 ± 0.054	–	–	0.488	0.132	−34.25
ASH	1/2H	14.729 ± 13.415	−0.019 ± 0.977	−0.0001 ± 0.034	−2.973 ± 12.046	−0.049 ± 4.855	0.519	2.499	12.92
	1/4H	11.076 ± 12.370	−0.025 ± 0.954	−0.0001 ± 0.035	−4.702 ± 11.201	−0.270 ± 4.568	0.584	2.362	10.66
	3/4H	13.016 ± 12.686	−0.014 ± 0.940	−0.0001 ± 0.034	−4.065 ± 11.385	−0.188 ± 4.605	0.500	2.531	14.54
	3/8H	8.336 ± 11.373	−0.006 ± 0.892	−0.0001 ± 0.036	−7.892 ± 10.256	−0.869 ± 4.181	0.571	2.153	11.24
	5/8H	12.246 ± 11.429	−0.012 ± 0.861	−0.00009 ± 0.031	−4.104 ± 10.296	−0.397 ± 4.179	0.566	2.083	10.56
	7/8H	13.473 ± 12.568	−0.035 ± 0.922	−0.00008 ± 0.032	−3.724 ± 11.290	−0.085 ± 4.563	0.519	2.523	14.89
	G	3.522 ± 14.959	−0.005 ± 1.005	−0.0001 ± 0.033	−10.842 ± 13.237	−1.236 ± 5.258	0.629	1.648	6.20
	H	13.316 ± 12.395	−0.0176 ± 0.870	−0.0001 ± 0.029	−10.020 ± 11.034	−0.890 ± 4.427	0.572	2.624	16.56

Notes: a, b, c, d, f: parameters that define scale and shape of curve in the model; WDnlin: Wood’s linear model and 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.

Except for 1/4H, 3/8H, 7/8H and H groups, the mean lactation curves estimated by WD_nlin, WL, WD_lin and ASH may be used to represent the real shape of the average lactation curve (Figure 1). For the 1/4H group, the mean individual curves estimated by WD_nlin, WL and ASH did not represent the observed curve, since the peak of lactation of observed curve was not fitted adequately. For the 3/8H group, only the curve obtained by the ASH model fitted the initial stage of lactation compared to observed curve. For the 7/8 group, the trajectory of the lactation curve estimated by WD_nlin did not follow the ascending stage. Similarly, for H group, when the curve estimated by WD_nlin was compared to observed curve, the trajectory did not follow either the ascending or descending stage of the lactation curve. In general, the WD_lin model estimated lower milk yield values than those in the observed lactation curve on all stages. The visual comparison of lactation curves estimated by WD_nlin, WL, ASH with the observed curve showed that the best fit was for 1/2H e G groups (Figure 1).

The estimate mean of PY and PT differed between models and genetic groups (Table 4 and 5). Different results were found for the estimates of TMY and persistency (P_2:1, P_3:1, P_3:2 and P_weller), which presented smaller differences between models and genetic groups.

Table 4. Estimated (mean ± SE) of peak yield (PY), peak time (PT), persistency (P_2:1; P_3:1, P_3:2 and P_weller) and 305-day milk yield (TMY) for different genetic groups estimated by nonlinear models.

		Breed (Composition)
Model		H	1/4H	1/2H	3/4H	3/8H	5/8H	7/8H	G
WDnlin	PY	27.87 ± 2.56^a	19.13 ± 2.70^d	21.65 ± 2.65^b	20.95 ± 2.51^c	18.56 ± 2.55^de	18.24 ± 2.53^e	21.23 ± 2.50^bc	16.64 ± 2.56^f
	PT	93.26 ± 17.37^a	63.64 ± 8.02^bcd	70.26 ± 8.30^bc	73.80 ± 8.27^b	62.40 ± 8.27^bcd	60.28 ± 8.55^d	70.76 ± 8.28^bc	68.69 ± 51.95^c
	P2:1	103.69 ± 4.17^a	96.85 ± 5.01^d	101.30 ± 4.72^b	102.92 ± 4.67^a	94.98 ± 4.75^d	96.85 ± 4.67^d	100.86 ± 4.52^bc	100.08 ± 5.08^c
	P3:1	89.66 ± 5.01^a	79.24 ± 6.18c^d	87.20 ± 7.88^ab	87.67 ± 5.60^ab	77.05 ± 5.67^d	82.11 ± 8.53^c	85.40 ± 5.49^b	78.81 ± 5.75^d
	P :2	85.21 ± 3.91^a	78.22 ± 4.97^d	83.13 ± 5.97^b	83.50 ± 4.33^b	78.48 ± 4.66^d	81.68 ± 6.02^c	82.65 ± 4.43^bc	76.73 ± 4.68^d
	Pweller	78.98 ± 4.75^a	70.53 ± 6.04^ab	79.09 ± 12.64^a	77.09 ± 5.29^a	69.64 ± 5.67^ab	76.32 ± 13.27^a	75.81 ± 5.31^a	68.03 ± 5.49^b
	TMY	7321.86 ± 1556.14^a	4782.31 ± 1939.43^c	5441.27 ± 1975.17^b	5384.15 ± 1803.25^b	4623.76 ± 1806.46^c	4644.33 ± 1794.35^c	5372.36 ± 1769.58^b	3737.18 ± 1584.82^d
WL	PY	21.55 ± 2.36^a	18.34 ± 2.61^c	20.78 ± 2.62^a	20.39 ± 2.48^b	18.19 ± 2.51^c	17.78 ± 2.48^c	20.59 ± 2.48^ab	15.61 ± 2.40^d
	PT	106.35 ± 10.62^d	81.19 ± 5.81^ab	82.27 ± 5.86^a	84.92 ± 6.39^b	85.01 ± 5.77^abc	88.77 ± 5.84^c	85.91 ± 5.85^bc	81.08 ± 6.09^a
	P2:1	98.58 ± 4.54^ab	94.94 ± 5.74^cd	97.27 ± 4.78^b	98.27 ± 4.31^a	91.10 ± 4.47^e	93.17 ± 4.38^de	96.97 ± 4.37^abc	87.63 ± 4.31^f
	P3:1	99.34 ± 6.53^a	82.20 ± 8.45^cd	86.41 ± 6.94^bc	87.59 ± 6.18^b	75.22 ± 6.54^ef	79.14 ± 6.24^de	85.00 ± 6.17^bc	73.20 ± 5.77^f
	P3:2	97.31 ± 4.76^ab	73.57 ± 10.55^ab	81.37 ± 10.94^b	90.81 ± 18.55^a	79.00 ± 6.70^ab	78.51 ± 8.30^b	85.73 ± 6.70^ab	78.58 ± 7.14^b
	Pweller	98.96 ± 5.95^a	73.50 ± 7.84^de	79.48 ± 7.24^bc	80.88 ± 6.15^b	70.40 ± 7.17^e	74.15 ± 6.39^e	77.50 ± 7.93^cd	72.68 ± 6.18^e
	TMY	7307.96 ± 1856.67^a	4794.28 ± 1923.09b^c	5449.82 ± 1969.00^b	5391.54 ± 1794.81^b	4627.30 ± 1800.01^c	4652.60 ± 1787.73^c	5383.58 ± 1759.36^b	3744.07 ± 1577.05^d

Notes: PY: peak yield; PT: peak time; P_2:1, P_3:1 and P_3:2: milk yield persistency measures proposed by Johansson and Hansson (1940); P_weller: milk yield persistency measure proposed by Weller et al. (2006); TMY: 305- day milk yield; WD_nlin: Wood’s nonlinear model; WL: Wilmink; 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. ^a-dEstimates of least square means with different letters in rows are significantly different (p < .05).

The means of PY and TMY reached for the genetic groups when estimated by the WD_nlin, WL, WD_lin and ASH models ranged from 16.64–27.87 kg and 3737–7321 kg, from 15.61–21.55 kg and 3744–7307 kg, from 13.01–23.10 kg and 3374–7021 kg, from 20.70–27.00 kg and 3745–7303 kg, respectively. Regardless of the fitted model, the highest mean estimates of PY and TMY were obtained for H group, followed by group 1/2H. The lowest average values of PY and TMY were estimated by the WD_nlin, WL and WD_lin models for the G group.

Unlike those models, the ASH model estimated the lowest mean of PY for 5/8H, but, similarly to the other models, the lowest mean of TMY was found for G group.

The mean values of PT estimated by the WD_nlin, WL, WD_lin and ASH for the crossbred groups ranged from 60.28–73.28 days, 81.08–88.77 days, 58.03–80.26 days and 90.22–100.01 days, respectively. The PT values estimated by WD_nlin and WL were 93.26 and 106.35 days for H and 68.69 and 81.08 for G, respectively (Table 4). The PT values estimated by linear models were 72.48 and 106.96 days for H and 68.64 and 107.55 days in milk for G, respectively (Table 5).

Table 5. Estimated (mean ± SE) of peak yield (PY), peak time (PT), persistency (P_2:1; P_3:1, P_3:2 and P_weller) and 305-day milk yield (TMY) for different genetic groups estimated by linear models.

		Breed (Composition)
Model		H	1/4H	1/2H	3/4H	3/8H	5/8H	7/8H	G
WDlin	PY	23.10 ± 3.37^a	15.08 ± 3.81^cef	18.86 ± 3.38^b	18.03 ± 3.31^bc	14.07 ± 3.72^def	15.14 ± 3.22^e	17.58 ± 3.14^bcd	13.01 ± 2.51^f
	PT	72.48 ± 7.55^de	80.26 ± 23.72^e	60.87 ± 6.97^ab	65.19 ± 6.92^bc	58.03 ± 7.70^ab	60.28 ± 4.76^a	66.91 ± 7.50^bd	68.64 ± 6.99^cd
	P2:1	133.40 ± 4.37^e	138.08 ± 5.18^bc	137.64 ± 4.85^cd	137.65 ± 4.72^cd	140.53 ± 4.70^b	138.79 ± 4.84^bc	136.52 ± 4.73^d	144.35 ± 4.79^a
	P3:1	154.90 ± 5.72^d	164.98 ± 6.79^bc	163.28 ± 6.41^c	162.96 ± 6.21^c	170.25 ± 6.02^b	166.09 ± 6.36^b	161.12 ± 6.23^c	175.59 ± 6.19^a
	P :2	114.69 ± 3.28^d	116.85 ± 3.81^bc	116.56 ± 3.66^c	116.51 ± 3.58^c	120.31 ± 3.99^a	117.88 ± 6.68^b	116.13 ± 3.58^c	120.27 ± 3.25^a
	Pweller	132.57 ± 4.62^d	138.63 ± 11.77^bc	137.61 ± 5.17^c	137.34 ± 5.01^c	145.99 ± 6.95^a	140.24 ± 5.28^b	136.37 ± 5.05^c	145.49 ± 4.89^a
	TMY	7021.06 ± 396.03^a	4384.89 ± 466.80^c	5048.44 ± 447.43^b	5026.83 ± 433.39^b	4256.89 ± 455.43^c	4275.31 ± 450.88^c	5031.83 ± 432.14^b	3374.78 ± 434.61^d
ASH	PY	27.00 ± 3.84^a	20.85 ± 3.76^e	22.53 ± 3.81^b	22.20 ± 3.83^c	20.88 ± 3.76^e	20.70 ± 3.77^e	22.09 ± 3.78^c	21.80 ± 3.97^d
	PT	106.96 ± 9.44^a	90.22 ± 9.09^e	95.27 ± 9.17^d	99.45 ± 9.26^b	93.85 ± 9.12^de	97.22 ± 9.23^c	100.01 ± 9.28^b	107.55 ± 9.27^a
	P2:1	159.77 ± 3.96^b	160.48 ± 4.84^ab	161.91 ± 3.54^ab	162.13 ± 2.80^ab	161.79 ± 2.90^ab	161.61 ± 6.53^ab	162.24 ± 2.81^ab	165.36 ± 13.52^a
	P3:1	189.58 ± 4.78^b	191.02 ± 5.09a^b	192.61 ± 4.04^ab	193.10 ± 3.29^ab	192.53 ± 3.47^ab	192.14 ± 6.42^b	193.16 ± 3.29^ab	196.56 ± 13.54^a
	P3:2	118.27 ± 2.12^e	118.58 ± 2.24d^e	118.82 ± 1.83^cd	119.02 ± 1.53^b	118.87 ± 1.84^abd	118.77 ± 1.89^d	118.99 ± 1.44^abc	119.14 ± 1.42^a
	Pweller	144.85 ± 3.37^e	145.73 ± 3.11^d	146.20 ± 2.75^cd	146.66 ± 2.34^b	146.31 ± 2.66^abd	146.06 ± 2.85^d	146.54 ± 2.25^bc	146.95 ± 2.11^a
	TMY	7303.04 ± 1977.90^a	4791.88 ± 2025.67^c	5439.55 ± 2073.17^b	5381.50 ± 1917.69^b	4611.45 ± 1883.70^d	4646.65 ± 1874.90^cd	5375.95 ± 1881.14^b	3745.81 ± 1637.25^e

Notes: PY: peak yield; PT: peak time; P_2:1, P_3:1 and P_3:2: milk yield persistency measures proposed by Johansson and Hansson (1940); P_weller: milk yield persistency measure proposed by Weller et al. (2006); TMY: 305-day milk yield; WD_lin: Wood’s linear model; ASH: Ali and Shaeffer’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. ^a-dEstimates of least square means with different letters in rows are significantly different (p < .05).

The values of persistency measures (P_2:1, P_3:1, P_3:2 and P_weller) were higher for nonlinear models (WD_nlin and WL) than those for linear models (WD_lin and ASH). In general, the most persistent lactation curves were found for H group, which presented the highest mean values of P_2:1 (103.69% and 98.58%), P_3:1 (89.66% and 99.34%), P_3:2 (85.21% and 97.31%) and P_weller (78.98% and 98.96%), when fitted by WD_nlin and WL, respectively. On the other hand, the linear models (WD_lin and ASH) indicated the most persistent lactation curves for G group, considering all the persistency measurements.

Considering the nonlinear models, the heterosis effect was evaluated only by the WL model, whose choise was based on the lowest R ² _Adj, RMSE and AIC values for the genetic groups with the largest number of animals, that is, 1/2H, 5/8H, 7/8H, G and H groups (Table 2). In addition, in spite of the lowest values of R ² _Adj, RMSE and AIC between the linear models, WD_lin did not present the best precision in the fit quality of the lactation curves for the different genetic groups when the trajectory of the curves was compared (Figure 1).

The effect of heterosis was significant (P < .001) for TMY, PY and the parameter ‘a’ of model WL (initial production) for cows of different genetic groups that comprise the Girolando breed (Table 6). The heterosis effect was also significant (P < .05) for PY, P_2:1 and the parameter ‘a’ of WL model. On the other hand, there were not significant (P > .05) heterosis effect for the different measures of persistency (P_3:1, P_3:2, P_weller) and for the ‘c’ parameter of WL model. The magnitude of heterosis effect of TMY and PY estimated for the different genetic groups were + 809 kg and + 4.59 kg larger than that in parental breeds. It indicates that the performance of Girolando cows were 14.64% (TMY) and 20.60% (PY) larger than the average of parental breeds.

Table 6. Heterosis effect (mean ± SE) for Wilmink parameters (‘a’, ‘b’ and ‘c’), peak yield (PY), peak time (PT), persistency measures (P_2:1, P_3:1, and P_3:2) and 305-day milk yield (TMY) estimated Wilmink nonlinear model.

Parameters, components and trait	Heterosis	Heterosis (%)	P-value
a	4.511 ± 0.412	23.05	<.0001
b	−1.861 ± 0.840	9.64	.0267
c	−0.003 ± 0.002	7.26	.1615
PY	4.594 ± 0.386	20.60	<.0001
PT	0.063 ± 0.029	0.007	.0333
P2:1	2.609 ± 1.184	2.56	.0276
P3:1	4.519 ± 2.487	5.36	.0693
P3:2	3.700 ± 2.872	4.56	.1977
Pweller	5.098 ± 2.758	6.93	.0645
TMY	809.77 ± 75.524	14.64	<.0001

Notes: PY: peak yield; PT: peak time; P_2:1, P_3:1 and P_3:2: milk yield persistency measures proposed by Johansson and Hansson (1940); P_weller: milk yield persistency measure proposed by Weller et al. (2006).

The largest heterosis effect was for ‘a’ parameter estimated by WL (associated to the initial stage of lactation), which presented an initial production of 23.05% higher than the average of parental pure breeds. The lowest heterosis effect was for PT (0.007%). Among the four persistency measures (P_2:1, P_3:1, P_3:2 and P_weller) only P_2:1 presented heterosis effect (2.56%).

4. Discussion

The Use of mathematical modeling of the lactation pattern provides an important tool for the management of milk production because the selection of animals can be based on the prediction of milk yield, considering the variation between and within genetic groups (Hossein-Zadeh 2017). Thus those models provide an indication of the nutritional management for each genetic group and the group with the highest probability of achievement of a desired level of milk yield (Hossein-Zadeh 2016). However, the studies about the modeling of lactation curve and its components should be used with caution, because the wrong choice may lead to economical losses (Pereira et al. 2016).

In this study, the individual lactation curves in different genetic groups of cows from Girolando cattle were described by linear and nonlinear models in order to study the heterosis effect on the components of the lactation curve. Those heterosis effects were obtained for the ‘a’, ‘b’ and ‘c’ parameters of mathematical model, peak yield, time to peak yield, different persistency measures. The values of the ‘a’, ‘b’, and ‘c’ parameters from WD_nlin and WD_lin were different in the genetic groups. In both models, the values of those parameters were positive, with ‘b’ and ‘c’ near to zero. Similar values were found for Holstein cattle by Jeretina et al. (2015) and Torshizi et al. (2011), using WD_nlin and WD_lin, respectively. The values of the parameters estimated by WD_lin and WD_nlin were associated to the typical pattern of the lactation curves in the different genetic groups. The typical pattern of lactation curves from Wood’s model was due to the positive values of ‘a’, ‘b’ and ‘c’, with ‘b’ between 0 and 1 (Vinay-Vadillo et al. 2012). The negative values of ‘b’ and/or ‘c’ in that model are considered a problem (Pollott and Gootwine 2000). Consequently, they are not indicated for calculating the estimates of PY, PT and TMY (Wood 1967).

In the WL model, the values of the parameter ‘a’ were positive but ‘b’ and ‘c’ were negative in genetic groups. Larger values of the ‘a’, ‘b’, and ‘c’ parameters estimated by WL for Holstein cattle were found by Torshizi et al. (2011). However, the results in this study using WL were close to those found by Pereira et al. (2016) for cows from Bos taurus taurus x Bos taurus indicus.

The ASH model presented the highest estimates for parameter ‘a’ compared to the other models, with ‘b’ and ‘d’ negative and ‘f’ close to zero. In the ASH model, the inferences based on the parameter values are not indicated, because their parameters do not present a biological meaning (Macciotta et al. 2011). The graphical presentation of lactation curves fitted by ASH model for the different groups was similar to the observed TMY (Table 1 and Figure 1).

The curves fitted by ASH model presented the typical shape of the lactation curve. Atypical shapes are characterized by the absence of the peak yield (Olori et al. 1999; Macciotta and Vicario 2005). The absence of a typical lactation curve shape may be indicative of the low quality of fit, but other factors may occur to change the shape of the lactation curve. Lactation curve shapes different from the standard may occur, especially when individual patterns are fitted (Macciotta et al. 2011). One of the possibilities is a continuously increasing or continuously decreasing curve with the absence of the peak yield, being called atypical shapes (Olori et al. 1999; Macciotta and Vicario 2005). Another possibility is a reversed shape, with an initial decreasing phase to a minimum followed by an increase, that is common for fat and protein contents (Macciotta et al. 2011). Another exception to the standard shape is represented by the existence of a second lactation peak in cows calving in autumn in pasture-based farming systems (García and Holmes, 2001).

The tests of the quality of fit indicated that the curves fitted by WD_nlin were best for the 1/4H, 3/4H, 3/8H groups while the 1/2H, 5/8H, 7/8, G and H groups presented the best fit when fitted by WL model (Table 2). Olori et al. (1999) estimated higher values of R ² _adj in WL than those in WD_nlin for Holstein cattle. Those same authors reported that R ² _adj values higher than 0.7 indicate the models with the best fit while values lower than 0.4 the worst. In regard to the linear models, WD_lin presented lower values of R ² _adj compared to ASH model for all the genetic groups. The quality of fit increased slightly with the increase in the number of parameters, which was similar to the results found by Jamrozik et al. (1997). However, WD_lin presented the lowest values of RMSE and AIC for all genetic groups (Table 3). The differences of quality of fit between models may be attributed to differences between breeds, mathematical functions, differences between test-day milk yields and the amount of available data (Khan et al. 2012). Besides the mathematical functions, the model that best fit the shape of lactation curves also depend on the calving order (Şeahin et al. 2015) and the biological nature of the parturition itself, which varies randomly between cows (Olori et al. 1999).

Despite the use of data of first parity cows from the same region, analyses were realized in data of each cow individually. It could be associated to the fact that each criterion indicated a different model as the best fit for the lactation curves of cows. Although the RMSE and AIC criteria indicated WD_lin as the best choice, that model did not fit the lactation curves with similar precision to WD_nlin, WL and ASH, when the fitted curves were compared to the observed curves (Figure 1). There is not a consensus in literature about the best model for a situation because different criteria may indicate different models (Cobuci et al. 2011). Consequently, such a choice can become a difficult task. In several studies, WD_nlin, WL and ASH models were successfully applied in the adjustment of individual lactation curves (Macciotta and Vicario 2005; Silvestre et al. 2006). The use of ASH, WL and WD_nlin models is commonly recommended in dairy cattle because they are suitable to describe the lactation curves, provided there is no limitation in the amount of available data as lack of test-day records in a certain stage of lactation or lower number of animals of groups of animals (Khan et al. 2012).

The lower number of test-day records of the 1/4H, 3/8H and 7/8H groups could have influenced the quality of fit of models, since none of the fitted curves followed perfectly the trajectory of the observed curves along days in milk (Table 1). However, despite the higher number of animals in the H group, the curves fitted by WD_nlin presented lower milk yields between 5 and 100 days in milk and higher milk yields between 101 and 305 days than those in the observed lactation curve. Torshizi et al. (2011) also found a similar result for Holsteins using Wood’s nonlinear model. Although there are limitations in the Wood’s model, one of its main advantage is the fact that it may fit the lactation curves of cows with atypical shapes (Dijkstra et al. 2010).

The estimate mean of PY and TMY differed between genetic models and groups (Table 5 and 6). Overall, the models that best estimated the PY were WD_nlin and ASH. The latter one estimated the values of TMY for 1/2H, 3/4H, 7/8H and G with higher precision, whose range were closer to the observed values (Table 1). It could be due to the fact that the models with higher number of parameters have a better performance in the quality of fit (Steri et al. 2009).

Regardless of the model, the highest mean values of PY and TMY were obtained for group H, followed by group 1/2H, with positive association between PY and TMY. Hossein-Zadeh (2016) reported that cows with high peak yields also present high 305-day milk yields compared to cows with low peak yields. However, such an association is not perfect, which could allow the selection for cows with lower peak yields and high 305-day milk yields. Thus the selection of cows could be based on the peak yields (Hossein-Zadeh, 2014). On the other hand, the selection should be combined in an index that could take into account peak yield and 305-day milk yields with other selected traits.

WD_lin model presented the lowest estimates of PY, which could be associated to the lowest values estimated for TMY. According to Prasad (2003), Wood’s nonlinear model tend to underestimate the peak yields. Similarly, Torshizi et al. (2011) indicated that Wood’s nonlinear model presented a better fit of the lactation curves of Holstein cows than the linear model. According to Pollott and Gootwine (2000), nonlinear functions present a better fit than their linear equivalent.

For WL and ASH models, the time to peak yield increased with the increase in the proportion of Holstein genes. Some studies showed the Holstein breed can reach the peak yield around 90 days in milk (Cobuci et al. 2004; Torshizi et al. 2011), while Gyr breed can reach the peak around 60 days in milk (Herrera et al. 2008). However, the time to peak estimated by WD_nlin and WD_lin, presented higher variation between genetic groups. Such a variation in the day of peak milk yield for the different genetic groups can be attributed to the choise of the model or function, and the smaller number of test-day records that influenced the shape of the lactation curve (Oliveira et al. 2007; Glória et al. 2010; Torshizi et al. 2011).

Persistency of milk yield was the parameter describing the course of lactation curve that presented the lowest variation between models and genetic groups. Persistency can be defined as the ability of the cow to maintain milk yield after achieving the maximum milk production (Hickson et al. 2006). Thus such an ability of maintaining a similar level of production along days in milk with low peak yields may dilute the costs of production along time (food and management costs) and decrease the health treatments caused by metabolic disorders (Hossein-Zadeh, 2014; Hossein-Zadeh, 2016). Additionally, the welfare of cows is improved as a result of the increase of persistency (Cole and VanRaden 2006).

In general, the persistency measures estimated by WD_nlin and WL showed that the most persistent lactation curves were for H group, which also presented the highest level of production among genetic groups. That result was corroborated by Gengler (1996), which reported that persistency was associated to the level of production. However, the estimates of persistency by linear models (WD_lin and ASH) were more persistent for G group. In the estimates of ASH model, the results for persistency may be occurred because of the highest estimates of PT for G group compared to H group as well as the high correlation between persistency and PT (Albarran-Portillo and Pollott 2011). Thus it suggests that peak yields that occur later are associated to more persistent lactation curves.

Many studies have revealed that there is an interest of heterosis effect on milk yield (Norberg et al. 2014; Lembeye et al. 2015). In Brazil, the effect of heterosis has been indicated as the most important aspect of milk yield in the different genetic groups of the Girolando cattle (Facó et al. 2005; Facó et al. 2008). However, there is a lack of studies evaluating the effect of heterosis in Girolando cows, especially in regard to aspects associated to the shape of the lactation curve.

Between the nonlinear models, WL model was chosen in order to estimate the effect of heterosis for cows of different Girolando breed groups, whose choice was based on the quality of fit criteria and the shape of lactation curves. Additionally, the goodness of fit of groups with the largest number of animals (1/2H, 5/8H, 7/8H, G e H) was higher for WL compared to WD_nlin.

Thus the heterosis effect was significant (P < .001) for PY, TMY and the ‘a’ parameter of WL model. The heterosis effect of ‘a’, PY and TMY were 4.51, 4.59, and 809.77 kg, which represent 23.05%, 20.60%, and 14.64% larger yield values than the mean between parental breeds, respectively. In fact, the magnitude of heterosis depends on the degree of genetic dominance of the trait, but it is also related to the genetic distance between the parental breeds, so that, in general, the higher the genetic distance, the higher the heterosis effects (Mäki-Tanila 2007).

Several studies investigated the effect of heterosis on milk yield in different crossbred animals that involved the Holstein breed. Akbas et al. (1993) and Boichard et al. (1993) reported heterosis effect of 135 and 104 kg of milk yield in Holstein x French Black and White cattle and Holstein and European Friesian crossbred cows, respectively. Recently, Penasa et al. (2010) and López-Villalobos et al. (2010) found heterosis effect of 477 and 496 kg more milk yield in Holstein x Jersey crossbred cows than the average of milk yield of pure breeds.

The results of heterosis effect for initial production (‘a’ parameter), PY, and TMY indicated that cows that had heterosis effect for initial production also presented higher 305-day milk yield. One of the factors that may have influenced the heterosis effect of TMY was the fact that ‘a’ parameter could be associated to the level of production (Wilmink 1987). Another factor is that PY and TMY were positively associated, once the results showed that cows with higher PY also presented higher TMY (Table 4). Several previous studies have reported that there is a genetic correlation between PY and TMY, ranging from 0.82–0.90 in Holstein cows (Shanks et al. 1981; Rekaya et al. 2000). Similarly, Buckley et al. (2003) found that cows with the highest peak milk yield were those with the highest 305-day milk yield. The results of the heterosis effect found for initial production (parameter ‘a’), PY and TMY may also indicate the adaptation of the Girolando cows to the environmental conditions of Brazil. In this case, the term adaptation should be considered in a broad sense because the genes expressed in the Girolando genotypes provided the best performance of those animals (McManus et al. 2008).

Many studies have reported that the environmental causes an effect on the expression of heterosis (Bryant et al. 2007; Penasa et al. 2010). That possibility could be a factor that may have influenced (or at least could explain) the higher heterosis expression of Girolando cows for initial production (parameter ‘a’), PY and TMY. On the other hand, the ‘c’ parameter from WL model did not present a significant effect of heterosis. It is known that ‘c’ parameter is associated to the increase in the milk yield up to the peak yield (Wilmink 1987). That result permitted to infer that the effect of heterosis in PY was independent of the effect of heterosis in parameter ‘c’, both components associated with the initial phase of the lactation curve. There was also no effect of heterosis for parameter ‘b’ in cows of different genetic groups of the Girolando breed. The ‘b’ parameter is associated to the decrease of milk yield after peak yield (Wilmink 1987), which may be directly associated to the persistency of milk yield.

Similarly, most of the persistency measures (P_3:1, P_3:2 and P_weller) did not show heterosis effect (P > .05), indicating that cows from the different Girolando genetic groups had less persistent lactations when compared to the mean persistency of the parental breeds (H and G). Thus, the absence of heterosis effect on ‘b’ parameter was associated to the absence of heterosis effect on the different persistency measures (P_3:1, P_3:2 and P_weller), because all represent the same component of the lactation curve (persistency in milk production). Another factor that may have influenced the absence of the heterosis effect in the different measurements of persistency was the small effect of heterosis in PT in this study.

Němečková et al. (2015) estimated a genetic correlation of 0.54 ± 0.07 between the day of peak and persistency in Holstein cows. In contrast, a significant heterosis effect (P < .05) was found for P_2:1 in Girolando cows, although with no expressive effect (2.56%) when compared to the mean persistency of parental breeds (H and G). Part of the explanation could be attributed to the fact that the initial production (‘a’ parameter) presented positive heterosis effect. It is known that the level of the initial stage of the lactation curve influences the cow’s ability to maintain milk yield after the peak of milk yield (Wood 1967). It is worth mentioning that P_2:1 is the measure more associated to the first two stages or parts of lactation because the lactation has been divided in three parts when persistency had been calculated. It could explain the heterosis effect on that persistency measure. Therefore, the higher heterosis effect was found for the initial stage of milk yield (parameter ‘a’ of WL model), followed by PY and TMY. The significant effect of heterosis on TMY is commonly expected due to the genetic distance between Holstein breed (Bos taurus taurus) and Gyr (Bos taurus indicus).

However, the most important benefit of the heterosis effect obtained in this type of crossing may be due to the effect of heterosis on the components of the curve which were associated mainly with initial stage of milk yield and milk yield at the peak and with a lesser extent with the two final stages (persistency). The heterosis effect in the first stages of lactation is the most important benefit of crossbreeding in the adaptation of dairy cattle in the tropical conditions of many countries as Brazil. Farmers with animals more adapted to their environment are more competitive due to the lower costs with labour, feed and higher productivity (López-Villalobos et al. 2000; McManus et al. 2008).

5. Conclusions

The quality of fit of the mathematical models was different for all the different models and genetic groups. Based on the criteria of goodness of fit, the results of this study showed that WL provided the best fit of lactation curve for most genetic groups between nonlinear models. In regard to the linear models, WD_lin provides a better adjustment of the lactation curve for the crossbred genetic groups of cows. The crossbreeding between Holstein and Gyr breeds may bring many benefits for the producers, since the cows from the different crosses between them present heterosis higher than 14% for 305-day milk yield. The heterosis effect is more associated to the initial stage of lactation curve and the peak yield.

Acknowledgements

We would like to thank Associação dos Criadores de Gado Holandês de Minas Gerais (ACGHMG), Associação Brasileira dos Criadores de Gir Leiteiro (ABCGIL), Associação Brasiliera de Girolando (ABCG) and EMBRAPA Gado de Leite for providing the valuable data.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

Darlene dos Santos Daltro http://orcid.org/0000-0003-2278-8129

Marcos Vinícius Gualberto Barbosa da Silva http://orcid.org/0000-0001-5449-1413

Luís Telo da Gama http://orcid.org/0000-0003-3894-3488

Jaime Araujo Cobuci http://orcid.org/0000-0002-1414-4967

Additional information

Funding

This study was financed in part the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) – Finance Code 001.