Breeding Dairy Cattle for Female Fertility and Production in the Age of Genomics

Simple Summary Long-term selection should lead to reduction in heritability, due to fixation of positive alleles. Selection on multiple traits should lead to negative genetic correlations because alleles with positive effects on both traits will respond first to selection, leaving alleles with opposing effects. Efficient selection for milk production traits began in Israel in the 1980s, and selection for female fertility in 2000. Introduction of genomic selection increased rates of genetic gain. Many studies have shown the negative relationships between milk production traits and female fertility. Phenotypic and genetic changes for female fertility and production traits in the Israeli dairy cattle population over the last three decades were studied in order to determine if long-term selection has resulted in reduced heritability and increased negative genetic correlations. Heritabilities were 0.4 for protein production and 0.05 for conception status. The genetic correlation between conception status and protein yield was −0.38. Heritabilities decreased with increase in parity for protein but remained the same for conception status. For milk, fat, and protein production and female fertility, heritabilities increased or stayed the same over the entire period of 30 years. There is no indication that a selection plateau is imminent for dairy cattle. Abstract Phenotypic and genetic changes for female fertility and production traits in the Israeli Holstein population over the last three decades were studied in order to determine if long term selection has resulted in reduced heritability and negative genetic correlations. Annual means for conception status, defined as the inverse of the number of inseminations to conception in percent, decreased from 55.6 for cows born in 1983 to 46.5 for cows born in 2018. Mean estimated breeding values increased by 1.8% for cow born in 1981 to cows born in 2018. Phenotypic records from 1988 to 2016 for the nine Israeli breeding index traits were divided into three time periods for multi-trait REML analysis by the individual animal model. For all traits, heritabilities increased or stayed the same for the later time periods. Heritability for conception status was 0.05. The first parity genetic correlation between conception status and protein yield was −0.38. Heritabilities decreased with the increase in parity for protein but remained the same for conception status. Realized genetic trends were greater than expected for cows born from 2008 through 2016 for somatic cell score, conception status and herd-life, and lower than expected for the production traits.


Introduction
In 1960 Falconer [1] noted that selection on a multiple trait index should lead to the generation of negative genetic correlations among the index traits. This will occur because those polymorphic genes for which the same alleles have positive effects on both traits will respond first to selection, leaving the genes with economically negative correlations. It is also conventional wisdom that selection should result in a reduction in genetic variance, as alleles with positive effects reach fixation. However, as shown by [2] by simulation, and by [3] based on theoretical considerations, in some cases selection can lead to temporary increases in genetic variance due to the increase in frequency of positive rare alleles.
Genomic selection based on genotypes for thousands of genetic markers has been incorporated into nearly all commercial dairy cattle breeding programs over the last decade [4]. In all cases, genomic evaluation is performed on each individual trait included in the selection index without consideration of possible genetic correlations [4].
Female fertility traits were first included in national breeding indexes by the Scandinavian countries in the 1970s. Female fertility was included in the total merit index for Norwegian dairy cattle since 1972 by considering the 56-d nonreturn rate in virgin heifers, that is, the fraction of cows with a first insemination that were not re-inseminated [5]. Israel incorporated "conception status" in the national index in 2000 [6], and the US included "daughter pregnancy rate" in 2003 [7]. Development of the Israeli breeding index since 1985 is given in Table 1. Since 1990, protein production has been given the greatest emphasis in the breeding index. This has been the case for nearly all dairy commercial breeding programs since 2000 [8]. Genomic selection began in Israel in 2015 [9]. 1 Somatic cell score. 2 Negative values are economically favorable. 3 Conception status, defined in the text.
Much has been written about both the phenotypic and genetic economically negative correlation between milk production traits and female fertility, reviewed by [10,11]. Numerous studies have attempted to estimate the negative genetic correlations between milk production traits and various measures of female fertility, reviewed by [12].
One major problem with these studies, as noted by [13], is that nearly all measures of female fertility are biased with respect to milk production traits. With respect to the non-return rate, the fact that the cow was not re-inseminated does not mean that the first insemination was successful. The cow could have died, or the farmer could have decided to cull, and farmers are more likely to cull low-producing cows. Similarly, other measures of fertility, such as the number of inseminations per lactation and days open will also be biased, because low production cows are more likely to be culled before a days-open record is generated, and cows with the lowest fertility do not become pregnant. Most studies that have attempted to estimate genetic correlations between fertility and production traits have used sire models and have considered only first parity cows [12].
Insemination data in Israel is unique in that nearly all cows that are inseminated are vet-checked for pregnancy after 60 days [13]. Routine genetic evaluations are computed in Israel for "conception status" (CS), calculated as the inverse of the number of inseminations to conception in percent. For cows that are culled prior to conception or cows for which conception has not yet been recorded, the expected number of inseminations to conception is estimated [14]. Thus, all cows that are inseminated at least once have a record for this trait, and biases with respect to production traits are likely to be minimal.
In order to determine if long-term selection for fertility and production has in fact resulted in reduced heritability and more negative genetic correlations, we estimated the phenotypic and genetic changes for CS over the last three decades, and the environmental and genetic correlations among the traits included in the Israeli breeding index by the individual animal model, including CS and milk, fat, and protein production. Genetic and environmental variance components were computed for first parity for the nine index traits over the last 30 years, and for protein and CS for the first three parities over the last eight years. Genetic trends over the last decade were computed for all traits included in the Israeli breeding index, and the realized values were compared to the expected values derived from the principles of selection index. Finally, genome-wide association studies (GWAS) were performed for both protein production and CS, and correlations were computed between the marker effects.

The Data Sets and Traits Analyzed
The seven data sets analyzed are described in Table 2. The first data set included phenotypic records of cows born from 1988 through 2016 with valid records for herd-life and valid first parity records for the other eight index traits listed in Table 1. For cows born prior to 1988, some of the index traits were recorded only on subsets of the data. Records for cows born after 2016 were incomplete, especially for herd-life. This data set was divided into three parts for estimation of variance components and genetic parameters by restricted maximum likelihood (REML) methodology, as described in Table 2.  4 A single multi-trait analysis was performed including parities 1-3 of protein and CS. 5 Regression of the cows' breeding values on their birth years. 6 Birth years and number of records refer to bulls with reliabilities >0.5.
The second data set included phenotypic records of cows born between 2008 and 2015 with valid first parity records for protein and CS. This data set was used to estimate variance components and genetic parameters by REML for protein and CS for parities 1 through 3. Genetic and environmental correlations were computed among the three parities for both traits. Both data sets were analyzed by the multi-trait individual animal model (IAM).
The third and fourth data sets were used to compute genetic evaluations for the three milk production traits and CS, respectively, for all cows with valid records born from 1983 through 2018 by the multi-trait IAM, including valid records for parities 1 through 5. Data set 3 included only parities with valid records for all three milk production traits. Data set 5 included estimated breeding values (EBV) for the nine index traits from cows born between 2009 and 2018 and was used to compute realized genetic trends in the Israeli Holstein population. Data sets 6 and 7 included bulls with reliabilities >0.5 for protein and CS respectively, from the analyses of data sets 3 and 4, and genotypes for >40,000 SNPs. These two data sets were used for the GWAS analyses.
Edits applied for all traits included in PD19 prior to genetic evaluation were as described previously [15,16]. All the traits included in PD19, except for herd-life, were analyzed by a multi-trait IAM, with each parity considered a separate trait. Except for the calving traits, all valid parities up to fifth were included in the analyses. In addition to the additive genetic effects, the analysis models included the effects of herd-year-season and parity. The single parity evaluations were then combined into a multi-parity index as described previously [15]. Herd-life was computed as the number of days from first calving to culling and analyzed by a single trait animal model. For cows that have not yet been culled, expected herd-life was computed as described previously [17]. First and second parity dystocia and rate of stillbirth were analyzed jointly by a multi-trait IAM including the effects of the cow calving and the sire of calf as described by [18]. Reliabilities for all traits were estimated as described previously [18,19].

Statistical Analyses
Variance components were estimated for data set 1 with the MTC program [20] and the AIREMLf90 program for data set 2 [21]. The AIREML90 could not be applied to data set 1, in which nine traits were analyzed, because of software limitations. AIREML90 can accept different numbers of records per trait, differing analysis models for each trait and estimates of standard errors for estimates of all variance components and genetic parameters. Both data sets were analyzed by the multi-trait IAM. In addition to the random residual and the random additive genetic effect, data sets 1 and 2 included fixed herd-yearseason (HYS) effects. Two seasons were defined for each herd-year, based on freshening month, beginning in April and October of each year. In the analysis of data set 1, the same HYS classes were defined for all nine traits relative to first parity. In the analysis of data set 2, different HYS classes were defined for each parity. In the analysis of data set 1 the same number of records were analyzed for each trait, as required by the MTC program. Therefore, cows with missing records for any of the 9 index traits were deleted. In the analysis of data set 2, the numbers of valid records decreased with an increase in parity. To avoid possible bias, later parity records were included only if there were valid records for the previous parities. In both data sets, all known parents and grandparents of cows with first parity records and of sires of cows with records were included to construct the relationship matrix among animals. The total numbers of animals included in each analysis are also given in Table 2. For both data sets, two genetic groups were defined for animals with unknown parents or animals of the first generation in the data set, one for males and one for females.
Data sets 3 and 4, which included all valid records for cows born between 1983 and 2018, were analyzed by the multi-trait IAM. These models included fixed parity and HYS effects, as defined previously, in addition to the random additive genetic and residual effects. Records were pre-adjusted for birth and calving month and days open as described previously [15,16]. The magnitudes of the residual and genetic variance components were assumed to be known. A total of 92 groups for the milk production traits and 80 groups for CS were defined, based on the sex of the animal with unknown parents, which parent was unknown, and the birth year. In addition, separate groups were defined for sire of cows of breeds other than Holstein. Although <2% of the cows were sired by bulls of other breeds, these bulls were a significant fraction of the total number of bulls, and an even larger fraction of the bulls with unknown parents. The genetic base of the evaluations for all traits was set as the mean EBV of calves born in 2015.
The contribution of each of the nine traits included in the PD19 selection index (c j ) was computed as proposed by [22]: where b j = the index coefficient for trait j, g j = the genetic standard deviation for trait j, J = 9, the total number of traits, and "abs" denotes "absolute value". The values of b j were derived from Table 1, and the values for g j were derived from the REML analysis of data set 1 for cows born between 2008 and 2016. Realized genetic gains for the nine index traits were estimated from data set 5 as the linear regressions of the cows' EBV for each trait on their birth dates. The vector of expected genetic changes over 10 years of selection on PD19 (Φ) were computed using the following equation [23]: where i = the selection intensity, b = the vector of breeding index coefficients for PD19, G = the genetic variance matrix among the 9 traits, and P = the phenotypic variance matrix, computed as the sum of the genetic and residual variance matrices. The G and P matrices were derived from the REML analysis of data set 1 for cows born between 2008 and 2016. The index coefficients for the Israeli breeding index are given in Table 1. The expected economic gain for PD19, TEG, was computed as follows: with all terms as defined previously. The realized gain for PD19 was computed in the same manner with Φ replaced by the vector of realized gains, as derived from the analysis of data set 5. The selection intensity in equation [2] is a function of the selection intensities along the four paths of selection and is only known approximately. Therefore, i was set to 3.02 so that the expected economic gain for PD19 should be equal to the realized gain.

Genomic Analysis
The genomic analysis included all Israeli Holstein bulls with genotypes born from 1991 with reliabilities >0.5 from the analyses of data sets 3 and 4. Of the 1749 Israeli Holstein bulls genotyped, 1663 were born since 1991 and had genetic evaluations for protein with reliabilities >0.5, and 1610 had genetic evaluations for CS. As genotyping of these bulls was performed using several SNP chip platforms, we included only those markers that were genotyped in >90% of the tested cohort. A total of 40,498 SNPs were retained. GWAS files were prepared and formatted as described [24,25] using the plink software [26]. The allele substitution effects and the nominal probabilities for the hypothesis of no effect were computed using EMMAX software [27]. To account for the effects of relationships among individuals, we generated a pseudo-relationship matrix based on the identity by state matrix calculated using the emmax-kin-intel64 algorithm and the -v -s -d 10 flags. The GWAS was computed using the EMMAX algorithm with the -v -d 10 -t flags, and the relationship matrix was included in the analysis using the -k argument. Experimentwise probabilities accounting for multiple testing were computed based on the Bonferroni correction. To assess the variance explained by all SNPs, we used the kinship matrix and the bull phenotypes list and calculated the EMMAX software genomic REML that provides pseudo-heritability estimates [27].

Results
First parity phenotypic annual means for CS for cows born from 1983 through 2018 and annual means of EBV for CS of cows born from 1981 through 2018 are given in Figure 1. In correspondence with the results presented by [10], annual phenotypic means for CS decreased from 55.6% for cows born in 1983 to 46.5% for cows born in 2018, while mean EBV for CS increased from −1.9 for cow born in 1981 to −0.1 for cows born in 2018, the last year with nearly complete data. As noted previously, the genetic base was set to the mean of cows born in 2015. CS was added to the Israeli breeding index in 2000 (Table 1). Mean annual EBV increased by 2.9% from 2000 to 2018. Annual genetic means for the EBV of the three milk production traits are given in Figure 2. The slopes for fat and protein EBV were nearly linear during the entire period from 1982 through 2018. The slope for milk is higher in the early years and declines after 1992. This corresponds to the major changes in the breeding index in 1990 and 1991 (Table 1). Genetic gains were nearly equal for fat and protein, despite the fact that the index coefficient for protein was more than double the coefficient for fat since 1991. This occurs because both the heritability and genetic variance for fat are greater than for protein. Genetic parameters derived from the REML analysis of data set 1 are given in Table 3. Heritabilities were generally highest for milk, followed by fat and protein. Heritability for CS was between 0.04 and 0.05 for the three time periods. For all traits, heritabilities increased or stayed the same for the later time periods, despite the "conventional wisdom" that long-term selection should reduce heritability. Genetic correlations between the 3 milk production traits and CS were all negative, and the absolute value was highest for the correlation between CS and protein, between −0.37 and −0.38. Similar to the heritabilities, genetic correlations did not decrease over time. Environmental correlations between the three milk production traits and CS were very close to zero, despite previous studies that indicate a negative phenotypic relationship between milk production traits and female fertility [9]. Genetic parameters for protein and CS ± standard errors in the first three parities as derived from the REML analysis of data set 2 are given in Table 4. Since this data set was analyzed by the AIREMLf90 program, approximate standard errors were computed. All standard errors of the heritabilities and genetic and environmental correlations were <0.05. First parity heritabilities for protein and CS was nearly equal to the first parity heritabilities in Table 3 for the most recent time period. For both traits, heritabilities decreased with increase in parity, but proportionally less for CS. Genetic correlations among parities were all <0.65 for protein, and >0.89 for CS. All genetic correlations between protein and CS were negative, but the absolute value of −0.34 was highest for the first parity. This absolute value is slightly lower than the absolute value of −0.38 in the analysis of data set 1, but not significantly different. The absolute values of genetic correlations between protein and CS decrease with the increase of parity to −0.12 for the third parity. Environmental correlations between protein and CS for the three parities were positive, but all were <0.07, as compared to nearly zero in data set 1. Genetic and environmental variance components from the analysis of data set 1 for cows born from 2008-2016 are given in Table 5. These values are similar to previous analyses of the Israeli dairy cattle population [28]. Genetic standard deviations (SD), fraction of the Israeli breeding index, PD19 as computed by equation [1], the expected genetic trends, as computed by equation [2], and realized genetic trends, as computed by the regression of EBV on the cows' birth dates for cows born from 2008 through 2016, are given in Table 6. The genetic SD and the genetic and phenotypic variances in equation [2] were derived from the genetic and environmental variances in Table 5. The difference between the realized and expected trends, divided by the genetic SD of each trait, is also given. For the traits in which negative values are economically favorable, the sign of realized-expected is reversed so that a positive value indicates that the realized genetic gain in the desired direction was greater than the expected gain. Realized trends were greater than expected for SCS, CS and herd-life, with the largest difference for herd-life. This is not too surprising, since natural selection also favors increased herd-life and fertility, and lower mastitis. Farmers also prefer high fertility bulls. In addition, the expected trends were based only on first parity records, while the realized trends were based on all parities up to the fifth. As shown in Table 4, heritabilities decreased with an increase in parity proportionately more for protein, as compared to CS, and genetic correlations among parities were higher for CS. Both factors should result in increased realized genetic gain for CS, as compared to protein. The absolute values of the discrepancy between the realized and expected genetic trends were lowest for protein, the trait with the largest fraction of the index, and the highest for milk. The large discrepancy for milk is probably due to the fact that farmers prefer bulls with high fat and protein concentrations, even though the concentration traits are not included in the index.
Production quotas are in fluid milk, while payment is for fat and protein. Thus, farmers can increase total income using high concentration bulls. The EMMAX G-REML results for CS and protein were 0.63 and 0.87. That is, the cumulative effects of all the markers explained 63% and 87% of the variance among the sire evaluations for these traits. After the Bonferroni correction for multiple comparisons, none of the markers for CS were significant at the 5% experiment-wise level, and only three markers met this criterion for protein.
The scatter plot of the allele substitution effects of 31,744 SNPs with minor allele frequencies (MAF) >0.10 on protein as a function of their effects on CS is given in Figure 3. As can be seen on the level of all markers, there was virtually no relationship between the marker effects on the two traits. The coefficient of determination was <10 −4 . The scatter plot of the 10 markers with the lowest probability values for protein as a function of the effects of these markers on CS is given in Figure 4. For these markers the coefficient of determination was 0.01, that is a correlation of~0.1, and the slope was not significantly different from zero. Thus, the overall genetic correlation between these traits was not reflected by the substitution effects of the individual markers.

Discussion
"Cow conception rate" in the US dairy cattle population is defined as the percentage of inseminated cows that become pregnant at each service, and "daughter pregnancy rate" is defined as the percentage of nonpregnant cows that become pregnant during each 21-day period [29]. Heritabilities for the daughter pregnancy rate and cow conception rate were estimated at 0.04 and 0.02 in 2014, and the genetic correlation between the two traits was estimated as 0.87. Genetic correlations with protein production were −0.18 and −0.15 [30], but details of the analysis model were not given. CS is apparently closer to the cow conception rate and has similar heritability. The number of services, pregnancy/conception to first service, and pregnancy within a given period were included among the fertility traits considered in the meta-analysis of [11]. Mean heritabilities for these traits among the studies analyzed ranged from 0.02 to 0.03. Thus, heritability of CS in the current study is at the high end for traits that estimate the probability of cow conception. The mean of genetic correlations of protein production with the three traits considered by [11] were 0.35, −0.37 and −0.17. Thus, correlations of protein with number of services (which was positive, because a high value indicates low fertility) and pregnancy/conception to first service were very similar to the correlation between protein and CS in first parity, despite the fact that CS should be less biased by the culling of cows with low production [13,14].
There is an apparent contradiction between the major phenotypic reduction in CS over the last 30 years, and the fact that the environmental correlations between CS and the milk production traits were all very low. It should be noted, however, that the HYS effect was considered fixed in the variance component analyses. Thus, changes in management factors detrimental to fertility would be excluded from the environmental correlations.
The genetic trend for CS in the Israeli Holstein population became positive in 2001, the year after this trait was included in the Israeli breeding index. Genetic trends in the UK and Ireland for calving interval were negative until 2004 and then positive since then, corresponding to the changes in the selection indices of these countries [11]. Genetic evaluations for the DPR of US Holsteins first became available to the industry in 2003. The genetic trend for DPR was negative until 2009 and has since been positive [31]. Thus, it is clear from the results for all four countries that positive selection for female fertility is possible without a major decrease in selection for production traits.
In the current study, changes in heritability and genetic correlations were investigated for close to 30 years, or approximately six generations. In long-term selection experiments, the selection response usually ends after 20 to 30 generations [32], although, in some cases, a significant response has continued for over 100 generations [33]. Despite these considerations, the analysis of several studies indicates that heritability of lactation milk yield in dairy cattle has risen from~25% in the 1950s to~35% in the past decade, although this may be largely due to improved management [34]. To the best of our knowledge, this is the first study that attempted to estimate changes in heritabilities and genetic correlations over time for economic traits in dairy cattle in a specific population. The slight increase in heritabilities corresponds to the results of [34]. Thus, as noted previously by [3], there is no reason to assume that these traits are approaching a selection plateau.
Current methods applied for genomic evaluation in nearly all populations are based on the evaluation of bull EBV separately for each trait [4]. These methods cannot be readily adapted to multi-trait genomic evaluation. The method of [35] is based on the analysis of the phenotypic records with the inclusion of the genomic relationship matrix, as derived from the marker genotypes. It should be theoretically possible to apply this method to multi-trait genomic evaluation. Results of the GWAS demonstrate that the overall negative genetic correlation between protein production and female fertility is not reflected in the effects of the individual markers.

Conclusions
Annual means for CS decreased from 55.6 for cows born in 1983 to 46.5 for cows born in 2018, while mean estimated breeding values increased from −1.9 for cows born in 1981 to −0.1 for cows born in 2018. Heritability for CS increased from 0.04 to 0.05 from the first through the third time periods. For all traits, heritabilities increased or stayed the same for the later time periods. The first parity genetic correlation between conception status and protein was −0.38. Environmental correlations between CS and the milk production traits were all close to zero. Heritabilities decreased with an increase in parity for protein but remained the same for CS. Realized trends were greater than expected for cows born from 2008 through 2016 for SCS, CS and herd-life, and were lower than expected for the production traits. No correlation was found for the effects on protein and conception status of the 31,744 SNPs with minor allele frequency >0.1 included in the GWAS analysis, or for the 10 markers with the lowest probability values for the substitution effect on protein.
Thus, there is no apparent reason to change the current procedures that compute genomic evaluations separately for each trait.