Mitigation of greenhouse gases in dairy cattle via genetic selection: 2. Incorporating methane emissions into the breeding goal.

The objective of this study was to analyze the impact of incorporating enteric methane into the breeding objective of dairy cattle in Spain, and to evaluate both genetic and economic response of traits in the selection index under 4 scenarios: (1) the current ICO (Spanish total merit index), used as benchmark; (2) a hypothetical penalization of methane emissions through a carbon tax; (3) considering methane as a net energy loss for the animal; and (4) desired genetic response to reduce methane production by 20% in 10 yr. A bio-economic model was developed to derive the economic values for production and methane traits in each scenario. The estimated economic values for methane were estimated at -€1.21/kg and -€0.32/kg for scenarios 2 and 3, respectively. When merged with other traits in the selection index, methane had less economic importance (1-5%) than milk protein yield (39-42%) or milk fat yield (27-28%). Under these scenarios, selection resulted in an unfavorable response in methane emissions when it was included with an economic weight, with an increase in methane estimated from 0.52 to 0.60 kg/cow per year. Small differences in total profit per cow per year were observed between indices. The incorporation of methane production into the breeding objective had a negligible effect on production, with minor reductions in the expected genetic gain for fat and protein yields and in total economic benefits. However, total methane emissions in the dairy industry in Spain were estimated to decrease between 2 and 5% in the next 10 yr due to positive genetic trends for milk yield and an expected decrease in the total number of dairy cows. Additionally, methane intensity per 1 billion liters of milk would decrease in all scenarios. The uncertainty in the genetic parameters of methane and in carbon prices were tested in a sensitivity analysis, resulting in small deviations from the benchmark scenario. A major effect was observed only under the desired genetic response scenario. In this case, it was possible to achieve a 20% reduction of methane production in 10 yr via selective breeding but at the expense of a larger ad hoc weight (33%) of methane in the selection index and decelerating the genetic gain for production traits from 6 to 18%. This study shows the potential of including environmental traits in the selection indices while retaining populations profitable for producers.

Mitigation of greenhouse gases in dairy cattle via genetic selection: 1. Genetic parameters of direct methane using noninvasive methods and proxies of methane.

Records of methane emissions from 1,501 cows on 14 commercial farms in 4 regions of Spain were collected from May 2018 to June 2019. Methane concentrations (MeC) were measured using a nondispersive infrared methane detector installed within the feed bin of the automatic milking system during 14- to 21-d periods. Rumination time (RT; min/d) was collected using collars with a tag that registered time (minutes) spent eating and ruminating. The means of MeC and methane production (MeP) were 1,254.28 ppm and 182.49 g/d, respectively; mean RT was 473.38 min/d. Variance components for MeC, MeP, and RT were estimated with REML using pedigree and genomic information in a single-step model. Heritabilities for MeC and MeP were 0.11 and 0.12, respectively. Rumination time showed a slightly larger heritability estimate (0.17). The genetic correlation between MeP and MeC was high (>0.95), suggesting that selection on either trait would lead to a positive correlated response on the other. Negative correlations were estimated between RT and MeC (-0.24 ± 0.38) and MeP (-0.43 ± 0.35). Methane concentration and MeP had slightly positive correlations with milk yield (0.17 ± 0.39 and 0.21 ± 0.36), protein percentage (0.08 ± 0.32 and 0.30 ± 0.45), protein yield (0.22 ± 0.41 and 0.31 ± 0.35), fat percentage (0.02 ± 0.40 and 0.27 ± 0.36), and fat yield (0.27 ± 0.28 and 0.29 ± 0.28) from bivariate analyses. Rumination time had positive correlations with milk yield (0.41 ± 0.75) and protein yield (0.26 ± 0.57) and negative correlations with fat yield (-0.45 ± 0.32), protein percentage (-0.15 ± 0.38), and fat percentage (-0.40 ± 0.47). A positive approximated genetic correlation was estimated between fertility and MeC (0.10 ± 0.05) and MeP (0.18 ± 0.05), resulting in slightly higher CH4 production when selecting for better fertility [days open estimated breeding values (EBV) are expressed with mean 100 and SD 10, inversely related to days from calving to conception; that is, greater days open EBV implies better fertility]. Positive correlations were also estimated for stature with MeC and MeP (0.30 ± 0.04 and 0.43 ± 0.04, respectively). Other type traits (chest width, udder depth, angularity, and capacity) were positively correlated with methane traits, possibly because of higher milk yield and higher feed intake from these animals. Rumination time showed positive EBV correlations with production traits and type traits, and negative correlations with somatic cell count and body condition score. Based on the genetic correlations and heritabilities estimated in this study, methane is measurable and heritable, and estimates of genetic correlations suggest no strong opposition to current breeding objectives in Spanish Holsteins.

Genomic sequences of five novel HLA class I alleles: A*30:129, B*08:195, B*51:01:62, C*01:147 and C*12:195:02.

Five new HLA class I alleles are described, A*30:129, B*08:195, B*51:01:62, C*01:147 and C*12:195:02.

Characterization of a novel HLA-DRB1*07 allele, DRB1*07:83.

A new DRB1*07 allele, DRB1*07:83, was described in a Caucasian Spanish donor.

HLA-A*30:125 shows a new antigen binding domain created from HLA-A*30:01 and HLA-A*30:02.

A new HLA-A allele, A*30:125, was characterized in a Spanish individual.

Characterization of two new HLA-B alleles, HLA-B*07:299 and HLA-B*35:350.

Two novel HLA-B alleles, B*07:299 and B*35:350, were characterized by genomic full-length sequencing.

Genomic sequences of two novel HLA-C alleles, HLA-C*15:143 and HLA-C*07:109:02.

Two novel HLA-C alleles, C*07:109:02 and C*15:143, were characterized in two Spanish individuals.

Exon 2 sequencing of the new HLA-DRB1 allele, DRB1*13:216.

A novel HLA-DRB1*13:216 allele was characterized in a Spanish bone marrow donor.

Natural killer cell subsets in cerebrospinal fluid of patients with multiple sclerosis.

Changes in blood natural killer (NK) cells, important players of the immune innate system, have been described in multiple sclerosis (MS). We studied percentages and total cell counts of different effector and regulatory NK cells in cerebrospinal fluid (CSF) of MS patients and other neurological diseases to gain clearer knowledge of the role of these cells in neuroinflammation. NK cell subsets were assessed by flow cytometry in CSF of 85 consecutive MS patients (33 with active disease and 52 with stable MS), 16 with other inflammatory diseases of the central nervous system (IND) and 17 with non-inflammatory neurological diseases (NIND). MS patients showed a decrease in percentages of different CSF NK subpopulations compared to the NIND group. However, absolute cell counts showed a significant increase of all NK subsets in MS and IND patients, revealing that the decrease in percentages does not reflect a real reduction of these immune cells. Remarkably, MS patients showed a significant increase of regulatory/effector (CD56(bright) /CD56(dim) ) NK ratio compared to IND and NIND groups. In addition, MS activity associated with an expansion of NK T cells. These data show that NK cell subsets do not increase uniformly in all inflammatory neurological disease and suggest strongly that regulatory CD56(bright) and NK T cells may arise in CSF of MS patients as an attempt to counteract the CNS immune activation characteristic of the disease.

Cerebrospinal fluid immunological biomarkers associated with axonal damage in multiple sclerosis.

BACKGROUND AND PURPOSE: Cerebrospinal fluid (CSF) neurofilament light protein (NFL) is a promising biomarker of axonal injury and neurodegeneration. Here CSF lymphocyte subpopulations and antibodies, potential players of neurodegeneration, are examined in relation to CSF NFL shedding in MS. METHODS: Cerebrospinal fluid NFL from 127 consecutive untreated MS patients was analysed. Samples from 37 age-matched patients with other central nervous system non-inflammatory neurological diseases (NIND) were also assessed. CD4+, CD8+, CD56+ and CD19+ cell subsets were studied by flow cytometry. Oligoclonal IgG and IgM bands (OCMB) against lipids were studied by isoelectric focusing and immunoblotting. These data were analysed in relation to clinical and magnetic resonance imaging features. RESULTS: A CSF NFL cut-off value of 900 ng/l (mean + 3 SD of NIND values) was calculated. MS patients with increased NFL values showed significantly higher Multiple Sclerosis Severity Score and magnetic resonance imaging lesion number. The presence of OCMB (P < 0.0001) and elevated T and B lymphocyte counts was associated with increased levels of CSF NFL. CONCLUSIONS: High CSF NFL levels are associated with elevated CSF lymphocyte cell counts and intrathecal synthesis of IgM against lipids. These findings support a role for OCMB in the axonal damage of MS offering a rationale for the association of these antibodies with disability and brain atrophy progression in MS.

Identification of the major HHV-6 antigen recognized by cerebrospinal fluid IgG in multiple sclerosis.

BACKGROUND AND PURPOSE: Different data show an association between human herpesvirus 6 (HHV-6) and multiple sclerosis (MS). Intrathecal anti-HHV-6 immunoglobulin G (IgG) was detected in MS patients, but the antigen recognized by cerebrospinal fluid (CSF) IgG has not been characterized yet. Our objective was to identify the HHV-6 antigens recognized by IgG present in the CSF of patients with MS. METHODS: Cerebrospinal fluid IgG of 15 MS patients and eight patients with other neurological diseases was purified on protein G Sepharose columns. Purified IgG from every patient was linked to a CNBr-activated Sepharose 4B column. Fifty micrograms of viral extract was applied to each column. Bound proteins were eluted and analysed by SDS-PAGE and silver staining. The viral protein was characterized by mass spectrometry. RESULTS: A protein of 150 kD was eluted from CSF IgG columns of three of eight patients with primary progressive MS and one of seven with relapsing-remitting MS. After digestion and mass spectrometry analysis 10 peptides were found with 100% homology with the major capsid protein of the HHV-6A. DISCUSSION: These findings confirm the presence of anti-HHV-6 IgG in CSF of MS patients, particularly in progressive forms, and identify major capsid protein as the major antigen recognized by CSF IgG from MS patients.

Assets of imputation to ultra-high density for productive and functional traits.

The aim of this study was to evaluate different-density genotyping panels for genotype imputation and genomic prediction. Genotypes from customized Golden Gate Bovine3K BeadChip [LD3K; low-density (LD) 3,000-marker (3K); Illumina Inc., San Diego, CA] and BovineLD BeadChip [LD6K; 6,000-marker (6K); Illumina Inc.] panels were imputed to the BovineSNP50v2 BeadChip [50K; 50,000-marker; Illumina Inc.]. In addition, LD3K, LD6K, and 50K genotypes were imputed to a BovineHD BeadChip [HD; high-density 800,000-marker (800K) panel], and with predictive ability evaluated and compared subsequently. Comparisons of prediction accuracy were carried out using Random boosting and genomic BLUP. Four traits under selection in the Spanish Holstein population were used: milk yield, fat percentage (FP), somatic cell count, and days open (DO). Training sets at 50K density for imputation and prediction included 1,632 genotypes. Testing sets for imputation from LD to 50K contained 834 genotypes and testing sets for genomic evaluation included 383 bulls. The reference population genotyped at HD included 192 bulls. Imputation using BEAGLE software (http://faculty.washington.edu/browning/beagle/beagle.html) was effective for reconstruction of dense 50K and HD genotypes, even when a small reference population was used, with 98.3% of SNP correctly imputed. Random boosting outperformed genomic BLUP in terms of prediction reliability, mean squared error, and selection effectiveness of top animals in the case of FP. For other traits, however, no clear differences existed between methods. No differences were found between imputed LD and 50K genotypes, whereas evaluation of genotypes imputed to HD was on average across data set, method, and trait, 4% more accurate than 50K prediction, and showed smaller (2%) mean squared error of predictions. Similar bias in regression coefficients was found across data sets but regressions were 0.32 units closer to unity for DO when genotypes were imputed to HD density. Imputation to HD genotypes might produce higher stability in the genomic proofs of young candidates. Regarding selection effectiveness of top animals, more (2%) top bulls were classified correctly with imputed LD6K genotypes than with LD3K. When the original 50K genotypes were used, correct classification of top bulls increased by 1%, and when those genotypes were imputed to HD, 3% more top bulls were detected. Selection effectiveness could be slightly enhanced for certain traits such as FP, somatic cell count, or DO when genotypes are imputed to HD. Genetic evaluation units may consider a trait-dependent strategy in terms of method and genotype density for use in the genome-enhanced evaluations.

Reaction norm of fertility traits adjusted for protein and fat production level across lactations in Holstein cattle.

A total of 304,001 artificial insemination outcomes in up to 7 lactations from 142,389 Holstein cows, daughters of 5,349 sires and 101,433 dams, calving between January 1995 and December 2007 in 1,347 herds were studied by a reaction norm model. The (co)variance components for days to first service (DFS), days open, nonreturn rate in the first service (NRFS), and number of services per conception were estimated by 6 models: 3 Legendre polynomial degrees for the genetic effects and adjustment or not for the level of fat plus protein (FP) production recorded at day closest to DFS. For all traits and type of FP adjustment, a second degree polynomial showed the best fit. The use of the adjusted FP model did not increase the level of genetic (co)variance components except for DFS. The heritability for each of the traits was low in general (0.03-0.10) and increased from the first to fourth calving; nevertheless, very important variability was found for the estimated breeding value (EBV) of the sires. The genetic correlations (rg) were close to unity between adjacent calvings, but decreased for most distant parities, ranging from rg=0.36 (for DFS) to rg=0.63 (for NRFS), confirming the existence of heterogeneous genetic (co)variance components and EBV across lactations. The results of the eigen decomposition of rg shows that the first eigenvalue explained between 82 to 92% and the second between 8 to 14% of the genetic variance for all traits; therefore, a deformation of the overall mean trajectory for reproductive performance across the trajectory of the different calving could be expected if selection favored these eigenfunctions. The results of EBV for the 50 best sires showed a substantial reranking and variation in the shape of response across lactations. The more important aspect to highlight, however, is the difference between the EBV of the same sires in different calvings, a characteristic known as plasticity, which is particularly important for DFS and NRFS. This component of fertility adds another dimension to selection for fertility that can be used to change the negative genetic progress of reproductive performance presented in this population of Holstein cows. The use of a reaction norm model should allow producers to obtain more robust cows for maintenance of fertility levels along the whole productive life of the cows.

The gradient boosting algorithm and random boosting for genome-assisted evaluation in large data sets.

In the next few years, with the advent of high-density single nucleotide polymorphism (SNP) arrays and genome sequencing, genomic evaluation methods will need to deal with a large number of genetic variants and an increasing sample size. The boosting algorithm is a machine-learning technique that may alleviate the drawbacks of dealing with such large data sets. This algorithm combines different predictors in a sequential manner with some shrinkage on them; each predictor is applied consecutively to the residuals from the committee formed by the previous ones to form a final prediction based on a subset of covariates. Here, a detailed description is provided and examples using a toy data set are included. A modification of the algorithm called "random boosting" was proposed to increase predictive ability and decrease computation time of genome-assisted evaluation in large data sets. Random boosting uses a random selection of markers to add a subsequent weak learner to the predictive model. These modifications were applied to a real data set composed of 1,797 bulls genotyped for 39,714 SNP. Deregressed proofs of 4 yield traits and 1 type trait from January 2009 routine evaluations were used as dependent variables. A 2-fold cross-validation scenario was implemented. Sires born before 2005 were used as a training sample (1,576 and 1,562 for production and type traits, respectively), whereas younger sires were used as a testing sample to evaluate predictive ability of the algorithm on yet-to-be-observed phenotypes. Comparison with the original algorithm was provided. The predictive ability of the algorithm was measured as Pearson correlations between observed and predicted responses. Further, estimated bias was computed as the average difference between observed and predicted phenotypes. The results showed that the modification of the original boosting algorithm could be run in 1% of the time used with the original algorithm and with negligible differences in accuracy and bias. This modification may be used to speed the calculus of genome-assisted evaluation in large data sets such us those obtained from consortiums.

Comparison of methods for the implementation of genome-assisted evaluation of Spanish dairy cattle.

The aim of this study was to evaluate methods for genomic evaluation of the Spanish Holstein population as an initial step toward the implementation of routine genomic evaluations. This study provides a description of the population structure of progeny tested bulls in Spain at the genomic level and compares different genomic evaluation methods with regard to accuracy and bias. Two bayesian linear regression models, Bayes-A and Bayesian-LASSO (B-LASSO), as well as a machine learning algorithm, Random-Boosting (R-Boost), and BLUP using a realized genomic relationship matrix (G-BLUP), were compared. Five traits that are currently under selection in the Spanish Holstein population were used: milk yield, fat yield, protein yield, fat percentage, and udder depth. In total, genotypes from 1859 progeny tested bulls were used. The training sets were composed of bulls born before 2005; including 1601 bulls for production and 1574 bulls for type, whereas the testing sets contained 258 and 235 bulls born in 2005 or later for production and type, respectively. Deregressed proofs (DRP) from January 2009 Interbull (Uppsala, Sweden) evaluation were used as the dependent variables for bulls in the training sets, whereas DRP from the December 2011 DRPs Interbull evaluation were used to compare genomic predictions with progeny test results for bulls in the testing set. Genomic predictions were more accurate than traditional pedigree indices for predicting future progeny test results of young bulls. The gain in accuracy, due to inclusion of genomic data varied by trait and ranged from 0.04 to 0.42 Pearson correlation units. Results averaged across traits showed that B-LASSO had the highest accuracy with an advantage of 0.01, 0.03 and 0.03 points in Pearson correlation compared with R-Boost, Bayes-A, and G-BLUP, respectively. The B-LASSO predictions also showed the least bias (0.02, 0.03 and 0.10 SD units less than Bayes-A, R-Boost and G-BLUP, respectively) as measured by mean difference between genomic predictions and progeny test results. The R-Boosting algorithm provided genomic predictions with regression coefficients closer to unity, which is an alternative measure of bias, for 4 out of 5 traits and also resulted in mean squared errors estimates that were 2%, 10%, and 12% smaller than B-LASSO, Bayes-A, and G-BLUP, respectively. The observed prediction accuracy obtained with these methods was within the range of values expected for a population of similar size, suggesting that the prediction method and reference population described herein are appropriate for implementation of routine genome-assisted evaluations in Spanish dairy cattle. R-Boost is a competitive marker regression methodology in terms of predictive ability that can accommodate large data sets.

Genotyping strategies for genomic selection in small dairy cattle populations.

This study evaluated different female-selective genotyping strategies to increase the predictive accuracy of genomic breeding values (GBVs) in populations that have a limited number of sires with a large number of progeny. A simulated dairy population was utilized to address the aims of the study. The following selection strategies were used: random selection, two-tailed selection by yield deviations, two-tailed selection by breeding value, top yield deviation selection and top breeding value selection. For comparison, two other strategies, genotyping of sires and pedigree indexes from traditional genetic evaluation, were included in the analysis. Two scenarios were simulated, low heritability (h 2 = 0.10) and medium heritability (h 2 = 0.30). GBVs were estimated using the Bayesian Lasso. The accuracy of predicted GBVs using the two-tailed strategies was better than the accuracy obtained using other strategies (0.50 and 0.63 for the two-tailed selection by yield deviations strategy and 0.48 and 0.63 for the two-tailed selection by breeding values strategy in low- and medium-heritability scenarios, respectively, using 1000 genotyped cows). When 996 genotyped bulls were used as the training population, the sire' strategy led to accuracies of 0.48 and 0.55 for low- and medium-heritability traits, respectively. The Random strategies required larger training populations to outperform the accuracies of the pedigree index; however, selecting females from the top of the yield deviations or breeding values of the population did not improve accuracy relative to that of the pedigree index. Bias was found for all genotyping strategies considered, although the Top strategies produced the most biased predictions. Strategies that involve genotyping cows can be implemented in breeding programs that have a limited number of sires with a reliable progeny test. The results of this study showed that females that exhibited upper and lower extreme values within the distribution of yield deviations may be included as training population to increase reliability in small reference populations. The strategies that selected only the females that had high estimated breeding values or yield deviations produced suboptimal results.

Genetic evaluation of susceptibility to clinical mastitis in Spanish Holstein cows.

This study had 3 objectives: to estimate genetic parameters and predict sires' transmitting abilities for clinical mastitis in a Spanish Holstein population, to propose a methodology for comparing models with different response variables by using a cost-based loss function, and to evaluate alternative genetic evaluation models by using this methodology. On-farm records for clinical mastitis from herds in 3 Spanish regions were analyzed as a binary trait (CM) and as number of episodes (NCM) per lactation. Linear and probit models were fitted for CM, whereas linear and Poisson models were used for NCM. Predictive ability of the models was evaluated by using the average predicted residual sum of squares from cross-validation and an alternative cost-based loss function. The loss function for model comparison was calculated by using average mastitis costs depending on the NCM and average cost per infected lactation. The average cost per infected lactation was $345.58, whereas the cost per lactation ranged from $204.86 to $985.44 for lactations with 1 to 5 cases, respectively. Management and hygiene practices on individual farms had a large impact on clinical mastitis because the herd-year variance was larger than that of other random effects considered. The sire variance was significantly different from zero, confirming that genetic variation exists for clinical mastitis. Estimates of heritability for CM using the linear and probit models were 0.07 and 0.10 on the underlying scale, respectively. For NCM, the estimate of heritability for the linear model was 0.10 and estimates for the Poisson model evaluated at the mean and the median of lambda on the underlying scale were 0.09 and 0.07, respectively. Regarding ranking of sires, the definition of response variable (CM or NCM) was of greater importance than the choice of statistical model. Cross-validation results indicated that models with the best fit for CM and NCM were the probit model and the linear model, respectively. However, a comparison across all models using the alternative cost-based loss function showed that using NCM as a response variable with a Poisson model provided the most accurate predictions of future costs associated with clinical mastitis.

Selection for female fertility using censored fertility traits and investigation of the relationship with milk production.

Bivariate models (censored linear-linear and censored threshold-linear) were used to estimate genetic parameters for production and fertility traits in the Spanish Holstein population. Records on 71,217 lactations from 41,515 cows were used: 30 and 36% of lactations were censored for days open (DO) and number of inseminations to conception (INS), respectively. Heritability estimates for production traits (milk, fat, protein) ranged between 0.18 and 0.25. Heritability of days to first service (DFS) and DO was 0.05; heritability of INS on the liability scale was 0.04. Genetic correlations between fertility traits were 0.41, 0.71, and 0.87 for DFS-INS, DO-INS, and DO-DFS, respectively. Days open had a larger genetic correlation (ranging from 0.63 to 0.76) with production traits than did DFS (0.47 to 0.59) or INS (0.16 to 0.23). Greater antagonism between production and DO may be due to voluntary management decisions for high-yielding cows, resulting in longer lactation lengths. Inseminations to conception appeared to be less correlated with milk production than were the other 2 female fertility traits. Including INS in a total merit index would be expected to increase genetic gain in terms of profit, but profit would decrease if either DO or DO and DFS were included in the index. Thus, INS is the trait to be preferred when selecting for female fertility. The genetic correlation between actual milk yield and 305-d standardized milk yield was 0.96 in the present study, suggesting that some reranking of sires could occur. Because the target of attaining a 12-mo calving interval, as implied by a 305-d standardized lactation length, is changing in the dairy industry, routine genetic evaluation of actual total lactation milk yield should be considered.

Genetic and phenotypic relationships among locomotion type traits, profit, production, longevity, and fertility in Spanish dairy cows.

The phenotypic and genetic relationships of 3 locomotion traits with profit, production, longevity, and fertility traits were studied to determine the importance of locomotion traits for dairy producers. Two data sets including official milk records and type classification scores of 62,293 cows, and reproductive records of 24,561 cows from the Basque and Navarra Autonomous Regions were analyzed. Higher scores for feet and legs (FL), foot angle (FA), and rear legs set (RLS) were positively related to production and functional traits, whereas fertility was not significantly affected. The cows that scored the highest for FL were $213/yr more profitable, produced 575 kg more milk per year, and remained in the herd for 307 more functional days than the cows scoring the lowest. Feet and legs was the trait most genetically correlated to profit, although a low value (0.10) was obtained, whereas RLS was the trait most correlated to milk production (0.12). Genetic correlations among FL, FA, RLS, and longevity traits (from -0.10 to 0.05) were low. Quadratic curves were the best fit for both profit and functional herd life for EBV of each of the 3 locomotion traits. Further studies dealing with profitability and lameness, instead of using conformation traits, could be performed directly if a larger data pool of lameness was routinely recorded.

Genetic parameters for female fertility traits and a fertility index in Spanish dairy cattle.

Genetic correlations among female fertility traits (linear and binary) were estimated using 225,085 artificial insemination records from 120,713 lactations on 63,160 Holstein cows. Fertility traits were: calving interval, days open, a linear transformation of days open, days to first insemination, interval between first and last insemination, number of inseminations per service period, pregnancy within 56 and 90 d after first insemination, and success in first insemination. A bivariate animal model was implemented using Bayesian methods in the case of binary traits. Low heritabilities (0.02 to 0.06) were estimated for these fertility traits. Strong genetic correlations (0.89 to 0.99) were found among traits, except for days to first service, where the genetic correlation with other fertility traits ranged from -0.52 to -0.18 for binary traits, and from 0.50 to 0.82 for days to first service, calving interval, and days open. Four fertility indices were proposed utilizing information from insemination records; these indices combined one indicator of the beginning of the service period and one indicator of conception rate. Two additional indices used information from the milk-recording scheme, including calving interval and a linear transformation of days open. The fertility index composed of days to first service and pregnancy within 56 d achieved the highest genetic gain for reducing fertility cost, reducing days to first service, and reducing the number of inseminations per lactation ($8.60, -1.31 d, and -0.03 AI, respectively). This index achieved at least 15% higher genetic gain than obtained from indices with information from the milk recording scheme only (calving interval and days open).