Multibreed genomic evaluation for production traits of dairy cattle in the United States using single-step genomic best linear unbiased predictor

Official multibreed genomic evaluations for dairy cattle in the United States are based on multibreed BLUP evaluation followed by single-breed estimation of SNP effects. Single-step genomic BLUP (ssGBLUP) allows the straight computation of genomic (G)EBV in a multibreed context. This work aimed to develop ssGBLUP multibreed genomic predictions for US dairy cattle using the algorithm for proven and young (APY) to compute the inverse of the genomic relationship matrix. Only purebred Ayrshire (AY), Brown Swiss (BS), Guernsey (GU), Holstein (HO), and Jersey (JE) animals were considered. A 3-trait model with milk (MY), fat (FY), and protein (PY) yields was applied using about 45 million phenotypes recorded from January 2000 to June 2020. The whole data set included about 29.5 million animals, of which almost 4 million were genotyped. All the effects in the model were breed specific, and breed was also considered as fixed unknown parent groups. Evaluations were done for (1) each single breed separately (single); (2) HO and JE together (HO_JE); (3) AY, BS, and GU together (AY_BS_GU); (4) all the 5 breeds together (5_BREEDS). Initially, 15k core animals were used in APY for AY_BS_GU and 5_BREEDS, but larger core sets with more animals from the least represented breeds were also tested. The HO_JE evaluation had a fixed set of 30k core animals, with an equal representation of the 2 breeds, whereas HO and JE single-breed analysis involved 15k core animals. Validation for cows was based on correlations between adjusted phenotypes and (G)EBV, whereas for bulls on the regression of daughter yield deviations on (G)EBV. Because breed was correctly considered in the model, BLUP results for single and multibreed analyses were the same. Under ssGBLUP, predictability and reliability for AY, BS, and GU were on average 7% and 2% lower in 5_BREEDS compared with single-breed evaluations, respectively. However, validation parameters for these 3 breeds became better than in the single-breed evaluations when 45k animals were included in the core set for 5_BREEDS. Evaluations for Holsteins were more stable across scenarios because of the greatest number of genotyped animals and amount of data. Combining AY, BS, and GU into one evaluation resulted in predictions similar to the ones from single breed, especially when using about 30k core animals in APY. The results showed that single-step large-scale multibreed evaluations are computationally feasible, but fine tuning is needed to avoid a reduction in reliability when numerically dominant breeds are combined. Having evaluations for AY, BS, and GU separated from HO and JE may reduce inflation of GEBV for the first 3 breeds.


INTRODUCTION
Official multibreed genomic evaluations for dairy cattle in the United States are based on multibreed BLUP evaluations followed by single-breed estimation of SNP effects.This is the so-called multistep method (VanRaden, 2008;VanRaden et al., 2009), where genomic PTA is obtained as the combination of direct genomic value based on markers and parent average.The effectiveness of genomic evaluations relies on the availability of a large reference population where the marker effects can be estimated (Goddard, 2009).The accuracy of this estimation strongly depends on the size of the reference population.VanRaden et al. (2009) demonstrated that, when the genotyping was performed almost exclusively on elite bulls, there was a linear increase in estimation reliability with the number Multibreed genomic evaluation for production traits of dairy cattle in the United States using single-step genomic best linear unbiased predictor A. Cesarani,1 * D. Lourenco, 1 S. Tsuruta, 1 A. Legarra,2 E. L. Nicolazzi, 3 P. M. VanRaden, 4 and I. Misztal 1  of bulls added to the reference population.Nowadays, the availability of genotypes drastically increased and the number of animals with genomic information recently hit 5 million only in the United States (https: / / queries .uscdcb.com/Genotype/ counts .html).
Single-step genomic BLUP (ssGBLUP) has been adopted in lieu of the multistep method for genetic evaluations in several farm animal species over the years.The ssGBLUP involves the inverse of a realized relationship matrix (H −1 ; Aguilar et al., 2010), which combines the pedigree relationship matrix (A) and the genomic relationship matrix (G).Despite past limits due to problems with unknown parent groups (UPG) and computational cost, recent studies demonstrated the validity of this method for evaluating breeding values in several different livestock species including dairy (Himmelbauer et al., 2021;Liu and Alkhoder, 2021;Pimentel et al., 2021) and beef cattle (Lourenco et al., 2015), buffalo (Aspilcueta-Borquis et al., 2015;Cesarani et al., 2021a), goats (Teissier et al., 2018), and sheep (Cesarani et al., 2019;Macedo et al., 2020).Recently, better prediction features of ssGBLUP compared with BLUP in US Holstein were reported (Cesarani et al., 2021b).Those authors used about 860k genotyped animals and demonstrated that old phenotypic records (i.e., before 2000) can be removed without reducing prediction accuracy of young selection candidates.Moreover, they concluded that the best way to consider UPG in ssGBLUP is to include them in both A and the relationship matrix among genotyped animals (A 22 ), as also described by Tsuruta et al. (2019) and Masuda et al. (2021).Another advantage of ssGBLUP is the consideration of genomic preselection (Masuda et al., 2018), that is, animals selected based on their GEBV before they have phenotypes available (Patry and Ducrocq, 2011;Jibrila et al., 2020).
Some studies have already investigated multibreed genomic evaluation models.Winkelman et al. (2015) applied the multibreed evaluation, using different blended or hybrid genomic models, to Holstein (HO), Jersey (JE), and their crosses from New Zealand; these authors reported high reliabilities for pure-and crossbreds, with the latter showing higher accuracy for some traits.Different types of SNP panels and genomic models (i.e., GBLUP and Bayesian methods) were tested in crossbreeds, and the highest accuracy for crossbreeds was achieved when their data were considered (Khansefid et al., 2020).Creating a genomic relationship matrix among all genotyped animals belonging to different breeds in ssGBLUP can be challenging; however, it may help overcome the limitation due to the small reference population size in some breeds.This research aimed to develop ssGBLUP multibreed genomic predictions for US dairy cattle using purebred Ayrshire (AY), Brown Swiss (BS), Guernsey (GU), HO, and JE animals.

Data
No animals were used in this study, and ethical approval for the use of animals was thus deemed unnecessary.
Data used in the official multibreed genomic evaluations for US dairy cattle breeds were provided by the Council on Dairy Cattle Breeding (Bowie, MD).Only purebred AY, BS, GU, HO, and JE animals were considered.It is important to note that animals are considered purebred if they have more than 90% pedigree-based ancestry from one breed.The data sets were analyzed in 4 different scenarios: (1) each single breed separately (single); (2) HO and JE together (HO_JE); (3) AY, BS, and GU together (AY_BS_GU); and (4) all the 5 breeds together (5_BREEDS).Number of phenotypic lactation records (and cows), genotyped animals, and total animals in the pedigree (traced back 3 generations) are shown in Table 1.Milk (MY), fat (FY), and protein (PY) yields recorded from January 2000 to June 2020 were considered.The 305-d yields for the first 5 lactations included projected records for the final lactation and for lactations not yet completed by June 2020.The projected records were expanded to have the same genetic variance as completed records but given less weight to account for their larger error variance (VanRaden, 1997).According to the real lactation length, different weights are applied to project the records to 305 d: the shorter the lactation, the lower the weight.Animals were genotyped with 48 different arrays ranging from <3,000 to >600,000 usable SNPs.Genotypes were imputed, within each breed, to a common set of 79,294 selected SNPs using Findhap version 3 (VanRaden, 2016).
Single-step GBLUP was used in all the scenarios.The algorithm for proven and young (APY; Misztal et al., 2014) was used when the number of genotyped animals exceeded 50k (i.e., in single for Holstein and Jersey evaluations with 15k random core animals, in HO_JE with 30k random core following an equal representation of the 2 breeds, in AY_BS_GU, and in 5_ BREEDS with 15k random core animals).According to Pocrnic et al. (2016), the number of largest eigenvalues explaining 98% of the variance of G (EIGEN98), which represents the number of core animals in APY, in the US cattle populations varies from 11k to 19k; however, this number can be higher in multibreed populations if there is less overlap among breeds (Pocrnic et al., 2019).

Cesarani et al.: MULTIBREED GENOMIC EVALUATION IN US DAIRY CATTLE
Because the EIGEN98 when combining AY, BS, and GU was 19,184, being 14,909 for BS alone, we tested one sub-scenario under AY_BS_GU with 15k BS and all AY and GU (AY_BS_GU_30k).Additionally, we compared AY_BS_GU and AY_BS_GU_30k against one analysis without APY, meaning G with almost 61k animals was directly inverted (AY_BS_GU_direct).The 5_BREEDS scenario had a sub-scenario where 45k core animals were used, being composed by 5k AY, 5k BS, 5k GU, 15k HO, and 15k JE (5_BREEDS_45k).

Model and Analysis
Two evaluation methods were considered: (1) traditional BLUP (BLUP) with UPG; (2) ssGBLUP with UPG for A and A 22 , implemented as in Tsuruta et al. (2019) and Cesarani et al. (2021b).The UPG were defined based on breed, sex, and year of birth.The latter consisted of different groups for animals born before 2000, from 2001 to 2005, from 2006 to 2010, and from 2011 to 2020.Therefore, a total of 8 UPG were assigned within each breed.
A 3-trait repeatability BLUP animal model was applied in the single scenario according to Cesarani et al. (2021b), whereas the model used in the present study for ssGBLUP is described as SS_UPG2 in the same paper.In particular, the following 3-trait repeatability animal model was applied for the traditional BLUP: where y was the vector of the phenotypic records (i.e., milk, fat, and protein lactation yields); b was the vector of the considered fixed effects of herd-management, age-parity, inbreeding coefficient covariates, and heterosis; h was the vector of the random effect of the herdsire combination; g was the vector of the fixed effect of UPG; a was the vector of animal additive genetic effect; p was the vector of permanent environmental effect; and e was the vector or residuals.Heterosis is included because purebred animals can have up to 10% of other breeds (VanRaden et al., 2007).A single heterosis regression coefficient per animal was computed, as 1 − Σs i d i , where s i and d i are fractions of the sire's and dam's genes from breed i (VanRaden, 1992); this results in 0 for a 100% purebred animal, but is different from 0 if, say, the bull was 100% Holstein and the dam was 95% Holstein and 5% Jersey.X was the incidence matrix relating each phenotypic record to the fixed effects in b; Q a was an incidence matrix relating animals in vector a to UPG in vector g a ; Z h , Z a , and Z p were incidence matrices relating phenotypic records in y to herd-sire, animal, and permanent environment effects, respectively.The h vector had mean zero and variance I ⊗ V h , where I was an identity matrix and V h was a 3 × 3 matrix of herd-sire variances and covariances, and ⊗ was the Kronecker product.The vector a had mean zero and variance A ⊗ V a , where A was the numerator relationship matrix and V a was a 3 × 3 matrix of additive genetic variances and covariances.Finally, the vector p had mean zero and variance I ⊗ V p , where I was an identity matrix and V p was a 3 × 3 matrix of permanent environment variances and covariances.The ssGBLUP model had the same effect structure of BLUP but the covariance structure for Z a a and Z a Q a g a also contained genomic relationships.In particular, the covariance matrix for the additive genetic effect in ssG-BLUP was given by , where H * was the inverse of the realized relationship matrix with UPG added to the pedigree relationship matrix A and to the matrix of pedigree relationships among genotyped animals A 22 ; A * was the inverse of the pedigree relationship matrix with UPG, modified with the QP transformation (Quaas, 1988); G APY −1 was the inverse of the genomic relationship matrix (G) built with or without the APY according to the scenario (see above).The G matrix was constructed according to the method 1 of VanRaden, 2008 (i.e., centered and scaled by current allele frequencies from all genotyped animals).The genomic matrix was combined with 5% of A 22 to avoid singularity; the latter matrix was built based on the algorithm proposed by Colleau (2002).Moreover, G was scaled based on A 22 so the 2 matrices had same average values for both diagonal and off-diagonal elements.
Variance components for MY, FY, and PY were obtained from VanRaden et al. ( 2014): heritability was 0.20 for all the 3 traits; genetic correlations were 0.45 (MY vs. FY), 0.81 (MY vs. PY), and 0.60 (FY vs. PY).The heterogeneous herd variance was adjusted according to Wiggans and VanRaden (1991).The pre-adjustments for herd-year variance also adjust all breeds to have equal genetic variance (VanRaden et al., 2007).Higher heritability in different breeds is accounted for in the model by higher weights obtained from the ratios of genetic to error variance.The same pre-adjustments and variance ratios are applied to all the 3 traits.
In AY_BS_GU, HO_JE, and 5_BREEDS scenarios, fixed effects were made breed-specific by concatenating a breed indicator to the effect levels and nesting covariables for inbreeding and heterosis within breed.
The mixed model equations (MME) were solved using iteration on data and a block preconditioning conjugate gradient with the software BLUP90IOD2 (ver. 3.113;Tsuruta et al., 2001;Tsuruta and Misztal, 2008).If the MME are Ax = b, the convergence criterion was b − Ax 2 /b 2 , with a termination criterion of 10 −15 .The computations were carried out on a Linux server (x86_64) equipped with Intel Xeon E5-2683 v4 2.10 GHz processor with 32 cores.

Validation of Breeding Values
The (G)EBV were validated within each tested scenario.Two data sets were considered: (1) complete, with phenotypes for milk, fat, and protein recorded in cows born from 1992 to 2018; (2) reduced, with phenotypic records from cows born up to 2014.Genotyped cows born between 2015 and 2018 with phenotypes in the complete but not in the reduced data set were in the validation set.The numbers of validation cows were 181 for AY, 2,423 for BS, 750 for GU, 577,340 for HO, and 90,666 for JE.Genotyped bulls with no daughters with phenotypes in the reduced but at least 10 daugh-ters in the complete data set were in the validation group for AY (n = 17), BS (n = 107), and GU (n = 28), whereas 50 daughters were requested for HO (n = 3,278) and JE (n = 471) validation bulls.
In this study, theoretical accuracies were not analyzed, and models were compared based on validation reliabilities.The validation for bulls was carried out using daughter yield deviations (DYD), estimated in the complete data sets for all bulls following the method by Liu et al. (2004) and the algorithm by Mrode and Swanson (2004).The validation reliability of predictions for bulls was measured using the coefficient of determination (R 2 ) from the regression of DYD on (G) EBV.The validation for cows was based on the estimates of phenotypes adjusted for fixed effects (Y adj ), computed based on the complete data set using PRE-DICTF90 (ver.1.5; Misztal et al., 2018).Note that the R 2 cannot be, strictly speaking, interpreted as a reliability (squared correlation between true and estimated BV) unless the number of daughters in the DYD is very high.The predictive ability [i.e., the correlation between Y adj and (G)EBV] was used as a measure of accuracy of predictions for validation cows.Note that this correlation cannot be interpreted as an accuracy (correlation between true and estimated BV).Finally, the inflation of predictions was measured through the regression coefficients (b 1 ) of the regression of DYD or Y PRED on (G)EBV.
There are several reasons to use one criterion for bulls and another for cows.First, in this manner, we can compare each of these criteria with preceding studies for bulls (e.g., VanRaden et al., 2009) and for cows (e.g., Bengtsson et al., 2020;VanRaden et al., 2020).Second, measures of goodness of fit for bulls and cows are not comparable, as bulls have been strongly selected, which results in a decrease in cross-validation accuracy (Bijma, 2012).Third, it is not easy to develop a simple statistic that can be used for bulls and cows simultaneously.

RESULTS AND DISCUSSION
Table 2 reports the comparison between (G)EBV estimated in single and 5_BREEDS_45k scenarios.Since animals were not connected by the pedigree and the breed was correctly accounted for by the model, BLUP results for single and all multibreed scenarios were the same: correlations between raw EBVs estimated in BLUP single and BLUP multibreed were 1.00 for all breeds and traits.In the genomic models, correlations between raw GEBVs estimated in ssGBLUP single and ssGBLUP 5_BREEDS_45k ranged from 0.95 (MY for AY and BS) to 1.The latter was the case for all correlations involving HO and JE, probably because of the greatest number of genotypes.Moreover, we also computed the slope of the regression of (G)EBV in 5_BREEDS on (G)EBV in single models.As expected, slopes for BLUP were always 1, confirming that breed was correctly considered, and animals were not connected by the pedigree.Acceptable slope values were also found for ssGBLUP as values ranged from 0.93 (MY and PY for BS) to 1.02 (MY for AY).
The coefficients of determination (R2 ) for bulls in single-breed and multibreed analyses are shown in Table 3.As expected, the highest values for BLUP were observed for the numerically dominant breeds (i.e., Jersey and Holstein).All the 7 tested genomic scenarios had greater values compared with BLUP; however, some differences were observed among the genomic scenarios.Within the single scenario, ssGBLUP showed values that were on average (±SD) 0.23 (±0.14) points greater than the ones from BLUP in all considered breeds.The greatest increases moving from BLUP to ssGBLUP were observed for the breeds with larger numbers of genotyped animals (i.e., BS, JE, and HO).Looking at the 5_BREEDS scenario, there was a considerable reduction in R 2 for BS compared with the single scenario.Contrasting, there was an overall increase of 16 points in R 2 for AY when moving from single to 5_BREEDS.A small drop in R 2 was observed for fat and protein in GU.The validation results for AY, BS, and GU should be taken with caution because of the limited numbers of validation bulls, which were 17, 107, and 28, respectively.When the number of validation bulls was large enough (i.e., for HO and JE), R 2 was similar between single and 5_BREEDS scenarios.A successful transition from single-breed to multibreed evaluations would produce R 2 and dispersion in the latter that are at least as good as in the former.
A factor contributing to the drop in R 2 in our study was the breed representation in the core group.When randomly selecting 15k animals among the 5 breeds to be in the core set, the breeds with less genotyped animals were not well represented.While there were 13,049 Holsteins in the core, there were only 1,720 JE, 182 BS, 32 AY, and 17 GU.After increasing the core size to 45k (5_BREEDS_45k), which included 15k HO, 15k JE, 5k AY, 5k BS, and 5k GU, R 2 increased to similar or larger levels as in the single scenario.
Multibreed analyses for the 3 breeds with fewer data were performed to investigate whether the genomic information from HO and JE was responsible for the loss in prediction accuracy and validation reliability.Issues in multibreed evaluations are common when some breeds have many more genotyped animals than others.In such a case, SNP effects and frequencies are driven by the numerically dominant breeds, which causes a drop in validation reliability and an increase in inflation of GEBV for smaller breeds (van den Berg et al., 2020).Validation reliabilities for the AY_BS_GU scenario were not homogeneous across breeds when comparing with the single scenario.We observed larger R 2 for AY, slightly larger for BS, and slightly lower for GU.Comparing to AY_BS_GU_30k, which had 15k BS and all AY and GU in the core, results were similar for BS and GU.When the AY_BS_GU evaluation involved the direct inverse of G (AY_BS_GU_direct), the average R 2 across traits remained 0.40 for BS and 0.26 for GU.The large oscillation in R 2 for AY confirms the number of validation bulls is too small to draw any conclusions.
Having a multibreed evaluation with combined data from the breeds with fewer data is a possibility; therefore, investigating an evaluation with data from the numerically dominant breeds becomes interesting.The evaluation combining HO and JE produced a similar R 2 for HO bulls across all traits and a slight increase in R 2 for milk in JE when compared with 5_BREEDS_45k.Small changes of 0.01 or 0.02 for some traits when changing the model to account for different breeds can be because of the scaling of G to match A 22 .In such a case, diagonals and off-diagonals of G are set to have the same average as in A 22 (Chen et al., 2011;Vitezica et al., 2011), which may be different by breed; however, the statistics for A 22 are computed across all the breeds included in the model.One way to possibly solve this issue is to scale G differently based on A 22 constructed separately for each breed.Another way would be to use metafounders (Legarra et al., 2015), which scales A 22 to match G and sets pseudoindividuals as proxies for the base animals.
Predictive abilities (i.e., validation reliability for cows) for validation cows in single-breed and multibreed analyses are shown in Table 4.The number of validation cows was large enough for most breeds, so the results for BS and GU were more stable for cows than bulls.Ayrshire still had a small number of validation cows (n = 181), which resulted in unexpected predictivities for some scenarios and a small gain of 0.02 (milk and protein) and 0.03 (fat) when switching from BLUP to single-breed ssGBLUP.Overall, the gain in predictivity by using genomic single-breed models was 0.16 for BS, 0.09 for GU, 0.21 for HO, and 0.17 for JE.When moving from single to 5_BREEDS scenarios, the average predictivity for HO and JE was similar, whereas the average predictivity for AY, BS, and GU cows decreased from 0.42 to 0.36.Again, this drop was mainly driven by the number and composition of the core animals in 5_BREEDS, which did not have a good representation of AY, BS, and GU.Using 45k core animals in APY for the multibreed evaluation (5_BREEDS_45k) increased the average predictivity for AY, BS, and GU cows to 0.45.Predictivity for HO and JE did not change.When evaluations for AY, BS, and GU were done separately from HO and JE (i.e., AY_BS_GU, AY_BS_GU_30k, and AY_BS_GU_direct), predictivities were overall similar to the single scenario.Likewise, predictivities for HO and JE did not considerably change from single to HO_JE scenarios.However, small fluctuations were observed across traits and breeds.
One of the concerns with multibreed evaluations is the increase in inflation of GEBV by combining breeds that have been selected differently.In this study, when breeds were combined in the same evaluation, we created breed-specific fixed effects and used UPG based on breed in addition to year of birth interval and sex.Thus, levels of effects were defined within breed (i.e., the effect of age by parity number 1 was assumed to be different in different breeds).Whereas UPG take care of the genetic differences among breeds, the breed-specific fixed effects account for other nongenetic variation, such as different effects of parity number or different levels of average within breed heterosis.An initial investigation showed that not considering breed-specific fixed effects in the model resulted in lower validation reliability and predictivity and increased inflation of EBV and GEBV (results not shown).
The inflation and deflation of (G)EBV were assessed in this study by investigating the slope (b 1 ) of the regression of DYD on (G)EBV for bulls (Table 5) and of adjusted phenotypes on (G)EBV for cows (Table 6).According to Tsuruta et al. (2011), estimates within 15% from 1 are acceptable, whereas estimates within 5% are considered good predictions.Across traits, EBV for bulls were inflated for most of the breeds.Average b 1 for BS, GU, and HO was 0.63, 0.71, and 0.88.A deflation happened for milk and protein in JE.Adding genomic information in single-breed evaluations helped to reduce inflation and deflation.Combining the 5 breeds (5_BREEDS) increased the average inflation of BS GEBV compared with the single scenario, and small fluctuations were observed across breeds and traits.When using 45k core animals (5_BREEDS_45k), inflation and deflation levels were, in general, similar to 5_BREEDS.The only exception was AY that showed more inconsistencies across scenarios because of the small number of bulls.Overall, the AY_BS_GU scenarios resulted in b 1 that were similar to the single scenario for those breeds.The b 1 for HO did not vary much across 5_BREEDS and HO_JE scenarios compared with the single scenario, whereas small changes were observed for JE.Inflation and deflation patterns were very similar between bulls and cows, with b 1 values from BS and GU cows in AY_BS_GU closer to the ones in the single scenario, stable b 1 across scenarios for HO, and small variation for JE.The exception in this comparison between bulls and cows was AY that had less variation in the evaluations using 3 breeds compared with the single scenario.
Validation statistics for AY and GU in multibreed evaluations in the United States have not been reported before mainly because of the small number of validation bulls available and the late onset of genomic selection for those breeds.Whereas genomic selection was officially implemented for HO, JE, and BS in 2009, it only started in 2013 for AY and 2016 for GU.Olson et al. (2012) investigated the use of multibreed reference populations under multistep methods in the United States, but only genotypes for HO, JE, and BS were used.The authors concluded that predictivity across breeds was limited but breeds with fewer data slightly benefited from a multibreed estimation of SNP effects.who included 5 breeds and many more animals.However, the estimation of SNP effects is still in a singlebreed fashion for routine US dairy evaluations because of the small gains.This is true also for the evaluations in other countries (Su et al., 2010;Shabalina et al., 2020).
Although the estimation of SNP effects is single-breed for most of the dairy evaluations, the deregressed EBV used to estimate SNP effects in the multistep process are computed based on multibreed BLUP evaluations.VanRaden et al. (2007) developed an all-breed model for routine BLUP evaluations in the United States that considered AY, BS, GU, HO, JE, and Milking Shorthorn, in addition to crossbreds.The advantage of combining different breeds is to have updated breed differences, and a reason to include data for crossbred animals is to have extra, accurate information for both parents, which can increase their EBV reliability.
In our study, we used ssGBLUP, which is an extension of BLUP to include relationships based on genomic information.We showed that having multibreed genomic evaluations in a single run is possible because validation reliabilities, predictabilities, and regression coefficients were not considerably different from those in single-breed evaluations, especially for the breeds with enough validation animals.Multibreed evaluations should deliver GEBV with at least the same reliability and inflation-deflation level as the single-breed evaluations.Although we did not include data on crossbred animals in this study, several studies have shown the usefulness of ssGBLUP in multibreed and crossbred evaluations and that accuracy of predictions for crossbred animals greatly improves when their genotypes are in the evaluation (Lourenco et al., 2016;Xiang et al., 2016;Pocrnic et al., 2019).The current genomic evaluation for crossbred dairy in the United States is based on purebred direct genomic values (DGV) weighted by breed proportion (VanRaden et al., 2020), and the genomic predictions from this method are more accurate than parent average.Additional accuracy could possibly be added by having a multibreed and crossbred joint ssGBLUP evaluation, which is currently under investigation.
The use of ssGBLUP in a multibreed context, however, implies that SNP effects are the same for all the breeds in the data set.Similar to the other results, SNP estimates can also be more conditioned to the more numerically dominant breeds and different from those estimated in single-breed models.Because the number of genotypes is constantly increasing, DGV for the newly genotyped animals can be estimated by using indirect predictions instead of solving the whole MME (Lourenco et al., 2015).The DGV estimated using SNP effects from multibreed models will be different from those of single-breed models, but they could be more accurate.In our laboratory, studies on indirect predictions and a new algorithm to approximate their theoretical reliabilities are under investigation.
Because of the large number of genotyped animals, APY was used to compute G −1 for BS, HO, and JE as single breed, AY_BS_GU, AY_BS_GU_30k, all 5_ BREEDS scenarios, and HO_JE.Pocrnic et al. (2016) found the number of core animals in purebred cattle populations would vary from 10 to 19k; however, this may not apply to multibreed populations.In such a case, the number of core animals should represent the dimensionality within each breed, which implies using the number of core animals for each breed as the number of eigenvalues explaining 98% of the variance of G.This was confirmed by the fact that using 15k core animals randomly sampled among all the genotyped animals caused some fluctuations in validation reliability, predictivity, and b 1 for bulls and cows in the multi-compared with single-breed evaluations.Therefore, if the multibreed evaluation involves the 5 breeds, 45k core animals may be enough, whereas 30k is enough for AY_BS_GU and 30k for HO_JE.Adding crossbred genotypes to the evaluation may not increase the number of core animals because of the sharing of independent chromosome segments with the purebred parents (Pocrnic et al., 2019).Independent chromosome segments can also be shared among different breeds, especially if there was recent introgression from the use of a few sires across breeds.
When combining multiple breeds in a joint ssGBLUP evaluation with APY, it is important to have a good representation of all breeds in the core (Mäntysaari et al., 2017;Vandenplas et al., 2018).Additionally, modeling the breed effects in a proper way guarantees more accurate predictions in any multibreed evaluations.In this study we added breed code to all the fixed effects to ensure animals from different breeds were not mistakenly in the same group.Breed effects were also added to the model as a fixed UPG.Therefore, GEBV resulting from the multibreed models were adjusted for breed-specific fixed effects but contained genetic breed effects combined with selection pathway by sex and year of birth (i.e., the original definition of UPG).This model ensured that predictions from multibreed and single breed were more similar.When the model does not correctly account for the breed and breed-specific effects, extremely high reliabilities from multi-compared with single-breed models may be observed due to the prediction of breed effects instead of GEBV (Bermann et al., 2021).However, this is not the case in this study because all validations are done within breed.
Table 7 shows number of rounds and computing time (i.e., only the preconditioning conjugate gradient with iteration on data) to reach convergence assuming a criterion of 10 −15 in all the scenarios.The computing cost for AY, BS, and GU was small in the single scenario and in the 3 AY_BS_GU scenarios.For HO, it took 2.7 h for BLUP and 7.4 h for ssGBLUP (3.4 million genotyped animals) in the single scenario.Larger computing time was recently reported by Cesarani et al. (2021b), who analyzed Holstein data up to 2018 with 861k genotypes.The larger computing costs reported in that previous work could be because of the use of selected genotypes only for cows with records and their parents, in addition to more UPG.When the 5 breeds were combined, the computing time increased 1.8 times for BLUP, 2.7 times for 5_BREEDS, and 8.6 times for 5_BREEDS_45k compared with single HO.The largest number of iterations suggests a possible lack of compatibility between genomic and pedigree information after combining different breeds.Additionally, because scaling of G is based on an average A 22 across breeds, the combination of the 5 breeds possibly resulted in an H −1 that has a larger condition number.As far as the time needed for the creation of the relationship matrices, it took about 2.5 d to create the matrices for the 3.4 million HO, about 3.5 d for 5_BREEDS, and about 5 d for 5_BREEDS_45k.About half of this time was used in the creation of A 22 for the scaling of G.A new algorithm is under investigation to reduce the computing time required for the construction of A 22 .This algorithm considers the APY structure and sparsity, meaning that coefficients of A 22 are computed separately for core and noncore animals.

CONCLUSIONS
Results of the present study demonstrated that large-scale ssGBLUP multibreed evaluations are computationally feasible, and that genomic EBV are as reliable as in single-breed evaluation when the number of genotyped animals in the APY core represent the dimensionality of the genomic information within the breeds.Accounting for genetic and nongenetic breed effects is important when combining several breeds in a joint evaluation; for instance, fixed effects need to be redefined within breed.Fine tuning is still needed to avoid a reduction in reliability for less represented breeds when numerically dominant breeds are in the evaluation.Having evaluations for Ayrshire, Brown Swiss, and Guernsey separate from Holstein and Jerseys may reduce inflation of GEBV for the first 3 breeds.Hokkaido, Japan) are greatly acknowledged.The authors have not stated any conflicts of interest.
This study was partially funded by Agriculture and Food Research Initiative Competitive Grant no.2020-67015-31030 from the US Department of Agriculture's National Institute of Food and Agriculture.The authors thank the Council on Dairy Cattle Breeding (CDCB, Bowie, MD) for providing access to the data.The contribution of dairy producers who supplied data through their participation in the Dairy Herd Improvement program and the Dairy Records Processing Centers that edited and relayed information to the Council of Dairy Cattle Breeding are also acknowledged.Discussions with Duane Norman (CDCB, Bowie, MD) and Yutaka Masuda (Rakuno Gakuen University, Ebetsu, Cesarani et al.: MULTIBREED GENOMIC EVALUATION IN US DAIRY CATTLE

Table 1 .
Number of records, cows, genotypes, and total animals in 5 single breeds and in the combined scenarios 1 Cesarani et al.: MULTIBREED GENOMIC EVALUATION IN US DAIRY CATTLE

Table 2 .
Cesarani et al.: MULTIBREED GENOMIC EVALUATION IN US DAIRY CATTLE Comparison between genomic EBV estimated in single and 5_BREEDS_45k scenarios 1

Table 7 .
Computational costs for all considered scenarios in terms of number of rounds and seconds per round needed to reach a convergence criterion of 10 −15 1 AY = Ayrshire; BS = Brown Swiss; GU = Guernsey; HO = Holstein; JE = Jersey; 5_BREEDS = all 5 breeds together.2Time to reach the convergence (i.e., only time for the preconditioning conjugate gradient with iteration on data).ssGBLUP = single-step genomic BLUP.