Estimates of genetic and crossbreeding parameters for 305-day milk yield of Girolando cows

Abstract

The objective of this work was to evaluate the effects of breed and heterosis on the estimates of genetic parameters of 305-day milk yield (305MY), by fitting milk yield data to the Wood, Mixed Log, Morgan, and Wilmink non-linear models. A total of 258,891 test-day milk yield records of 37,965 Holstein (H), Gyr (G), and Girolando (1/2 H, 1/4 H, 3/4 H, 3/8 H, 5/8 H and 7/8 H genetic groups) cows from 1840 herds were collected in the state of Minas Gerais, Brazil, in the period of 1998 to 2014. The pedigree file included 36,640 animals, 3677 bulls and 24,472 cows. The genetic parameters were obtained by one-character analysis using the AIREMLF90 software. The heritability estimates of 305MY, and of 305MY fitted to the different nonlinear models ranged from 0.14 to 0.20. Most genetic groups had heterosis effect, with the greatest effects for 305MY on 1/2 H (1/2 Holstein × 1/2 Gyr) cows, presenting a milk yield of 1112.73 kg higher than the average of their parents. The results indicate that most genetic groups of Girolando present heterosis effect for 305MY estimated through nonlinear models, and the possibility of obtaining moderate genetic gains when these animals are subjected to selection procedures.

HIGHLIGHTS

• The use Girolando cows in production systems in Brazil has shown positive results of production.

• This greater productive efficiency may be because of heterosis on the production and adaptation of these animals.

• The effect of breed and heterosis can be included in traditional models for the improvement of selection.

Keywords:

• Dairy cows
• heritability
• heterosis
• mathematical model

Introduction

Strategies for crossbreeding Bos taurus taurus with Bos taurus indicus are important in Brazil to achieve the balance between yield and adaptability of animals to adverse conditions in the country (Santos et al. 2013) and explore the benefits of heterosis and complementarity originated from the crossbreeding between different breeds (Canaza-Cayo et al. 2014).

The first crosses between Gyr and Holstein cattle intending to form the Girolando breed started in the 1940s (Facó et al. 2002; Silva et al. 2015). The use of this breed in production systems in Brazil has shown positive results of production. Silva et al. (2016) reported an increase of 41% in milk production of Girolando cows between 2000 and 2014. Although Girolando cows present lower milk production than Holstein cows, they present greater resistance to ectoparasites and hot weather and greater productive efficiency when farmed in medium to low technological production systems (Facó et al. 2002; McManus et al. 2008). This greater productive efficiency may be because of heterosis on the production and adaptation of these animals. This shows the need for conducting studies on the effect of heterosis in milk production of Girolando cows.

The use of mathematical models to describe milk production during lactation in genetic breeding programmes (Hossein-Zedeh 2017) allows the establishment of strategies to optimising the selection, and the search for efficient and profitable genotypes (Oliveira et al. 2007). Moreover, these models allow the evaluation of different components of the lactation curve of the animals, such as milk yield peak, yield peak time, and lactation persistency. Knowing the mathematical model that best fits the lactation curve enables the estimation of genetic parameters and prediction of genetic values. Moreover, the effect of breed and heterosis can be included in traditional models for the improvement of selection procedures for the traits of economic importance that are affected by these variables.

In the literature, there are several models widely used to describe lactation curves of different breeds. There are also those models that have not yet been used for this purpose, especially in the Girolando breed. Therefore, several of those models were used to study heterosis in this breed. Therefore, the objective of this work was to evaluate the effects of breed and heterosis on the estimates of genetic parameters of 305-day milk yield (305MY), by fitting test day milk yield records to the Wood, Mixed Log, Morgan, and Wilmink non-linear models.

Material and methods

Data

Milk production data were obtained from the Association of Holstein Breeders of Minas Gerais, Brazilian Association of Dairy Gyr, and Brazilian Association of Girolando. The dataset consisted of 258,891 test day milk yield records of 37,965 cows from 1840 herds, which were collected in the state of Minas Gerais, Brazil, in the period of 1998 to 2014. These cows were from the Holstein (H), Gyr (G), and Girolando breeds. Six genetic crossbreedings of Holstein (H) × Gyr (G) crossbreeds (1/4 H, 3/4 G (1/4 H); 3/8 H, 5/8 G (3/8 H); 1/2 H, 1/2 G (1/2 H); 5/8 H, 3/8 G (5/8 H); 3/4 H, 1/4 G (3/4 H); 7/8 H, 1/8 G (7/8 H)) were used; animals from these groups are officially called Girolando in Brazil.

The effect of breed and heterosis on the 305-day milk yield (305MY) was evaluated by fitting the test day milk yieldsrecords to the Wood (305MY-WD), Mixed Log (305MY-ML), Morgan (305MY-MG), and Wilmink (305MY-WL) models.

Only first lactation cows were maintained for the analyses, considering complete lactations (305 days in milk) when the cow had a minimum of 4 and maximum of 10 test days. The test-day milk yield, and 305-day milk yield records were not considered when the daily milk yield was not between 3 kg and 45 kg. The production data were edited and the breed group, heterosis coefficient, genetic group of the animal, genetic group of the sire, genetic group of the dam, and number of cows were determined (Table 1).

Table 1. Number of cows per genetic group according to heterosis coefficient.

Breed group	Heterosis coefficient	Genetic group	Genetic group (sire)	Genetic group (dam)	Number of cows	Number of sires per cow genetic group
1	0	H	H	H	10704	545
		G	G	G		622
2	0.250	7/8H	H	3/4H	770	48
			3/4H	H
3	0.375	3/4H	3/4H	7/8H	357	202
4	0.406	3/4H	5/8H	7/8H	278
5	0.438	5/8H	5/8H	3/4H	616	140
6	0.469	5/8H	5/8H	5/8H	455
7	0.500	3/4H	H	1/2H	3130	202
		5/8H	3/4H	1/2H		140
		1/4H	1/2H	1/2H		15
8	0.625	5/8H	H	3/8H	385	140
		1/2H	3/4H	1/4H		157
9	0.750	5/8H	H	1/4H	1009	140
		3/8H	G	3/4H		14
10	1	1/2H	G	H	2484	157
			H	G

H: Holstein; G: Gyr.

The contemporary groups (CG) were defined by herd-year of calving in the 305MY analysis. Only CG with at least three controls of cows from at least two bulls were considered. Two calving periods were considered, October to March and April to September. The pedigree file included 36,640 animals, 3677 bulls and 24,472 cows. The birth period of animals in the pedigree file was 1979 to 2009, and included the Holstein, Gyr, and Girolando breeds.

Lactation curve models

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

1. Wood gamma model (WD) (Wood 1967): $$Y_t={at}^b{e}^{-ct}$$

2. Mixed Log (ML) (Guo and Swalve 1995): $$Y_t={a+bt}^{0.5}+c\log t$$

3. Morgan model (MG) (Thornley and France 2007): $$Y_t=a{b}^{c}c\frac{t^{c-1}}{(t^c+b^c{)}^2}$$

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

wherein $Y_t$ is the average daily milk yield in the t^th test day of lactation, a is the initial milk yield just after calving, and b and c are the parameters of the ascending phase to the production peak, and descending phase after the production peak, respectively. The constant k in the WL model was determined in a previous analysis as 0.05.

Statistical analyses

The normal distribution test of milk production levels for different genetic groups was applied. Data were tested for non-normality by the Shapiro–Wilk and Kolmorgorov–Smirnov test. The non-linear models were adjusted to the milk yield records as the iteration method of Gauss–Newton.

The models were tested using root mean square error (RMSE), Akaike’s information criterion (AIC), and Bayesian information criterion (BIC).

RMSE (generalised standard deviation) was calculated using the equation: $$\mathit{RMSE}=\sqrt{\frac{\mathit{RSS}}{n-p-1}},$$ wherein RSS is residual sum of squares, n is the number of observations (data points), and p is the number of parameters of the equation.

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

BIC was calculated using the equation: $$\mathit{BIC}=n\mathrm{ }ln\left(\frac{\mathit{RSS}}{n}\right)+p\ln \left(n\right).$$

The lowest numerical values of RMSE, AIC and BIC indicated the better fit when comparing the different models.

Breed and heterosis effects

The proportion of genes was calculated for each cow using the equation: $${\alpha }_i^p=\mathrm{ }\frac{\left({\alpha }_i^s+\mathrm{ }{\alpha }_i^d\right)}{2}$$ where ${\alpha }_i^p$ is the porportion 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.

These specific effects of heterosis were used for the six genetic groups of Girolando, because the distribution of cows across 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 previously calculated coefficients of specific heterosis.

Genetic parameters

The 305MY was evaluated by fitting the test day milk yields records to the 305MY-WD, 305MY-ML, 305MY-MG, and 305MY-WL models, later they were analysed by means of the animal model, in one-character analysis. The components of variance were estimated by the maximum likelihood-restricted method, using a derivative-free algorithm, using the software AIREMLF90 (Misztal et al. 2014). The convergence criterion was 1 × 10⁻¹².

Model and analysis

The model included the direct additive genetic effect and the residual effect as random, the systematic effects of GC, breed, heterosis, and the covariant cow age (linear and quadratic effect). The model can be represented in matrix form by the equation: $$y=X\beta +Zu+e.$$ wherein $y$ is the vector of records for 305MY; $\beta$ is the vector of the fixed effects, $u$ is the vector of random animal effects, including animals without records; $e$ is the vector of random residual effects, e ∼ NID (0, ${\sigma }_e^2$); and $X$ and $Z$ are incidence matrices assigning observations to fixed and random animal effects, respectively.

Results

Considering all models, the lowest mean RMSE (Table 2) was found for the G group, followed by the 7/8 H, and 5/8 H groups. The lowest mean AIC and BIC were found in the 1/4 H group.

Table 2. Quality-of-fit of 305-day milk yield means to Wood (WD), Mixed Log (ML), Morgan (MG) and Wilmink (WL) models.

Genetic group	Model	RMSE	AIC	BIC
H	305MY-WD	5.90	774531	122726
	305MY-ML	6.23	774628	195018
	305MY-MG	6.47	778537	198927
	305MY-WL	5.91	774700	195089
G	305MY-WD	5.72	477960	122734
	305MY-ML	6.06	477952	122726
	305MY-MG	6.12	478609	123382
	305MY-WL	5.72	477953	122726
1/2H	305MY-WD	6.49	259045	72604
	305MY-ML	6.79	259045	72604
	305MY-MG	6.90	259605	73164
	305MY-WL	6.49	259011	72570
1/4H	305MY-WD	6.41	18359	6260
	305MY-ML	6.73	18359	6260
	305MY-MG	6.82	18404	6305
	305MY-WL	6.41	18357	6258
3/4H	305MY-WD	5.90	386369	101832
	305MY-ML	6.24	386405	101869
	305MY-MG	6.36	387516	102979
	305MY-WL	5.91	386426	101889
3/8H	305MY-WD	6.55	18384	6310
	305MY-ML	6.86	18385	6312
	305MY-MG	6.93	18420	6346
	305MY-WL	5.91	386426	101889
5/8H	305MY-WD	5.88	259710	70240
	305MY-ML	6.21	25712	70242
	305MY-MG	6.30	260227	70757
	305MY-WL	5.88	259712	70604
7/8H	305MY-WD	5.74	80476	23478
	305MY-ML	6.08	80484	23486
	305MY-MG	6.23	80788	23790
	305MY-WL	5.74	80484	23486

AIC: Akaike's information criterion; BIC: Bayesian information criterion; H: Holstein; G: Gyr; RMSE: root mean square error; 305MY-WD: 305-day milk yield estimated by Wood model; 305MY-ML: 305-day milk yield estimated by Mixed Log; 305MY-MG: 305-day milk yield estimated by Morgan model; 305MY-WL: 305-day milk yield estimated by Wilmink model.

The highest milk production was found for cows of the H group (Table 3), followed by the 7/8 H, 1/2 H, and 3/4 H groups, regardless of the 305MY fitting. Contrastingly, cows in the G group had the lowest performance for 305MY, 305MY-WD, 305MY-ML, 305MY-MG and 305MY-WL.

Table 3. Least mean squares and standard errors for 305-day milk yield (kg) in different non-linear models for Holstein cows (H), Gyr (G), and different genetic groups of the Girolando breed.

Genetic group^1	305MY	305MY-WD	305MY-ML	305MY-MG	305MY-WL
H	6982.98^a ± 22.21	6988.05^a ± 22.24	6971.19^a ± 22.15	7016.90^a ± 22.92	6970.62^a ± 22.22
G	4256.04^d ± 26.11	4473.07^f ± 26.15	4396.01^f ± 26.04	4610.26^f ± 26.94	4380.73^f ± 26.12
1/2H	5118.23^b ± 34.99	5631.76^b ± 35.04	5592.60^b ± 34.90	5727.87^b ± 36.11	5593.45^b ± 35.01
1/4H	4307.01^c,d ± 118.30	4865.70^d,e ± 118.69	4819.26^d,e ± 117.99	4876.04^d,e ± 122.06	4784.10^d,e ± 118.36
3/4H	5076.57^b ± 30.12	5414.94^c ± 30.18	5388.22^c ± 30.05	5509.51^c ± 31.08	5392.67^c ± 30.14
3/8H	4569.97^c ± 120.52	4999.65^d ± 120.65	4933.99^d ± 120.20	5116.16^d ± 124.35	4935.03^d ± 120.58
5/8H	4388.13^c ± 35.53	4721.97^e ± 35.58	4683.66^e ± 35.44	4771.54^e ± 36.66	4684.29^e ± 35.55
7/8H	5211.12^b ± 63.08	5502.03^b,c ± 63.19	5470.28^b,c ± 62.92	5578.15^c ± 65.09	5465.48^b,c ± 63.12

¹H: Holstein, G: Gyr.

^a–fEstimates of least mean square with different letters in the columns are significantly different (p < 0.05). 305MY: 305-day milk yield; 305MY-WD: 305-day milk yield estimated by Wood model; 305MY-ML: 305-day milk yield estimated by Mixed Log; 305MY-MG: 305-day milk yield estimated by Morgan model; 305MY-WL: 305-day milk yield estimated by Wilmink model.

The genetic variance estimates for 305MY were higher than those obtained by the 305MY-WD, 305MY-ML, and 305MY-WL (Table 4). Similar results were found for estimates of residual variance. Heritability estimates for 305MY (0.14 ± 0.01), and similar estimates were found for 305MY-WD (0.18 ± 0.02), and 305MY-WL (0.17 ± 0.02). The highest estimates of heritability were found for 305MY-MG (0.20 ± 0.02).

Table 4. Estimates of variance components and heritability (h²) of 305-day milk yield (kg) estimated by different nonlinear models.

	${\sigma }_{\text{a}}^2$	${\sigma }_{\text{e}}^2$	${\sigma }_{\text{total}}^2$	h^2
305MY	264980	1654300	1919280	0.14 ± 0.01
305MY-WD	249620	1088300	1337920	0.18 ± 0.02
305MY-ML	241010	1097100	1097100	0.18 ± 0.02
305MY-MG	269590	1052300	1321890	0.20 ± 0.02
305MY-WL	237540	1103400	1340940	0.17 ± 0.02

305MY: 305-day milk yield; 305MY-WD: 305-day milk yield estimated by Wood model; 305MY-ML: 305-day milk yield estimated by Mixed Log; 305MY-MG: 305-day milk yield estimated by Morgan model; 305MY-WL: 305-day milk yield estimated by Wilmink model; σ²_a: animal additive genetic variance; σ²_e: residual error variance; σ²_total: sum of all variances.

The effect of breed (Table 5) showed that only cows of group 7/8 H presented positive performance (440.95 kg to 305MY). The effects of breed found in the different genetic groups for 305MY-ML, 305MY-MG, and 305MY-WL were similar to those found for 305MY. The 305MY-WD estimate was the one that most differed from the 305MY.

Table 5. Breed effect with standard error estimated by the Wood (WD), Mixed Log (ML), Morgan (MG), and Wilmink (WL) models.

Genetic group (sire)	Genetic group (dam)	Genetic group	305MY	305MY-WD	305MY-ML	305MY-MG	305MY-WL
1/2H	1/2H	1/4H	−1312.89 (140.81)	219.46 (1802.85)	−1243.98 (137.04)	−1250.43 (143.73)	−1258.80 (138.13)
G	3/4H	3/8H	−1595.99 (243.29)	−28.13 (1813.64)	−1462.73 (241.21)	−1293.47 (253.22)	−1427.64 (243.24)
H	G	1/2H	−974.95 (936.99)	150.68 (1890.29)	−1342.91 (949.37)	−1252.08 (997.43)	−1352.99 (957.80)
5/8H	5/8H	5/8H	−521.96 (224.89)	985.17 (1814.07)	−479.14 (225.56)	−550.86 (236.88)	−489.85 (227.51)
3/4H	1/2H	5/8H	−646.08 (185.47)	892.55 (1808.24)	−488.35 (184.57)	−485.94 (193.80)	−477.04 (186.15)
H	1/4H	5/8H	−476.25 (584.83)	803.98 (1702.79)	−645.58 (582.92)	−709.42 (612.12)	−636.27 (587.93)
H	1/2H	3/4H	0 (0)	1447.84 (1799.05)	0 (0)	0 (0)	0 (0)
H	3/4H	7/8H	440.95 (129.17)	1292.57 (1702.19)	333.49 (121.67)	314.36 (127.49)	323.87 (122.57)

H: Holstein; G: Gyr; 305MY: 305-day milk yield; 305MY-WD: 305-day milk yield estimated by Wood model; 305MY-ML: 305-day milk yield estimated by Mixed Log; 305MY-MG: 305-day milk yield estimated by Morgan model; 305MY-WL: 305-day milk yield estimated by Wilmink model.

The effect of heterosis (Figure 1) was higher for animals with heterosis coefficients of 1.0, represented by the 1/2 H group, with superior milk yields to their parents 1112.73 kg, 875.09 kg, 1421, 46 kg, 1335.83 kg, and 1433.35 kg for 305MY, 305MY-WD, 305MY-ML, 305MY-MG, and 305MY-WL, respectively. The estimated heterosis decreased with decreasing the heterosis coefficient, which was zero for the heterosis coefficient of 0.250 (7/8 H). The animals with heterosis coefficients of 0.375 and 0.406 were from the same genetic group (3/4 H), but from different crosses (Table 1) and did not present heterosis effect for 305MY, 305MY-WD, 305MY-ML, 305MY-MG, and 305MY-WL. However, the animals with heterosis coefficients of 0.438, 0.469, and 0.500 represented by the 5/8 H group had positive heterosis effect for 305MY, 305MY-ML, 305MY-MG, and 305MY-WL.

Figure 1. Estimation of heterosis (kg) for 305-day milk yield (305MY) according to the heterosis coefficient of the different genetic groups of Girolando cows, obtained by the Wood’s model (305MY-WD), Mixed Log function (305MY-ML), Morgan model (305MY-MG) and Wilmink model (305MY-WL).

Different from other models, the 305MY-WD indicated that only animals with heterosis coefficients of 0.625 (566.86 kg), 0.750 (404.49 kg), and 1.0 (875.09 kg) were affected by positive heterosis.

The linear regression of the heterosis curves of the animals with heterosis coefficients of 0.250 to 1.0 presented trends of 1837.70 kg, 1800.10 kg, 2173.3 kg, 2080.10 kg, and 2203.20 kg for 305MY, 305MY-WD, 305MY-ML, 305MY-MG, and 305MY-WL, respectively (Table 6). The R² values ranged from 0.64 to 0.81. The WD model underestimated the effect of heterosis for all animals, regardless of their heterosis coefficients (Figure 1).

Table 6. Linear regression of the heterosis curve.

	Model	R^2
305MY	y = 1837.7x − 518.18	0.818
305MY-WD	y = 1800.1x − 939.65	0.646
305MY-ML	y = 2173.3x − 705.35	0.815
305MY-MG	y = 2080.1x − 667.37	0.745
305MY-WL	y = 2203.2x − 731.10	0.812

305MY: 305-day milk yield; 305MY-WD: 305-day milk yield estimated by Wood model; 305MY-ML: 305-day milk yield estimated by Mixed Log; 305MY-MG: 305-day milk yield estimated by Morgan model; 305MY-WL: 305-day milk yield estimated by Wilmink model.

Discussion

The breed composition of cows affected the 305-day milk yield (305MY) (Table 3); cows with higher proportion of Holstein genes (H, 7/8 H, 3/4 H, and 1/2 H) had higher milk yield than cows with higher proportion of Gyr genes (1/4 H and 3/8 H). The mean 305MY of Holstein cows was 2726.94 kg higher than that of Gyr cows.

The effect of genetic groups on 305MY in Girolando cows in Brazil was also observed by McManus et al. (2008), and Balancin Júnior et al. (2014). Balancin Júnior et al. (2014) identified the genetic groups 7/8 H, 3/4 H and 1/2 H as the ones with the highest milk production. McManus et al. (2008) found higher mean 305MY for the 3/4 H group, followed by the 1/2 H, 3/8 H and 1/4 H.

When test day milk yield records were fitted to the different nonlinear models (WD, ML, MG, and WL), the highest averages were found for the H group, followed by the 7/8 H, 1/2 H, and 3/4 H; and the lowest were found for the G group. The different models (WD, ML, MG, and WL) followed the same trend of the 305MY. However, when evaluating the quality-of-fit criteria, larger differences for 305MY estimated by nonlinear models were found, with differences between fitting models varying according to the genetic group.

The WD model presented the lowest differences between mean 305MY estimates for the H and 7/8 H groups. This result agrees with the RMSE, AIC and BIC found, which were lower in these groups (H and 7/8 H) through the WD model. Torshizi et al. (2011) evaluated different non-linear models in primiparous Holstein cows in Irãn and found a lower RMSE for the WD model when compared to the WL model, showing a better prediction of milk production. The Wood model is one of the best and most popular mathematical models to describe lactation curves in dairy cattle (Gradiz et al. 2009; Macciotta et al. 2011).

The WL model was the one that best described the 305MY for the G, 1/2 H, and 1/4 H groups; it presented the lowest differences between the mean 305MY estimates and lowest RMSE, AIC, and BIC. Contrastingly, Bangar and Verma (2017) evaluated different non-linear models for Gyr cows in India and found better prediction using the WD model, in terms of RMSE and AIC, when compared to the WL model. WD and WL models have been successfully used for fitting individual lactation curves of dairy cattle (Macciotta et al. 2005; Torshizi et al. 2011).

The ML and WL models were the ones that best described the milk production of cows of the 3/4 H, 3/8 H, and 5/8 H groups, regarding the difference between the production estimated by the models and the observed production (305MY). However, they presented a high quality-of-fit variation for these genetic groups when tested by the RMSE, AIC and BIC criteria. This variation in the quality-of-fit of the models for the lactation curve depends not only on the mathematical functions, but on the calving order of the cows (Şeahin et al. 2015), and biological nature of the lactation itself, which varies randomly between the cows (Gantner et al. 2010). Despite only primiparous cows from the same region were used, the analyses were performed with different genetic groups and individual lactations; this may have caused these quality-of-fit variations.

The 305MY and its variation was estimated by different mathematical models; thus, it probably affected the variation of the heritability estimates, which were moderated. These moderate and expected heritability indicates that the additive genetic action had limited effect, and genetic and environmental factors were important in its expression (Dangar and Vataliya 2017). Similarly, Kim et al. (2009) found heritability of 0.14 to 305MY in Holstein cows. Low heritability for 305MY were also found by Singh et al. (2001) (0.11) and Ulmek (1990) (0.12) for the Gyr breed. Higher heritability for 305MY in Girolando cows were reported by Canaza-Cayo et al. (2018) (0.27) and Facó et al. (2008) (0.21). The differences in heritability estimates found in these studies can be attributed to several factors, such as production level, population size, analysis model, variable assessed (fitted, or total milk production), and environmental effects (Canaza-Cayo et al. 2018). Regardless of the variations between heritability estimates by the mathematical models for 305MY, all estimates showed moderate heritability, that is, these mathematical models can be used to estimate genetic parameters of milk production at 305 days. Therefore, effects of heterosis on other characteristics can be studied, such as the peak and persistence of lactation, in further genetic evaluations fitted to nonlinear models.

The effect of heterosis is important in different Girolando genetic groups in Brazil, when evaluating 305-day milk yield (Facó et al. 2002, 2008). Animals with heterosis coefficient of 1.0 (1/2 H) presented a greater heterosis effect to 305MY, as well as when fitted to the different nonlinear models (305MY-WD, 305MY-ML, 305MY-MG and 305MY -WL) (Figure 1). This result was expected because heterosis is more intense in this group because the parental breeds were more genetically distant (Wentworth 1927). When parental breeds have different alleles, or frequency of different alleles, offspring will show greater heterogeneity and heterosis compared to crosses between parental breeds with similar allele frequencies (Sorensen et al. 2008). This is fully expressed in animals 1/2 H, in which all pairs of genes have one gene from each breed.

Heterosis effects for 305MY, 305MY-ML, 305MY-MG and 305MY-WL were observed in all heterosis coefficients, except for animals with coefficients of 0.250, 0.375 and 0.406 (Figure 1). This indicates that most genetic groups had higher average milk production at 305 days than their parents. These results demonstrate the importance of heterosis in the expression of milk production of cows in populations composed of animals of different genetic groups, such as the Girolando breed in Brazil. Thus, it is noteworthy that the inclusion of the heterosis effect in the models used in the genetic evaluation of the breed to predict the genetic values of bulls is of fundamental importance, possibly increasing the prediction accuracy.

In practical terms, such results reveal to breeders and industry that it is possible to increase performance gains through simple processes of exploiting the effect of heterosis among zebu. If this effect can be obtained by careful combination of genetic material (animals of different breed proportions), this gain can be achieved without raising production costs.

According to McManus et al. (2008), animals that showed heterosis effects may have better adaptation to the environment to which they were subjected. The term adaptation in this case should be understood in a broad sense. It means that the types of genes acting on the genotypes of Girolando animals make them to have better performances.

Contrastingly, animals from the 3/4 H group from different crosses (Table 1) with heterosis coefficients of 0.375 and 0.406 did not show positive heterosis effect for 305MY, or for the 305MY fitted to the different non-linear models. This may have occurred due to the loss of positive epistatic effects resulting from the gene recombination that occurs by crossings (Facó et al. 2008). Breakdown of favourable epistatic interactions may occur in the recombination process during meiosis when the breeds involved in the crossing are selected in different directions for a long time (Dickerson 1973). According to Facó et al. (2008), the genetic recombination observed in some types of mating produces detrimental effects to milk production.

Canaza-Cayo et al. (2018) found better productive performance for the 3/4 H and 7/8 H groups and attributed this result to the improved environment of these herds, which favoured the performance of animals with higher proportion of Holstein genes. However, animals with heterosis coefficients of 0.438 and 0.469 (5/8 H group) from different crosses showed a positive heterosis effect (Table 1). The 5/8 H group is composed of pure animals of the Girolando breed and come from different parents. Thus, heterosis does not guarantee that crosses can bring favourable epistatic interactions. In this case, a subsequent selection programme is important because the large genetic compositions of the Girolando breed can guarantee moderate genetic gains for milk production and contribute to the improvement of the productive efficiency of the herd when properly used in breeding programmes (Canaza-Cayo et al. 2018).

The evaluation of the heterosis effect showed that only the Wood model differed from the other nonlinear models (ML, MG, and WL) when estimating milk production at 305 days (Figure 1). This denotes that the estimates obtained by the other models (ML, MG and WL) were more similar to 305MY than that obtained by the WD model. The quality-of-fit of the models is shown by several factors, such as the differences between races (Khan et al. 2012), and individual production, which varies between lactations of the same animal (Gantner et al. 2010). Therefore, not all cows in the same herd have the same lactation curves, especially when there is a genetic distance such as that between the H (Bos taurus taurus) and G (Bos taurus indicus) groups.

The lactation curve found for Holstein cows can be divided into three phases: ascending, peak, and descending (Cruz et al. 2009). However, it is different in zebu breeds and their crosses, the lactation curve is generally composed of two segments, the initial production, and the rate of decline (Cobuci et al. 2001). Thus, WD model presented the best fit for milk production records on the day of the control of animals of the H and 7/8 H group, that is, for the animals with a higher proportion of Holstein genes. This may have caused the WD model to have a poor estimating of the heterosis effect of the different genetic groups that make up the Girolando breed when compared to the other models studied.

Conclusions

The estimated lactation curve for cows of different Girolando genetic groups is dependent on the mathematical model chosen. It is utopian to expect a single model to properly adjust the lactation curves of various genetic groups. The effect of heterosis on milk yield was important for most of the breed genetic groups, and 1/2 H 1/2 G animals presented higher heterosis among the groups. A careful combination of genetic material (animals of different racial proportions) enables breeders and the dairy industry to benefit from the performance gains generated by heterosis without rising production costs. Considering the observed importance of heterosis it is essential that the genetic evaluation processes of the breed consider the heterosis effect to predict the genetic values of bulls and cows in order to increase their prediction accuracy.

Ethical approval

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

Acknowledgments

The authors thank the Associação dos Criadores de Gado Holandês de Minas Gerais (ACGHMG); the Associação Brasileira dos Criadores de Gir Leiteiro (ABCGIL); the Associação Brasileira de Girolando (ABCG); the Empresa Brasileira de Pesquisa Agropecuária – EMBRAPA Gado de Leite for providing the data and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES), for funding part of this study – Finance Code 001.

Disclosure statement

No potential conflict of interest was reported by the authors.