Variance components and genetic parameters for milk production and lactation pattern in an ethiopian multibreed dairy cattle population.

The objective of this study was to estimate variance components and genetic parameters for lactation milk yield (LY), lactation length (LL), average milk yield per day (YD), initial milk yield (IY), peak milk yield (PY), days to peak (DP) and parameters (ln(a) and c) of the modified incomplete gamma function (MIG) in an Ethiopian multibreed dairy cattle population. The dataset was composed of 5,507 lactation records collected from 1,639 cows in three locations (Bako, Debre Zeit and Holetta) in Ethiopia from 1977 to 2010. Parameters for MIG were obtained from regression analysis of monthly test-day milk data on days in milk. The cows were purebred (Bos indicus) Boran (B) and Horro (H) and their crosses with different fractions of Friesian (F), Jersey (J) and Simmental (S). There were 23 breed groups (B, H, and their crossbreds with F, J, and S) in the population. Fixed and mixed models were used to analyse the data. The fixed model considered herd-year-season, parity and breed group as fixed effects, and residual as random. The single and two-traits mixed animal repeatability models, considered the fixed effects of herd-year-season and parity subclasses, breed as a function of cow H, F, J, and S breed fractions and general heterosis as a function of heterozygosity, and the random additive animal, permanent environment, and residual effects. For the analysis of LY, LL was added as a fixed covariate to all models. Variance components and genetic parameters were estimated using average information restricted maximum likelihood procedures. The results indicated that all traits were affected (p<0.001) by the considered fixed effects. High grade B×F cows (3/16B 13/16F) had the highest least squares means (LSM) for LY (2,490±178.9 kg), IY (10.5±0.8 kg), PY (12.7±0.9 kg), YD (7.6±0.55 kg) and LL (361.4±31.2 d), while B cows had the lowest LSM values for these traits. The LSM of LY, IY, YD, and PY tended to increase from the first to the fifth parity. Single-trait analyses yielded low heritability (0.03±0.03 and 0.08±0.02) and repeatability (0.14±0.01 to 0.24±0.02) estimates for LL, DP and parameter c. Medium heritability (0.21±0.03 to 0.33±0.04) and repeatability (0.27±0.02 to 0.53±0.01) estimates were obtained for LY, IY, PY, YD and ln(a). Genetic correlations between LY, IY, PY, YD, ln(a), and LL ranged from 0.59 to 0.99. Spearman's rank correlations between sire estimated breeding values for LY, LL, IY, PY, YD, ln(a) and c were positive (0.67 to 0.99, p<0.001). These results suggested that selection for IY, PY, YD, or LY would genetically improve lactation milk yield in this Ethiopian dairy cattle population.

INTRODUCTION

The majority of Ethiopian cattle (e.g., Boran and Horro) are low producers of milk with high individual variability (Yohannes et al., 2002) which is attributed to natural selection to the tropical environment and management. Milk yield could be increased through within and between breed selection and crossbreeding with cattle from improved dairy breeds. The strategy currently employed in Ethiopia to increase milk yield is by crossbreeding locally adapted Bos indicus breeds capable of withstanding heat stress, diseases, low feed supply and low management level with Bos taurus dairy breeds of high genetic ability for milk yield. This mating scheme has resulted in a variety of crossbred groups with different levels of additive and non-additive genetic abilities.

Variance components and genetic parameters are needed for genetic improvement programs to predict the breeding values of candidates for genetic selection, to choose among mating plans and to predict selection response (Montaldo et al., 2012). Knowledge of the type and amount of genetic variation and distribution of animals for traits considered for selection in the population can help design optimum breeding programs (Willham and Pollak, 1985). Additive genetic variation has been instrumental for the genetic improvement of traits of economic importance in dairy cattle populations (Willham and Pollak, 1985).

Numerous estimates of genetic parameters for production traits have been reported (Kahi et al., 2000; Demeke et al., 2004; Espinoza et al., 2007). These parameter estimates are specific to particular populations. Estimates of genetic parameters could change due to selection, migration of genes from one population to another, or changing environmental conditions (Van Der Werf and De Boer, 1989). Parameter estimates for one population or breed may not be appropriate for use in another population. Limited research has been done to estimate variance components and genetic parameters for lactation pattern and milk production traits and parameters of a lactation curve function in Ethiopian dairy cattle populations. Thus, the objective of this study was to estimate variance components and genetic parameters (i.e., heritability, repeatability, genetic and phenotypic correlations) for lactation pattern and milk production traits and parameters of a lactation curve function in a multibreed Ethiopian dairy cattle population.

MATERIALS AND METHODS

Description of the study area, animals and data

The study was based on dairy cattle milk data from Bako, Debre Zeit and Holetta Research Centers, Ethiopia. The Bako Agricultural Research Center is located 250 km West of Addis Ababa at an altitude of 1,650 m above sea level. The area receives a mean annual rainfall of 1,200 mm in a bimodal distribution, 80% of which falls from May to September. In this area, the average relative humidity is 59%, and the average minimum and maximum temperatures are 13.5°C and 27°C, respectively (Gebreyohannes et al., 2013).

The Debre Zeit Research Station of the International Livestock Research Institute (ILRI) is located approximately 50 km south East of Addis Ababa, in the Ethiopian highlands (9°N and 39°E), at an altitude of approximately 1,850 metres above sea level. The average annual rainfall in the Debre Zeit area is approximately 866 mm. The annual average temperature is 18.7°C and the average monthly relative humidity is 52.4% (Haile et al., 2011).

The Holetta Agricultural Research Centre is located 45 km west of Addis Ababa at 38.5°E longitude and 9.8°N latitude, and elevation of 2,400 metres above sea level. It is situated in the central highlands of Ethiopia. Average annual rainfall is approximately 1,200 mm. The annual average temperature is 18°C and the average monthly relative humidity is 60% (Demeke et al., 2004; Haile et al., 2011).

The herds in these three centers were genetically linked by the use of common sires supplied from the National Artificial Insemination Centre (NAIC). All of them used at least one common indigenous breed as dams of the replacement cows. Mating took place throughout the year using both artificial insemination and natural service. The bulls for natural services were those selected from the available males in the herd by visual appraisal. The Boran (B) and Horro (H) cattle (Bos indicus) used as a foundation stock for crossbreeding were bought from farmers in the Southern and Western Ethiopia, respectively. The herd from Debre Zeit research center was B and their crosses with Friesian (F).

Animals in this population were purebred B and H (Alberro and Hailemariam, 1982), and their crosses with F, Jersey (J) and Simmental (S) recorded by the Bako and Holetta research centers, and B and F×B crosses in Debre Zeit research station. This mating strategy generated 23 straightbred and crossbred groups (Table 1).

Table 1

Number of observations in each breed group

Breed group number	Breed group1	Number of records
1	Boran (B)	222
2	1/4F 3/4B	38
3	1/2F 1/2B	1,857
4	9/16F 7/16B	15
5	5/8F 3/8B	183
6	3/4F 1/4B	375
7	13/16F 3/16B	10
8	7/8F 1/8B	43
9	1/2J 1/2B	642
10	9/16J 7/16B	16
11	5/8J 3/8B	37
12	3/4J 1/4B	115
13	1/2S 1/2B	310
14	Horro (H)	362
15	1/4F 3/4H	84
16	3/8F 5/8H	17
17	1/2F 1/2H	409
18	3/4F 1/4H	64
19	1/4J 3/4H	11
20	1/2J 1/2H	372
21	3/4J 1/4H	64
22	1/2S 1/2H	240
23	3/4S 1/4H	21
All		5,507

F = Friesian, J = Jersey, and S = Simmental.

Lactation curve function and parameters

The modified incomplete gamma function (MIG; Papajcsik and Bodero, 1988) was chosen to be used in this study for its goodness of fit to the monthly test-day milk data and predictive ability of lactation milk yields (Gebreyohannes et al., 2013). The incomplete gamma (IG) is represented by y = ate where t is the length of time since calving. The MIG assumes the parameter b of the IG function to be equal to one, i.e., y = ate. Here, the log-transformed linear form of MIG ((ln(y/t) = ln(a)+(−ct)) was used to fit to the monthly test-day milk data using the regression procedure (SAS, 2003). The y is milk yield (kg) at time t (d) after calving, and a, b and c are parameters of the functions. The parameter a is a scaling factor associated with the average yield, b is related to pre-peak curvature and c is related to post-peak curvature. Parameters of the MIG function were obtained as an output from the regression analysis fitting the MIG function to each lactation of each cow. Cows with lactation lengths shorter than 90 d (three monthly test-day records) were excluded from the analysis and longer lactations were truncated to 305-d lactation length. A total of 59,413 monthly test-day milk records were used to estimate the parameters of the MIG function.

Data analysis

The data covered a period from 1977 to 2010 for Bako and Holetta, and 1989 to 2006 for Debre Zeit research center. The analysis used 5,507 lactation milk records of 1,639 cows from 283 sires. The traits considered were lactation milk yield (LY, kg), initial milk yield (IY, kg), peak milk yield (PY, kg), average daily milk yield (YD, milk yield per day of lactation, kg), days to peak (DP, d), lactation length (LL, d), and parameters of the MIG function (i.e., ln(a) and c). Phenotypic means by breed group and parity were compared using least squares means (LSM). The model used to compute LSM included the effects of herd-year-season of calving, parity and breed group subclasses for all traits (LY, IY, PY, YD, DP, LL, ln(a) and c). In addition, the model for LY included LL as a covariate to account for LL differences among cows. Computations were done with the GLM procedure of the Statistical Analysis System (SAS, 2003). The model for IY, PY, YD, DP, LL, ln(a) and c was:

Where

Y = Observation on lth cow that calved in ith herd-year-season of calving, jth parity, kth breed group subclasses,

μ = Overall mean,

HYS = ith calving herd-year-season subclasses (i = 1 to 325),

P = jth parity (j = 1, 2, 3, 4, 5, 6 and ≥7),

BG= kth breed group (k = 1 to 23; Table 1), and

e = residual associated with y.

In addition to all the effects in the model above, the model for LY included the covariate term b LL, where b = regression coefficient of LY on LL, and LL = lactation length of the lth cow that calved in ith herd-yr-season of calving, jth parity, and from kth breed group. The LSM for breed groups and parities were compared using Bonferroni t-tests (SAS, 2003).

Estimation of variance components and genetic parameters

Variance components were estimated using an Average Information Restricted Maximum Likelihood (AI-REML) procedure with the ASREML software (Gilmour et al., 2009). Single and two-trait animal repeatability models composed of fixed (i.e., herd-year-season, parity and regressions on cow H, F, J, and S breed fraction and general heterosis) and random (i.e., animal additive genetic, permanent environment and residual) effects were used to estimate variance components and to predict breeding values. Lactation length was considered as a covariate only for LY analysis. The general heterosis of cow involved interbreed interactions between alleles of any two different breeds among the 5 breeds (H, B, F, J, and S) present in the population. The H, F, J, and S breed effects were computed as a deviation from the Boran. Eight single-trait (LY, IY, PY, YD, DP, LL, ln(a) and c) analyses were performed to estimate variances, heritability and repeatability, and predict breeding values. The single trait animal repeatability model in matrix notation could be described as: where:

y = the vector of observations (i.e., LY, IY, PY, YD, DP, LL, ln(a) and c),

b = vector of fixed effects of herd-year-season and parity subclasses, and LL covariate for LY only,

g = vector of fixed breed effects,

a = vector of random animal additive genetic effects,

pe = vector of random permanent environmental effects,

e = vector of random residual effects,

X = incidence matrices relating records to fixed herd-year-season and parity, and LL covariate for LY only,

Q = matrix relating observations to breed effect (through H, F, J, and S breed fractions) and general heterosis (through heterozygosity fraction),

Z = incidence matrices relating records to animal additive genetic effects, and

W = incidence matrices relating records to permanent environmental effects.

The model assumed an expected value of y to be Xb+ Qg. The vector of direct animal additive genetic effects was assumed to have a normal distribution with mean zero and variance , where A is the additive numerator relationship matrix among animals in the population, and is the additive genetic variance. The vectors of permanent environmental and random residual effects were assumed to have independent normal distributions with mean zero and variances and where I is the identity matrix, and and are the permanent environment and residual variances, respectively. The phenotypic variance was . Heritabilities ( ) and repeatabilities (( )/ ) were calculated from the animal additive genetic, permanent environmental and residual variances, which were estimated in the single-trait analyses. Genetic and phenotypic correlations between traits (LY, IY, PY, YD, DP, LL, ln(a) and c) were computed using two-trait animal repeatability models. The two-trait animal repeatability model in matrix notation could be described as:

where y is vector of observation for the ith trait; b is vector of fixed effects (i.e., HYS, parity, and LL covariate for LY only), gi is vector of cow breed and general heterosis effects for the ith trait; a is vector of animal direct genetic effects for the ith trait; X, Z, and W are incidence matrices that relate fixed (herd-year-season and parity subclasses and covariate), animal direct genetic and permanent environmental effects to observations for the ith trait, respectively and Qi is matrix relating observations to breed effect (through H, F, J, and S breed fractions) and general heterosis (through heterozygosity fraction). The model assumed the expected value for trait i (y) to be X+Q. The vectors of direct animal additive genetic effects, permanent environment and residual effects for each trait were assumed to be normally distributed with mean zero. The variance-covariance structure of the random effects for a bivariate animal repeatability model could be described as:

where and are variances due to additive genetic effects for traits 1 and 2, and are variances due to permanent environmental effects for traits 1 and 2, and are residual variances for traits 1 and 2, σa1a2, σpe1pe2 and σe1e2 are additive, permanent environment and residual covariances between traits 1 and 2, A is the numerator relationship matrix, and I is an identity matrix.

Estimation of breeding values and sire rank correlations

Sire estimated breeding values (EBV) were computed from the Best Linear Unbiased Prediction (BLUP) solutions obtained from ASREML (Gilmour et al., 2009). The EBV for an animal was calculated as the sum of the generalized least squares (GLS) estimate of the difference between its breed group and B plus its predicted random additive genetic value. In matrix notation the EBV for each animal could be computed as:

Where û is the EBV of an individual for a particular trait, ĝ is a vector of GLS estimates of differences between breeds H, F, J, and S and breed B, p' is the transpose of the vector of fractions of H, F, J, and S breeds in an individual, and â is the predicted random animal additive genetic value (i.e., from p' ĝ ; Elzo and Famula, 1985; Arnold et al., 1992; Koonawootrittriron et al., 2002). The ranks of the sires EBV for each trait were tested using Spearman’s Rank correlations for significant association between the ranks of the sires for pairs of traits using the CORR procedure of the SAS (SAS, 2003). Genetic trends across years were evaluated using a regression analysis of mean yearly EBV for each trait calculated from solutions for each herd-year-season contemporary group and regressed on year. The year data collection started, i.e., 1977, was considered as the base year.

RESULTS AND DISCUSSION

Factors affecting the lactation pattern and milk production traits

Herd-year-season, parity and breed group subclasses were important factors (p<0.001) for all traits (LY, IY, PY, YD, DP, LL, ln(a) and c). Crossbred cows had higher least squares means (LSM) for LY, IY, PY, and YD, and longer LSM for LL than B and H breeds (Table 2). The 13/16F 3/16B crossbred cows had the highest LSM for LY (2,490±178.9 kg), IY (10.5±0.8 kg), PY (12.7±0.9 kg) and YD (7.6±0.55 kg) and the longest LL (361.4±31.2 days). Upgrading crossbred F×B or F×H from 50% to 87.5% F showed no significant increase in LSM for LY, IY, PY, and YD. However, the H cows yielded higher LSM for LY (1,210±37.9 kg vs 947±42.3 kg), IY (3.9±0.2 kg vs 1.0±0.2 kg), PY (5.5±0.2 kg vs 3.5±0.2 kg) and YD (2.8±0.12 kg vs 1.8±0.13 kg) and milked longer (234.9±6.5 kg vs 211.1±7.2 kg) than B cows. These differences among breed groups agreed with previous reports (Kahi et al., 2000; Demeke et al., 2004; Wolf et al., 2005) suggesting that biological differences among indigenous (B and H) and exotic parental populations (F, J, and S) may help explain differences in additive and non-additive genetic effects among breed groups in this study.

Table 2

Least squares means (±standard errors) for lactation pattern and milk production traits and lactation curve parameters

Breed group2	LY (kg)	IY (kg)	PY (kg)	YD (kg)	DP (d)	LL (d)	ln(a)	c
Pure breed
Boran (B)	947±42.3e	1.0±0.2d	3.5±0.2e	1.8±0.13e	53.9±1.9a	211.1±7.2c	−2.25±0.03d	0.016±0.0003c
Horro (H)	1,210±37.9de	3.9±0.2c	5.5±0.2d	2.8±0.12d	29.1±1.7b	234.9±6.5bc	−1.95±0.03cd	0.017±0.0003c
Two breed cross
1/4 F 3/4B	1,476±95.2cd	5.9±0.4b	8.2±0.5bc	4.3±0.29bc	33.1±4.4b	279.0±16.6abc	−1.51±0.07ab	0.015±0.0007bc
1/2F 1/2B	2,031±20.9a	8.4±0.1a	11.0±0.1a	6.4±0.06a	32.6±1.0b	337.2±3.6a	−1.32±0.02a	0.012±0.0001a
9/16F 7/16B	1,682±147.4bcd	7.1±0.7b	9.8±0.7ab	5.6±0.46ab	26.0±6.8b	323.3±25.7a	−1.33±0.11a	0.013±0.0011abc
5/8F 3/8B	2,134±52.1a	8.8±0.2a	11.4±0.2a	6.6±0.16a	31.5±2.4b	329.1±9.1a	−1.28±0.04a	0.012±0.0004ab
3/4F 1/4B	2,240±35.9a	8.9±0.2a	11.6±0.2a	7.0±0.11a	32.3±1.6b	343.2±6.3a	−1.25±0.03a	0.011±0.0003a
13/16F 3/16B	2,490±178.9a	10.5±0.8a	12.7±0.9a	7.6±0.55a	34.8±8.2b	361.4±31.2a	−1.10±0.13a	0.012±0.0013ab
7/8F 1/8B	2,410±92.4a	9.8±0.4a	12.1±0.4a	7.6±0.29a	24.8±4.2b	299.0±16.1ab	−1.20±0.07a	0.012±0.0007ab
1/2J 1/2B	1,788±26.5bc	7.0±0.1b	9.4±0.1b	5.6±0.08ab	33.3±1.2b	315.3±.6a	−1.51±0.02ab	0.011±0.0002a
9/16J 7/16B	1,496±141.7cd	5.4±0.7bc	7.8±0.7bcd	4.7±0.44bc	46.8±6.5ab	312.4±24.7a	−1.69±0.10abc	0.011±0.0010a
5/8J 3/8B	1,707±96.7bc	6.8±0.5b	9.0±0.5bc	5.4±0.30ab	26.3±4.4b	320.6±16.9a	−1.54±0.07ab	0.012±0.0007ab
3/4J 1/4B	1,832±56.0abc	6.8±0.3b	9.2±0.3b	5.7±0.17ab	37.5±2.6b	302.8±9.8a	−1.57±0.04ab	0.011±0.0004a
1/2S 1/2B	1,900±35.0ab	7.3±0.2ab	10.4±0.2a	6.0±0.11a	38.6±1.6b	339.1±6.1a	−1.48±0.03a	0.011±0.0003a
1/4F 3/4H	1,346±68.2 cde	5.2±0.3bc	7.0±0.3cd	3.7±0.21cd	30.4±3.1b	271.0±11.9abc	−1.70±0.05bc	0.015±0.0005bc
3/8F 5/8H	1,809±152.9abc	7.9±0.7ab	10.2±0.7ab	5.8±0.47ab	27.0±7.0b	324.5±26.7a	−1.33±0.11a	0.013±0.0011abc
1/2F 1/2H	1,836±31.6abc	7.1±0.1b	9.8±0.2ab	5.7±0.10ab	36.1±1.5b	321.0±5.5a	−1.46±0.02a	0.012±0.0002ab
3/4F 1/4H	2,184±72.8a	8.7±0.3a	11.5±0.3a	6.8±0.23a	38.2±3.3b	360.7±12.7a	−1.31±0.05a	0.011±0.0005a
1/4J 3/4H	1,404±174.3cde	5.4±0.8bc	7.8±0.8bcd	4.1±0.54bcd	34.6±8.0b	285.0±30.4abc	−1.78±0.13bcd	0.012±0.0012ab
1/2J 1/2H	1,621±33.1cd	6.5±0.2b	8.5±0.2bc	4.9±0.10b	32.5±1.5b	303.8±5.8a	−1.60±0.02ab	0.012±0.0002ab
3/4J 1/4H	1,724±73.9bc	7.4±0.3ab	9.7±0.4ab	5.5±0.23ab	29.4±3.4b	329.0±12.9a	−1.53±0.05ab	0.011±0.0005a
1/2S 1/2H	1,662±40.0cd	6.5±0.2b	8.8±0.2bc	5.2±0.12b	34.3±1.8b	299.2±7.0a	−1.57±0.03ab	0.012±0.0003ab
3/4S 1/4H	1,876±125.9abc	6.6±0.6b	10.1±0.6ab	5.9±0.39ab	47.0±5.8ab	309.3±22.0a	−1.42±0.09a	0.012±0.0009ab

LY = Lactation milk yield; IY = Initial milk yield; PY = Peak milk yield; YD = Average yield per d; DP = Days to peak; LL = Lactation length; ln(a) and c = Parameters of the MIG function.

F = Friesian, J = Jersey and S = Simmental.

Least squares means within a column group with different letters differ significantly (p<0.001).

Upgrading from 50% to higher F fractions for both F×B and F×H crosses showed no significant differences for milk yields (Table 2). This outcome was different from Koonawootrittriron et al. (2012), who found that the LSM for milk yield tended to increase as Holstein fraction increased from less than 0.625 (4,075±109 kg) to 0.875 (4,285±47 kg) and then it decreased towards purebred Holstein (4,120±76 kg). This difference may be attributed to limitations in the provision of the necessary management that high percentage F cows require, and differences in adaptation to the tropical environment and feed availability. In addition, the preferential treatment given to improved breed groups relative to low milk yield indigenous cows may also have contributed to the higher performance observed in crossbred cows than that of B and H cows. The milk yield and lactation length for B and H cows reported in this study were higher than values reported for the same breeds by Demeke et al. (2004). This could be due to inclusion of cows with lactation lengths above 90 d in the present study compared to the 60 d considered in the previous study (Demeke et al., 2004). Among the three exotic sire breeds (F, J, and S), the F were superior followed by S and J for the considered traits.

First parity cows had longer DP (39.6±1.1 d) and LL (322.6±4.2 d) compared to the second and later (≥3) parity cows (p<0.001; Table 3) while LSM for LY, IY, PY, and YD were lower for the first parity cows followed by cows in the second parity compared to cows in later parity (≥3). Other authors have observed increased daily milk yields with parity number (Epaphras et al., 2004). The difference among parities could be associated with the maturity of the cows, high feed intake of mature cows than heifers, differences in partitioning of more feed to maintenance and growth at the expense of milk production in heifers than mature cows and early culling of cow with lower milk yield (Gradiz et al., 2009). Longer DP and LL observed in first parity cows compared to later parities could be associated to the better persistency of first lactation cows (Stanton et al., 1992).

Table 3

Least squares means (±standard errors) for lactation pattern and milk production traits and lactation curve parameters

Parity	N	LY (kg)	IY (kg)	PY (kg)	YD (kg)	DP (d)	LL (d)	ln(a)	c
1	1,361	1,443±24.2c	4.9±0.1d	7.5±0.1d	4.4±0.07d	39.6±1.1a	322.6±4.2a	−1.79±0.02d	0.012±0.0002a
2	1,138	1,680±25.1b	6.4±0.1c	8.9±0.1c	5.1±0.08c	35.4±1.2b	317.2±4.4ab	−1.59±0.02c	0.012±0.0002ab
3	917	1,805±26.9a	7.0±0.1b	9.5±0.1b	5.5±0.08b	35.5±1.2b	311.2±4.7ab	−1.49±0.02b	0.013±0.0002bc
4	739	1,878±28.7a	7.5±0.1a	9.9±0.1ab	5.7±0.09ab	33.6±1.3bc	314.6±5.0ab	−1.45±0.02ab	0.013±0.0002bc
5	540	1,892±31.8a	7.7±0.1a	10.1±0.2a	5.8±0.10a	31.0±1.5c	305.9±5.5bc	−1.40±0.02a	0.013±0.0002cd
6	380	1,837±35.8a	7.3±0.2ab	9.7±0.2ab	5.6±0.11ab	34.9±1.6bc	302.3±6.2bc	−1.45±0.03ab	0.013±0.0003cd
≥7	432	1,884±35.6a	7.6±0.2a	9.8±0.2ab	5.8±0.11ab	30.4±1.6c	290.8±6.2c	−1.38±0.03a	0.014±0.0003d

LY = Lactation milk yield; IY = Initial milk yield; PY = Peak milk yield; YD = Average yield per d; DP = Days to peak; LL = Lactation length; ln(a) and c = parameters of the MIG function.

Least squares means within a column group with different letters differ significantly (p<0.001).

Covariance component and genetic parameter estimates

Variance component estimates (Table 4) from single-trait analysis for LY, IY, PY, YD, DP, LL, ln(a), and c produced low heritability (0.03±0.03 to 0.08±0.02) estimates for LL, DP, and parameter c of MIG. Heritability estimates for LY, IY, PY, YD, LL, and ln(a) were medium (0.21±0.03 to 0.33±0.04). Repeatability estimates ranged from 0.21±0.02 to 0.53±0.01 for LY, IY, PY, YD, LL, ln(a), and c, and 0.14±0.01 for DP (Table 4).

Table 4

Variance component and genetic parameter estimates for lactation pattern and milk production traits and lactation curve parameters

Trait1	Variance components2				Genetic parameters3
	σa2	σpe2	σe2	σp2	h2	r
LY (kg)	84,166	71,259	139,090	294,514	0.29±0.04	0.53±0.01
IY (kg)	2.2	0.5	4.1	6.8	0.33±0.04	0.39±0.02
PY (kg)	2.1	1.2	3.4	6.8	0.32±0.04	0.49±0.02
YD (kg)	0.9	0.6	1.4	2.9	0.30±0.04	0.51±0.02
DP (d)	52	38	557	647	0.08±0.02	0.14±0.01
LL (d)	547	1,395	7,358	9,300	0.06±0.03	0.21±0.02
ln(a) (kg)	0.04	0.01	0.12	0.16	0.21±0.03	0.27±0.02
c (kg/d)	0.0000004	0.0000034	0.0000120	0.0000157	0.03±0.03	0.24±0.02

LY = Lactation milk yield; IY = Initial milk yield; PY = Peak milk yield; YD = Average yield per d; DP = Days to peak; LL = Lactation length; ln(a) and c = parameters of the MIG function.

= Additive genetic variance, = Residual variance, = Permanent environmental variance, = Phenotypic variance.

h2 = Heritability, r = Repeatability.

The heritability estimate for LY found here was comparable to values reported by Albuquerque et al. (1996) and Mashhadi et al. (2008) for Holstein cows, but higher than estimates (0.24±0.04) from Demeke et al. (2004) for B and F×B crosses, Yilmaz et al. (2011) for Brown Swiss cattle (0.25) and Ayied et al. (2011) for Holstein cattle (range 0.05 to 0.21). However, it was lower than values obtained in a Thai multibreed dairy population (0.43±0.24; Seangjun et al., 2009) and in a Dutch Black and White dairy cattle population (0.48; Hoekstra et al., 1994). Further, the heritabilities estimated in this study were higher for IY and PY, but lower for LL and DP than the values reported for a multibreed Thai dairy population (0.23±0.18 for IY; 0.20± 0.18 for PY and 0.52±0.28 for DP; Seangjun et al., 2009). Heritability estimates are specific to a particular population (Van Der Werf and De Boer, 1989). Heritability estimates vary with lactation (Albuquerque et al., 1996; Cilek and Sahin, 2009). Cilek and Sahin (2009) attributed the lower heritability for all lactation milk yields than the heritabilities of first and second lactations to culling and selection that may have reduced the phenotypic variance of older cows compared to first-lactation cows. Hoekstra et al. (1994), for instance, attributed the high heritability estimate in their study to the high milk production. Ayied et al. (2011) compared heritability estimates for milk yield and lactation length using four methods of estimation (MINQUE, ANOVA, ML, and REML) and reported heritability estimates ranging from 0.07 (REML) to 0.21 (MINQUE) for milk yield and 0.02 to 0.09 for lactation length. Espinoza et al. (2007) compared heritability estimated for the first, second, third and fourth parities lactation milk yield using univariate analysis with animal repeatability model and reported lower heritability values (ranged 0.10 to 0.15) from the univariate analysis than from the animal repeatability model (0.16).

The heritability estimate for YD (0.33±0.04, Table 4) was comparable to the heritability estimate of 0.28 reported by Lopez-Villalobos and Spelman (2010). Atili and Khattab (2005) reported lower heritability estimates for lactation milk yield and lactation length for Holstein Friesian cows using animal repeatability model compared to heritability estimated from sire model. In general, moderate heritability estimates indicate sufficient additive genetic variance that would allow selection within the population. Thus, considering heritability estimates, LY, IY, PY, YD, and ln(a) could be improved through selection of either one or a combination of these traits in this dairy cattle population. The low heritability and repeatability estimates for LL, DP, and c indicate that these traits were largely affected by environmental factors, perhaps largely related to the management and feeding of the animals. The medium repeatability obtained for LY, IY, PY, and YD (0.39 to 0.53) in this study suggested that information from the first parity could be used for early prediction of EBV and cow selection. These early EBV would later be updated as data from later lactations became available.

The genetic correlations between LY, IY, PY, YD, LL, and ln(a) were high (range 0.59±0.15 to 0.99±0.03) and positive (Table 5). Conversely, the genetic correlations coefficients between DP and LY (−0.27±0.15), IY (−0.57± 0.12), PY (−0.41±0.15), YD (−0.31±0.16), ln(a)(−0.42±0.16) were negative. Cows with higher IY and PY tended to reach their peak yield earlier than cows with lower IY and PY. The high and positive genetic correlations obtained here between LY, IY, PY, YD, ln(a), and LL indicated that cows with high genetic potential for milk production traits (LY, IY, PY, YD, and ln(a)) tended to show longer LL. Thus, selection for longer LL would tend to improve lactation milk yield in this Ethiopian population. To improve LY in this Ethiopian population, selection should consider sires and dams with higher IY, PY, and YD and longer LL. Similar correlations were also reported in different studies (Batra et al., 1987; Varona et al., 1998; Tekerli et al., 2000; Gradiz et al., 2009). Cows producing at higher levels in the beginning of lactation tended to have the highest peak yield and consequently the highest 305-d yield (Ali and Schaeffer, 1987; Tekerli et al., 2000).

Table 5

Genetic correlations (above diagonal) and phenotypic correlations (below diagonal) for lactation pattern and milk production traits

Traits	Traits1
	LY	IY	PY	YD	DP	LL	ln(a)	c
LY		0.91±0.03	0.96±0.02	0.98±0.01	−0.27±0.15	0.73±0.12	0.93±0.05	−0.02±0.29
IY	0.63±0.01		0.96±0.02	0.92±0.03	−0.57±0.12	0.59±0.15	0.90±0.05	−0.09±0.28
PY	0.78±0.01	0.80±0.01		0.97±0.02	−0.41±0.15	0.61±0.16	0.99±0.03	−0.04±0.28
YD	0.88±0.01	0.64±0.01	0.81±0.01		−0.31±0.16	0.67±0.16	0.96±0.04	0.13±0.25
DP	0.04±0.02	−0.30±0.01	−0.10±0.02	0.04±0.02		**	−0.42±0.16	**
LL	0.50±0.01	0.19±0.02	0.21±0.02	0.08±0.02	0.04±0.01		0.67±0.24	−0.09±0.77
ln(a)	0.57±0.01	0.53±0.01	0.65±0.01	0.63±0.01	0.09±0.01	−0.13±0.02		0.001±0.27
c	0.15±0.02	0.08±0.02	0.11±0.02	0.16±0.02	0.15±0.01	0.54±0.01	−0.43±0.01

LY = Lactation milk yield; IY = Initial milk yield; PY = Peak milk yield; YD = Average yield per d; DP = Days to peak; LL = Lactation length; ln(a) and c = parameters of the MIG function.

Genetic correlation between DP and LL and DP and c could not be estimated.

Estimated breeding values and genetic trends

Estimates for general heterosis were 489±44.9 kg for LY, 2.6±0.2 kg for IY, 3.0±0.2 kg for PY, 1.7±0.1 kg for YD, −4.7±1.7 d for DP, 50.9±6.8 d for LL, 0.33±0.03 kg for ln(a), and 0.002±0.0003 kg/d for c. Table 6 show estimates of breed differences between H, F, J, and S relative to Boran. These values were used to estimate breeding values. The estimates of breed effects for F, S, and J relative to B were large and influenced the ranking of sires across breed groups for all traits. Differences between F, S, and J breeds and B were positive for all traits, except for DP that were negative. Sire EBV ranged from −392.5 to 1,962 kg for LY, −2.9 to 10.9 kg for IY, −2.5 to 10.6 kg for PY, −1.8 to 7.4 kg for YD, −26.3 to 139.4 d for LL, −28.7 to 10.1 d for DP, −0.495 to 1.301 kg for ln(a), and −0.002 to 0.01 kg/d for c, respectively (Table 7). Comparison of EBV across sire breed groups showed that the highest EBV were observed for straightbred exotic sires (Friesian, Jersey and Simmental) for LY, IY, PY, YD, LL, and ln(a). Spearman’s rank correlations between sire EBV were significant for all pairs of traits (p<0.001) except for the correlation between DP and parameter c (Table 8). The rank correlation for pairs of traits involving LY, LL, IY, PY, YD, ln(a), and c, were positive and significant which indicated that sires with higher EBV for LY would likely also have higher EBV for the other positively correlated traits. Conversely, the rank correlation between DP with LY, IY, PY, and YD was negative which could be explained by the low negative genetic correlation (range −0.27±0.15 to −0.57±0.12) observed in this study (Table 5).

Table 6

Estimates (±standard errors) for breed differences from Boran for lactation pattern and milk production traits

Trait	Breed
	Horro	Friesian	Jersey	Simmental
LY (kg)	10.9±69.5	1,134±89.3	616±98.3	707±128.9
IY (kg)	0.6±0.3	6.3±0.4	3.9±0.5	3.8±0.6
PY (kg)	−0.1±0.3	6.5±0.4	3.5±0.5	4.3±0.6
YD (kg)	0.1±0.2	4.5±0.3	2.9±0.3	3.4±0.4
DP (d)	−10.0±2.5	−11.5±3.2	−9.8±3.5	−0.9±4.5
LL (d)	−13.7±9.8	115.7±12.6	71.4±13.9	91.4±17.8
ln(a) (kg)	0.01±0.05	0.81±0.06	0.42±0.07	0.48±0.08
c (kg/d)	−0.001±0.0004	0.01±0.001	0.01±0.001	0.01±0.001

LY = Lactation milk yield; IY = Initial milk yield; PY = Peak milk yield; YD = Average yield per d; DP = Days to peak; LL = Lactation length; ln(a) and c = Parameters of the MIG function.

Table 7

Number of sires (N) and ranges (minimum and maximum) of estimated breeding values of sires for lactation pattern and milk production traits across different sire breed groups

Sire breed group	N	Traits
		LY		IY		PY		YD		DP		LL		ln(a)		c
		Min	Max	Min	Max	Min	Max	Min	Max	Min	Max	Min	Max	Min	Max	Min	Max
B	21	−393	196	−2.9	0.9	−2.5	1.0	−1.8	0.6	−4.5	10.1	−17.4	16.2	−0.50	0.17	−0.0003	0.001
H	28	−43	277	0.0	2.9	−0.4	1.6	−0.3	1.1	−22.0	−3.5	−26.3	5.5	−0.10	0.25	−0.0017	−0.001
F	77	729	1,962	3.5	10.9	3.7	10.6	3.0	7.4	−28.7	0.2	94.5	139.4	0.52	1.30	0.0046	0.005
J	32	505	1,069	2.7	5.4	2.5	6.2	2.3	4.3	−15.1	−4.5	47.4	114.8	0.32	0.74	0.0047	0.005
S	37	503	1,287	2.8	5.1	3.2	6.5	2.8	4.9	−6.9	2.8	62.4	108.5	0.34	0.67	0.0056	0.006
1/2F 1/2B	31	150	709	1.3	4.5	1.5	4.1	0.9	2.9	−13.4	−0.7	36.9	77.2	0.26	0.53	0.0022	0.003
3/4F 1/4B	4	718	1,006	4.3	5.6	4.4	5.7	2.9	4.1	−11.2	−6.6	72.0	90.5	0.56	0.70	0.0036	0.004
1/2J 1/2B	17	144	333	0.6	2.6	1.0	2.5	1.1	1.7	−8.8	−2.3	26.1	61.7	0.07	0.34	0.0023	0.003
3/4J 1/4B	2	432	564	3.1	3.1	2.6	3.1	2.2	2.5	−8.6	−4.4	52.5	55.6	0.29	0.40	0.0037	0.004
1/2S 1/2B	3	236	385	1.7	3.3	1.6	3.4	1.3	2.3	−1.7	0.9	38.3	57.1	0.14	0.46	0.0029	0.003
1/2F 1/2H	11	213	694	2.4	4.1	1.2	4.1	0.9	2.9	−16.8	−8.1	38.4	61.3	0.14	0.52	0.0016	0.002
3/4F 1/4H	2	899	967	5.3	5.3	5.1	5.4	3.6	3.8	−12.4	−9.8	77.9	94.5	0.68	0.69	0.0034	0.004
1/2J 1/2H	11	113	567	1.3	3.2	0.3	2.7	0.7	2.0	−12.5	−4.9	18.8	41.6	−0.01	0.36	0.0019	0.002
3/4J 1/4H	1	496	496	3.2	3.2	2.7	2.7	2.5	2.5	−8.4	−8.4	59.8	59.8	0.32	0.32	0.0036	0.004
1/2S 1/2H	6	37	463	1.1	2.3	0.3	2.6	1.1	2.2	−13.0	−0.8	25.8	50.7	0.06	0.33	0.0022	0.003

LY = Lactation milk yield; IY = Initial milk yield; PY = Peak milk yield; YD = Average yield per d; DP = Days to peak; LL = Lactation length; ln(a) and c = Parameters of the MIG function.

Table 8

Spearman’s Rank Correlation coefficients of sires estimated breeding values for lactation pattern and milk production traits and lactation curve function parameters

Traits	Trait1
	IY	PY	YD	DP	LL	ln(a)	c
LY	0.97*	0.99*	0.99*	−0.33*	0.95*	0.96*	0.77*
IY		0.98*	0.96*	−0.44*	0.93*	0.94*	0.73*
PY			0.99*	−0.34*	0.95*	0.97*	0.76*
YD				−0.30*	0.96*	0.95*	0.81*
DP					−0.26*	−0.37*	0.07ns
LL						0.92*	0.82*
ln(a)							0.67*

LY = Lactation milk yield; IY = Initial milk yield; PY = Peak milk yield; YD = Average yield per d; DP = Days to peak; LL = Lactation length; ln(a) and c = parameters of the MIG function.

p<0.001, ns = Not significant.

The only trait with a significant regression of EBV on calving year was LL. The EBV for LL declined across calving years (b = −0.91; p = 0.02). In contrast, the EBV for LY (corrected for LL) across years was highly variable (Figure 1) perhaps preventing the upward trend of EBV yearly means becoming significant. The high degree of variation of EBV for LY across years may have been influenced by genetic variability among sires used across years and by management variability across years within and across the three locations involved in this study.

Figure 1

Yearly mean estimated breeding value (EBV) of sires for lactation milk yield (LY, kg).

IMPLICATIONS

Crossbreeding and upgrading could help increase milk production in this Ethiopian dairy population. The high additive genetic variance, heritability and repeatability estimates obtained here suggested the possibility of improving LY, IY, PY, and YD by genetic selection. High genetic and phenotypic correlation estimates among the traits (LY, IY, PY, and YD) also suggested that genetic improvement in one of these traits would result in improvement of other genetic correlated traits. The medium to high repeatability estimates suggested that first lactations could be used for early prediction of animal genetic abilities in this population. These early predictions will later be updated as information from later lactations becomes available.