Research| Volume 104, ISSUE 11, P11779-11789, November 01, 2021

# Genomic prediction of crossbred dairy cattle in Tanzania: A route to productivity gains in smallholder dairy systems

Open AccessPublished:August 04, 2021

## ABSTRACT

Selection based on genomic predictions has become the method of choice for genetic improvement in dairy cattle. This offers huge opportunity for developing countries with little or no pedigree data, and preliminary studies have shown promising results. The African Dairy Genetic Gains (ADGG) project initiated a digital system of dairy performance data collection, accompanied by genotyping in Tanzania in 2016. Currently, ADGG has the largest body of dairy performance data generated in East Africa from a smallholder dairy system. This study examines the use of genomic best linear unbiased prediction (GBLUP) and single-step (ss)GBLUP for the estimation of genetic parameters and accuracy of genomic prediction for daily milk yield and body weight in Tanzania. The estimates of heritability for daily milk yield from GBLUP and ssGBLUP were essentially the same, at 0.12 ± 0.03. The heritability estimates for daily milk yield averaged over the whole lactation from random regression model (RRM) GBLUP or ssGBLUP were 0.22 and 0.24, respectively. The heritability of body weight from GBLUP was 0.24 ± 04 but was 0.22 ± 04 from the ssGBLUP analysis. Accuracy of genomic prediction for milk yield from a forward validation was 0.57 for GBLUP based on fixed regression model or 0.55 from an RRM. Corresponding estimates from ssGBLUP were 0.59 and 0.53, respectively. Accuracy for body weight, however, was much higher at 0.83 from GBLUP and 0.77 for ssGBLUP. The moderate to high levels of accuracy of genomic prediction (0.53–0.83) obtained for milk yield and body weight indicate that selection on the basis of genomic prediction is feasible in smallholder dairy systems and most probably the only initial possible pathway to implementing sustained genetic improvement programs in such systems.

## INTRODUCTION

Genomic selection has become the method of choice for improvement of dairy cattle in developed countries due to the accelerated rate of genetic progress, mostly because of reduction in generation interval and increased prediction accuracy.
• García-Ruiz A.
• Cole J.B.
• Wiggans G.R.
• Ruiz-López F.J.
• Van Tassell C.P.
Changes in genetic selection differentials and generation intervals in US Holstein dairy cattle as a result of genomic selection.
reported that genetic gains for milk, fat, and protein yields for registered cows in the United States were 50, 2.2, and 1.6 kg/yr before genomic selection was adopted but increased by more than 2-fold to 109, 6.0, and 4.1 kg/yr after genomic selection was introduced. In developing countries, with limited or no pedigree data, the use of genomic information offers a huge opportunity for prediction of the genetic merit of animals using the G matrix. In East Africa, preliminary application of genomic prediction for the estimation of genetic parameters for milk yield have been reported based on Dairy Genetics East Africa (DGEA) data, using various models such as fixed repeatability and random regression genomic best linear unbiased prediction (GBLUP;
• Brown A.
• Ojango J.
• Gibson J.
• Coffey M.
• Okeyo M.
• Mrode R.
Short communication: Genomic selection in a crossbred cattle population using data from the Dairy Genetics East Africa project.
;
• Ojango J.M.K.
• Mrode R.
• Rege E.O.
• Mujibi D.
• Strucken E.M.
• Gibson J.
• Mwai O.
Genetic evaluation of test-day milk yields from smallholder dairy production systems in Kenya using genomic relationships.
). The data were rather limited, with low validation accuracies reported (0.28–0.41).
In general, the calibration slopes for genomic prediction from small cattle populations are widely varied, either less than unity (0.29–0.97) or much higher (1.16–1.24;
• Mrode R.
• Ojango J.M.K.
• Okeyo A.M.
• Mwacharo J.M.
Genomic selection and use of molecular tools in breeding programs for indigenous and crossbred cattle in developing countries: Current status and future prospects.
), indicating huge underpredictions for the study, with a regression coefficient of 0.29. However, the African Dairy Genetic Gains project (ADGG; https://www.ilri.org/research/projects/african-dairy-genetic-gains), which began initiating a digital system of dairy performance data collection in Tanzania in 2016, accompanied by genotyping, is currently the largest body of data generated in East Africa from a smallholder dairy system setting. This paper examines several models for the estimation of genetic parameters of daily milk yield and body weight, and accuracy of genomic prediction using GBLUP and single-step (ss)GBLUP in Tanzania.

## MATERIALS AND METHODS

### Genotypes

Genotypic data for this study were from 5,268 cows and bulls genotyped with the GeneSeek Genomic Profiler Bovine 50K chip. About 47,843 SNP were returned from the laboratory, and, after the usual edits, 40,581 SNP were available for analysis. These were imputed to the Illumina HD chip using a reference population consisting of crossbred cattle from a previous East Africa Dairy Genetics Gain project and several European Holstein-Friesian, Jersey, Guernsey, and Ayrshire purebred animals (
• Aliloo H.
• Mrode R.
• Okeyo A.M.
• Ni G.
• Goddard M.E.
• Gibson J.P.
The feasibility of using low density marker panels for genotype imputation and genomic prediction of crossbred dairy cattle of East Africa.
). The sampled cows were crosses between indigenous East African breeds, most of the latter being ancient admixtures of African Bos taurus and Bos indicus (
• Kim K.
• Kwon T.
• Dessie T.
• Yoo D.A.
• Mwai O.A.
• Jang J.
• Sung S.
• Lee S.B.
• Salim B.
• Jung J.
• Jeong H.
• Tarekegn G.M.
• Tijjani A.
• Lim D.
• Cho S.
• Oh S.J.
• Lee H.-K.
• Kim J.
• Jeong C.
• Kemp S.
• Hanotte O.
• Kim H.
The mosaic genome of indigenous African cattle as a unique genetic resource for African pastoralism.
), represented by N'dama and Small Eastern Africa Zebu, respectively, and 5 exotic Bos taurus dairy breeds. The percentages of indigenous and exotic genes for each cow were estimated by an admixture analysis (
• Aliloo H.
• Mrode R.
• Okeyo A.M.
• Ni G.
• Goddard M.E.
• Gibson J.P.
The feasibility of using low density marker panels for genotype imputation and genomic prediction of crossbred dairy cattle of East Africa.
). The exotic dairy percentage in each cow was then computed as the sum of the estimated percentage contributions of each of the 5 exotic dairy breeds. Four classes of animals were then created, based on the percentage of exotic genes, as follows: cows with >87.5%, 61−87.5%, 36−60%, and <36% exotic genes. These will account for effects of breed composition in all subsequent analyses.

### Performance Data

Performance data consisting of test day milk yield and body weight were extracted from the International Livestock Research Institute database (https://www.adgg.ilri.org/uat/auth/auth/login) for cows registered in the ADGG project from November 2016 until May 2020 in Tanzania. The data consisted of 134,987 test day records for 14,741 cows from 8,735 farms or herds, with pedigree information comprising 63,889 records. The herds were nested within 673 villages, within 263 wards, within 34 districts, and within 12 regions of Tanzania. Wards in the study area are administrative units comprising several villages. The data were subjected to several edits, including the identification of either sire or dam; daily milk yields values being in the range of >1 kg/d to ≤45 kg/d; age at first calving of at least 18 mo; and DIM restricted to between 4 and 500 d. The upper limit was based on a preliminary analysis that showed a reduction in heritability for milk yield when records with DIM greater than 500 were included. Also, a lower limit of 305 d, for instance, does not reflect the real production pattern in the population due to tendency of farmers to milk cow longer for various reasons: to continue to provide a source of income, inability to buy a replacement cow on time, or failure of cow to become pregnant, due to fertility issues. After edits, 3,473 cows with 36,376 test day milk records were available for analysis, with lack of either sire or dam and genotypes accounting for 73% of lost records. No requirement for the presence of either parent was imposed in the editing procedure for cows that had genotypes. Body weight was captured at the time of milk recording using a heart girth tape, and these were converted to body weight using the formulas published for smallholder systems by
• Lukuyu M.N.
• Gibson J.P.
• Savage D.B.
• Duncan A.J.
• Mujibi F.D.N.
• Okeyo A.M.
Use of body linear measurements to estimate liveweight of crossbred dairy cattle in smallholder farms in Kenya.
.
For the purposes of parameter estimation, 2 data files were created: Data1, consisting of 1,906 cows and 19,015 test day records with genotypes and minimum of 3 test days in first recorded lactation, and Data2, which was similar to Data1 but included cows with at least 1 parent known, resulting in a total of 2,716 cows with 25,474 test day records and a pedigree with 6,550 animals.
A fixed regression model (FRM) was implemented for both data sets for daily milk yield and body weight using ASRemL (
• Gilmour A.R.
• Gogel B.J.
• Cullis B.R.
• Thompson R.
ASReml User Guide Release 3.0.
). The FRM fitted for each trait, in matrix form, was as follows:
y = 1 + Φβ2 + Z1hd + Z2a + Z3pe + e,
[1]

where y is the vector of observations (milk yield or body weight); β1 is the vector of fixed effects, consisting of wards, calving year-season, test year-month, the 4 classes for the proportion of exotic genes, and age effects nested within lactation class; β2 are vectors of parameters for the fixed lactation curves modeled by Legendre polynomials of order 2, which were nested within lactation class; hd, a, and pe represent the random effects of herd, animal, and permanent environmental effects, respectively. The design matrices X, Z1, Z2, and Z3 related observations to fixed effects, random herd, pe, and animal effects, respectively; and Φ was the matrix of Legendre polynomials of order 2 to model fixed lactation curves. The choice of Legendre polynomial of order 2 was based on a preliminary analysis that examined several orders (1 to 5) to discover which gave the best fit using the Akaike and Schwarz's Bayesian information criteria (
• Wolfinger R.
Covariance structure selection in general mixed models.
). Although 5 parities were represented in the data, parity effects were fitted as only 3 classes: parities 1 and 2 separately, and parities 3 to 5 pooled into the third class. Broadly 2 season subclasses (dry or wet) were created per year, with January to June regarded as the dry seasons, but with some variations by regions. For Data1, model [1] was referred to as FRM-1, and the assumptions for the random effects were as follows: var(herd) = $Iσhd2$; var(pe) = $Iσpe2$; and var(a) = $Gσa2$ where G was the genomic relationship matrix. The matrix G was computed using method 1 of
Efficient methods to compute genomic predictions.
, as $G=MM′/2∑jn[pj(1-pj)],$ where M is the matrix of SNP genotypes with elements centralized as zij − 2pj, where pj is the allele frequency of the minor allele of the jth SNP. For Data2, model [1] was also fitted, and it is referred to as FRM-2, but with var(a) = $Hσa2$, where the H matrix was computed from the G matrix and the pedigree of animals (
• Misztal I.
• Legarra A.
• Aguilar I.
Computing procedures for genetic evaluation including phenotypic, full pedigree and genomic information.
). The inverse of H was computed as follows:
$H-1=A-1+(000G-1+A22-1).$

In FRM-2, 6 classes of breed composition were fitted instead of 4 classes as in FRM-1. The 2 additional classes were created for the 810 cows in Data2 that had no genotypes but had at least 1 parent identified. The lack of genotypes implies that breed compositions computed from the admixture analysis were not available for these animals; therefore, the breed types assigned by farmers were used. The breed types assigned by farmers indicated that 635 of these cows were from crosses of indigenous cattle with exotic breeds such as Ayrshires, Holstein-Friesians, Jerseys, and Guernsey, and these were assigned to one group, whereas 175 were termed indigenous or crosses among the indigenous, such as Small Eastern African Zebu, Maasai Zebu, or unknown, and were assigned to a second additional group.
In the various analyses detailed previously, herd effect was fitted as a random effect to accommodate the small sizes of herds (
• Powell O.
• Gaynor R.C.
• Jenko J.
• Gorjanc G.
• Mrode R.
• Hickey J.M.
Can genomic data enable genetic evaluation with phenotypes recorded on smallholder farms? Interbull Meeting, Auckland, New Zealand.
). However,
• Ojango J.M.K.
• Mrode R.
• Rege E.O.
• Mujibi D.
• Strucken E.M.
• Gibson J.
• Mwai O.
Genetic evaluation of test-day milk yields from smallholder dairy production systems in Kenya using genomic relationships.
demonstrated that the effects of herds are modeled either as fixed or random effects on estimates of heritability. To briefly investigate this influence, the FRM-1 model was implemented with herd effects fitted as a fixed effect nested within wards rather than a random effect.
In general, test day models for milk yield have been associated with higher estimates of heritabilities, due to the better correction of fixed and environmental factors relevant for each test day record and accounting for the shape of the lactation curve at the phenotypic and genetic levels. In an attempt to investigate the possible benefits of a test day model for milk yield in Tanzanian data, Data1 and Data2 were analyzed with a random regression model (RRM), which, in matrix form, was as follows:
y = 1 + Φβ2 + hd +Φ1a + Φ2pe + e,
[2]

where vectors of solutions β1, β2, and hd and design matrices X and Φ were as defined for model [1]; a and pe are vectors of animal and permanent environmental effects modeled with Legendre polynomials of orders 2 (Φ1) and 1 (Φ2), respectively. Attempts to fit Legendre polynomials of order 2 for both random effects did not converge, probably due to the limited data size. However, fitting Legendre polynomials of order 1 for both random effects gave a poorer fit based on Akaike or Schwarz's Bayesian information criteria. For Data1, the RRM was termed RRM-1, and the assumptions were as follows: var(pe) = P = IP, where P is a matrix of order 2 and ⊗ is the Kronecker product; and var(animal) = var(a) = GD, where D is a matrix of genetic covariance between random regression coefficients and a constant residual variance $(Iσe2)$ was assumed for all DIM. In the case of Data2, the only difference in terms of the assumptions was that var(a) = HM. In all the RRM models, daily estimates of genetic variance (vai) and permanent environmental variance (vpei) for day i were computed as vai = ϕ1i1i′ and vpei = ϕ2i2i′, respectively, where ϕji is the vector of Legendre polynomials corresponding to the ith day for the jth variance component. Daily estimates of heritability for the ith day were then obtained as $vai/(vai+vpei+σhd2+σe2)$. Genomic breeding values (gEBVik) for day i for the kth animal was computed as gEBVik = ϕ1iak, with ak = vector of random regression coefficients for the kth animal.
The models for the analyses of Data1 (FRM-1 or RRM-1), as previously, have involved use of the G matrix and therefore are classified as GBLUP subsequently in the paper. Also, the models for the analyses of Data2 (FRM-2 and RRM-2) have been termed ssGBLUP, as the inverse of H is used in these analyses. These ssGBLUP analyses involved the use of pedigree information; therefore, attempts were made to verify the accuracy of the pedigree information collected under a smallholder system. Using genotype data, parent-offspring relationships identified in the pedigree were verified in terms of Mendelian inconsistencies in the genotypes of offspring given the genotypes of parents.

### Validation of Genomic Prediction

Two types of validation were carried out for the genomic predictions. The first validation was an attempt to answer the question, “How accurate is the genomic system for predicting the genetic merit of candidate animals of different breed proportions, as classified by the 4 classes of animals by exotic genes?” A cross-validation was implemented that entails the exclusion of records for each class of cows of certain proportions of exotic genes and then undertaking genomic prediction. For the purposes of this cross-validation, yield deviations for all cows were computed as daily milk yield or body weight, corrected for all fixed and pe effects solutions from either model [1] or model [2] using the full data (no data excluded) and then averaged per cow. The accuracy of genomic predictions for the category of cows whose records were excluded from the prediction analysis was computed as the correlation of their gEBV with their yield deviations. Also, the degree of under- or overprediction was estimated by regressing their gEBV on their yield deviations. The cross-validation was implemented only for the fixed regression models FRM-1 and FRM-2 for milk yield and body weight.
The second validation was a forward validation, attempting to address the accuracy of predicting young or future candidate animals with only genotypes but no performance records. The forward validation was implemented by excluding the milk or body weight records of 270 cows born in 2014 and afterward in the genomic prediction system. Then the accuracy of genomic prediction and degree of under- or overprediction were computed as described for the cross-validation, using the gEBV and yield deviations for these 270 cows.

## RESULTS AND DISCUSSION

### General Data Summary

The overall means (±SD) for daily milk yield and body weight from Data1 were 8.3 ± 4.3 L and 299.9 ± 56.4 kg. The mean for daily milk yield for cows with >87.5%, 61−87.5%, 36−60%, and <36% exotic genes were 9.6 ± 4.2, 8.2 ± 4.2, 6.3 ± 3.6, and 4.5 ± 2.6 L, respectively, indicating a trend in milk production with increasing proportion of exotic genes. A similar trend was observed for body weight. Summaries of data for some variables in Data1 fitted in the various models are presented in Table 1. Also, simple plots of daily milk yield by weeks in milk for the 3 parity classes are shown in Figure 1. The plots for body weight have not been presented, as these showed very little variation by DIM.
Table 1Summary statistics of the data of 1,906 genotyped cows (Data1) analyzed for milk yield and BW
VariableNumber of cowsNumber of test daysMeans ± SD
Milk yield (L)BW (kg)
Parity
11,0265,8467.4 ± 4.1286.3 ± 60.0
21,2176,9448.3 ± 4.2298.1 ± 55.6
3–57876,2259.2 ± 4.4314.5 ± 49.9
Percentage of exotic genes
>87.57407,3599.6 ± 4.2316.7 ± 50.6
61–87.57727,9568.2 ± 4.2298.0 ± 58.0
36–603103,0056.3 ± 3.6273.7 ± 58.0
<36846954.5 ± 2.6255.3 ± 52.2
As expected, mean milk yield and body weight showed an increasing trend as the number of parities increased. The trends for mean milk yield and body weight for cows confirm the expected improvement in productivity with higher levels of exotic genes, but this was more pronounced for milk yield. The mean milk yield in cows with >87.5% exotic genes was twice that for cows with less than 36%. However, for body weight, this was only 25% higher in cows with >87.5% exotic genes compared with those with the least amount of exotic genes (<0.36%).

### Heritability Estimates for Milk

The estimates of heritability for Data1 and Data2 from the fixed regression models (FRM-1 and FRM-2, respectively) are presented in Table 2. The heritability estimates for daily milk from FRM-1 and FRM-2 were essentially the same at 0.12 ± 0.03, indicating that GBLUP and ssGBLUP gave similar results for this data set. The estimates of total variance due to herds followed the same pattern but the total variance due to permanent environment (PE) was slightly higher for the ssGBLUP. Working with various fixed regression models on the East Africa Dairy Genetics Gain data in Kenya,
• Ojango J.M.K.
• Mrode R.
• Rege E.O.
• Mujibi D.
• Strucken E.M.
• Gibson J.
• Mwai O.
Genetic evaluation of test-day milk yields from smallholder dairy production systems in Kenya using genomic relationships.
reported heritabilities that ranged from 0.05 ± 0.03 to 0.27 ± 0.05, but the estimate of heritability from the closest model to that fitted in this study was 0.09 ± 0.05, which is generally similar to the estimate obtained in this study. Generally, heritability estimates for daily test day milk from a fixed regression model are lacking in literature, because most estimates are based on 305-d milk yield or a random regression. It was expected that the use of ssGBLUP, which included more information, would have yielded higher heritability compared with GBLUP in this study, but similar estimates were obtained instead. Attempts were therefore made to understand the relationship between the genotyped and ungenotyped animals in the H matrix. First, examining the 1,906 cows that had genotypes and were used to construct the G matrix, the pedigree indicated that 12, 8, 206, and 680 cows had, respectively, both parents, only sire, only dam, and no parents identified. Moreover, nonzero additive genetic relationships from the pedigree were observed only among 81 cows of all genotyped cows, with an average additive genetic relationship of 0.457 and a range of 0.25 to 0.5. This clearly underlines the importance of the availability of genotypic information in enabling prediction of the genetic merit in smallholder systems, as the pedigree relationships are clearly inadequate. Of the 801 ungenotyped cows included in the ssGBLUP analysis, 62, 43, and 705 had, respectively, both parents, only sire, and only dam identified. However, nonzero additive genetic relationships occurred only among 251 of these cows, with an average additive relationship of 0.255, with a range of 0.0625 to 0.5. In addition, only 56 nonzero additive genetic relationships occurred between the 1,906 genotyped and 801 of the ungenotyped cows, with an average relationship of 0.337, but with additive genetic relationships ranging from 0.125 to 0.5. It appears that the limited pedigree relationships in the H matrix in general could explain the lack of any improvements in genetic parameters in ssGBLUP compared with GBLUP in this study.
Table 2Genetic parameters for daily milk yield (L) and BW (kg) from various fixed and random regression models (±SE)
TraitModel
FRM-1 and FRM-2 = fixed regression model (model [1]; genomic BLUP) on 1,906 genotyped cows (Data1) and on 1,906 genotyped plus 810 ungenotyped cows (Data2), respectively. RRM-1 and RRM-2 = random regression (model [2]; single-step genomic BLUP) on Data1 and Data2, respectively.
HeritabilityVariance due to PE
PE = permanent environment.
Variance due to herdPhenotypic variance
Milk yieldFRM-10.12 ± 0.030.10 ± 0.030.23 ± 0.029.73 ± 0.19
FRM-20.12 ± 0.030.12 ± 0.030.22 ± 0.029.68 ± 0.16
RRM-10.220.140.219.76
RRM-20.240.150.219.72
BWFRM-10.24 ± 040.20 ± 0.040.22 ± 0.031,287.6 ± 33.2
FRM-20.22 ± 040.22 ± 040.26 ± 031,338.4 ± 29.9
1 FRM-1 and FRM-2 = fixed regression model (model 1; genomic BLUP) on 1,906 genotyped cows (Data1) and on 1,906 genotyped plus 810 ungenotyped cows (Data2), respectively. RRM-1 and RRM-2 = random regression (model 2; single-step genomic BLUP) on Data1 and Data2, respectively.
2 PE = permanent environment.
Although the pedigree information was scanty, as previously indicated, the recorded pedigrees were generally accurate. Using all available genotypes and pedigree information to investigate the accuracy of the pedigree relationships indicated that, for 7 progeny who were genotyped and whose parents both were known in the pedigree and also genotyped, Mendelian inconsistency varied from 0 to 0.3%. The same trend was observed for 4 genotyped progeny with only their sires known and genotyped. However, for 290 progeny that were genotyped as well as only their dams, 85% had Mendelian inconsistency of 0, and 4% and 11% had Mendelian inconsistencies of about 5 to 10% and 11 to 14% respectively. This level of accuracy is encouraging and indicates the need to support and prompt more farmers to record pedigree information, to avoid reliance on genotypic data, which is rather expensive for the prediction of genetic merit of animals, especially under the less commercially competitive smallholder systems. As farmers are both encouraged and supported to record more pedigree data over time, the availability of more pedigree information will help in designing effective strategic genotyping to help achieve increases in the accuracy of genetic predictions.
The corresponding heritability estimates for daily milk yield averaged over the whole lactation from RRM (Table 2) were about twice the estimates from the fixed regression model, but estimates of total variance due to PE or herds were generally similar for both models. The estimates for daily genetic, PE variances, and heritabilities for daily milk yield are shown in Figure 2. Apart from the heritability of 0.42 at 4 DIM, the estimates varied from 0.17 to 0.33, which is similar to trends reported for smallholder Holstein-Friesian cows in Ethiopia with a similar model (
• Meseret S.
• Tamir B.
• Negussie E.
Estimation of genetic parameters for first-lactation test day milk yield in Holstein Friesian cows fitting random regression models.
). In another study in Ethiopia, involving Horro, Boran, Friesian, Jersey, and Simmental breeds and their crosses from 3 research centers,
• Gebreyohannes G.
• Koonawootrittriron S.
• Elzo M.A.
• Suwanasopee T.
Estimation of genetic parameters using a random regression monthly test-day model in an Ethiopian dairy cattle population.
reported increasing heritability estimates from 0.17 at 30 DIM to 0.42 at 300 DIM from an RRM using the modified, incomplete, gamma function. However, the estimates reported by
• Ojango J.M.K.
• Mrode R.
• Rege E.O.
• Mujibi D.
• Strucken E.M.
• Gibson J.
• Mwai O.
Genetic evaluation of test-day milk yields from smallholder dairy production systems in Kenya using genomic relationships.
for DGEA Kenya crossbred dairy cattle were higher at most of stages of lactations, with minimum of 0.25 at 60 DIM and maximum of 0.41 at 150 DIM. The last 2 studies did not model PE effects with any function, and this could partly account for differences between their results and those obtained in this study.
Estimates of genetic parameters for highly admixed populations, as in the population studied, are uncommon, but the trend in estimates of genetic variances and heritabilities are similar to those reported for small population of Holstein cows in Brazil (
• Cobuci J.R.
• Euclydes R.F.
• Lopes P.S.
• Costa C.N.
• de Almeida Torres R.
• Pereira C.S.
Estimation of genetic parameters for test-day milk yield in Holstein cows using a random regression model.
). The higher estimates of genetic variances or heritabilities at the extremes of the stage of lactation has been referred to by
• Schaeffer L.R.
• Jamrozik J.
Random regression models: A longitudinal perspective.
as an artifact of using Legendre polynomials in RRM. They concluded that the use of other functions such as Splines have not removed this artifact and indicated that the best method for removing the artifact would be to have a genetic effect for each day in milk over the lactation period. However, such a model will not be feasible in practice, as daily milk records for the whole lactation is not available on a national scale.

### Heritability Estimates for Body Weight

The estimates of genetic parameters for body weight from the various models and data sets are presented in Table 2. The heritability of 0.24 ± 04 for FRM-1 from GBLUP is slightly higher than the estimates based on ssGBLUP (FRM-2). However, the reverse occurred in terms of total variance due to PE and herd effects, with estimates from ssGBLUP about 10% higher than those from GBLUP. This could be due to the limited genetic links between genotyped and ungenotyped animals in the H matrix used for ssGBLUP in this study.
A limited number of studies exist in the literature on body weight of crossbred dairy cows in smallholder systems. The heritability estimates obtained in this study (0.22–0.24) are lower than those estimated for Holsteins in developed countries, with the latter values varying from 0.39 to 0.60 (
• Berry D.P.
• Buckley F.
• Dillon P.
• Evans R.D.
• Rath M.
• Veerkamp R.F.
Genetic parameters for level and change of body condition score and body weight in dairy cows.
;
• Vallimont J.E.
• Dechow C.D.
• Daubert J.M.
• Dekleva M.W.
• Blum J.W.
• Barlieb C.M.
• Liu W.
• Varga G.A.
• Heinrichs A.J.
• Baumrucker C.R.
Genetic parameters of feed intake, production, body weight, body condition score, and selected type traits of Holstein cows in commercial tie-stall barns.
), but these estimates are usually from random regression models, given the large variation in body weight at different stages of lactation. This is due to the strong selection over many generations on stature and increased milk production in the Holstein-Friesian in most developed countries, which is accompanied by increasing body weight. This is in sharp contrast to the rather flat lactation curve (Figure 1) with no peak observed in our data, due to the lack of selection for milk yield in the studied populations, and may also explain the lower heritability estimates.

### Modeling Herd Effect

The modeling of herd effects as fixed effects rather than random effects in the FRM-1 model for milk yield gave a higher heritability estimate of 0.23 ± 0.08, which is 92% higher than the estimate from model [1]; a slightly lower total variance due to PE effects of 0.08 ± 0.06, but total phenotypic variance was much reduced to 7.7 ± 0.24 compared with phenotypic variance of 9.73 ± 0.19 from model [1]. In the case of body weight, the estimate of heritability from modeling of herd as a fixed effect in the model FRM-1 was 0.33 ± 0.10, which is 67% higher than the estimate from the model with random herd effect; whereas variance due to PE increased by 50% to 0.24 ± 0.09, but total phenotypic variance was lower at 1,011.5 ± 46.9, compared with a phenotypic variance of 1,287.6 ± 33.2 from model [1]. In general, the sampling variances associated with the estimates for FRM-1 with fixed herd effects for both traits were larger compared with estimates from models with random herd effects. The increased heritability estimates from models with herd modeled as fixed rather than random, are due mainly to increases in genetic variances. This could indicate the difficulty of separating animal and herd effects when herd sizes are small (i.e., 1–3 cows per herd), particularly in cases where herds have only 1 cow and both herd and animal effects are therefore fully confounded. It seems that a random herd effect provides a better approach to modeling herd effects with genetic links across herds, due to the G matrix helping separate both effects (
• Powell O.
• Gaynor R.C.
• Jenko J.
• Gorjanc G.
• Mrode R.
• Hickey J.M.
Can genomic data enable genetic evaluation with phenotypes recorded on smallholder farms? Interbull Meeting, Auckland, New Zealand.
)

### Genetic Correlation Estimates

The genetic and phenotypic correlations among different DIM from the RRM are presented in Table 3. As expected, estimates of genetic correlation among different stages of lactation were high between adjacent DIM but decreased as intervals increased. The phenotypic correlations were low to medium in value, ranging from 0.07 to 0.50, and tended to follow the same pattern. The genetic correlations between early (4–84 DIM) and late (324–484 DIM) stages of lactation were slightly negative, but these correlations were generally very low (−0.17 to −0.05) and could essentially be regarded as zero.
• Ojango J.M.K.
• Mrode R.
• Rege E.O.
• Mujibi D.
• Strucken E.M.
• Gibson J.
• Mwai O.
Genetic evaluation of test-day milk yields from smallholder dairy production systems in Kenya using genomic relationships.
reported similar results for the DGEA data in Kenya, but the range of estimates was larger, varying from −0.57 to −0.05. In their study, the limited data size implied that Legendre polynomials of order 2 and 0 were fitted for the additive genetic and PE effects, respectively. However,
• Meseret S.
• Tamir B.
• Negussie E.
Estimation of genetic parameters for first-lactation test day milk yield in Holstein Friesian cows fitting random regression models.
also reported the problem of negative genetic correlations for Holstein-Friesian data in Ethiopia, from a model with Legendre polynomials of order 1 and 2 for additive genetic effects and PE effects, respectively; the problem was resolved by fitting an order of 2 for both additive genetic effects and PE effects. In this study Legendre polynomials of order 2 and 1 were fitted for additive genetic effects and PE effects, respectively, and problems of convergence occurred when higher orders were fitted. This illustrates the difficulty of choosing the best model for the analysis of the limited data observed in smallholder dairy systems. Therefore, the final model fitted may be based on what is feasible given the data structure, rather than what might be considered as theoretically optimum. In general, negative genetic correlation estimates between early and later milk test day records have been commonly reported for dairy animals raised in tropical environments (
• Rekaya R.
• Carabano M.J.
• Toro M.A.
Use of test day yields for the genetic evaluation of production traits in Holstein-Friesian cattle.
;
• Bignardi A.B.
• El Faro L.
• Cardoso V.L.
• de Albuquerque L.G.
Random regression models to estimate testday milk yield genetic parameters Holstein cows in South eastern Brazil.
;
• Singh A.
• Singh A.
• Singh M.
• Prakash V.
• Ambhore G.S.
• Sahoo S.K.
• Dash S.
Estimation of genetic parameters for first lactation monthly test-day milk yields using random regression test day model in Karan Fries cattle.
). Although limited data structure could partly explain this, it has also been attributed to the rapid decline in milk (low persistency) generally associated with Zebu cows. As more data accumulates, the influence of levels of exotic dairy genes on these genetic correlations could be examined by splitting the data based on the level of breed proportion. However, the genetic correlations obtained from early test day yield could be used to predict yield up to about 320 DIM.
Table 3Genetic (below diagonal) and phenotypic (above diagonal) correlations among several DIM with heritabilities (bold) on the diagonal from the random regression model (model 2) using 1,906 genotyped cows (Data1)
DIMDIM
44484124164204244284324364404444484
40.420.500.410.330.250.180.120.090.070.070.090.120.17
440.980.330.410.350.280.220.180.150.130.120.130.150.18
840.900.970.260.350.300.260.220.200.180.170.180.180.20
1240.740.860.960.220.320.290.260.240.230.220.210.210.21
1640.530.690.850.960.210.300.290.270.260.250.250.240.23
2040.330.510.710.870.970.200.300.290.290.280.270.260.24
2440.150.340.560.760.910.980.200.300.300.300.290.280.26
2840.010.200.430.650.820.930.980.190.310.310.310.300.29
324−0.090.090.300.530.720.850.940.980.180.320.320.320.32
364−0.15−0.000.190.400.590.730.840.920.970.170.330.340.35
404−0.17−0.070.070.240.410.550.660.770.870.960.160.360.38
444−0.14−0.10−0.030.070.180.290.410.530.670.810.950.180.41
484−0.06−0.08−0.09−0.08−0.040.030.130.250.400.590.800.950.22
The genetic correlation and phenotypic correlations between milk yield and body weight from the bivariate model were 0.20 ± 0.01 and 0.34 ± 0.11. Most studies examining the relationship between body weight and milk yield have been reported for herds in developed countries, usually from an RRM, as explained earlier. Therefore, such results are not really comparable with the estimates obtained in this study. However,
• Berry D.P.
• Buckley F.
• Dillon P.
• Evans R.D.
• Rath M.
• Veerkamp R.F.
Genetic parameters for level and change of body condition score and body weight in dairy cows.
reported a genetic correlation of 0.32 between body weight and milk yield at 60 DIM in Holstein-Friesian cows, but negative estimates were reported between milk yield and body weight for subsequent stages of lactation. The estimates they reported between cumulative milk yield at d 240 and body weight at various stages varied between 0.07 to 0.18, which is slightly lower than the estimates in this study.
In smallholder systems, feed is a major limiting factor; therefore, an important selection goal could be efficient cows that produce as much milk as possible with less feed, without a corresponding increase in body weight. Such a goal would favor highly producing but small to medium-sized cows. Although the aim in the long term is to develop a proper index that takes into consideration economic factors, the results from this study have been used to develop a restricted selection index for Tanzania using milk yield and body weight, with the goal of improving the rate of milk production but keeping body weight constant. Selection for milk yield will result in more efficient animals in feed utilization (
• Li B.
• Fikse W.F.
• Løvendahl P.
• Lassen J.
• Lidauer M.H.
• Mäntysaari P.
• Berglund B.
Genetic heterogeneity of feed intake, energy-corrected milk, and body weight across lactation in primiparous Holstein, Nordic Red, and Jersey cows.
), and the restriction in body weight will reduce energy requirements for body maintenance and some foot and leg diseases associated with heavy body weight. However, restricting body weight while selecting for production could also produce unfavorable genetic responses on cows' fertility and health, due to body tissue mobilization for high yield. Therefore, the long-term goal is to include fertility and health traits in the index as well, to ensure the selection of more robust cows. Currently, ADGG is piloting the collection of somatic cell counts, and data on treatments of cows for various diseases and fertility are gradually accumulating.

### Validation

The results of the cross-validation for milk yield and body weight based on the percentage of exotic genes are presented in Table 4. As expected, the accuracy of genomic prediction was higher for body weight than for milk yield, because body weight has a higher heritability. The accuracy estimates for body weight varied from 0.47 to 0.71 for the GBLUP and from 0.48 to 0.66 for ssGBLUP, indicating higher accuracies from the former, most likely due to the limited additional pedigree information used in ssGBLUP. The estimates of regression coefficients were generally close to unity, indicating very good calibration for the prediction of future animals of different breed proportions. However, the regression coefficients varied from 0.71 to 0.78 for the groups of cows with the highest amounts of exotic genes, indicating some degree of overprediction for these groups. The cross-validation accuracy for daily milk yield varied from 0.32 to 0.49 for GBLUP, and estimates from ssGBLUP were generally similar. The regression coefficients were only unity for the 2 group of cows with the highest proportion of exotic genes (>87.5% and 61–87.5%), but a high degree of underprediction occurred for cows with 36 to 60% exotic genes, and slight overprediction was observed for the cows with the least amount of exotic genes, again possibly due to the limited data size given the low heritability observed in this study.
Table 4Cross-validation results for daily milk yield (L) and BW (kg)
TraitModel and parameters
FRM-1 and FRM-2 = fixed regression model (model [1]; genomic BLUP) on 1,906 genotyped cows (Data1) and on 1,906 genotyped plus 810 ungenotyped cows (Data2), respectively. Corr. = correlations. Reg. = regression coefficients.
Percentage of exotic genes
>0.875 (n = 740)0.61–0.875 (n = 772)0.36–0.60 (n = 310)<0.36 (n = 84)
Daily milk yieldFRM-1
Corr.0.370.370.490.32
Reg.1.11.11.40.85
FRM-2
Corr.0.390.390.490.31
Reg.1.01.01.30.78
BWFRM-1
Corr.0.470.600.710.68
Reg.0.711.31.10.95
FRM-2
Corr.0.480.560.660.57
Reg.0.781.21.10.88
1 FRM-1 and FRM-2 = fixed regression model (model 1; genomic BLUP) on 1,906 genotyped cows (Data1) and on 1,906 genotyped plus 810 ungenotyped cows (Data2), respectively. Corr. = correlations. Reg. = regression coefficients.
The average percentage of exotic genes for the 270 cows born in 2014 and afterward, which were excluded in the forward validation data, was 0.80, indicating a high degree of exotic genes in younger cows in Tanzania. The accuracies of genomic prediction for the forward validation for daily milk yield and body weight are presented in Table 5. Again, higher accuracies were obtained for body weight compared with milk yield, with values of 0.83 and 0.77 from GBLUP (FRM-1) and ssGBLUP (FRM-2), respectively, for body weight. Corresponding values for milk yield were 0.57 and 0.59, respectively, indicating an inconsistent pattern in terms of trend of estimates from GBLUP versus ssGBLUP. However, the estimates for accuracies from both methods are generally similar. Compared with the FRM, the accuracies for milk yield from the RRM were lower, at 0.55 (GBLUP) and 0.53 (ssGBLUP), but the difference was more pronounced for the latter. These slightly lower accuracies from the RRM could be due to the limited data structure, as indicated earlier, and this will be investigated further as more data accumulate in the project. The estimates of regression coefficients from the RRM varied from 0.92 to 1.0, indicating a little or no overprediction.
Table 5Forward validation results for daily milk yield (L) and BW (kg)
TraitMethod
FRM-1 and FRM-2 = fixed regression model (model [1]; genomic BLUP) on 1,906 genotyped cows (Data1) and on 1,906 genotyped plus 810 ungenotyped cows (Data2), respectively. RRM-1 and RRM-2 = random regression (model [2]; single-step genomic BLUP) on Data1 and Data2, respectively.
CorrelationRegression
Milk yieldFRM-10.571.1
FRM-20.591.0
RRM-10.551.0
RRM-20.530.92
BWFRM-10.831.0
FRM-20.771.1
1 FRM-1 and FRM-2 = fixed regression model (model 1; genomic BLUP) on 1,906 genotyped cows (Data1) and on 1,906 genotyped plus 810 ungenotyped cows (Data2), respectively. RRM-1 and RRM-2 = random regression (model 2; single-step genomic BLUP) on Data1 and Data2, respectively.
The small number of genotyped animals in genomic prediction studies in developing countries makes it difficult to clearly define separate reference and validation populations. Consequently, most validations are undertaken in test data sets created by either random or structured sampling from all genotyped animals, usually using a combination of both bulls and cows (
• Mrode R.
• Ojango J.M.K.
• Okeyo A.M.
• Mwacharo J.M.
Genomic selection and use of molecular tools in breeding programs for indigenous and crossbred cattle in developing countries: Current status and future prospects.
).
The cross-classified accuracies reported in this study for milk yield are slightly higher than the estimates of 0.32 to 0.41 reported by
• Brown A.
• Ojango J.
• Gibson J.
• Coffey M.
• Okeyo M.
• Mrode R.
Short communication: Genomic selection in a crossbred cattle population using data from the Dairy Genetics East Africa project.
for the DGEA data in Kenya but similar to the estimates for Gyr cattle in Brazil (
• Boison S.A.
• Utsunomiya A.T.H.
• Santos D.J.A.
• Neves H.H.R.
• Carvalheiro R.
• Meszaros G.
• Utsunomiya Y.T.
• do Carmo A.S.
• Verneque R.S.
• Panetto J.C.C.
• Garcia J.F.
• Solkner J.
• da Silva M.V.B.
Accuracy of genomic predictions in Gyr (Bos indicus) dairy cattle.
). The intermediate accuracies (0.60–0.71) for body weight for cows of 0.36 to 0.875 exotic gene composition indicate that candidate cows of these breed compositions could be selected with a good degree of confidence. The accuracies of 0.32 to 0.37 for milk yield from FRM-1 for cows with 61 to >87.5% and <36% exotic genes indicate that these can be used as guides to screen candidate individuals or select a team of bulls for breeding.
The estimates of accuracy from the forward validation are similarly higher than most of the estimates summarized by
• Mrode R.
• Ojango J.M.K.
• Okeyo A.M.
• Mwacharo J.M.
Genomic selection and use of molecular tools in breeding programs for indigenous and crossbred cattle in developing countries: Current status and future prospects.
for several traits from genomic studies in developing countries. Although these estimates are lower than estimates from developed countries for dairy production traits (0.70–0.90) with medium to high heritability (
• Moser G.
• Khatkar M.S.
• Hayes B.J.
Accuracy of direct genomic values in Holstein bulls and cows using subsets of SNP markers.
;
• Wiggans G.R.
• Cole J.B.
• Hubbard S.M.
• Sonstegard T.S.
Genomic selection in dairy cattle: The USA experience.
), the estimates of 0.57 (milk yield) to 0.83 (body weight) from this study are encouraging and are expected to improve with time and provide the opportunity for selecting candidates animals in smallholder dairy systems, where currently no genetic improvement program is in place.

## CONCLUSIONS

From these results, we can conclude that levels of additive genetic variation and, hence, heritability for both body weight and milk yield for crossbred populations in countries such as Tanzania are high enough to enable significant gains to be achieved through selection. This study provides useful parameters that can be used to inform the design and implementation of future breeding programs, in similar smallholder crossbred-based dairy systems. The importance of initially genotyping animals at the early stages of the breeding improvement program, to provide the basis of genetic prediction given the lack of pedigree in a smallholder dairy systems, followed by systematic procedures to subsequently record pedigree information, is clearly elucidated. The moderate to high levels of accuracy obtained indicate that selection on the basis of genomic prediction is not only feasible in smallholder dairy systems but is the most realistic pathway to realizing sustained genetic gains in these systems.

## ACKNOWLEDGMENTS

The authors thank the Bill and Melinda Gates Foundation (Seattle, WA) for funding the African Dairy Genetic Gains project (ADGG, hosted at the International Livestock Research Institute, Nairobi, Kenya). We also acknowledge contributions by the Tanzania Livestock Research Institute (TALIRI; Dodoma, Tanzania), the University of New England (Armidale, NSW, Australia), and Scotland's Rural College, (Edinburgh, Scotland) for cofinancing. We gratefully acknowledge contributions by TALIRI colleagues, as well as contributions from the other national partners and all the smallholder farmers in Tanzania who participated in the ADGG project and made their animals available for monitoring. The authors have not stated any conflicts of interest.

## REFERENCES

• Aliloo H.
• Mrode R.
• Okeyo A.M.
• Ni G.
• Goddard M.E.
• Gibson J.P.
The feasibility of using low density marker panels for genotype imputation and genomic prediction of crossbred dairy cattle of East Africa.
J. Dairy Sci. 2018; 101 (30077450): 9108-9127
• Berry D.P.
• Buckley F.
• Dillon P.
• Evans R.D.
• Rath M.
• Veerkamp R.F.
Genetic parameters for level and change of body condition score and body weight in dairy cows.
J. Dairy Sci. 2002; 85 (12214996): 2030-2039
• Bignardi A.B.
• El Faro L.
• Cardoso V.L.
• de Albuquerque L.G.
Random regression models to estimate testday milk yield genetic parameters Holstein cows in South eastern Brazil.
Livest. Sci. 2009; 123: 1-7
• Boison S.A.
• Utsunomiya A.T.H.
• Santos D.J.A.
• Neves H.H.R.
• Carvalheiro R.
• Meszaros G.
• Utsunomiya Y.T.
• do Carmo A.S.
• Verneque R.S.
• Panetto J.C.C.
• Garcia J.F.
• Solkner J.
• da Silva M.V.B.
Accuracy of genomic predictions in Gyr (Bos indicus) dairy cattle.
J. Dairy Sci. 2017; 100: 5479-5490
• Brown A.
• Ojango J.
• Gibson J.
• Coffey M.
• Okeyo M.
• Mrode R.
Short communication: Genomic selection in a crossbred cattle population using data from the Dairy Genetics East Africa project.
J. Dairy Sci. 2016; 99 (27423951): 7308-7312
• Cobuci J.R.
• Euclydes R.F.
• Lopes P.S.
• Costa C.N.
• de Almeida Torres R.
• Pereira C.S.
Estimation of genetic parameters for test-day milk yield in Holstein cows using a random regression model.
Genet. Mol. Biol. 2005; 28: 75-83
• García-Ruiz A.
• Cole J.B.
• Wiggans G.R.
• Ruiz-López F.J.
• Van Tassell C.P.
Changes in genetic selection differentials and generation intervals in US Holstein dairy cattle as a result of genomic selection.
Proc. Natl. Acad. Sci. USA. 2016; 113 (27354521): E3995-E4004
• Gebreyohannes G.
• Koonawootrittriron S.
• Elzo M.A.
• Suwanasopee T.
Estimation of genetic parameters using a random regression monthly test-day model in an Ethiopian dairy cattle population.
Agric. Nat. Resour. (Bangk.). 2016; 50: 64-70
• Gilmour A.R.
• Gogel B.J.
• Cullis B.R.
• Thompson R.
ASReml User Guide Release 3.0.
VSN International, 2009
• Kim K.
• Kwon T.
• Dessie T.
• Yoo D.A.
• Mwai O.A.
• Jang J.
• Sung S.
• Lee S.B.
• Salim B.
• Jung J.
• Jeong H.
• Tarekegn G.M.
• Tijjani A.
• Lim D.
• Cho S.
• Oh S.J.
• Lee H.-K.
• Kim J.
• Jeong C.
• Kemp S.
• Hanotte O.
• Kim H.
The mosaic genome of indigenous African cattle as a unique genetic resource for African pastoralism.
Nat. Genet. 2020; 52 (32989325): 1099-1110
• Li B.
• Fikse W.F.
• Løvendahl P.
• Lassen J.
• Lidauer M.H.
• Mäntysaari P.
• Berglund B.
Genetic heterogeneity of feed intake, energy-corrected milk, and body weight across lactation in primiparous Holstein, Nordic Red, and Jersey cows.
J. Dairy Sci. 2018; 101 (30146279): 10011-10021
• Lukuyu M.N.
• Gibson J.P.
• Savage D.B.
• Duncan A.J.
• Mujibi F.D.N.
• Okeyo A.M.
Use of body linear measurements to estimate liveweight of crossbred dairy cattle in smallholder farms in Kenya.
Springerplus. 2016; 5 (26839756): 63
• Meseret S.
• Tamir B.
• Negussie E.
Estimation of genetic parameters for first-lactation test day milk yield in Holstein Friesian cows fitting random regression models.
Int. J. Genet. Mol. Biol. 2015; 7: 25-31
• Misztal I.
• Legarra A.
• Aguilar I.
Computing procedures for genetic evaluation including phenotypic, full pedigree and genomic information.
J. Dairy Sci. 2009; 92: 4648-4655
• Moser G.
• Khatkar M.S.
• Hayes B.J.
Accuracy of direct genomic values in Holstein bulls and cows using subsets of SNP markers.
Genet. Sel. Evol. 2010; 42: 37
• Mrode R.
• Ojango J.M.K.
• Okeyo A.M.
• Mwacharo J.M.
Genomic selection and use of molecular tools in breeding programs for indigenous and crossbred cattle in developing countries: Current status and future prospects.
Front. Genet. 2019; 9 (30687382): 694
• Ojango J.M.K.
• Mrode R.
• Rege E.O.
• Mujibi D.
• Strucken E.M.
• Gibson J.
• Mwai O.
Genetic evaluation of test-day milk yields from smallholder dairy production systems in Kenya using genomic relationships.
J. Dairy Sci. 2019; 102: 5266-5278
• Powell O.
• Gaynor R.C.
• Jenko J.
• Gorjanc G.
• Mrode R.
• Hickey J.M.
Can genomic data enable genetic evaluation with phenotypes recorded on smallholder farms? Interbull Meeting, Auckland, New Zealand.
• Rekaya R.
• Carabano M.J.
• Toro M.A.
Use of test day yields for the genetic evaluation of production traits in Holstein-Friesian cattle.
Livest. Prod. Sci. 1999; 57: 203-217
• Schaeffer L.R.
• Jamrozik J.
Random regression models: A longitudinal perspective.
J. Anim. Breed. Genet. 2008; 125 (18479264): 145-146
• Singh A.
• Singh A.
• Singh M.
• Prakash V.
• Ambhore G.S.
• Sahoo S.K.
• Dash S.
Estimation of genetic parameters for first lactation monthly test-day milk yields using random regression test day model in Karan Fries cattle.
Asian-Australas. J. Anim. Sci. 2016; 29 (26954137): 775-781
• Vallimont J.E.
• Dechow C.D.
• Daubert J.M.
• Dekleva M.W.
• Blum J.W.
• Barlieb C.M.
• Liu W.
• Varga G.A.
• Heinrichs A.J.
• Baumrucker C.R.
Genetic parameters of feed intake, production, body weight, body condition score, and selected type traits of Holstein cows in commercial tie-stall barns.
J. Dairy Sci. 2010; 93 (20855024): 4892-4901