Data from 1,567 animals of Holstein Friesian born to eight sire origins and 136 sires, 1,431 cows that calved between 1979 and 2013 were used. The data was edited in such a way that incomplete lactation record on lactation milk yield less than 1,000 kg and lactation lengths less than 220 days and greater than 450 days were excluded to avoid bias Ayied et al [3]. In addition, all lactation length records between 21 and 320 days were standardized to 305-days milk yield using the projection procedures of Rege [4]. Lactations greater than 321 days were adjusted by a regression model that considered coefficient of determination and variance inflation factor for each independent variable in the model as described by Khan et al [5]. All animals with AFC between 18 to 60 months were considered as normal expected age range and those animals with records greater than 60 months age for AFC were considered in the terminal category, similarly calving interval between 300 and 900 days and days open between 21 and 500 days were included in the final analysis [6]. Finally a total of 3,552, 3,733, 2,938, 1,125,2,764, and 2,773 records on LMY, 305 days milk yield (305-d MY), LL, AFC, CI, and DO, respectively were used in the final analysis of the mixed models fitted to generate the estimates (Table 1).

Seasons of calving was classified into three groups based on weather and climatic conditions of the area. These included June to September as long rainy season, March to May as short rainy season and October to February as dry season [7]. Further, lactations were classified into three parities as 1, 2, and 3⁺. Party three and above were all pooled together as parity three due to very few number of observations in later lactations as well as due to the higher correlation between 3rd and later lactations. Sires used in the farm were assigned into eight groups based on their source or country of origin. i.e. those imported from Cuba, Finland, Kenya, Israel, Italy, United States of America, and the remaining recruited in Holetta and unknown group from Ethiopia.

Estimation of variance and (co)variance components were made using the multivariate mixed model program package (DMU) using average information restricted maximum likelihood algorithm [8]. To estimate variance and (co)variance components and correlations, single and multi-trait animal models were used for the different traits depending on the nature of the data. A full multi-trait animal model for six traits was initially tried, but it had convergence problems. Other multiple-trait analyses were tried by reducing the number of traits step wisely one at a time from 6 but only the 2 trait and 3 trait multi-trait models resulted in logical estimates. Preliminary analysis was performed to identify the fixed effects which should be considered in the genetic analysis. The presence of any significant differences was checked by using Tukey Kramer multiple comparison tests. Fixed effects which were significant (p<0.05) were fitted in to the model to estimate the genetic parameters. Two different models, model 1 and 2 as shown below were used for estimation of variance components. In general, the fixed effects included in these models were season of calving, year of calving/birth, origin of sire and parity whilst permanent environment, animal genetic and residuals as random effects. The description of the models in matrix notation was:

Model 1 is used for univariate analysis for AFC

Where, Y is the vector of observation; b includes vector of fixed effects; a is the vector of solutions for the coefficients of direct animal (additive) genetic random effects; e is the vector of residual effects; and X and Z are incidence matrices of the fixed effects and additive genetic random effects, respectively.

Model 2 was used for univariate and multivariate (two trait and three trait) analysis of productive and reproductive traits (LMY, LL, 305-d MY, CI, and DO) fitting permanent environmental effect due to repeated records per cow.

Where, Y is the vector of observations; b includes vector of fixed effects; a is the vector of solutions for the coefficients of direct animal (additive) genetic random effects; pe is the vector of solution for permanent environmental effects; e is the vector of residual effects; and X, Z, W, are the correspondent incidence matrices of the fixed effects and additive genetic and permanent environmental random effects, respectively.

In the two trait analysis we considered 305-d MY as a major trait to evaluate the milk yield of dairy cows as in most modern dairy farms, a lactation length of 305 days commonly accepted as a standard. This standard allows for ideal calving every 12 months with a 60-day dry period. Therefore, we fitted lactation length and reproductive traits with 305-d MY and presented the heritability estimates of all traits in the two trait model. We expected a different h² estimate for 305-d MY. However, we found that all the estimates were similar and we just picked the estimates from one of the runs. Likewise, in the three trait model analysis we considered 305-d MY as a major trait and AFC as second trait based on its h² estimate (genetic variability) from the univariate analysis and economic value of the traits as indicated by Emana [9]. So we have made a combined analysis of 305-d MY and AFC with LL, CI, and DO and we got h² estimates for all traits just like two trait analysis. In both two and three trait we used LMY instead of 305-d MY merely to test its association with production and fertility traits.

Heritability was estimated as the ratio of the additive genetic variance to total phenotypic variance; and repeatability, as the ratio of the sum of the additive genetic variance and permanent environmental variance to phenotypic variance [10].

Where, h² is heritability; r is repeatability; σa2 is direct additive genetic variance; σpe2 is permanent environmental variance related to repeated records and σe2 is residual variance.

The genetic correlations were estimated as the ratio of additive covariance of the two traits (trait i and j) to the products of additive genetic standard deviations of the two traits whilst phenotypic correlations were estimated similarly as ratio of the sum of the genetic and environmental covariances (phenotypic covariances) to the products of phenotypic standard deviations of the two traits.

Where, r_g is genetic correlations; r_p is phenotypic correlation; σai2 is additive genetic variance for trait i; σaj2 is additive genetic variance for trait j; σpi2 is phenotypic variance for trait i, σpj2 is phenotypic variance for trait j; σ_aij is additive genetic covariance between trait i and j; and σ_pij is phenotypic covariance between trait i and j.

The heritability estimates for LMY in the present study from single and three trait models were similar (0.17±0.04) whilst in the two trait model was 0.25±0.02 (Table 2). These values of heritability for LMY indicated presence of modest genetic variability among the small population in Holetta. The increment in the estimates in two trait analyses for LMY and 305-d MY record of an animal under the study indicated the presence of high correlation among the traits that amplified the heritability estimates. The current estimate from the univariate model is comparable to those estimates reported for Holstein Friesian dairy cattle [11–13]. On the other hand, it is higher than those reported for first lactation Holstein Friesians in Egypt [14].

The heritability estimates for adjusted 305-d MY from univariate, 2 trait and 3 trait analysis were 0.15±0.04, 0.16±0.04, and 0.15±0.04, respectively (Table 2). These values of heritability for 305-d MY indicate presence of modest genetic variability among the small population in the study. However, unlike the estimates for LMY the heritability estimates in univraite, two and three traits models did not show much change, this could be due to adjustment made on lactation days that triggers weak correlation to other traits. The current estimates from all models were comparable with 0.11 for Egyptian Holstein Friesian dairy cattle [14]. Similarly, closely heritability estimate of 0.19 reported for the same breed in a temperate environment [15]. However, higher estimates were reported by several other authors from different production environments [4,16,17]. Such differences in the heritability estimates are expected and are mainly results of size of dataset, genetic variation within population, management and environmental conditions and the methods used for parameter estimation. In addition, we believed that a long term biased phenotypic selection and involuntary culling of animals from the existing herd might have contributed to the decreased hereditary differences and genetic variability within the population. Although estimate of heritability for LMY and adjusted 305-d MY indicates that still there is an opportunity for improving these traits through selection, future improvement of the herd could be meaningful if the current phenotypic selection is changed to selection based on genetic merit combined with improved herd management as well as introducing of new animals from other source to increase genetic variability within the herd.

The heritability estimates for lactation length ranged from 0.03±0.03 to 0.08±0.03 (Table 2) and these estimates correspond to those reported by Ojango and Pollott [17]. It is slightly higher than heritability estimates of 0.003±0.08 in Sudan using paternal half-sib correlation [18]. These low heritability estimates obtained in this study could be attributable to the high residual variance. In the study the heritability estimates of LL shows slight improvement from univariate model to 3 traits model. This enhanced accuracy is due to the ability of the model to account for relationship between the traits and better connections in the data due to residual covariance between traits [19]. Moreover, multivariate models gave better reliabilities than univariate models in the prediction of genetic merit which is mainly a result of the model’s ability to use extra information from correlated traits [20].

Heritability estimates for AFC were 0.47, 0.38, and 0.38 from univariate, 2 trait and 3 trait models, respectively (Table 2). These high heritability estimates for AFC from this study may be explained by weak selection made on this trait in contrast to the intense selection practiced for milk production and presence of genetic variability in the small herd as the result of the different sources of sires being used in the study period. The current estimates were comparable with estimates of Rege [4] and Ojango and Pollott [17] who reported 0.47 and 0.38, respectively for the same breed in Kenya. Differences between the estimates of heritability obtained in the study and estimates from other countries are most likely caused by management and climate differences affecting genetic and environmental variances and methods of estimation [21]. Our study suggests that improvement of this trait through selection in line with several other authors who have come to the same conclusion [4,22].

Heritability estimates for calving interval were 0.11±0.04, 0.11±0.04, and 0.13±0.04 from univariate, 2 trait and 3 trait models, respectively (Table 2). These estimates for CI is still within the literature estimates for the traits, indicating the existence of certain level of genetic diversity. The present result is lower than estimates of 0.17 for Iraqi Friesian cows [13]. However, it is higher than the estimates reported for the same breed in other tropical countries [17,23,24]. This relatively low heritability estimates for CI could be explained by large environmental variance. Therefore, improvements in nutrition and reproductive management should lead to a considerable decrease in length of CI than making purely genetic selection alone. However, in the study, multi-trait analysis with milk production trait and AFC indicates slight improvement in the estimates of heritability of calving interval. Incorporating information on correlated traits may enhance the accuracy of selection and reduce selection bias made by univariate analysis [22]. Therefore, a multiple trait model selection including milk production with fertility traits could improve the selection accuracy for selecting the right animals with good fertility traits [15,25].

Heritability estimates for days open in the same study were 0.09±0.03, 0.09±0.03, and 0.10±0.04 from univariate, 2 trait and 3 trait models, respectively (Table 2). The heritability estimates from all models were similar and coincide with most literatures. The heritability estimates for DO is also indicative of the presence of genetic basis that deserve selection or culling of an animal. The current estimates were close to 0.082 reported by Zambrano and Echeverri [24], but higher than 0.042 that was reported in temperate zone [15]. On the contrary, higher estimates of days open have been reported in Egypt and Iraq for the same breed [26,13].

In the current study the repeatability estimate for LMY and adjusted 305-d MY is moderate (Table 2). For LMY the repeatability estimates were 0.39 and 0.45 form univariate and 2 trait models, respectively, which are comparable with that reported by Ojango and Pollott [17], but higher than the reported value of 0.26±0.08 for the same breed in Pakistan [12]. As for adjusted 305-d MY the repeatability estimate were 0.42 and 0.41, these values were in close agreement with 0.43 reported by Dematawewa and Berger [15]. Nevertheless, Rege [4] and Kadarmideen et al [16] have reported higher repeatability estimates of 0.55 and 0.58, respectively. These moderate repeatability estimates in the present study indicated that selection of dairy cattle for milk yield based on early performance records is reliable.

Regarding lactation length the repeatability estimate was very low which correspond with Ojango and Pollott [17]. However, slightly higher estimate of 0.19 was reported for the same breed [12]. The low repeatability estimate implied that variation in lactation length was more a result of variation in management and feeding in the given environment of a particular lactation rather than of factors associated genetic makeup of the cow [17,12].

The repeatability estimates for calving interval and days open from all models were very low (Table 2). Estimates for CI and DO were comparable with Rege [4] who reported 0.25 for the same breed in Kenya, but higher than estimates reported by Zambrano and Echeverri [24]. The low repeatability estimates obtained in this study suggest that CI and DO are strongly influenced by temporary environmental factors, this is mainly because of the complex nature of reproductive traits, difficulties in detection of oestrus, and various other managerial and nutritional factors.

Genetic correlations between milk production traits were high. The highest genetic correlation was 0.99±0.01 between LMY and 305-d MY and the lowest was 0.73±0.19 between 305-d MY and LL (Table 3). This is comparable with Ahmad et al [27]. The highest genetic correlation between milk productions traits were result from pleiotropy [10]. This means that genetic improvement of one milk production trait could result in a correlated response in the correlated trait. Generally the very high correlation between these two traits is indicative of the fact that using 305-d MY can be sufficient to meet selection for lactation milk yield.

The genetic correlation between 305-d MY and AFC was −0.24±0.11 (Table 3). This is in close agreement with −0.22±0.08 reported in Ethiopia [28]. On the contrary, a relatively higher and positive (0.54) genetic correlation was reported for Holstein Friesian dairy cattle in Kenya [17]. The negative genetic correlation between 305-d MY and AFC indicates that genes which affect the 305-d MY trait positively will also be responsible for an early onset of puberty among the current population found at Holetta bull dam herd.

Similarly genetic correlation between 305-d MY and CI found to be −0.10±0.15. This result is accords with −0.03 correlation estimate of Friesian cattle in Kenya [4]. However, higher and positive genetic correlations of 0.54±0.01 was reported for the same breed in Ethiopia [28]. Regarding genetic correlation between 305-d MY and DO, it was close to zero (−0.02±0.16). Positive and lower estimates were reported in Turkey [29] and Egypt [14]. On the other hand, higher and positive estimates of genetic correlation was observed for Ethiopian Holstein [28]. The negative genetic correlation of 305-d MY and fertility traits indicates that improvement of fertility is possible to certain degree in this herd as milk yield increase. However, since the standard error of the estimates were high, the results needs to be interpreted with caution. The higher standard errors obtained from pair wise combinations of CI and DO with 305-d MY could be due to the smaller proportion of CI and DO records in the dataset.

In the current study, positive genetic correlations between various reproductive traits were observed. This positive genetic correlation between reproductive traits indicates that improvement of one trait has positive impact on the other trait. The genetic correlation between CI and DO was 0.99±0.01 (Table 3). This is in close agreement to a nearly perfect genetic correlation of 0.99 [30] and 0.98 [24]. The positive and strong genetic correlation in the present study indicate the evidence of common genetic and physiological mechanisms controlling those traits. In the tropics in general and in the Holetta bull dam herd in particular the recording system is still at its early stage of development. As a result, recording of calving interval is relatively easier than days open. Therefore, evaluation based on calving interval would be easy to handle by the available human resource and recording infrastructure than evaluation based on days open.

In this study phenotypic correlations between milk productions traits were high (Table 3). The highest phenotypic correlation was 0.84±0.01 between 305-d MY and LMY and the lowest was 0.05±0.02 between 305-d MY and LL. This explains that these traits are influenced by more or less same genetic and environmental factors. The highest phenotypic correlation among LMY and 305-d MY was in agreement with the estimate reported for crossbreeds dairy cattle in Egypt [26,27].

Phenotypic correlations between 305-d MY with fertility traits (AFC, CI, and DO) were in general negative (Table 3). The phenotypic correlation between 305-d MY and AFC in present study was higher than Ahmad et al [27] for crossbreed cattle (−0.02) and −0.02 for first parity Holstein Friesian [17]. Conversely, positive phenotypic correlation was reportedfor Jersey cattle in Ethiopia [28]. In the current study the phenotypic correlation between CI and DO is 0.98±0.01 (Table 3). This strong and positive phenotypic correlations is in agreement with results of Kadarmideen et al [25], Ghiasi et al [30], and Zambrano and Echeverri [24] who reported 0.95, 0.95, and 1.0 for Holstein Friesian cattle in different production environment, respectively.