Dataset

The phenotypic dataset consisted of 2,314AFC records of Murrah buffalo cows born between 1995 and 2017. Contemporary groups (CG) were formed considering animals born on the same farm, year and season of birth, which was divided into two (dry and rainy season). Animals with records outside the range between ±3 standard deviations from their CG averages and animals belonging to GCs with less than five individuals were removed.

From a pedigree file containing 3320 buffaloes, 553 animals (539 dams and 14 sires) were genotyped with the Axiom Buffalo Genotyping Array 90 K (Iamartino et al., Reference Iamartino, Nicolazzi, Van Tassell, Reecy, Fritz-Waters, Koltes, Biffani, Sonstegard, Schroeder, Ajmone-Marsan, Negrini, Pasquariello, Ramelli, Coletta, Garcia, Ali, Ramunno, Cosenza, Oliveira, Drummond, Bastianetto, Davassi, Pirani, Brew and Williams2017). Only SNPs present on autosomal chromosomes (BBU1-BBU24 referenced by UOA_WB_1 genome assembly), as well as those with call rate > 95%, MAF > 3% and significance level for Hardy Weinberg equilibrium test was 10⁻⁵ , were remained in the analysis. All samples had a call rate >90%. The database description (phenotypic and genotypic information) after the consistency step is in Table 1.

Table 1. Description of age at first calving dataset used in this study

PBLUP, ssGBLUP and WssGBLUP models

This set of methodologies is based on the use of relationship matrices between individuals in mixed model equations, with the differences consisting of the type of information used (pedigree or genomic) and the way these matrices are constructed. The model using pedigree based BLUP (PBLUP) can be described as:

$$y = {\boldsymbol X}\beta + {\boldsymbol Z}\alpha + e$$

where, y is the vector of phenotypic records (AFC); β is the fixed effect vector (CG), X is the incidence matrix that associates β with y; α is the vector of additive genetic effects, Z is an incidence matrix associating α with y and, e is the vector of residuals. The following assumptions were made:

$$E( y ) = {\boldsymbol X}\beta \quad E( \alpha ) = 0\quad E( e ) = 0$$

$${\rm Var}( \alpha ) = A\sigma _\alpha ^2 \quad {\rm Var}( e ) = {\boldsymbol I}_n\sigma _\varepsilon ^2 $$

where $\sigma _\alpha ^2$ and $\sigma _\varepsilon ^2$ represent the additive genetic and residual variances; A and I_n are the pedigree relationship and an identity matrices, respectively.

The single step GBLUP (ssGBLUP) method is an extension of the PBLUP model, in which the pedigree-based (A) and genomic (G) relationship matrices are combined into a single matrix (H), as described by Aguilar et al. (Reference Aguilar, Misztal, Johnson, Legarra, Tsuruta and Lawlor2010). So, its inverse can be obtained as:

$${\boldsymbol H}^{{-}1} = {\boldsymbol A}^{{-}1} + \left[{\matrix{ 0 \hfill & 0 \hfill \cr 0 \hfill & {{\boldsymbol G}^{{-}1}-{\boldsymbol A}_{22}^{{-}1} } \hfill \cr } } \right]$$

where A⁻¹ and G⁻¹ are the inverse matrices of A and G respectively, and ${\boldsymbol A}_{22}^{{-}1}$ is the inverse of the section of A related to the genotyped animals. The G matrix was obtained according to VanRaden (Reference VanRaden2008) as:

$${\boldsymbol G} = \displaystyle{{{\boldsymbol ZD{Z}^{\prime}}} \over {2\;\sum {\,p_i} ( 1-p_i{) }^{\prime}}}$$

where p _i is the minor allele frequency (MAF) of SNP i, Z is a matrix relating genotypes of each locus centred by allele frequencies, and D is a diagonal matrix of weights for SNP variances, (with elements d _i = SNP i  weight). In ssGBLUP model is set D = I, so that the weight of all SNPs is equal to 1.

In WssGBLUP method (Wang et al., Reference Wang, Misztal, Aguilar, Legarra and Muir2012) markers are assigned with different weights using an iterative process described with the following steps:

1. 1. Set D ^(t) = I, when t = 1

2. 2. The genomic relationship matrix is setup for t as ${\boldsymbol G}^{( t ) } = {\boldsymbol Z}{\boldsymbol D}^{( {\boldsymbol t} ) }{\boldsymbol {Z}^{\prime}}/( 2\;\sum p_i( {1-p_i} )$

3. 3. The Genomic Breeding Value (GEBV or $\hat{\alpha }$) for t are obtained;

4. 4. 4.The GEBV $\hat{\alpha }^{( t ) }$are converted to SNP effects ($\hat{u}$) as $\hat{u}^t = {\boldsymbol D}^{( t ) }{\boldsymbol {Z}^{\prime}}( {\boldsymbol G}^{( t ) }) ^{{-}1}\hat{\alpha }^{( t ) }$

5. 5. The weight of the ith SNP (the ith element of D or d _i) is calculated;

6. 6. SNP weights are normalized in D^(t+1) to have constant genetic variances of SNP effects;

7. 7. Loop to step 2 three times and exit.

For the methods mentioned above, the program BLUPF90+ (Misztal et al., Reference Misztal, Tsuruta, Lourenço, Masuda, Aguilar, Legarra and Vitezica2014) was used for the estimation of the variance components and for the validation step. The estimation of the effects of the markers and the calculation of the weighting (WssGBLUP) were performed with the program POSTGSF90 (Misztal et al., Reference Misztal, Tsuruta, Lourenço, Masuda, Aguilar, Legarra and Vitezica2014). The heritability estimates for AFC were obtained from these estimated variance components.

Single-step Bayesian regression models (ssBR)

The ssBR models combine all available data from genotyped and non-genotyped animals and use imputed marker covariates for animals that are not genotyped (Fernando et al., Reference Fernando, Dekkers and Garrick2014), can be described as:

$$\left[{\matrix{ {y_n} \cr {y_g} \cr } } \right] = \left[{\matrix{ {{\boldsymbol X}_{\boldsymbol n}} \cr {{\boldsymbol X}_g} \cr } } \right]\beta + \left[{\matrix{ {{\boldsymbol Z}_n} & 0 \cr 0 & {{\boldsymbol Z}_g} \cr } } \right]\left[{\matrix{ {{\boldsymbol M}_n\alpha + \epsilon } \cr {{\boldsymbol M}_g\alpha \;\;\;\;\;\;\;} \cr } } \right] + e$$

where vectors and matrices for non-genotyped animals are denoted with the subscript n and for genotyped animals with the subscript g. Thus, y _n and y _g represent the vectors of phenotypic observations; β is the vector of systematic effects of GC (equivalent to the fixed effect in frequentist methods); X_n and X_g design matrices for the fixed effects; Z_n and Z_g represent incidence matrices associated with the genomic values of the animals; α, $\epsilon$ and e, represent in this order, the vectors of marker effects, imputation residual and model residual; M_g is the marker matrix for the genotyped animals and M_n is the imputed marker matrix for the non-genotyped animals. The matrix M_n can be obtained as:

$${\boldsymbol M}_n = {\boldsymbol A}_{ng}{\boldsymbol A}_{22}^{{-}1} {\boldsymbol M}_g$$

where A_ng represents the pedigree-based relationship matrix between genotyped and non-genotyped animals and; ${\boldsymbol A}_{2g2g}^{{-}1}$ represents the inverse of the pedigree-based relationship matrix between genotyped animals. The following ssBR methods were used: BayesA (ssBA- Meuwissen et al., Reference Meuwissen, Hayes and Goddard2001), BayesBπ (ssBBπ- Meuwissen et al., Reference Meuwissen, Hayes and Goddard2001), BayesCπ (ssBCπ- Habier et al., Reference Habier, Fernando, Kizilkaya and Garrick2011), Lasso Bayesian (ssBL- Yi and Xu, Reference Yi and Xu2008) and Bayesian ridge regression (ssBRR- Campos et al., Reference Campos, Hickey, Pong-Wong, Daetwyler and Calus2013), where the assumptions for markers effects are shown in Table 2. The hibayes package (Yin et al., Reference Yin, Zhang, Li, Zhao and Liu2022), available for the R program (R Core Team, 2021), was used to perform the ssBR analyses. A total of 350 000 samples and a burn-in period of 150 000 samples were generated. The convergence was evaluated using graphical analysis. The heritability estimates for AFC were also obtained using these models.

Table 2. Description of different prior for markers effects used in Bayesian single-step regression methods used for genomic analysis of age at first calving in dairy buffaloes

Validation of (G)EBVs

The predictive ability of the methods was assessed using accuracy and dispersion statistics proposed by Legarra and Reverter (Reference Legarra and Reverter2018). This methodology is based on predicting (G)EBV using full ($\hat{u} _w$) and a partial $( \hat{u} _p)$ datasets. The partial dataset (reference group) was obtained by omitting phenotypic records of genotyped animals born after 2010 (validation group). Thus, to assess the prediction ability of the models, only genotyped animals from the validation group were used to calculate the following statistics:

$$\matrix{ {{\rm Accuracy}} \hfill & {{\rm Dispersion}} \hfill \cr {r_{A\hat{A}} = \sqrt {\displaystyle{{{\rm cov}( {\hat{u}_w, \;\hat{u}_p} ) } \over {( {1-\bar{F}} ) \hat{\sigma }_g^2 }}} } \hfill & {b_{w, p} = \displaystyle{{{\rm cov}( {\hat{u}_w, \;\hat{u}_p} ) } \over {{\rm var}( {\hat{u}_p} ) }}} \hfill \cr } $$

where $\bar{F}$ is the average inbreeding coefficient of the validation individuals and $\hat{\sigma }_{g, i}^2$ is the estimate of additive genetic variance.