This study used 4,118 BoD crossbred calves sired by Belgian Blue (WBB) or Angus (ANG) and born to Holstein (HOL) dams. Between December 2021 to December 2023, the weight of calves was recorded once at the age of 0–40 days on 225 dairy farms in Schleswig-Holstein, Germany. Information regarding the insemination date and parity number of the dam, birth date and sex of the calf, type of birth (singleton or twin) and calving difficulty (CD) recorded on the official German scale (Arbeitsgemeinschaft Deutscher Rinderzüchter; organization of cattle production in Germany) and converted into binary scale (coded, 0 = no difficulty, 1 = difficult calving) was available. Gestation length (GL) was calculated as the duration (in days) between the date of insemination and the date of calving. Recorded values for weight and gestation length deviating ±4 SD from the mean were removed to filter for outliers. Additionally, twin calvings and farms with less than two observations were excluded from the analysis. After editing, 285 animals were excluded from the analysis and the final dataset contained 3,833 calves from 116 farms.

3,530 crossbred calves were genotyped using EuroG MD BeadChip (Illumina Inc). Quality control of genotypes was performed separately within ANG and WBB crossbreds, respectively, using PLINK 1.9 (Chang et al., 2015). SNPs with a minor allele frequency of <1%, a call rate lower than 90% and SNPs that deviate from the Hardy–Weinberg equilibrium at threshold of 1 × 10⁻⁶ along with variants located on sex chromosomes were excluded from the analysis (Supplementary Table S1).

In order to approximate the birth weight (BW) of the calves recorded not at the day of birth, a correction of measured weight was performed within each breed using the following linear regression model (Equation 1): y g = μ + AGE g + e g (1) where y_g was the measured weight (kg) of a calf within the sire breed WBB or ANG, AGE_g represented the regression of calf’s age at the day of measurement (g = 0, … ,40) and e_.g., was the residual term. The model was fitted in R (Becker et al., 1988) using the package lme4 (Bates et al., 2015). BW of the calf was calculated by subtracting the average predicted weight gain over the period of time from the measured weight of the calf.

Population stratification of the animals under study was examined using principal component analysis (PCA) in the individual breeds but also in the combined population (COM) via PLINK 1.9 (Chang et al., 2015). The first three components of PCA were used for visualization.

Haplotype phasing was conducted using SHAPEIT v2 (Delaneau et al., 2013) for each autosomal chromosome with default parameters for MCMC iteration combined with pedigree information using–duoHMM option to apply pedigree based post hoc haplotypes correction (O’connell et al., 2014). Haplotypes for Bos taurus autosomes (BTA) 1 to 29 were combined and converted into the binary format using the R package GHap (Utsunomiya et al., 2016). For further analysis, haplotype blocks were generated with a sliding window of five consecutive SNPs as this approach has been found to have the highest power to detect associated QTLs (Braz et al., 2019).

Genetic parameters for BW, GL and CD were estimated by using both univariate additive models and dominance models. SNP-based variance components were estimated using restricted estimation of maximum likelihood (REML) using the software GCTA (Yang et al., 2011) while narrow sense heritability for the trait was estimated as h 2 = σ a 2 / σ p 2 , where σ a 2 represents the additive genetic variance and σ p 2 is the phenotypic variance.

GWAS were performed for each trait using univariate single SNP regression mixed linear models in GCTA (Yang et al., 2011). The additive model was (Equation 2): y = Xb + Wa + Zu + e (2) where y is the vector of phenotypes (BW, GL or CD); b is a vector of fixed effects (herd-year-season, sex, parity number, first three principal components as covariates), a is the vector of polygenic effects with a ∼ N 0 , G σ a 2 , where G is the genomic relationship matrix (grm) (Yang et al., 2010) and σ a 2 is the additive genetic variance, u is the vector for SNPs for homozygous major and minor (0,2) and heterozygous alleles (1) and e is the random residual effect with e ∼ N 0 , I σ e 2 , where σ e 2 is the residual variance and I, X, W and Z represent incidence matrices for b, a, u and e respectively. For the GWAS model based on haplotypes, a represents the vector of polygenic effects, where the estimation of G is based on haplotypes.

For the dominance model (Equation 3), a modified version of Equation 2 was used: y = Xb + Wa + Wd + Zu + e (3) where d is the vector of dominance effects with d ∼ N 0 , D σ d 2 , where D is the dominance grm (Equation 4) defined as D j k = 1 m ∑ i w D i j w D i k (4) where m is the number of SNPs, w _D(ij) and w _D(ik) dominance encoded genotypes, σ d 2 is the dominance variance and u is the vector with the modified version of SNP coding (homozygous (0) and heterozygous (1) alleles). All the remaining variables are the same as in the Equation 2. Dominance variation at all SNPs was defined as δ 2 = σ d 2 /( σ a 2 + σ d 2 + σ e 2 ). The significance threshold for the SNPs was determined using Bonferroni correction, where a p-value of P < (0.05/N), with N representing the number of SNPs or pseudo-SNPs analyzed in the GWAS, indicates statistically significant associations. A suggestive association level was set at P < (1/N).

Linkage disequilibrium (LD) for each crossbred population was measured as the correlation coefficient r² (Hill and Robertson, 1968) by using PLINK v1.9 (Chang et al., 2015). LD decay between adjacent markers up to the distance of 500 kb was visualized in R using package ggplot2 (Wickham, 2011). For the functional annotation of SNPs with significant associations, base pair position was checked within 150 kb up- and downstream of the respective SNPs against the genome assembly Bos taurus UMD 3.1.1 to find relevant genes in the region. This distance was defined based on the LD decay (Qanbari, 2019), where the highest drop in first 150 kb of distance was observed in both crossbred populations (Additional file 1: F1). For haplotype blocks, screening for genes was limited to the significantly associated blocks in the analysis. The origin of significant haplotype blocks was determined by comparing them with parental haplotypes. This was achieved by matching the blocks with the corresponding dam and sire haplotypes with same block positions and allelic patterns to confirm their parental origins (Eiríksson et al., 2021a).

Descriptive statistics for the crossbred calves are depicted in Table 1. The distribution of the calves per parity number and age group at weight measurement is shown in Additional file 1: F2 & F3. The majority of calves belonged to WBB (Supplementary Table S2). These crossbreds had slightly higher daily weight gain (0.77 kg) compared to ANG crossbreds (0.66 kg). Moreover, calves from WBB had higher birth weight (48.5 kg ± 7.2) and longer gestation length (281.4d ± 4.7) compared to ANG sired calves (BW: 45.2 kg ± 7.2; GL: 280.1d ± 4.7). Incidences of calving difficulty were higher for ANG crossbreds, where 12.5% of calvings were labeled as difficult calving compared to WBB crossbreds with 11.0% difficult calvings.

Based on the PCA (Figure 1), clear clusters were identified. PC1 and PC2 collectively explained 48% of the variation in the combined data. Slight overlapping was observed for crossbred population based on PC1 and PC2 but with the addition of PC3 (58% variance), distinct clusters were observed across two studied populations (Additional file 1: F3, F4, F5)

Descriptive statistics for haplotype blocks are shown in Table 2. In total 11,160 haplotype blocks were created. In the combined population, highest number of haplotype blocks were observed on BTA1 (n = 726) while the lowest number of blocks were found on BTA28 (n = 210). Further details regarding haplotype blocks and number of markers per chromosome are presented in Additional file 1: F6, F7.

Estimates of heritability (s.e.) for BW, GL and CD in the combined population using the additive model were 0.29 (0.03), 0.36 (0.04) and 0.09 (0.03), respectively. For WBB, Dominance variation at all SNPs δ 2 for BW and GL was close to zero (Table 3), while for CD it was estimated to be 0.08 (0.07) (Table 3; Supplementary Table S3).

The study used two distinctive models to estimate heritability for BW, CD and GL. The effect of crossbreeding on heritability estimates depends upon the population and varies from trait to trait, but generally slightly lower heritability is expected for crossbred populations (Wientjes and Calus, 2017). This trend is also observed in the current study where BW heritability for WBB (0.26) crossbreds is lower compared to the purebred WBB (0.38) calves (Mota et al., 2017). However, heritability for ANG crossbred is slightly higher (0.36) compared to purebred ANG (0.33) (Torres-Vázquez et al., 2018). Non-additive genetic variance is expected to be higher in crossbred populations resulting in lower estimates for heritability (Fuerst and Solkneri, 1994). In the current study heritability estimates for CD varied from 0.11 (0.03) in the combined population to 0.07 for ANG and WBB BoD calves. These estimates are closer to the one reported in British herds, where heritability for CD in BoD crossbreds is reported to be 0.09 (0.01) (McGuirk et al., 1998), but lower to the once reported for purebred ANG (0.21) and WBB (0.25–0.34) breed (Cubas et al., 1991; Vanderick et al., 2017).

Heritability for GL varied from medium 0.36 (0.04) in combined crossbreds to high in ANG BoD 0.50 (0.11). The estimates for WBB and ANG BoD are closer to the purebred counterparts with 0.33 (0.04) for WBB and 0.59 (0.01) for ANG calves (Mota et al., 2017; Gilleland et al., 2021). A strong influence of the sire breed on gestation length is expected and can be the reason for the higher heritability estimates of ANG BoD calves (Haile-Mariam and Pryce, 2019). Heritability estimates can vary across the populations based on traits under observation, but generally traits with higher h ² in purebred animals will also have higher h ² in crossbred population (Wientjes and Calus, 2017). For the traits with dominant genes, crossbred animals may have higher additive genetic variance compared to the purebred animals (Wei et al., 1991). Comparison across purebred and crossbred populations is challenging due to the higher environmental variance in crossbred populations due to the scale effect (Habier et al., 2007), although no general trend is observed in this regard (Wientjes and Calus, 2017). Compared to purebred animals, the less controlled environmental conditions and greater variability in farming practices associated with crossbred animals may contribute to the observed differences in this regard (Wei and Van der Werf, 1995).

Utilization of GWAS can help to understand the genetic architecture of relevant traits. In this study, GWAS was performed using three different models for BW, CD and GL in BoD crossbreds. GWAS based on haplotypes outperformed the SNP-based GWAS in terms of finding the associated genes for BW, but no significant differences were observed for GL and CD. The addition of dominance effects in the GWAS model was beneficial for BW only in terms of finding genomic regions with significant or suggestive associations. Based on SNP and haplotype GWAS, in total five genes with significant association and seven genes with suggestive association were identified for the BW trait. Role of these genes varied from ion-channels regulation to regulation of mineralization process of the bone. GABRG1 gene, identified through significant SNP association is responsible for the regulation of the ion channel and receptor activity along with the influence on puberty development (Tahir et al., 2021). The association of this gene with milk yield has also been previously reported in the Holstein population (Pedrosa et al., 2021). SPP1 gene identified based on the dominance model encodes for multifunctional cytokines (Matsumoto et al., 2019) and is vital for bone mineralization process (Rodriguez et al., 2014). Association between SPP1 gene and growth traits including yearling weight, live weight and average daily gain has been suggested in American beef crossbred herds (White et al., 2007). Moreover, the influence of this gene on the birth weight (Allan et al., 2007) and carcass weight has been established in beef cattle (Matsumoto et al., 2019). For the genes identified through haplotypes-based GWAS, the influence of CCSER1, KCNIP4, and LCORL genes has been reported on BW and growth traits in the Angus breed (Xia et al., 2017; Smith et al., 2022). The CCSER1 gene, controlling mitosis has association with milk yield in Holstein cows (Teng et al., 2023). Similarly, PPARGC1A gene, influencing fat deposition and energy metabolism has been associated with BW (Pasandideh, 2020). Previously FAM13A gene with regulatory role in metabolism has already been described in the Holstein breed for association with bone structure (Niu et al., 2021; Dominguez-Castaño et al., 2024) Moreover, associations between PPARGC1A, a gene with significant role in fat metabolism (Komisarek and Walendowska, 2012) and yearling weight along with carcass traits have been observed in beef cattle (Fonseca et al., 2015). In this study, higher statistical power was observed for Haplotype-based GWAS in identifying associations and candidate regions. Similar results have been previously described in various species, including cattle (Pryce et al., 2010; Barendse and Schneider, 2011) and pigs (Sato et al., 2016). To some extent these differences can be attributed to the window size for haplotypes and the number of markers per haplotype can potentially influence the accuracy of QTL mapping (Calus et al., 2009). In contrast to aforementioned studies, SNP-based GWAS was found to be more efficient for carcass traits in Simmental cattle (Wu et al., 2014).

Interest in the approaches in the partitioning of crossbred genome has gained significant interest in the recent years for the crossbred evaluation. Various methods with varying degree of success have been developed including breed base representation (BBR) and breed of origin of the allele (BOA) (VanRaden et al., 2020; Eiríksson et al., 2021b; Guillenea et al., 2023). Despite the advancements in the statistical models and the increased number of genotyped animals, 100% precise estimation of the share of the genome of each parent in the crossbred population is difficult due to the shared DNA content between breeds (VanRaden et al., 2020). In this study, despite the partial success of assigning haplotypes to the parental breeds, complete segregation of the haplotypes of the crossbred calves was not successful. For most of the haplotypes with significant associations, the origin of the haplotypes was traced back to both parental haplotypes, though six breed-specific haplotypes with similar pattern of alleles were identified for ANG (n = 3) and WBB (n = 3) breeds. All of these six haplotypes were also present in HOL population (Supplementary Table S5).

A key limitation of the current study was the insufficient number of genotyped dams as this may reduce the accuracy of haplotype phasing. Moreover, for calving difficulty and gestation length, maternal influence is significant, and unavailability of sufficient number of maternal genotypes and phenotype records can be a limiting factor in this regard. Furthermore, the unequal number of crossbreds from two beef breeds can have impact on the accuracy of genetic parameters and association studies.