Investigation of β‐hydroxybutyrate in early lactation of Simmental cows: Genetic parameters and genomic predictions

Genomic information allows for a more accurate calculation of relationships among animals than the pedigree information, leading to an increase in accuracy of breeding values. Here, we used pedigree‐based and single‐step genomic approaches to estimate variance components and breeding values for β‐hydroxybutyrate milk content (BHB). Additionally, we performed a genome‐wide association study (GWAS) to depict its genetic architecture. BHB concentrations within the first 90 days of lactation, estimated from milk medium infrared spectra, were available for 30,461 cows (70,984 records). Genotypes at 42,152 loci were available for 9,123 animals. Low heritabilities were found for BHB using pedigree‐based (0.09 ± 0.01) and genomic (0.10 ± 0.01) approaches. Genetic correlation between BHB and milk traits ranged from −0.27 ± 0.06 (BHB and protein percentage) to 0.13 ± 0.07 (BHB and fat‐to‐protein ratio) using pedigree and from −0.26 ± 0.05 (BHB and protein percentage) to 0.13 ± 0.06 (BHB and fat‐to‐protein ratio) using genomics. Breeding values were validated for 344 genotyped cows using linear regression method. The genomic EBV (GEBV) had greater accuracy (0.51 vs. 0.45) and regression coefficient (0.98 vs. 0.95) compared to EBV. The correlation between two subsequent evaluations, without and with phenotypes for validation cows, was 0.85 for GEBV and 0.82 for EBV. Predictive ability (correlation between (G)EBV and adjusted phenotypes) was greater when genomic information was used (0.38) than in the pedigree‐based approach (0.31). Validation statistics in the pairwise two‐trait models (milk yield, fat and protein percentage, urea, fat/protein ratio, lactose and logarithmic transformation of somatic cells count) were very similar to the ones highlighted for the single‐trait model. The GWAS allowed discovering four significant markers located on BTA20 (57.5–58.2 Mb), where the ANKH gene is mapped. This gene has been associated with lactose, alpha‐lactalbumin and BHB. Results of this study confirmed the usefulness of genomic information to provide more accurate variance components and breeding values, and important insights about the genomic determination of BHB milk content.

Keywords: cow validation; dairy cattle; health traits; ketosis; single‐step GBLUP

Breeding programmes in livestock populations have been traditionally based on BLUP that provide unbiased predictions when all available information is included in the analysis (Henderson, 1975). The development of genome‐wide dense marker maps generated a shift from the traditional pedigree BLUP to genomic methods (Hayes et al., 2009; Meuwissen et al., 2001). It is now well‐established that genomic information allows for a more accurate estimation of relationships among all genotyped animals, even if they are not related through the pedigree (Hayes et al., 2009). Single‐step genomic BLUP (ssGBLUP, Aguilar et al., 2010; Christensen & Lund, 2010) has become the approach of choice for genomic evaluation because it is simple to implement (Legarra et al., 2014), and it accounts for preselection (Patry & Ducrocq, 2011). Greater accuracy of breeding values from ssGBLUP, compared to genetic evaluation, was also reported in minor species, such as French dairy goat (Teissier et al., 2018), French Lacaune (Baloche et al., 2014) and Sarda (Cesarani, Gaspa, et al., 2019) dairy sheep, but no advantage was found for Latxa dairy sheep (Granado‐Tajada et al., 2020).

According to Garcia‐Ruiz et al. (2016), genomic selection is particularly efficient for traits with low heritability and difficult to measure, as in the case of many functional traits. An example is the resistance to ketosis. This metabolic disease occurs in cattle mainly during the transition period from late gestation to early lactation. The high energy demand for fetus growth and milk production (Ranaraja et al., 2018), together with a low feed intake, results in a negative energy balance (NEB). Consequently, a significant mobilization of adipose tissue occurs, with a subsequent abnormal increase of blood and milk concentrations of non‐esterified fatty acids (NEFA) and ketone bodies (hyperketonaemia). Among the latter, the β‐hydroxybutyrate (BHB) is the most stable circulating ketone body in milk, blood and other fluids (Duffield et al., 2009). This negative energy balance in the transition period is more severe in high‐yielding dairy cows (McCarthy et al., 2010), because the energy request for milk production is higher and so more time is needed to recover from NEB in the postpartum period when compared to cows with lower milk production. High‐yielding dairy cattle breeds have larger production requirements, and higher incidence of NEB, compared to breeds with lower milk production (i.e., dual‐purpose breeds). Moreover, high‐yielding cows show more intense, easier and faster mobilization of adipose tissues as a consequence of the selection for milk production and, especially, for the lactation peak. For the same body weight, cows belonging to specialized breeds have very higher lactation peaks, which are possible because of the higher mobilization of adipose tissue.

Hyperketonaemia has negative consequences on health and immune system and on milk production (McArt et al., 2013). Animals with ketosis also have a greater probability of developing other diseases such as metritis and mastitis (Suthar et al., 2013) and they have lower probability to get pregnant (Ospina et al., 2010). High BHB concentrations postpartum (>10 mg/dl) lead to a reduction in milk yield (Ospina et al., 2010), to an increase in milk fat because of the greater availability of BHB and fatty acids, and to a decrease in milk protein due to the reduced glucose availability for the rumen microflora (Ranaraja et al., 2018). For this reason, the milk of cows with hyperketonaemia is often characterized by a higher fat‐to‐protein ratio.

Hyperketonaemia is strongly related to ketosis, and BHB concentration in the blood is a useful parameter for the diagnosis of this metabolic disease. In fact, a BHB concentration higher than 1.2 mM is considered an indicator of a possible status of ketosis (McArt et al., 2012). Weigel et al., (2017) reported that values of BHB blood concentration between 1.2 and 2.9 mM and higher than 3.0 mM can be associated to subclinical and clinical ketosis, respectively. However, phenotyping large populations for this trait is rather complex and expensive (Nayeri et al., 2019). Milk and blood BHB concentrations are strongly correlated (>0.85; Denis‐Robichaud et al., 2014), but the former is cheaper to record. This study aimed to investigate the feasibility of implementing genomic evaluation for BHB (as ketosis indicator) in Italian Simmental cattle by using a single‐step genomic approach. Hence, the genetics of BHB milk concentration was investigated. In particular, variance components and genetic correlations between BHB and milk production traits were estimated to evaluate possible consequences of including BHB in the breeding scheme to improve ketosis resistance. Moreover, a genome‐wide association study was performed to identify genomic regions involved in the genetic determination of BHB milk concentration.

Animal Care and Use Committee approval was not needed as data were obtained from pre‐existing databases.

The data set provided by the Italian Simmental Association (Associazione Nazionale Allevatori Bovini di Razza Pezzata Rossa Italiana, ANAPRI) contained 70,984 BHB milk concentration phenotypes from 30,461 cows born from 2008 to 2018, which were recorded within the first 90 days of lactation. Test‐day phenotypes were recorded according to the official national recording scheme of ANAPRI, which is a mixed recording scheme AT5 and AT4, where test‐day records are collected every 4 or 5 weeks with alternated sampling during morning or afternoon milking (ICAR, https://www.icar.org/Guidelines/02‐Overview‐Cattle‐Milk‐Recording.pdf). BHB milk concentrations (mmol per litre) were estimated using Fourier transform infrared spectroscopy from the medium infrared milk spectra (Milkoscan, FOSS). The distribution of records and the average per year of birth for BHB concentration and the other traits considered in this study are reported in Table 1. On average, 2.33 records were available per cow.

1 TABLEDistribution of records for β‐hydroxybutyrate (mM) in milk and other phenotypes used in this study, across year of birth and descriptive statistics (average ± SD)

1 Note

2 Abbrewviations: BHB, β‐hydroxybutyrate in milk (mM); F:P, fat‐to‐protein ratio; FP, fat percentage; LC, lactose (%); MY, milk yield (kg/day); PP, protein percentage; SCSLOG, logarithmic transformation of somatic cells countl UR, urea (mg/100 ml).

A total of 9,123 animals were genotyped for 42,152 SNPs (mapped to the ARS‐UCD 1.2 assembly). Among those, 2,146 were females with BHB records and 6,977 were their relatives (4,250 males and 2,727 females). Genotypes were already preprocessed from ANAPRI. Animals were genotyped with different chips and SNPs were selected based on a list used for the official evaluation. Quality control was applied to the selected SNPs to discard the ones with minor allele frequency (MAF) <2%, call rate <97.5%, and deviation from the Hardy–Weinberg Equilibrium (HWE) (p > 0.01). Pedigree‐based imputation of missing SNPs was then carried out using the pedimpute software (Nicolazzi et al., 2012). Genotypic, phenotypic and pedigree information were provided by ANAPRI.

Data were analysed with the following repeatability animal model:

$$\mathbf{y}=\mathbf{Xb}+\mathbf{Zu}+\mathbf{Wp}+\mathbf{e},$$

where y is the vector of phenotypic records; b is the vector of the fixed effects of herd‐test‐day (17,400 levels), calving season (four levels), age at parity (24 levels) and days in milk (considered as both linear and quadratic covariates); X is the incidence matrix associating phenotypic records to fixed effects; u and p are the vectors of random direct additive genetic and permanent environmental effects, respectively; Z and W the incidence matrices relating animals to phenotypic records; and e is the vector of random residuals.

Two different approaches were used to estimate variance components and to predict breeding values: (a) the pedigree‐based approach, where the (co)variance structure of the random animal effect was modelled as $\mathbf{u}\sim \mathrm{N}(0,\mathbf{A}{\mathrm{\sigma }}_{\mathrm{a}}^2)$ where A is the pedigree relationship matrix and ${\sigma }_a^2$ is the direct additive genetic variance; (b) the single‐step genomic approach, where the random animal effect was modelled as $\mathbf{u}\sim \mathrm{N}(0,\mathbf{H}{\mathrm{\sigma }}_{\mathrm{a}}^2)$ where H is a matrix that combines pedigree and genomic relationships (Legarra et al., 2009). The vector of permanent environmental effect was modelled as $\mathbf{p}\sim \mathrm{N}(0,\mathbf{I}{\mathrm{\sigma }}_{\mathrm{PE}}^2)$ , where I is an identity matrix and ${\sigma }_{\mathrm{PE}}^2$ is the corresponding variance component.

The pedigree was traced back for three generations from animals with phenotypes and/or genotypes, leading to a total of 94,698 animals. BHB was analysed with a single‐trait linear mixed model without any previous transformation in its normal scale with a single‐trait model. Genetic correlations between BHB and milk production traits were estimated by running a series of bivariate models with BHB and one of the following traits: milk yield (MY), fat (FP) and protein percentage (PP), fat‐to‐protein ratio (F:P), urea (UR; mg/100 ml), lactose (LC; %), and logarithmic transformation of somatic cells count (SCSLOG) based on Ali and Shook (1980). Variance components, heritability (h2), repeatability (r2) and genetic correlations were estimated via the Gibbs sampler (GIBBS2F90; Misztal et al., 2014), sampling 50,000 rounds and storing every 5th sample. After discarding 5,000 samples as burn‐in, posterior means for all the parameters were calculated. The approach used to estimate variance components based on pedigree will be termed simply as GIBBS, whereas the one based on genomic information will be the single‐step genomic GIBBS (ssGGIBBS). Breeding values for the pedigree‐based approach were computed using BLUP and for the genomic approach were computed using ssGBLUP. Both used the BLUPF90 software suite (Misztal et al., 2014). In order to highlight only differences due to the approach, the same variance components (the ones estimated with GIBBS) were used to estimate breeding values with BLUP and ssGBLUP.

Genomic EBV from ssGBLUP was backsolved into SNP effects for the genome‐wide association study (GWAS) as described in Wang et al. (2012). Following Aguilar et al. (2019), p‐values were computed based on prediction error variance of SNP effects using POSTGSF90 (Misztal et al., 2014). Genes mapped (ARS‐UCD 1.2 cow genome assembly) within an interval of ±0.25 Mb from any significant SNPs were flagged as significant (Cesarani, Sechi, et al., 2019; Manca et al., 2020). Significant SNPs were the ones with p‐values smaller than a threshold based on a significance level of 0.05 with a Bonferroni correction for multiple testing.

Cows born from 2016 to 2018 with both genotypes and phenotypes (n = 344) were identified as validation animals. Phenotypes of these cows were removed from the whole data set to create a reduced data set. Adjusted phenotypes in the whole data set were computed using PREDICTF90 (Misztal et al., 2014).

The validation process was carried out using two strategies: (a) predictive ability (pred), which is the correlation between adjusted phenotypes in the whole data set and breeding values estimated in the reduced data set; (b) the linear regression (LR) method proposed by Legarra and Reverter (2018). In the latter strategy, BLUP and ssGBLUP evaluations were compared through the following statistics:

$$acc=\sqrt{cov({\widehat{\mathit{u}}}_w,{\widehat{\mathit{u}}}_r)/(1‐\overline{F}){\widehat{\sigma }}_a^2}$$

where $\mathit{acc}$ is the accuracy of (G)EBV, ${\widehat{\mathit{u}}}_w$ and ${\widehat{\mathit{u}}}_r$ are the (G)EBV of candidate animals in the whole and reduced data set, respectively; $\overline{F}$ is the average inbreeding coefficient for the validation animals. The correlation between the two sets of (G)EBV was computed to assess the consistency between two subsequent evaluations:

$${\rho }_{w,r}=\mathrm{cor}({\widehat{\mathit{u}}}_w,{\widehat{\mathit{u}}}_r)$$

Finally, the prediction dispersion (i.e., slope of the regression of ${\widehat{\mathit{u}}}_w$ on ${\widehat{\mathit{u}}}_r$ ) of the reduced estimates was computed as:

$$b_{w,r}=\mathrm{cov}({\widehat{\mathit{u}}}_w,{\widehat{\mathit{u}}}_r)/var\left({\widehat{\mathit{u}}}_r\right)$$

Genetic parameters for BHB and for the other studied traits, estimated using the single‐trait BLUP and ssGBLUP models, are reported in Table 2. The heritabilities for the other traits are in agreement with, or slightly lower than, those estimated in ANAPRI (personal information). The difference was mainly because only the first 90 days of lactation were considered in this study. The h2 estimates with genomic information (ssGGIBBS) were, in general, slightly higher than those from GIBBS, mainly because of larger additive genetic variance (Table 2). Moreover, ssGGIBBS provided lower standard errors for the variance components than GIBBS. Low heritability and moderate repeatability values were observed in both approaches for BHB. However, a small increase in additive genetic variance and a decrease in PE variance were observed for ssGGIBBS compared to GIBBS. Since residual variance was almost the same in the two approaches, the inclusion of genomic information resulted in a small shift of variance from PE to the additive genetic effect. Thus, the more accurate estimation of relationships in ssGBLUP allowed a slightly better distinction between the genetic and the permanent environmental components.

2 TABLEVariance components, heritability and repeatability (±SD) estimated for BHB and for all the other phenotypes using single‐trait models

• 3 Abbreviations: BHB, β‐hydroxybutyrate in milk (mM); F:P, fat‐to‐protein ratio; FP, fat percentage; LC, lactose; MY, milk yield (kg/day); PP, protein percentage; SCSLOG, logarithmic transformation of somatic cells count; UR, urea (mg/100 ml).
• 4 a Variance components for BHB were multiplied by 10,000, variance components for and F:P and LC were multiplied by 1,000; variance components for FP, PP and SCSLOG were multiplied by 100. All numbers were rounded to the second digits.

BHB heritability of present study is in agreement with previous estimates based on pedigree analysis. For instance, h2 in Italian Holsteins for BHB measured in the first 100 days of lactation was 0.08 ± 0.01 (Benedet et al., 2018). Koeck et al. (2014) reported h2 values ranging from 0.14 to 0.29 for Canadian Holsteins. Similar h2 was reported for Korean Holsteins (Ranaraja et al. 2018) ranging from 0.14 to 0.09 in the first and fourth parity, respectively. Weigel et al. (2017) estimated heritability for three BHB phenotypes (maximum concentration, square root scale, binary scale) in Holsteins using both BLUP and ssGBLUP (0.06–0.09). In Weigel et al. (2017), heritabilities were marginally higher when genomic information was used, for example h2 for BHBMAX was 0.059 ± 0.045 for the pedigree‐based and 0.074 ± 0.042 for the genomic‐based approach. These authors reported large standard deviations for heritability (~0.04) in both approaches, which are larger than those estimated in the present study (~0.01). Luke et al. (2019) reported genomic‐based h2 of 0.09 ± 0.04 for BHBLOG10 measured within 30 days after calving in 1,393 Holstein‐Friesian Australian cows. Benedet et al. (2018) observed a repeatability of 0.20 for BHB, whereas Klein et al. (2020) reported values lower than 10% due to small PE variances (genomic‐based). Repeatability estimates closer to the ones in the present study were found in the literature for blood BHB: Benedet et al. (2020) reported r2 of 0.26 ± 0.01 in early‐lactation (up to 35 days in milk) Holstein cows (pedigree‐based). In the present study, PE variances for BHB were greater than the additive genetic variances with GIBBS and ssGGIBBS. Because of that, r2 for BHB was on average twice as large as h2.

Table 3 shows phenotypic and genetic correlations between milk BHB and the other milk traits used in this study. Heritability of BHB in the two‐trait models was always between 0.09 and 0.10, which is similar to the values in the single‐trait analysis. Repeatability ranged from 0.23 ± 0.01 to 0.25 ± 0.01 (the largest value was obtained in the bivariate analysis with fat‐to‐protein ratio). Heritability for the other traits ranged from 0.07 ± 0.01 (SCSLOG) to 0.28 ± 0.02 (PP), whereas repeatability ranged from 0.22 ± 0.01 (FP) to 0.49 ± 0.01 (LC). These values were very similar to those highlighted in the single‐trait models. Estimates of h2 and r2 with GIBBS and ssGGIBBS were almost the same; however, some slight differences were observed for genetic correlation.

3 TABLEHeritability, repeatability and genetic correlation (±SD) estimated using the two‐traits models. Each model had BHB and one other important milk trait

• 5 Abbreviations:
• 6 BHB, β‐hydroxybutyrate in milk (mM); F:P, fat‐to‐protein ratio; FP, fat percentage; LC, lactose; SCSLOG, logarithmic transformation of somatic cells count; MY, milk yield (kg/dAY); PP, protein percentage; UR, urea (mg/100 ml).
• 7 a Genetic correlation in bold showed standard deviation larger than the estimates.
• 8 b Phenotypic correlation were all significantly ≠0 (p < 0.001).

Phenotypic and genetic correlations between BHB and milk traits (Table 3) were generally weak. Additionally, they were very similar between GIBBS and ssGGIBBS, with the latter showing slightly lower standard deviations due to the larger amount of information. Genetic correlation between BHB and UR and BHB and SCSLOG had standard deviations greater than the estimates, and therefore, could be not considered different from zero (for both considered approaches). BHB showed negative correlations with FP, PP and LC, and positive with MY and F:P.

Disregarding the sign, the weakest and strongest phenotypic correlations were found with MY (−0.01) and with the F:P (0.20), respectively. The largest genetic correlation was estimated between BHB and PP (−0.27 ± 0.06). The positive, albeit weak, genetic correlation of BHB with MY (0.12 ± 0.09) could be explained by the fact that NEB occurs more frequently in high‐yielding dairy cows (McCarthy et al., 2010). Several studies have reported that cows with high genetic merit for milk yield have NEB (Buckley et al., 2000; Gordon et al., 1995), and higher level of NEFA and BHB (Hart et al., 1978; Oldenbroek et al., 1997). Consistently, BHB showed weak, negative genetic correlations with fat and protein contents. The negative genetic correlation between milk BHB and PP in the present study is in agreement with previous studies on blood BHB in Norwegian Red cows (Belay et al., 2017).

The genetic correlation between BHB and F:P (0.13 ± 0.07) was of particular interest because F:P is an indicator of energy balance and could be used to identify subclinical ketosis. This finding is supported by the strong genetic correlation observed between ketosis and MIR predicted milk BHB in dairy cattle (Jamrozik et al., 2016; Koeck et al. 2015). In the present work, milk BHB was negatively correlated with the concentration of lactose (−0.14 ± 0.06). Costa, Egger‐Danner, et al. (2019) found a negative genetic correlation between lactose % and ketosis in Austrian Fleckvieh cows. The negative, albeit weak, genetic correlation observed between milk BHB and lactose could be explained by the unfavourable genetic correlation found between ketosis and mastitis (Costa, Egger‐Danner, et al., 2019). In fact, BHB and lactose are considered as indicators of ketosis and mastitis, respectively. A decrease in the milk lactose concentration is often associated to subclinical or clinical mastitis, whereas a positive concentration of BHB is an indicator of ketosis. High‐yielding cows, which are more susceptible to mastitis, and show lower milk lactose concentrations, could be more susceptible to ketosis, and therefore could exhibit higher concentrations of milk BHB. This result suggests that there may be an indirect (desired) selection for udder health when selecting for resistance to metabolic diseases (Costa, Egger‐Danner, et al., 2019).

Table 4 shows the validation statistics of the breeding values for BHB in the single‐trait and two‐trait models. As expected, in both cases, the use of genomic information led to more accurate and precise breeding values. ssGBLUP showed also a larger correlation between breeding values estimated in the reduced and whole data sets ( ${\rho }_{w,r}$ ), suggesting that this approach is more consistent between subsequent evaluations. The use of genomic information in both single and two‐trait models led to an increase of 0.06 in prediction accuracy. Moreover, ssGBLUP showed a lower inflation (i.e., a value closer to 1, which is desirable) compared to BLUP in both single‐trait and two‐trait models. A comparison between BLUP and ssGBLUP for milkability in Italian Simmental reported greater accuracies and correlations for the ssGBLUP, whereas the inflation was similar (Cesarani et al., 2021).

4 TABLELinear regression‐based statistics for breeding values of validation cows (n = 344)

• 9 Abbreviations: BHB, β‐hydroxybutyrate in milk (mM); F:P, fat‐to‐protein ratio; FP, fat percentage; LC, lactose; SCSLOG, logarithmic transformation of somatic cells count; MY, milk yield (kg/day); PP, protein percentage; UR, urea (mg/100 ml).
• 10 a Statistics: $\widehat{\mathit{acc}}$  = prediction accuracy; $b_{w,r}$  = prediction dispersion; ${\rho }_{w,r}$  = correlation between (G)EBV in the full and reduced data set.
• 11 b Predictive ability: correlation between adjusted phenotypes in the full data set and (G)EBV in the reduced data set.

The inclusion of genomic information in the single‐trait model led to an increase of 0.07 (Table 4) in predictive ability. Predictive ability was largely used as validation strategy in livestock species (Bengtsson et al., 2020; Lourenco et al., 2015; Magalhães et al., 2019; Manzanilla‐Pech et al., 2020).

Regarding the two‐trait models, validation statistics were similar to those obtained using the single‐trait model. No great advantages could be found using the two‐trait models because the genetic correlations between BHB and the other traits are generally weak (Table 3) and therefore, BHB does not gain much information from other phenotypes. However, some slight differences were observed. For both BLUP and ssGBLUP, higher prediction accuracies and stronger correlations were found when BHB was analysed with F:P. This is very interesting because both traits are considered indirect predictors of ketosis. BHB showed the strongest positive genetic correlation with F:P (0.13 ± 0.07) and this could have had an impact on the higher accuracies. On the contrary, the lowest dispersion and the weakest correlations were highlighted when BHB was analysed jointly with PP: this could be due to the negative genetic correlation between these two traits (see Table 3). Predictive ability for BHB was lower in three two‐trait models (BHB‐FP, BHB‐LC, BHB‐SCSLOG) compared to the single‐trait model: lowest values were 0.29 (BHB‐LC) and 0.36 (BHB‐LC and BHB‐FP) for BLUP and ssGBLUP, respectively. For all the other two‐trait models, we observed the same predictive abilities (0.31 for BLUP and 0.38 for ssGBLUP) as in the single‐trait model. These figures are in agreement with a study of Guo et al. (2014) that found the same reliabilities in single or multi‐trait models. These authors found higher reliabilities for multi‐trait in comparison with single‐trait models only when phenotypic records were missing for some animals. Additionally, VanRaden et al. (2014) analysed the differences in correlations of truncated with current evaluations and reported similar, but slightly better, predictive ability for multi‐trait than for single‐trait models.

A GWAS was performed to elucidate the genetic background of milk BHB (Figure 1). The four significant SNPs, according to the Bonferroni‐corrected multiple test (red line in Figure 1), were all located between 57.5 and 58.2 Mb on BTA20. This region was also significant on a GWAS for milk BHB in Holsteins, where deregressed proofs for bulls were used as pseudo‐phenotypes (Nayeri et al., 2019). In the same study, significant regions associated with BHB in the second and later lactations were found on BTAs 3, 6, 11, 14 and 25. The Inorganic Pyrophosphate Transport Regulator (ANKH) gene, involved with skeletal development, was mapped at 58.3–58.5 Mb on BTA20 (Nayeri et al., 2019): this gene was close (less than 100kb and therefore within the chosen boundaries) to the SNPs that passed the Bonferroni threshold in the present study. Particularly, the ANKH gene was associated with inorganic pyrophosphate transport regulator that helps to prevent the deposition of Ca and P in the bones (Sanchez et al., 2018). This gene was found to be highly expressed in the mammary gland compared to other 17 tissues in Holstein and Jersey cows (Kemper et al., 2015), and it was associated with lactose percentage in Holsteins, Jerseys (Lopdell et al., 2017) and in Fleckvieh cattle, which is closely related to Simmental (Costa, Schwarzenbacher, et al., 2019). Although intronic variants of ANKH were mainly reported to be associated with milk alpha‐lactalbumin in Montbéliarde (Sanchez et al., 2018), Normande and Holstein cattle (Sanchez et al., 2017; Zaalberg et al., 2020), this gene also explains a relatively large portion of the phenotypic variance in mineral content (Mg) and has effects on protein composition (Sanchez et al., 2019). Finally, the ANKH gene was found in a window explaining around 1.14% of the additive genetic variance for birth weight in Simmental beef cattle (Zhuang et al., 2020).

Identifying SNPs associated with traits of interest can be important for fine mapping purposes, that is the location and impact of important SNPs can be identified and used to guide other studies. Furthermore, the identification can possibly help to improve prediction accuracy if the SNPs explain a large proportion of variance on the trait, that is the important SNPs may receive more weights in the genomic prediction machinery (Wang et al., 2012). This is particularly true for traits with just few markers or QTLs contributing to the majority of the genetic variance. Unfortunately, only a few traits show high percentage of variance explained by just one or a few SNPs (such as DGAT1, ABCG2, CSN, GHRH). Most of the economically important traits are more polygenic, including functional traits such as milkability (Chen et al., 2020) and fertility (Ma et al., 2019).

Remarkably, the significant genomic region we found on BTA20 for milk BHB in Italian Simmental was also associated with this trait in Holstein (Nayeri et al., 2019). This implies that regardless of the sample size, this region could be truly associated with milk BHB, and therefore, further studies are needed to confirm this association. The effect of the ANKH gene on BHB milk content could be explained by the genetic correlation between the latter and other milk production traits (such as lactose and protein). This gene has been associated with lactose and protein in milk (Costa, Schwarzenbacher, et al., 2019; Lopdell et al., 2017; Sanchez et al., 2019). Fat‐to‐protein ratio and protein yield are indicators of dysmetabolism as well as BHB. Thus, the ANKH gene could play an indirect role on the expression of BHB.

In the present study, the genetic background of the β‐hydroxybutyrate milk content as an indicator of ketosis was investigated in Simmental dual‐purpose cattle. Although the limited genetic variability, it may be possible to consider this trait in the Italian Simmental selection scheme. In particular, the inclusion of this trait was not detrimental in the estimation of variance components for other traits already included in the selection programme because of the low genetic correlations. The use of genomic information did not considerably improve the estimates of genetic parameters; however, the accuracy of estimated breeding values increased. Therefore, ssGBLUP should be preferred to BLUP for selecting animals for milk BHB, and therefore for the risk of ketosis. The strong signal on chromosome 20 that was identified by the genome‐wide association study suggests the presence of a QTL and requires further investigation.

We gratefully acknowledge the Associazione Nazionale Allevatori Pezzata Rossa Italiana (ANAPRI) for providing data. The present research was conducted thank to the grant "Fondo di Ateneo per la ricerca 2019—Una tantum per la ricerca", founded by the University of Sassari.

The authors declare that they have no competing interest.

Data subject to third party restrictions. Data supporting this study are available from ANAPRI. Restrictions apply to the availability of these data, which were used under license for this study.