1st Page

Journal of Animal Science - Animal Genetics

Genetic principal components for reproductive and productive traits in dual-purpose buffaloes in Colombia1


This article in JAS

  1. Vol. 93 No. 8, p. 3801-3809
    Received: Jan 20, 2015
    Accepted: June 17, 2015
    Published: August 3, 2015

    2 Corresponding author(s):

  1. D. A. Agudelo-Gómez 2*†,
  2. R. Pelicioni Savegnago,
  3. M. E. Buzanskas,
  4. A. S. Ferraudo,
  5. D. Prado Munari and
  6. M. F. Cerón-Muñoz
  1. * Facultad de Ciencias Administrativas y Agropecuarias, Corporación Universitaria Lasallista, 055440 Caldas, Colombia
     Facultad de Ciencias Agrarias, Grupo de Investigación en Genética, Mejoramiento y Modelación Animal (GaMMA), Universidad de Antioquia, 050036 Medellín, Colombia
     Departamento de Ciências Exatas, Faculdade de Ciências Agrárias e Veterinárias, UNESP, 14884-900 Jaboticabal, Brazil


A multitrait model (MC) and 5 reduced-rank models with principal component structure (components PC1, PC1:2, PC1:3, PC1:4, and PC1:5) were compared. The objectives were to determine the most appropriate model for estimating genetic parameters and to evaluate the genetic progress of dual-purpose buffaloes in Colombia using that model. The traits evaluated were weaning weight (WW), yearling weight (W12), weight at 18 mo of age (W18), weight at 2 yr of age (W24), age at first calving (AFC), and milk yield at 270 d of first lactation (MY270). Genealogy and productive information from 34,326 buffaloes born in Colombia between 1997 and 2014 were used. Colombian Association of Buffalo Breeders (ACB) provided the data. Direct additive genetic and residual random effects were included for all the traits. In addition, the maternal additive genetic effect and permanent environmental random effect were included for WW, while a maternal additive genetic effect was included for W12. The fixed effects were contemporary group (farm, year, and calving season: January to April, May to August, or September to December; for all traits) and sex (for WW, W12, W18, and W24). Additionally, parity was included as a fixed effect for WW and W12. Age at weighing was used as a covariate for WW, W12, W18, and W24. Genetic progress of all traits was analyzed using a generalized smooth model (GAM). According to the Akaike information criteria (AIC), the best model was the one with reduced rank and first 3 principal components (PC1:3). This model maintained 100% of the original variance. Genetic parameters estimated with this model were similar to those estimated by MC, but with smaller standard errors. Heritability for weight-related traits ranged between 0.23 and 0.44. Heritabilities for AFC and MY270 were 0.14 and 0.24, respectively. The genetic correlations obtained between all weights (WW, W12, W18, and W24) were positive and high. Correlations between all weights with AFC were negative and moderate. Correlations between all weights with MY270 were positive and moderate, and between MY270 with AFC were negative and low.


Buffalo herds in Colombia are managed as dual-purpose systems (Bolívar et al., 2012). It is important to know the genetic parameters and genetic correlations between traits associated with breeding and production of milk and meat because of their economic relevance. Traits associated with live weights can be measured several times over the life of the animals, and the associated milk yield can be assessed throughout several calving events. Multitrait models (Arango et al. 2002; Boligon et al., 2009) and repeatability models are suitable for studying these traits (Boligon et al., 2009; Nephawe, 2004).

In multitrait models, the (co)variances and genetic correlations between measurements taken at different ages might vary, and overparameterization of the model occurs when the number of traits is very large (Nobre et al., 2003). In that context, the principal component analysis (PCA) is a statistical approach that allows decreasing the rank of the covariance matrix (Kirkpatrick and Meyer, 2004), reducing the dimensionality of a set of correlated traits, and generating a new reduced set of variables called principal components (PC). PC are orthogonal to each other and are linear combinations of the original traits, preserving the original maximum variability (Johnson and Wichern, 1992).

The PCA has been used to evaluate these models in genetic evaluation studies for beef cattle (Boligon et al., 2013) to estimate genetic parameters of BW and reproduction traits in female Canchim cattle and to explore the relationships among estimated breeding values (Buzanskas et al., 2013). It has also been used in Angus cattle to estimate genetic parameters (Meyer, 2005), as well as in poultry to estimate genetic parameters for traits associated with egg production and their relationships (Savegnago et al., 2011). It is also used in dairy cattle to estimate genetic parameters (Bignardi et al., 2012). In the aforementioned studies, PCA has allowed reduction of dimensionality, offering an appropriate estimate of genetic parameters.

The aims of this study were to select from a multitrait and 5 reduced-rank models the one that offers the best estimation of genetic parameters and to evaluate the genetic progress for the analyzed traits using that model.



The Colombian Association of Buffalo Breeders (ACB) provided the data used in this study. Information of 12 dual-purpose herds with once-a-day milking and weaning between 6 and 10 mo of age was included. The herds are located in the Colombian Caribbean region in a rainforest life-zone (height above sea level: 80 m, temperature: 28°C, precipitation: 2000 mm/year; Holdridge, 1996). The first mating occurs after females reached 330 kg. Males were kept whole and culled for slaughtered after reaching 430 kg (Ramírez Toro et al., 2011). The animals were managed on a pasture and received mineral supplementation. The breeding system consisted of controlled natural mating, except for few cases in which artificial insemination was used. Genetic connectivity between years and herds was verified using the Genetic Drift Variance methodology proposed by Kennedy and Trus (1993), who suggested that a small value would indicate a high degree of connectedness. Average connectivity between years was 0.01, and the maximum value was 0.06. Average connectivity between herds was 0.16, and the maximum value was 0.75. Genetic connectivity between herds occurs because herds exchange breeding males. Animals were predominantly crossbreds with a tendency toward high proportion of Murrah breed.

Records were obtained between 1997 and 2014. The relationship matrix included 34,326 animals with 8,572 mothers and 684 fathers. The number of animals between 1997 and 2014 was 42, 161, 297, 725, 1,042, 1,490, 1,937, 2,041, 2,619, 1,989, 4,148, 4,558, 4,855, 5,417, 5,371, 4,194, 1,897, and 347, respectively.

The traits evaluated were weaning weight (WW) from 12,429 males and 11,518 females, yearling weight (W12) from 3,040 males and 4,154 females, weight at 18 mo of age (W18) from 1,305 males and 2,624 females, weight at 24 mo of age (W24) from 451 males and 2,282 females, age at first calving (AFC) from 4,138 records, and milk yield from first calving to 270 d in milk (MY270) from 2,229 females. The age ranges (in days) for each trait were WW (180 to 300), W12 (330 to 390), W18 (450 to 510), W24 (680 to 760), and AFC (760 to 1500). MY270 was estimated following the guidelines of the International Committee of Animal Recording (ICAR, 2012). Analyses performed by means of the generalized linear model helped to define the fixed effects that were considered in the contemporary groups. Observations with standardized residuals of 3.5 standard deviations above and below the mean were excluded. Additionally, contemporary groups with fewer than 3 observations were also excluded. The general structure of the data is presented in Table 1.

View Full Table | Close Full ViewTable 1.

Description of data for weaning weight (WW), yearling weight (W12), weight at 18 mo (W18), weight at 2 yr (W24), age at first calving (AFC), and milk yield at 270 d of first lactation (MY270)1

Trait SX No. Observations Mean No. Observations Total Mean CV Min Max NCG
WW (kg) M 12429 208.70 23937 207.35 0.23 80 384 216
H 11508 205.89
W12 (kg) M 3040 213.00 7194 210.75 0.21 100 490 192
H 4154 209.88
W18 (kg) M 1305 263.48 3929 256.09 0.20 130 500 155
H 2624 253.25
W24 (kg) M 451 385.34 2733 358.48 0.15 250 759 143
H 2282 353.65
AFC (days) H 4138 1110.35 0.11 761 1499 162
MY270 (kg) H 2229 993.69 0.23 401 2192 149
1SX = sex, CV = coefficient of variation, Min = minimum value, Max = maximum value, NCG = number of contemporary groups.


A multitrait model was fitted for the traits. Several fixed and random effects were fitted for the different traits. An overview is provided in Table 2. The contemporary group for WW, W12, W18, and W24 was defined as farm, year, and calving season (January to April, May to August, and September to December). The contemporary group for AFC and MY270 was farm, year, and first calving season (January to April, May to August and September to December).

View Full Table | Close Full ViewTable 2.

Fixed, random, and covariate effects for weaning weight (WW), yearling weight (W12), weight at 18 mo (W18), weight at 2 yr (W24), age at first calving (AFC), and milk yield at 270 d of first lactation (MY270) included in the models1

Random effects
Fixed effects
Covariable (age)
Trait a m pe ε SX NC CG
W12 X X X X X X X
W18 X X X X X X
W24 X X X X X X
MY270 X X X
1a = additive genetic effect, m = maternal genetic effect, pe = maternal permanent environmental effect,ε = residual effect, SX = sex (male, female), CN = number of calving (1 to 14), and CG = contemporary group.

A multitrait model (MC) and 5 reduced-rank models were analyzed, considering a structure of principal components from 1 to 5 components (PC1, PC1:2, PC1:3, PC1:4, and PC1:5).

The matrix representation of the MC model was

The model used for PC analysis was obtained by reparametrizing MC, generating an equivalent model that fits the principal components rather than conforming to the original traits (Meyer, 2005).

The matrix representation of the PC model iswhere y is vector of observations, β is vector of fixed effects, a is vector of direct additive genetic random effects, m is vector of maternal genetic effect, pe is vector of maternal permanent environmental effect, and ε is random residual vector. X, Z1, Z2, and W are incidence matrices relating the fixed effects, direct additive genetic effects, maternal genetic effects, and permanent environmental effects, respectively.

Where , and E = matrix of the eigenvectors ei, Λ = diagonal matrix of eigenvalues λi and = direct product. Then, the direct additive genetic (co)variance matrix () can be decomposed in terms of E and Λ, with , and λi and ei are ordered from highest to lowest. Considering only the largest mPC and replacing E by Em, the matrix comprises only the first m columns of E(e1, …, em), and Z1 has a number of columns proportional to the range of m to 6. The number of equations is correspondingly smaller (replacing Λ by submatrix Λm formed by the first m rows and columns), and a* contains m elements for each individual (Meyer and Kirkpatrick, 2005).

Wombat software (Meyer, 2007a) was used, applying the restricted maximum likelihood (REML) information criterion. The convergence criterion was set at 10−9.

The best model was selected by comparing the Akaike information criterion (AIC), which allows comparisons between nonhierarchical models, penalizing those with more parameters (Schwarz, 1978), and can be represented aswhere p is the number of model parameters and log L is the logarithm of restricted maximum likelihood. The best model is the one with the lowest AIC.

Genetic Progress

A generalized model with smoothed estimates (GAM) was used to describe the genetic progress of WW, W12, W18, W24, AFC, and MY270. The model included year of birth as a smoothed variable. Statistical analyses were performed separately for males and females. The “gam” library (Wood, 2011) of R Project software was used for models and graphic construction (R Core Team, 2014).

The model was as follows:where Yi is the ith estimated breeding value (EBV) record of the animal in the ith year of birth (yi), b0 is the intercept, s(yi) is the nonparameterized smooth function of the ith year of birth, and εi is the residual effect. GAM models allow data analysis without assuming a functional form known (Wood, 2006).


According to the AIC, the best model was the PC1:3 because of the following reasons: it kept all the original variance in the 3 eigenvalues, was more parsimonious because it estimated 40 parameters instead of the 46 estimated with MC (Table 3), achieved the convergence criterion faster-decreasing computational requirements, and adequately described the covariance structure between 6 traits, keeping 100% of the genetic variance of the traits. Several studies indicate that PCA is efficient enough to decrease the range of the model, maintaining a high percentage of the variation of the original traits in several production systems.

View Full Table | Close Full ViewTable 3.

Number of estimated parameters (NP), Akaike information criterion (AIC), and percentage of direct additive genetic variance explained by each of the eigenvalues in each model

Model1 NP AIC λ1 λ2 λ3 λ4 λ5 λ6
PC1 30 354077,96 100.00
PC1:2 36 354036,87 79.82 20.18
PC1:3 40 354033,18 74.77 22.31 2.92
PC1:4 43 354142,18 76.66 19.63 3.19 0.52
PC1:5 45 354045,87 77.60 18.30 3.33 0.16 0.16
MC 46 354143,17 77.56 17.59 3.82 0.65 0.37 0.01
1PCn = reduced-analysis model for n principal components, MC = multi-trait model.

Similar results to those found in this study were reported by Meyer (2005), who found that the first 3 or 4 PC summarize most of the genetic variation of 8 traits measured in Angus cattle. The author also concluded that 7 PC can explain most of the variations of the 14 traits associated with beef production (Meyer, 2007a,b). Similarly, Boligon et al. (2013) reported that 3 PC are sufficient to explain 100% of the variation in 8 traits related to body weight in Nellore cattle and found genetic correlations of 0.47 to 0.86, which allow for reducing the analysis range to a few PC, thus avoiding redundant information of highly correlated traits.

In Brazil, Savegnago et al. (2011) reported that 4 PC were sufficient to account for 80.04% of the variation of 13 traits associated with egg production of laying hens. Buzanskas et al. (2013) evaluated 4 reproductive traits in Canchim Cattle and reported that the first 2 PC explained 73.37% of the total variation.

A Holstein study in Brazil analyzed 10 traits associated with milk production. According to the authors, a 2-PC model kept 94.25% of the genetic variation of all the traits, having to estimate 74 instead of the 110 parameters required by the complete range model. The authors concluded that the low number of PC required to estimate most of the variance was due to the high genetic correlations between the analyzed traits (Bignardi et al., 2012).

The PCA technique was used to estimate genetic parameters of 9 traits (6 associated with milk production, 2 with reproduction, and 1 with health status) in buffaloes, concluding that a 3-PC model kept 94.25% of the total variation, without losing estimation efficiency (Oliveira et al., 2014).

Live weights of animals for WW, W12, W18, and W24 are presented in Table 1. Males were heavier than females (2.81, 4.35, 9.65, and 30.18 kg, respectively). The study by Bolívar-Vergara et al. (2012) in Colombia reported 182 ± 42.46, 201.8 ± 38.70, 50.89 ± 278, and 363.4 ± 54.32 live weights for the aforementioned traits, respectively. Their WW and W12 were lower, while W18 and W24 were higher than those observed in this work.

The MY270 estimated in this work is lower than that reported in Brazil by Tonhati et al. (2008) for Murrah buffaloes at 305 d of lactation (993.69 vs. 1,495.08 ± 617.08 kg, respectively; Table 1). It is also lower than that reported in Italy by Rosati and Van Vleck (2002) for buffaloes up to 270 d of lactation (2,286.8 ± 492.1 kg). In Colombia, Hurtado-Lugo et al. (2006) reported 1,064.59 kg yield at 270 d. These 3 mentioned studies measured yield of animals with several parities.

The estimated AFC in this work (1,110.35 d; Table 1) was higher than that reported for buffaloes in Brazil and Colombia, which were 1064.59 ± 171.83 d, n = 2,436; and 1,140 ± 283.5 d, n = 216 (Bolívar Vergara et al., 2010; Oliveira et al. (2014), respectively). This difference may be due to lower number of records in relation to the amount analyzed in the present study.

Heritabilities of the direct additive genetic component for W12, W18, W24, and MY270 estimated by model PC1:3 were similar to heritabilities estimated by the MC model. Heritabilities for WW and AFC were the same in both models (Table 4).

View Full Table | Close Full ViewTable 4.

Heritabilities (on the diagonal bold), genetic correlations (below the diagonal) and phenotypic correlations (above the diagonal) estimated for dual-purpose buffaloes in Colombia using a reduced model with three principal components (PC1:3) and a complete multi-trait model (MC; standard errors are in parentheses)

Trait1 WW W12 W18 W24 AFC MY270
PC1:3 model
    WW 0.24 (0.02) 0.78 (0.02) 0.72 (0.03) 0.39 (0.04) -0.19 (0.03) 0.09 (0.93)
    W12 0.99 (0.02) 0.26 (0.02) 0.81 (0.01) 0.49 (0.01) -0.23 (0.10) 0.05 (0.04)
    W18 0.94 (0.01) 0.95 (0.01) 0.44 (0.03) 0.66 (0.03) -0.33 (0.02) 0.15 (0.10)
    W24 0.80 (0.04) 0.79 (0.04) 0.89 (0.03) 0.30 (0.04) -0.32 (0.08) 0.29 (0.04)
    AFC -0.51 (0.09) -0.52 (0.10) -0.51 (0.09) -0.76 (0.08) 0.14 (0.03) 0.11 (0.03)
    MY270 0.86 (0.10) 0.88 (0.11) 0.76 (0.09) 0.41(0.12) -0.13 (0.08) 0.24 (0.05)
MC model
    WW 0.24 (0.01) 0.75 (0.00) 0.72 (0.00) 0.40 (0.01) -0.19 (0.01) 0.08 (0.03)
    W12 0.91 (0.02) 0.27 (0.02) 0.79 (0.00) 0.49 (0.02) -0.22 (0.02) 0.56 (0.04)
    W18 0.95 (0.02) 0.86 (0.03) 0.45 (0.03) 0.66 (0.01) -0.33 (0.02) 0.12 (0.04)
    W24 0.78 (0.06) 0.75 (0.07) 0.85 (0.05) 0.30 (0.04) -0.33 (0.02) 0.17 (0.05)
    AFC -0.57 (0.09) -0.44 (0.11) -0.53 (0.09) -0.76 (0.11) 0.14 (0.03) 0.12 (0.03)
    MY270 0.43 (0.22) 0.60 (0.12) 0.53 (0.10) 0.51 (0.12) -0.14 (0.9) 0.26 (0.06)
1WW = weaning weight, W12 = yearling weigh, W18 = weight at 18 mo of age, W24 = weight at 24 mo, AFC = age at first calving, MY270 = milk yield at 270 d of first lactation.
Bold values indicate heritabilities.

Heritabilities estimated by model PC1:3 were lower for WW, W12, and W24, and they were higher for W18 when compared with values reported for the same traits in Colombia, which were 0.33, 0.42, 0.42, and 0.41, respectively (Bolívar-Vergara et al., 2012).

The heritability estimate for AFC was 0.14. Other researchers reported that heritability was between 0.08 and 0.16 for the same trait in Nellore heifers (Dias et al., 2004). These values are considered low to moderate, suggesting that AFC would respond better to environmental improvement programs than to selection processes.

The heritability for MY270 estimated in this study was similar to that reported in Murrah buffaloes in Brazil by Aspilcueta-Borquis et al. (2010), Mendes Malhado et al. (2013), and Tonhati et al. (2008), 0.26, 0.25, and 0.28, respectively. A heritability estimate equal to 0.22 was obtained for the same trait in Colombia by other researchers (Hurtado-Lugo et al., 2006), while 0.24 was reported for milk production in Italy (Rosati and Van Vleck, 2002).

Heritability of the maternal genetic component for W12 was 0.05 and 0.06 in models PC1:3 and MC, respectively, and it was 0.05 for WW in both models. Those values indicate that no major changes would be expected in the short term for the maternal genetic effect as a result of selection based on WW or W12. Boligon et al. (2013) reported similar estimates for these traits in the Nellore breed. Also in Nellore, the maternal genetic effect was low, increasing from 0.03 at birth to 0.14 at 240 d, and then it decreased to zero at 600 d (Galvão de Albuquerque and Meyer, 2001).

The genetic correlation between AFC and body weight traits was negative in both models. This correlation was between -0.76 and -0.44 in model PC1:3 and between -0.76 and -0.51 in model MC (Table 4), indicating that selected animals at heavier weights and earlier age will be more precocious, starting to breed earlier. Caetano et al. (2013) reported a genetic correlation of -0.34 between age at first calving and weight gain in Nellore between 210 and 365 d, while Buzanskas et al. (2013) reported -0.15 between AFC and weight at 420 d in Canchim cattle; both correlations are also negative but less intense than those estimated in this work.

The genetic correlations between AFC and MY270 were low in PC1:3 and MC (-0.13 and -0.14, respectively). In Brazil, Seno et al. (2010) reported a similar genetic correlation of -0.12 between AFC and MY270, suggesting that selection of animals with less AFC may slightly increase milk yield in the first lactation.

Genetic correlations between WW, W12, W18, and W24 were positive and high in the PC1:3 and MC models (Table 4). This can be explained by the existing part-whole component because the first weight records are correlated with weight data taken in subsequent periods. In model PC1:3, the highest genetic correlation was in W12 with W18, with a value of 0.99, and the lowest was between WW and W12 with W24, with a value of 0.79. This lower value may be because they are the most distant measurements over time.

Another study in buffaloes estimated 0.85, 0.54, and 0.91 genetic correlations between W12-W18, W12-W24, and W18-W24, respectively (Bolívar-Vergara et al., 2012). Genetic correlations in Bosmara cattle were 1.00, 0.97, and 0.97 between WW-W12, WW-W18, and W12-W18, respectively (Bignardi et al., 2014). Similar data were reported in Nellore (Boligon et al., 2009). The genetic correlation values estimated in this work suggest that the genes responsible for WW are also highly responsible for later weights, which would allow for conducting the selection process at an early age.

The genetic correlations of MY270 with WW, W12, W18, and W24 estimated by the PC1:3 model were greater than 0.41 (Table 4), suggesting that animals selected at high weight and early age may have increased milk production during first lactation. Some studies have reported negative genetic correlations between live-weight traits and milk yield (Vercesi Filho et al., 2007). Other studies have reported positive genetic correlations, in agreement with those reported in this work (Moore et al., 1991; Van der Waaij, Galesloot, and Garrick, 1997). According to the data obtained in this work, selection of heavier animals at weaning would increase weight at 18 and 24 mo, increase milk production in first lactation, and lower the age at first calving.

Genetic Progress

The genetic trend of EBV and the standard error (indicated by the shaded area) for traits associated with meat production, WW, W12, W18, and W24 are presented in Fig. 1. Genetic progress in males and females was nonexistent until 2010, possibly because of a lack of clear selection criteria. Between 2010 and 2012, males showed a decrease in EBV and females showed a slight increase. During 2000 to 2014, W18 was the trait with the greatest increase: 0.77 and 0.59 kg/yr for males and females, respectively, while WW had the lowest genetic progress: 0.43 and 0.36 kg/yr for males and females, respectively. The greatest standard error occurred during the 1998–2000 period because of fewer production records and pedigree information.

Figure 1.
Figure 1.

Genetic trend of estimated breeding values (EBV) for WW, W12, W18, and W24 in dual-purpose buffalo cattle in Colombia by year of calving.


Genetic progress of males and females for W12 was 0.54 and 0.46 kg/yr, and for W24 it was 0.63 and 0.57 kg/yr, respectively. In Creole cattle, the genetic trend for WW was reported in Thailand as 0.32 kg/yr (Intaratham et al., 2008). In Nellore, the tendency for weight at 205, 365, and 505 d was 0.27, 0.45, and 0.74 kg/yr, respectively, over a 26 yr period (de Assunção Sousa et al., 2013).

The genetic trend for AFC and MY270 is shown in Fig. 2. In the 2012–2014 period, the genetic trend for AFC was -0.24 and -0.18 d/yr for males and females, respectively. MY270 showed a genetic tendency of 3.27 and 2.61 kg/yr in males and females, respectively. In Chile, 1.91 kg/yr genetic trend increase for milk yield was reported for Overo Colorado cattle during 1978–1999 (Uribe and Smulders, 2004). In the United States, genetic trend increases for milk production in Holstein were 79, 102, and 116 kg/yr in the 1970s, 1980s, and 1990s, respectively (Hansen, 2000).

Figure 2.
Figure 2.

Genetic trend of estimated breeding values (EBV) for AFC and MY270 by year of calving in dual-purpose buffaloes in Colombia


During 2012 to 2014, an increasing tendency for EBV of genetic traits associated with meat production and MY270 and a decrease for EBV associated with AFC were observed (Fig. 1 and 2). This coincides with the estimated genetic correlations among all traits. There are young animals in this period that had not yet been used for reproduction, suggesting that if selection is conducted according to EBV, major changes in the genetic trend might be expected for future generations, resulting in further progress and increased productivity of buffalo herds in Colombia.


The trait that could have a greater response to the selection process is W18. Nevertheless, MY270 was the trait with the greatest genetic progress during the evaluation period, probably due to the empirical selection conducted for this trait. The maternal genetic component for WW, W12, and AFC would have a lower response to the selection process due to low heritabilities. Reduced range models allow for a good estimation of both genetic parameters and genetic evaluation for economically important traits in Colombian buffaloes.




Be the first to comment.

Please log in to post a comment.
*Society members, certified professionals, and authors are permitted to comment.