Genetic parameters and association between longevity and milk production in buffaloes using the ssGBLUP method

ABSTRACT The objectives of this work were to estimate the genetic parameters for the traits longevity (LG) and accumulated milk yield at 305 days (MY305) using a bitrait animal model and the single-step GBLUP method and estimate the genetic gain for LG through direct and indirect selection for MY305. We used 4,057 records of first lactations of Murrah dairy buffaloes, collected between 1987 and 2020, belonging to six Brazilian herds located in the states Ceará, Rio Grande do Norte, and São Paulo and 960 animals genotyped using the 90K Axiom Buffalo Genotyping (Thermo Fisher Scientific, Santa Clara, CA) to estimate the genetic parameters. The heritability estimate was 0.25 for MY305 and 0.13 for LG. The genetic gain for LG was 0.13 months under direct selection, and 0.14 months under indirect selection, which results in a relative selection efficiency of 11% under selection for MY305 compared with the direct selection. The genetic correlation between the two traits was 0.77, indicating that animals with genetic potential for high MY305 tend to live longer. The genetic trends for MY305 and LG were 0.22 kg/year and 5.20 days/year, respectively, indicating a positive response, which reaffirms its relationship with the high genetic correlation between the two traits.


Introduction
Buffalo milk is characterized by its high industrial yield and its constituents, especially due to the high levels of total solids (15.4 to 17.4%; Bernardes, 2013), being used in milk products such as mozzarella and yogurt (Michelizzi et al., 2010;Liu et al., 2018).Buffaloes also stand out for their high fertility rates and longevity (LG), which allows their exploration for both meat and milk (Leite et al., 2020).
Longevity (permanence of cows in the herd) is an economically important trait and has been related to profitability (Stefani et al., 2018), reducing its costs on the purchase of replacement females, as longer-lived animals generate more offspring for a longer period and have higher milk yield.It is known that buffaloes have a longer productive life than cattle (Gabr et al., 2015), reaching about 20 years of age (Marques et al., 2020).However, selection for LG is limited by its low heritability estimates (0.09 to 0.14) in cattle and buffaloes (Pander et al., 2002;Jenko et al., 2015).Thus, it becomes interesting to study LG from correlated traits, which may present higher heritability estimates and are easy to measure, such as milk yield (0.22 to 0.31; Aspilcueta-Borquis et al., 2010;Barros et al., 2016;Guzman et al., 2020).Accumulated milk yield at 305 days (MY305) is the most important selection criterion in buffalo breeding programs, and therefore selection for MY305 can be used to favor the indirect genetic gain for correlated traits such as LG.
Genomic evaluation approaches are a reality in animal breeding, as they allow the direct use of information from the individual's DNA in the genetic evaluation process, increasing selection accuracy, optimizing genetic gain, and reducing costs (Almeida et al., 2016).The Single-Step Genomic Best Linear Unbiased Predictor (ssGBLUP; Legarra et al., 2009;Misztal et al., 2009;Christensen and Lund, 2010) combines pedigree and genomic information to predict the genomic estimated breeding value (GEBV) of individuals.
This approach has been successfully employed for the genetic evaluation of different species; however, studies that genetically correlate LG and milk yield in buffaloes are still scarce in Brazil (Bashir et al., 2007;Aspilcueta-Borquis et al., 2022).Therefore, the objectives of this study were to estimate the genetic parameters LG and MY305, using a bitrait animal model by the ssGBLUP approach, and estimate the genetic gain for LG through direct and indirect selection for milk yield.

Phenotypic and genotypic data
We used 323,141 records of MY305 and LG from the first lactations of 4,057 and 2,783 Murrah buffaloes, respectively.Data were collected between 1987 and 2020, belonging to six herds located in three Brazilian states (Ceará, Rio Grande do Norte, and São Paulo).The traits analyzed were MY305, defined as the total amount of milk produced between the 5th and 305th days in milk, and LG, defined as the total amount of time that the animal remained in the herd.Age at first calving ranged between 22 and 60 months and was included as a covariate (quadratic for MY305).
Data consistency was checked for LG and MY305.Only animals that had already been discarded were evaluated for LG.For MY305, only lactations of at least 90-day duration were evaluated, with the first test day (TD) before 70 days postpartum, and animals that suffered at least 3 TD were selected for further analysis.Phenotypic records that deviated ±3 standard deviation from the contemporary group (CG) mean and with less than four animals per CG were removed (Table 1).The pedigree file included 93,633 animals.Descriptive statistics for the LG and MY305 traits are described in Table 1.

Genotype quality control
The quality control (QC) for the genotypes was performed using Plink v.1.9software (Purcell et al., 2007;Chang et al., 2015).A total of 978 animals were genotyped with the 90K panel (123,040 SNPs) Axiom ® Buffalo Genotyping (Affymetrix, Santa Clara, CA).The following criteria were considered Genetic parameters and association between longevity and milk production in buffaloes using... Carvalho et al.
3 in the quality control: minimum call rate of 0.95, P-value for Hardy-Weinberg equilibrium test less than 10 −6 , and a minimum allele frequency (MAF) greater than 0.03.Only samples with a call rate greater than 0.90 were kept.After the QC, 45,690 autosomal markers and 960 animals remained (212 males and 748 females).

Statistical analysis
The (co)variance components and genetic parameters for the MY305 and LG were estimated by the single-step Genomic Restricted Maximum Likelihood (ssGREML) method, and the breeding values were predicted by the ssGBLUP method.Bitrait analyses were performed under an animal model, using the AIREMLF90 software (Misztal et al., 2015).The bitrait animal model is described below: in which y is the vector of observations, X is the incidence matrix associating the fixed effects (CG, covariate) to the vector b, Z is the incidence matrix associating the random additive genetic effect to the vector g, and e is the random residuals vector.The contemporary groups (CG) were defined based on herd-year-calving season for MY305 and herd-year-birth season for LG, and the seasons for both traits were defined as dry season (April -September) and rainy season (October -March).Age at first calving was included as a quadratic covariate for MY305 and ranged between 22 and 60 months.
In addition, the following assumptions were made: E[y] = Xb; the additive genetic effect and the residual effect are normally distributed with mean zero and Var(g) = H⨂S g ; Var(e) = I⨂S e ; in which S g is the random additive genetic covariance matrix, S e is the random residual covariance matrix, and I is an identity matrix.
In the single-step approach, the inverse of the relationship matrix A is replaced by the inverse of relationship matrix H, that includes both pedigree-based relationships and differences between pedigree-based and genomic-based relationships (Aguilar et al., 2010) as follows: in which A is the matrix of numerator of Wright's coefficient of relationship; A 22 is the section of A related to the genotyped animals; and G is the genomic relationship matrix, as proposed by the "first method" by VanRaden (2008).

Accuracy and bias
To assess the accuracy and bias of the GEBV, a reduced subset was created from the complete data, excluding information from the last four years (records of animals genotyped up to 2010).Thus, 150 genotyped animals constituted the validation population.
To validate the genomic predictions, the estimated breeding values (EBV) predicted in the full dataset were regressed on the GEBV or on the parental breeding values (PA) in the reduced dataset, according to the following equation: in which y c refers to the EBV vector of the animals from the validation population; b 0 and b 1 refer to the intercept and the slope of the regression, respectively; â r refers to the GEBV/PA of the animals; and e refers to the residual component.b 1 was used as a bias indicator in the genomic prediction and in the traditional genetic evaluation.

Genetic gain and correlated response
The genetic gain and the expected correlated response were calculated for the studied traits, considering the same selection intensity (equal to one) for both, based on the following equations: in which ΔG LG is the genetic gain based on direct selection for LG; i is the selection intensity; σ p is the phenotypic standard deviation of the trait under selection; h 2 is the heritability of the trait under selection; RC LG-MY305 is the correlated response obtained for LG, through direct selection for MY305; r aLG-MY305 is the genetic correlation between MY305 and LG; h MY305 and h LG are the square roots of the estimates of heritabilities (h 2 ) for MY305 and LG, respectively; σ pLG is the phenotypic standard deviation of LG; ER LG-MY305 refers to the relative efficiency of indirect selection for LG relative to the direct response for MY305.
Additionally, the estimates of genetic trends were obtained by linear regression of the EBV and GEBV means on year of birth.

Variance components and genetic parameters
The heritability estimate for MY305 was 0.25.For the LG trait, the heritability presented low magnitude (0.13).The genetic correlation obtained between LG and MY305 was 0.77, being positive and of high magnitude (Table 2).

Accuracy, regression coefficient, and validation
The accuracies were similar between the pedigree and genomic approaches (Table 3).However, the approach based on genomic information were less biased for both traits.Genetic parameters and association between longevity and milk production in buffaloes using... Carvalho et al. 5

Genetic gain and genetic trend
The genetic gain for LG was 0.13 months through direct selection; however, considering the correlated response through direct selection for MY305, i.e., indirect selection, we obtained an estimate of 0.14 months, which resulted in a higher relative selection efficiency of 7% under selection for MY305.
The genetic trends estimated by the regression of breeding values on the year of birth of the animals were significant (P<0.05) and equal to 0.22 kg/year (y = 0.22x -3.117) and 5.20 days/year (y = 5,228x -59,575) for MY305 and LG, respectively, both positive.The genetic trend remained close to zero for many years (1969 to 1972) and began to increase after the 2000s (Figure 1).

Variance components and genetic parameters
The heritability estimate obtained in this study for MY305 (0.25, Table 2) was similar to those reported by several studies with buffalo species, whose estimates ranged from 0.16 to 0.35 (Breda et al., 2010;Aspilcueta-Borquis et al., 2012;Camargo et al., 2015;Liu et al., 2018;Abdel-Shafy et al., 2020;Guzman et al., 2020;Lázaro et al., 2021).These heritability values demonstrate that there is additive genetic variance for this trait; consequently, genetic gains can be achieved when using it as a selection criterion.
On the other hand, despite the importance of LG in the herd, this trait has still been little explored in the literature.Longevity is a trait that involves several parameters that influence the animals' lifespan, such as milk production, health, reproductive efficiency, type traits, production system, and manager's decision regarding slaughter.Therefore, it is often difficult to make comparisons between LG studies due to the different measurement methods and terms applied to this trait (Galeazzi et al., 2010;Schuster et al., 2020;Dallago et al., 2021).Despite this, some studies GEBV -genomic estimated breeding value.GEBV LG (days) GEBV MY305 (kg)

MY305
LG Linear (MY305) Linear (LG) have reported low-magnitude estimates of heritability for LG.Tamboli et al. (2021) reported a heritability of 0.06 in Nili-Ravi buffaloes.A higher estimate (0.20) was observed by Kamaldeep et al. (2015), when defining LG in days in Murrah buffaloes.Galeazzi et al. (2010), evaluating the ability to stay in the herd (stayability) between one and six years, reported heritability estimates ranging from 0.11 to 0.23 in Brazilian Murrah buffaloes.Similarly, estimates ranging from 0.05 to 0.18 have also been reported in cattle in several European countries (Forabosco et al., 2009;van Pelt et al., 2015).The low heritability estimate for LG indicates that the trait variation may be more associated with environmental influence and/or non-additive genetic effects than the direct additive genetic component.Therefore, direct selection would result in a slow response in the progress of the trait.
There is a lack of studies reporting the genetic correlation between MY305 and LG in buffaloes, but some studies have reported genetic correlation estimates for these traits in dairy cattle (0.04 to 0.35; Wasana et al., 2015).Silva et al. (2016) estimated a positive genetic correlation (0.61) between the permanence in the herd at 60 months and MY305 in a herd of dairy Gyr cattle, concluding that selection for MY305 could be used to improve the animals' ability to remain in the herd.Stefani et al. (2018) found a negative genetic correlation (−0.25) between permanence in the herd at 60 months of age and MY305 in Holstein cattle.In general, the negative genetic correlations can be attributed to the intense selection for productivity in dairy cows, in which high-performance animals may present greater chances of metabolic disorders and reduced body condition, which causes the involuntary culling of these animals, therefore negatively affecting their LG (Lucy, 2001;Valsalan et al., 2014;De Vries and Marcondes, 2020;Schuster et al., 2020).Hence, we can infer that the high and positive genetic correlation obtained in the present study (0.77) may be due to the low intensity of selection for milk production in buffaloes in the past, differently from what happened with cattle.

Accuracy, regression coefficient, and validation
According to Oliveira et al. (2019), the number of animals required to provide an additional increase in prediction accuracy needs to be greater for traits with low heritability (such as LG) than for traits with higher heritability (MY305).In general, the bias for GEBV and PA were less than 1.0, indicating that the predictions were inflated, i.e., they do not correspond exactly to the genetic reality of these animals, and the predicted differences in the genetic merit are biased upwards compared with observed future performance (Oliveira et al., 2019).However, it is important to mention that approaches based on genomic information were less biased than approaches based on pedigree, for both traits, which agrees with other works with buffaloes (Tonhati et al., 2008;El-Bramony et al., 2010).

Genetic gain and genetic trend
The direct genetic gain for LG was 0.13 indicating the increase in life span per generation by selecting LG.The correlated response for LG would be 0.14 if selection was based on MY305, equivalent to a relative selection efficiency of 7%.Thus, the indirect selection for LG based on MY305 would be more efficient than the direct selection for LG.This can be explained mainly due to the low heritability estimate observed for LG and the high heritability estimate of MY305, in addition to the high and positive genetic correlation (0.77) between LG and MY305.Furthermore, our results agree with some results reported in the literature, for instance, Bashir et al. (2007) in their work with Nili-Ravi buffaloes (lifelong milk production = 15 kg/year and productive life = four days/year).
The fact that the genetic trend remained close to zero from 1969 to 1972 and increased after the 2000s can be explained by the increase in information in the database maintained by the genetic breeding program of the Departamento de Zootecnia of the Faculdade de Ciências Agrárias e Veterinárias, Universidade Estadual Paulista (FCAV-UNESP), Jaboticabal, SP, Brazil, which allowed the selection of animals according to their EBV for MY305.Furthermore, selection to increase MY305 generated a correlated response in LG, due to the high genetic correlation observed.It is Genetic parameters and association between longevity and milk production in buffaloes using... Carvalho et al. 7 crucial to constantly evaluate genetic correlation and genetic and phenotypic trends for monitoring genetic progress over time.

Conclusions
The heritability value obtained for MY305 can be considered as a selection criterion in breeding programs, while LG is not a good criterion for direct selection.Thus, selecting LG through MY305 tends to enable superior genetic gains for LG, in addition to the fact that MY305 can be measured earlier in the animals' lives.Genetic trends for both traits were positive and significant, which can be explained by the high genetic correlation between MY305 and LG.At the same time, the inclusion of genomic information helped to obtain more accurate breeding values and, consequently, increased the reliability for milk yield and longevity in this study with dairy buffaloes.

Table 3 -
Accuracy measured by Pearson's correlation coefficient between estimated breeding values (EBV) MY305 -accumulated milk yield at 305 days; LG -longevity.Regression coefficient (b 1 ) for GEBV with EBV (b GEBV ) and for EBV and PA (b PA ) for MY305 and LG in Murrah buffaloes.