Genetic Analysis of Persistency for Milk Fat Yield in Iranian Buffaloes (Bubalus bubalis)

This study aimed to estimate heritabilities and genetic trends for different persistency measures for milk fat yield and their genetic correlations with 270-day milk yield in Iranian buffaloes. The records of test-day milk fat yield belonging to the first three lactations of buffaloes within 523 herds consisting of 43,818 records were got from the Animal Breeding Center and Promotion of Animal Products of Iran from 1996 to 2012. To fit the lactation curves based on a random regression test-day model, different orders of Legendre polynomial (LP) functions were selected. Three persistency measures were altered according to the specific condition of the lactation curve in buffaloes: (1) The average of estimated breeding values (EBVs) for test day fat yield from day 226 to day 270 as a deviation from the average of EBVs from day 44 to day 62 (PM1), (2) A summation of contribution for each day from day 53 to day 247 as a deviation from day 248 (PM2), and (3) The difference between EBVs for day 257 and day 80 (PM3). The estimates of heritability for PM1, PM2, and PM3 ranged from 0.20 to 0.48, from 0.36 to 0.47, and from 0.19 to 0.35 over the first three lactations, respectively. The estimate of genetic trends for different persistency measures of milk fat yield was not significant over the lactations (P > 0.05). Genetic correlation estimates between various measures of persistency were generally high over the first three lactations. Also, genetic correlations estimates between persistency measures and 270-day milk yield were mostly low and varied from 0.00 to 0.24 (between PM1 and 270-day milk yield), from −0.19 to 0.13 (between PM2 and 270-day milk yield), and from −0.02 to 0.00 (between PM1 and 270-day milk yield) over the first three lactations, respectively. Persistency measures that showed low genetic correlations with milk fat yield were considered the most suitable measures in selection schemes. Besides, medium to high heritability estimates for different persistency measures for milk fat yield indicated that relevant genetic variations detected for these characters could be regarded in outlining later genetic improvement programs of Iranian buffaloes.

INTRODUCTION One important step for reaching self-sufficiency in any country is to identify the productive potential of native breeds of animals. The great adaptability of native animals to harsh conditions such as high environmental humidity and temperature, irregular rainfall, the incidence of different diseases, weak management practices, and low quality of feeds causes native buffaloes of Iran to play an important role in supplying milk and meat as major protein sources. Currently, many Asian countries depend mainly on buffalo as a source of milk and dairy products, especially in rural areas (Safari et al., 2018).
One of the main factors in determining the total milk production over a lactation period is persistency (Muir, 2004). Persistency is defined as the potential of an animal to maintain milk yield at a high extent after reaching the peak of production. The other definition of persistency is the gradual decrease of daily milk production after reaching the peak of the lactation curve (Togashi and Lin, 2004). The major cause for the worth of buffaloes with more persistent curves is that they can relatively satisfy most parts of their nutrient requirements from roughages (Sölkner and Fuchs, 1987). Therefore, not only metabolic problems, reproduction disorders, and diseases are lower in cows with more persistent lactations, but also production costs would be lower (Dekkers et al., 1998). Determining the method of measuring persistency is a critical point in estimating genetic progress for this trait. However, no general agreement is existent on the most appropriate method to describe the persistency of lactation (Cole and VanRaden, 2006). Various measures were suggested for calculating persistency (Gengler, 1996): measures based on the functions describing persistency; measures based on a fraction of total yield, peak yield, or parts of lactation; and those based on the breeding value of animals derived from analyzing random regression models.
The method used for defining persistency measures would determine the genetic parameter estimates for these measures and their genetic relationship with milk production (Swalve and Gengler, 1999;Jakobsen et al., 2002;Khorshidie et al., 2012). A measure of persistency must have two characteristics: association with lactation curve flatness, and independent explanation from production level. The latter item implies that the genetic correlation between milk yield and persistency measures to be decreased because milk production explains some genetic variance of persistency measures under study (Muir, 2004;Cole and VanRaden, 2006;Khorshidie et al., 2012). The independence of these two traits causes genetic selection for persistency of lactation and total yield to avoid unfavorable consequences of peak yield stress in high-yielding cows. Also, the incidence of metabolic diseases and reproductive disorders would be minimized while high milk production is maintained (Elmaghraby, 2012;Ghavi Hossein-Zadeh et al., 2017).
Previous studies carried out on dairy cattle indicated that lactation persistency positively correlated with favorable reproductive performance and health status (Jakobsen et al., 2002;Muir, 2004). Such favorable correlations along with the positive economic value for persistency would support including lactation persistency in the genetic improvement programs of cattle and buffalo (Dekkers et al., 1998;Khorshidie et al., 2012;Ghavi Hossein-Zadeh et al., 2017).
The random regression models enable fitting random genetic and environmental effects at different stages of lactation, which results in higher accuracy of estimated breeding values (EBVs) compared with other statistical models (Li et al., 2020). These models provide insights about the temporal variation of biological processes and physiological implications underlying the studied traits. Therefore, random regression models generate relevant information to be exploited in breeding programs (Oliveira et al., 2019). The functions generally used to model the lactation curve include Wood's model (Wood, 1967), Wilmink's function (Wilmink, 1987), spline function (White et al., 1999), and Legendre polynomial (LP) function (Kirkpatrick et al., 1990). Because of variations in production environments and management systems, optimal functions for test-day models in various countries may be different (Mrode et al., 2003). But several studies have indicated that LPs performed well in random regression test-day models (Li et al., 2020).
Milk constituents can be used as a simple indicator of the nutritional status of the lactating animals. Because of the dilution effect, milk fat percentage shows the opposite direction of the lactation curve for milk yield (Eicher, 2004;Ghavi Hossein-Zadeh, 2016), but fat yield follows a variation trend similar to milk yield over the lactation. When trying to apply milk composition as a nutritional evaluation tool, these fluctuations should be noticed. Although several researchers have studied the genetic analysis of the persistency for milk yield and components in dairy cattle (Cole and Null, 2009;Khorshidie et al., 2012;Canaza-Cayo et al., 2015), limited studies have been performed to estimate genetic parameters of persistency for milk production traits in buffaloes (Ghavi Hossein-Zadeh et al., 2017). Therefore, the objective of the present study was to estimate the heritability and genetic trend of distinct persistency measures for milk fat yield and their genetic correlations with 270-day milk yield in Iranian buffaloes using random regression test-day models.

Data
Records of test-day milk fat yield belonging to the first three lactations of Iranian native buffaloes in 523 herds consisting of 43,818 records were provided by the Animal Breeding Center and Promotion of Animal Products of Iran during 1996-2012. According to climatic conditions, Iranian native buffaloes can be grouped into three main classes: Azari ecotype, Kuhzestani ecotype, and Mazandarani or North ecotype . Borghese (2005) and Ghavi Hossein-Zadeh (2015a,b) described the overall management practices and population structure of buffaloes in Iran. Outliers that appeared to deviate markedly from other observations in the original data set were discarded. Therefore, the subsequent analyses included only production records corresponding to the first three lactations in which days in milk (DIM) were between 5 and 270. Calving ages ranged between 24-60, 39-76, and 54-100 months for the first, second, and third lactations, respectively. The total number of test-day records per animal was from 4 to 9. Summary statistics of the edited data set are presented in Table 1. The number of animals, sires, and dams in the pedigree of Iranian buffaloes was 42,285, 549, and 6,376, respectively.

Statistical and Genetic Analysis
Legendre polynomial functions were chosen to fit the lactation curves in the framework of a random regression test-day model for estimating (co)variance components. Model specification and the choice of fixed effects to be included in the model were based on the backward elimination method and variables which were significant at P < 0.05 were considered in the model. To obtain the appropriate random regression test-day model for the genetic analysis of test day fat yield, with the minimum number of parameters, different orders of fit for random regression coefficients of additive genetic and permanent environmental effects were evaluated. Also, the optimum set of polynomials was selected according to the logarithm of the likelihood function at the point of conversion and the total number of parameters to be estimated. The difference of these models was based on the LPs applied to fit the covariance functions for additive genetic and permanent environmental effects. The maximum logarithm likelihood of the models was compared and models with the lowest values of this criterion were selected for further analysis. Test day records were analyzed using the following random regression model: Where, Y ijmnptv : test day record i obtained at DIM t of cow p calved at the n th age in herd-test date m, G i : fixed effect of i th breed or ecotype, YS j : fixed effect of j th calving year-season, HTD m : fixed effect of m th herd-test date, c f : the f th fixed regression coefficient for calving age, age n : the n th calving age, k: the order of fit for fixed regression coefficients (k = 2), β r : the r th fixed regression coefficient, k a : the order of fit for additive genetic random regression coefficients, k p : the order of fit for permanent environmental random regression coefficients, α pr : the r th random regression coefficient of additive genetic value for p th cow, γ pr : the r th random regression coefficient of permanent environmental effect for p th cow, ∅ r dim t : the r th coefficient of LPs evaluated at days in milk t, e ijmnptv : the random residual error. All random regression analyses were conducted using the Average Information Restricted Maximum Likelihood (AIREML) algorithm of the WOMBAT program (Meyer, 2006).

Lactation Persistency Measures
The following measures were used to describe lactation persistency in this study. These measures were modified based on the lactation curve conditions of buffaloes and adapted for 270 days lactation period: 1. The average of EBVs for test day fat yield from day 226 to day 270 as a deviation from the average of EBVs from day 44 to day 62 [adapted from Kistemaker (2003)]: A summation of contribution for each day from day 53 to day 247 as a deviation from day 248 [adapted from Cobuci et al. (2007) and Jakobsen et al. (2002)]: 3. The difference between EBVs for day 257 and day 80 [adapted from Cobuci et al. (2004Cobuci et al. ( , 2007]: Small absolute values of the abovementioned measures indicate a high lactation persistency. Ifα i was a (k a ×1) vector of the estimates of additive genetic random regression coefficients specific to the animal i, and Z t was a k a × 1 vector of LP coefficients evaluated at day t, the EBV of animal i for day t was calculated as follows: Therefore, the EBV of animal i for 270-day production was obtained by summing the EBVs from day 5 to day 270: Where, Z c270 is a vector of the summations of LPs corresponding to total lactation yield. In addition to the 270-day yield, we could estimate a Z c corresponding to each persistency measures used in the current study as follows: For the first lactation fat yield: For second lactation fat yield:

Estimation of Genetic Parameters and Genetic Trends
The following formulas were applied to estimate additive genetic, permanent environmental and residual variances and heritabilities for different measures of persistency for fat yield and 270-day milk yield: are genetic covariances and correlations between persistency measures and 270-day milk yield, respectively. Estimates of genetic trends for persistency measures were obtained by regressing the average EBVs on the calving year of animals.

RESULTS
The orders of fit for different random regression testday models of milk fat production are given in Table 2.
The maximum log-likelihood values of test-day models 1, 10, and 10 differed significantly (P < 0.05) from the other models for fat yield in the first three lactations, respectively. Thus, models 1, 10, and 10 were chosen to fit the additive genetic and permanent environmental effects for the analysis of fat production in the first three lactations of buffaloes, respectively. Heritability estimates of persistency measures for fat production and estimates of genetic correlation among distinct fat yield persistency measures with each other and with 270-day milk production in Iranian buffaloes are presented in Table 3. Heritability estimates for PM 1 , PM 2 , and PM 3 ranged between 0.20-0.48, 0.36-0.47, and 0.19-0.35 for the first, second, and third lactations, respectively. In general, heritability estimates fluctuated largely among lactations and persistency measures. The highest estimate of heritability was observed for PM 1 in the third lactation (0.48), while the lowest one was recorded for PM 3 also in the third lactation.
Genetic correlation estimates among various measures of persistency were generally high and ranged from 0.98 to 0.99 (between PM 1 and PM 2 ), from −0.98 to −0.87 (between PM 1 and PM 3 ), and from −0.99 to −0.95 (between PM 2 and PM 3 ) over the first three lactations, respectively. Also, genetic correlation estimates between persistency measures and milk yield were mostly low and varied from 0.00 to 0.24 (between PM 1 and 270-day milk yield), from −0.19 to 0.13 (between PM 2 and 270-day milk yield), and from −0.02 to 0.00 (between PM 3 and 270-day milk yield) across the first three lactations, respectively ( Table 3).
Variation of milk yield and milk fat yield across the first three lactations of Iranian buffaloes are depicted in Figures 1, 2. The trend of observed milk yield and milk fat yield for all lactations  increased from day 5 of lactation to a peak several weeks later, declining thereafter until day 270. Genetic trends of persistency measures for milk fat yield are illustrated in Table 4. In general, all estimates are very low and not significant (P > 0.05). Therefore, they would not be considered different from zero. Changes in EBVs of buffaloes for three persistency measures of milk fat yield according to calving year and lactations are illustrated in

DISCUSSION
For many years, the breeding objectives of dairy animals emphasized increased milk yield. But negative genetic associations were observed between numerous functional characters with production traits, and decreases in genetic excellence for fitness and health have been detected in dairy farms (Egger-Danner et al., 2015). The management practices must be directed toward the compensation of these effects and to equalize reproduction performance, metabolic diseases, and udder health vs. enhanced production to maximize profit without any negative impact on animal welfare. Because concerns on animal welfare and consumers' appeal for natural and health products are increasing, the functional traits have received greater importance in animal breeding programs (Egger-Danner et al., 2015). In this regard, it is required to have valid genetic parameter estimates for outstanding traits related to the farm profit, including functional traits, in the animal breeding programs (Fleming et al., 2018). Interest to include new traits in the current animal breeding programs is extending to improve simultaneously the production and reproduction performance along with animal health and well-being in dairy farms. Although, for the inclusion of a specific trait into a genetic selection program, it would be inheritable, profitable, quantifiable, and changeable (Wood et al., 2003). Although, there were some reports on the genetic analysis of persistency measures for milk components in dairy cows (Cole and Null, 2009;de Oliveira Biassus et al., 2010), to the knowledge of authors, this is the first report on the genetic analysis of fat production persistency measures in buffaloes. In general, medium to high heritability estimates for three milk fat persistency measures in this study could be due to the reasonable additive genetic variations for these traits indicating that improvement in these traits could be attained by genetic selection. Regardless of the simpler estimation of PM 3 in contrast to other measures of persistency, the estimate of heritability for this measure was between the estimates of heritability for PM 1 and PM 2 measures for fat yield in the first lactation and had the smallest estimate in second and third parities. If a measure of persistency had higher heritability compared with other measures, this measure would be an appropriate measure to be considered in the selection objective (Ghavi Hossein-Zadeh et al., 2017). Respecting this explanation, the PM 2 measure would be regarded as the selection criterion in the first lactation, but the measure of PM 1 would be included as a selection objective in the second-, and third parities. Although there is no report of genetic parameters for persistency measures of fat yield in buffaloes, Cole and Null (2009) reported the estimates of heritability for fat production persistency measure varied from 0.07 to 0.12 in five breeds of dairy cows. Also, de Oliveira Biassus et al. (2010) reported the heritability estimates for fat yield persistency measures ranged from 0.00 to 0.23 in primiparous Holstein cows. Besides, Gengler (1995) estimated the heritability for fat yield persistency measure would be equal to 0.06 in dairy cows. In general, several factors could influence the variation in heritability estimates for milk fat yield persistency obtained in different studies, including the breed of the animal, within-population genetic diversity, management procedures, environmental conditions, and methods used for estimating genetic parameters. According to de Oliveira Biassus et al. (2010), different factors would influence the variations of heritability estimates for persistency measures among studies: the definition of persistency measure as absolute or relative terms, the statistical adequacy of the specific measure of persistency for under study population, the lactation period used to calculate the measure of persistency, and the method or model used to calculate a specific persistency measure. Compared with the first and third lactations, less variation in heritability estimates between the three persistency measures in the second lactation would be due to the differences in lactation curves, yield persistencies, and variation of records across the first three lactations. The production difference in two different parts of the lactation would be evaluated by the PM 1 and PM 3 persistency measures (Ghavi Hossein-Zadeh et al., 2017). Compared with PM 1 and PM 3 measures, the PM 2 measure displayed a domain below the lactation curve at a definite time that has been adjusted for yield at the end section of that period (Khorshidie et al., 2012). The procedure for defining the PM 2 and PM 3 persistency measures resulted in a high and negative genetic association between them. The positive genetic correlations between PM 1 and PM 2 are proof for the same genetic and physiological systems managing these persistency measures and would cause the same ranking of buffaloes according to these criteria in breeding and genetic schemes (Ghavi Hossein-Zadeh et al., 2017). Contrarily, high negative genetic correlations between PM 3 with two other persistency measures implied the existence of various mechanisms to govern them. In general, low genetic correlations between different persistency measures for fat yield with milk production point out that selection for a persistency measure for milk fat yield would slightly affect milk yield. In a selection program, it would be favorable to have persistency measures that had low genetic correlations with milk yield (Dekkers et al., 1998;Ghavi Hossein-Zadeh et al., 2017). According to this explanation and regarding the low genetic correlations of persistency measures for fat yield with milk production in the present study, all three measures would be considered as selection criteria that were relatively independent of production level in buffaloes. This finding indicates that a buffalo cow with the highest EBV for 270-day milk yield does not necessarily has the highest EBV for fat yield persistency and vice versa. In the other words, low estimates of genetic association between fat yield persistency measures with milk production signified that buffaloes with the identical quantity of 270-day milk yield could have a distinct extent of persistency across the lactation period (Jamrozik et al., 1998;Cobuci et al., 2007;Ghavi Hossein-Zadeh et al., 2017). The appropriateness of genetic correlation between a specific persistency measure for milk fat yield and milk production depends on the positive or negative mean of the persistency measure in the population under study (Khorshidie et al., 2012). Generally similar to the results of the present study, Cole and Null (2009) observed the estimates of genetic associations between persistency measure of fat yield with 305-day milk production varied from 0.07 to 0.29 in five breeds of dairy cows.
Predicting accurately the animals' breeding value is an appropriate way to increase the genetic gain in a specific breeding scheme (Ghavi Hossein-Zadeh, 2012). The successfulness of a selection scheme would be assessed by testing the actual alteration in breeding value indicated as a fraction of the expected theoretical modification in the average breeding value of the character under study (Jurado et al., 1994;. Non-significant genetic progress estimated for all fat production persistency measures in the present study and irregular changes in average EBVs of animals over the years demonstrated the non-presence of a clear breeding design for making better the lactation persistency for fat yield in Iranian buffaloes until now. A possible reason for the non-significant genetic trends of milk fat persistency measures would be the low and close to zero estimates of genetic correlation between fat yield persistency measures and 270-day milk yield in the population under study.

CONCLUSION
The persistency measures of fat yield proposed in the present study had favorable low genetic correlations with 270-day milk production. These low correlations would be a benefit in designing a selection program to enhance the milk yield in Iranian buffaloes because buffaloes with the identical quantity of 270-day milk yield could have a distinct extent of persistency across the lactation period. The PM 2 measure had the highest heritability estimate for the first lactation buffaloes, but the PM 1 measure had the highest estimate in the second-and -third lactations. Therefore, the PM 2 measure would be regarded as the selection criterion in the first lactation, but the measure of PM 1 could be suggested as a selection objective in the secondand third parities. Based on the results of this study, it would be necessary to consider the persistency of fat yield in the selection objective of buffaloes in Iran together with main characters such as production and reproduction traits, and persistency for milk production.

DATA AVAILABILITY STATEMENT
The data analyzed in this study is subject to the following licenses/restrictions: The data analyzed in this study was obtained from the Animal Breeding Center and Promotion of Animal Products of Iran. Requests to access these datasets should be directed to http://abc.org.ir.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because datasets used in this study were obtained from pre-existing databases based on routine animal recording procedures.

AUTHOR CONTRIBUTIONS
MN participated in the acquisition of data, statistical, and genetic analyses of data. NH-Z designed and conceived this study and contributed to the statistical and genetic analyses of data, and prepared the manuscript. AS contributed to the conception of the study and assisted with the interpretation of the outputs. DK assisted with the interpretation of data. All authors read and approved the final manuscript.