Until recently, the effectiveness of dairy cattle breeding programs was inferred through genetic gains observed in the production traits (Stevenson and Britt, Reference Stevenson and Britt2017). However, due to unfavorable genetic correlations, continuous selection for increased milk production led to a significant decrease in herd fertility, resulting in negative economic impacts attributed to the increase in reproductive failures, involuntary culling of cows and increased extra costs with reproductive biotechnologies (Sun et al., Reference Sun, Madsen, Lund, Zhang, Nielsen and Su2010; Tiezzi et al., Reference Tiezzi, Maltecca, Cecchinato, Penasa and Bittante2012; Chiumia et al., Reference Chiumia, Chagunda, Macrae and Roberts2013). In this context, the number of studies evaluating different classes of reproductive traits (Averill et al., Reference Averill, Rekaya and Weigel2006; Berry et al., Reference Berry, Wall and Pryce2014) has increased with special emphasis to interval traits, since they address the efficiency of a cow in becoming pregnant or the start of a new reproduction cycle after calving, besides providing early information for genetic evaluations (Jorjani, Reference Jorjani2007; Gernand and König, Reference Gernand and König2017).
Genetic evaluations of reproduction interval traits have been widely based on repeatability models (González-Recio et al., Reference González-Recio, Pérez-Cabal and Alenda2004; VanRaden et al., Reference VanRaden, Sanders, Tooker, Miller, Norman, Kuhn and Wiggans2004), which assume that the same groups of genes are influencing the expression and/or regulation of these traits during the animals' productive live. Under these models, the environmental components between longitudinal measurements are assumed equally correlated, which is questionable since events closer in time should be more associated than distant ones. A suitable approach to overcome this problem may be fitting an autoregressive covariance structure to describe the time decaying correlation between consecutive events over time. This methodology has been successfully applied on genetic evaluations of milk yield and its components, where its efficiency in disentangling the genetic components resulted in higher heritabilities, higher reliabilities for estimated breeding values (EBV) and lower residual variances (Carvalheira et al., Reference Carvalheira, Blake, Pollak, Quaas and Duran-Castro1998, Reference Carvalheira, Pollak, Quaas and Blake2002; Silva et al., Reference Silva, Costa, Silva, Silva, Lopes, Santos, Thompson and Carvalheira2019b). Despite these positive results obtained in longitudinal production traits, to the best of our knowledge there are no reports of its application on genetic modeling of longitudinal reproductive traits.
Including reproductive traits in designing selection indexes is essential for the advancement of dairy breeding programs in order to mitigate losses due to diminishing fertility. In this study, we hypothesized that the autoregressive repeatability model (AR) will provide more accurate genetic evaluations to achieve that goal. Therefore, we aimed to evaluate the efficiency of AR models for genetic evaluation of longitudinal reproductive traits in Portuguese Holstein cattle and to compare those results with conventional repeatability models (REP).
Material and methods
Data
The data was provided by the Portuguese Dairy Cattle Breeders Association and comprised pedigree, management and reproduction records from the first three calvings of cows that occurred between 1995 and 2018 in 3457 herds. Preliminary analysis indicated that the voluntary post-partum waiting period (VWP) corresponded to approximately two estrous cycles (42 d) and for the purpose of this study, we assumed a standard gestation length of 280 ± 15 d.
The interval reproductive traits were calculated as the difference in days between the date of calving and the subsequent event under consideration: interval from calving to first service (ICF), expressed as number of days between the calving and the following first insemination; days open (DO), defined as number of days between calving and conception; calving interval (CI), computed as number of days between two consecutive calvings. Daughter pregnancy rate (DPR) was also evaluated and computed as the linear transformation of the percentage of estrous cycles (of 21 d) that the cow needed to get pregnant. DPR was calculated as DPR = 0.25(226-DO). This formula is an adaptation from the method presented by VanRaden et al. (Reference VanRaden, Sanders, Tooker, Miller, Norman, Kuhn and Wiggans2004), taking into account the specific VWP estimated for the Portuguese dairy population, while maintaining the relationship in which each increase of 1% in EBV for DPR equals a decrease of 4 d in EBV for DO.
The total of 495 516 records for ICF, 906 038 for DO and 907 014 for CI, were edited according to trait specifications ensuring that the records of each trait stayed within physiological limits and reliable management criteria. The following minimum and maximum limits were established: 42 to 168 d for ICF, 42 to 252 d for DO and 307 to 730 d for CI. Additionally, ICF and DO values between 21 and 41 d were set to 42 d in order to avoid penalize early-breeder cows, which corresponded to 3.8% and 2.4% of total data for ICF and DO, respectively. All fertile insemination dates were validated by requiring that the number of days to the following calving should be in the range of the assumed gestation length. When the insemination date between two consecutive calving was not available, the DO was assumed equal to CI minus 280 d, the mean gestation period of the population. Finally, it was also required that the cows' data had no missing information between consecutive calving orders. The descriptive statistics for all reproductive traits are presented in Table 1.
Phenotypic descriptive statistics for interval between calving to first service (ICF), days open (DO), calving interval (CI) and daughter pregnancy rate (DPR)
Statistical analysis
In a preliminary analysis, we investigated if the observations obtained from each calving order for each trait (ICF, DO, CI and DPR) could be considered as repeated records of the same trait. To perform this study, we applied a multiple trait model considering each calving order as a different trait within ICF, DO, CI and DPR, respectively. As expected, the genetic correlations between calving orders were high for all reproductive traits (online Supplementary Table S1), justifying the use of repeatability models for genetic evaluations.
Analysis based on AR and REP models were carried out using the following general model:
$${\bf y} = {\bf Xb} + {\bf Kc} + {\bf Za} + {\bf Wpe} + {\bf e}\comma \;$$where 
${\bf y}$ is the vector of observations for the first, second and third calving orders; 
${\bf b}$ is the vector of fixed effects considering twelve month levels and sixteen age-month levels associated to each calving order, which ranged from 18 to 44, 28 to 60 and 40 to 72 months for the first, second and third calving orders, respectively; 
${\bf c}$ is a vector of random herd-year-season (HYS) of calving effects, corresponding to 73 472 levels for ICF, 137 433 levels for DO and DPR and 141 131 levels for CI; 
${\bf a}$ is the random additive genetic effect; 
${\bf pe}$ is a vector of random permanent environmental (PE) effects; 
${\bf e}$ is a vector of residuals. 
${\bf X}\comma \;{\rm \;}{\bf K}\comma \;{\rm \;}{\bf Z}$ and 
${\bf W}$ are the incidence matrices relating records to fixed effects and to random effects of HYS, additive genetic and PE, respectively. The expectations and (co)variance for these models were:
$$\left[{\matrix{ {\bf y} \hfill \cr {\matrix{ {\bf c} \cr {\matrix{ {\bf a} \cr {\matrix{ {{\bf pe}} \cr {\bf e} \cr } } \cr } } \cr } } \hfill \cr } } \right]\sim N\left({\left[{\matrix{ {{\bf Xb}} \hfill \cr {\matrix{ {\matrix{ 0 \cr 0 \cr } } \cr 0 \cr } } \hfill \cr 0 \hfill \cr } } \right]\semicolon \;\left[{\matrix{ {\matrix{ {\bf V} \cr {\matrix{ {} \cr {\matrix{ {} \cr {\matrix{ {} \cr {} \cr } } \cr } } \cr } } \cr } } \hfill & {\matrix{ {\matrix{ {{\bf JK}} \cr {\matrix{ {\bf J} \cr {\matrix{ {} \cr {\matrix{ {{\rm Sym}} \cr {} \cr } } \cr } } \cr } } \cr } } & {\matrix{ {\matrix{ {{\bf GZ}} \cr {\matrix{ 0 \cr {\matrix{ {\bf G} \cr {\matrix{ {} \cr {} \cr } } \cr } } \cr } } \cr } } & {\matrix{ {\matrix{ {{\bf PW}} \cr {\matrix{ 0 \cr {\matrix{ 0 \cr {\matrix{ {\bf P} \cr {} \cr } } \cr } } \cr } } \cr } } & {\matrix{ {\bf R} \cr {\matrix{ 0 \cr {\matrix{ 0 \cr {\matrix{ 0 \cr {\bf R} \cr } } \cr } } \cr } } \cr } } \cr } } \cr } } \cr } } \hfill \cr } } \right]} \right)\comma \;$$ The differences between AR and REP models are on the modeling of the HYS and PE (
${\bf J}$ and 
${\bf P}$ matrices, respectively) where, as the name implies, first-order autocorrelation structures were fitted on the AR model.
For AR model,
$${\bf J} = {\rm \sigma }_{\rm c}^2 \left[{\matrix{ 1 \hfill & {{\rm \rho }_{\rm c}} \hfill & {{\rm \rho }_{\rm c}^2 } \hfill \cr {} \hfill & 1 \hfill & {{\rm \rho }_{\rm c}} \hfill \cr {{\rm Sym}} \hfill & {} \hfill & 1 \hfill \cr } } \right]\otimes {\bf I}\comma \;$$is a block diagonal matrix where 
${\rm \sigma }_{\rm c}^2 $ and ρc are the variance and the autocorrelation coefficient for HYS, respectively (the size of the autocorrelation matrix depends on the number of year-season levels within each herd – 3 in this example). 
${\bf I}$ is an identity matrix with size equal to the number of herds.
$${\bf P} = {\rm \sigma }_{{\rm pe}}^2 \left[{\matrix{ 1 \hfill & {{\rm \rho }_{{\rm pe}}} \hfill & {{\rm \rho }_{{\rm pe}}^2 } \hfill \cr {} \hfill & 1 \hfill & {{\rm \rho }_{{\rm pe}}} \hfill \cr {{\rm Sym}} \hfill & {} \hfill & 1 \hfill \cr } } \right]\otimes {\bf I}\comma \;$$is also a block diagonal matrix where 
${\rm \sigma }_{{\rm pe}}^2 $ and ρpe are the variance and the autocorrelation coefficient for PE, respectively (in this case, the size of the autocorrelation matrix depends on the number of repeated observations for each cow – 3 in this example). 
${\bf I}$ is an identity matrix with size equal to the number of animals with records.
On the other hand, for REP model, 
${\bf J} = {\bf I}{\rm \sigma }_{\rm c}^{2{\rm \;}} $, where 
${\bf I}$ is an identity matrix with size equal to number of HYS levels; and 
${\bf P} = {\bf I}{\rm \sigma }_{{\rm pe}}^2 $, where 
${\bf I}$ is an identity matrix with size equal to number of animals with records.
For both models, 
${\bf G} = {\bf A}{\rm \sigma }_{\rm a}^2 $, where 
${\rm \sigma }_{\rm a}^2 $ is the additive genetic variance component and 
${\bf A}$ is the numerator relationship matrix; and 
${\bf R} = {\bf I}{\rm \sigma }_{\rm e}^2 $, where 
${\bf I}$ is an identity matrix with size equal to the total number of observations and 
${\rm \sigma }_{\rm e}^2 $ is the residual variance component. Finally, 
${\bf V} = {\bf KJK}{\rm ^{\prime}} + {\bf ZGZ}{\rm ^{\prime}} + {\bf WPW}{\rm ^{\prime}} + {\bf R}$, represents the phenotypic (co)variance.
To estimate the variance components and autocorrelations, we used six sub-data sets extracted from the complete data for each trait (ICF, DO, CI and DPR). Each sub-data set consisted of 30 randomly sampled herds with records spanning across the all period of study and representing all regions of Portugal (online Supplementary Table S2). All parameters were estimated applying DFREML methodology (Smith and Graser, Reference Smith and Graser1986). Likelihood functions were maximized by the multivariate simplex algorithm (Nelder and Mead, Reference Nelder and Mead1965) programed in MATLAB (Math Works Inc., Natick MA). To achieve convergence, the variance among the simplex −2 log likelihoods had to be less than 10−8. The occurrence of local maxima was checked by running five consecutive cold starts without significant changes in the log-likelihood (up to four decimal places). This process was repeated for each sub-data set and the respective means were used as the variance components and autocorrelations estimates used in all subsequent analyses of each reproductive trait.
Model comparisons
Comparisons between AR and REP models were carried out on each of the six sub data for each trait. The best fit was evaluated by Akaike Information Criterion (AIC), computed as AIC = −2log (L) + 2p, where L is the maximum of likelihood function, and p is the number of effective parameters. Since, AIC provides only a qualitative comparison between models (smallest AIC is better), we also calculated the Akaike Weights (AW), presented by Burnham and Anderson (Reference Burnham and Anderson2004), using the formula 
${\rm A}{\rm W}_i = \exp \lpar {-\Delta_i/2} \rpar /\sum _{i = 1}^2 \exp \lpar {-\Delta_i/2} \rpar $, for i = 1 (AR) and 2 (REP). In this formula, Δi is the AIC difference between model i and the model that presented the smallest AIC value (the best model). The Δi for the best model equals zero (since i and the best model will coincide). With this approach, AWi represents the probability of model i to be better than the other model evaluated. Additionally, we also evaluated the Mean Square Error (MSE) for each model, calculated as 
${\rm MSE} = {\rm \;}n^{{-}1}\mathop \sum \nolimits \lpar \,{y-\hat{y}} \rpar ^2$, where 
$y\;\;{\rm and}\;\hat{y}$ are the observed and predicted vectors for each trait respectively, and n is the number of observations.
The EBVs obtained from AR and REP models were used to evaluate the changes in the ranking of sires using the Spearman's rank correlation among the top 100 and the total sire population. Complementarily, the coincidence percentage between the best 1, 10 and 50% of bulls was also used for comparisons. The significance of the difference between EBV reliabilities of these bulls, obtained using the Student's t-test at the 5% level of significance, was another criterion of comparison.
Additionally, plots of genetic trends based on the EBV means at year of birth and the respective genetic gain, derived by regressing yearly means for EBV, were used for comparisons of the potential genetic progress provided by the two models. The average EBV of cows born in 2010 was set as the genetic base-year for each trait. Only bulls with at least 10 daughters recorded in at least five distinct herds were considered in these analyses.
Results and discussion
Model-fitting criteria
The results for model-fitting criteria in the analysis of ICF, DO, CI and DPR are presented in Tables 2 and 3. The AIC estimates for AR model were lower for all reproductive traits, indicating a better fitting when compared to REP model. The AW estimates for AR model were all close to one while for the REP models these values were closed to zero. The MSE estimates were lower for AR, corresponding to the least amount of bias, and therefore, higher predictive ability.
Akaike Information Criterion (AIC) and Akaike Weights (AW) for the Autoregressive repeatability model (AR) and the traditional repeatability model (REP) from the analyses of six sub data of interval between calving to first service (ICF), days open (DO), calving interval (CI) and daughter pregnancy rate (DPR)
Mean Square Error (MSE) for the autoregressive repeatability model (AR) and the traditional repeatability model (REP) from the analyses of interval between calving to first service (ICF), days open (DO), calving interval (CI) and daughter pregnancy rate (DPR)
Lower AIC values, AW probability close to unity and the lower MSE estimates are strong indications of the superiority of AR model in fitting (time) reproductive traits. Modeling the HYS and the PE effects as autoregressive variables did contribute to a better explanation of the variation observed in ICF, DO, CI and DPR. These findings are in line with the reports of Quaas et al. (Reference Quaas, Anderson and Gilmour1984) and Carvalheira et al. (Reference Carvalheira, Pollak, Quaas and Blake2002), where they suggested relaxing the usual assumption of unitary correlation between repeated events by implementing first-order autoregressive structures, which may result in larger additive genetic variances and heritabilities and estimates of individual genetic merit more accurate for longitudinal traits.
Rank correlation, percentage of coincidence and reliabilities of EBVs for bulls
Re-ranking between the breeding values predicted by AR and REP models were relatively small when considering the total population of bulls in all traits (rank correlation varied between 0.95 and 0.97 as shown in Fig. 1). On the other hand, for the top-100 bulls the rank correlation were significantly lower (no overlapping of 95% confidence intervals was observed) and ranged from 0.71 to 0.88 (Fig. 1). These results are consistent with those presented by Sawalha et al. (Reference Sawalha, Keown, Kachman and Van Vleck2005), which reported rank correlation estimates ranging from 0.71 to 0.87 for top animals of population when evaluating the inclusion of autoregressive covariance structure in the test-day models for analyses of productive traits in dairy cattle. The moderate rank correlation estimates for the top 100 bulls may indicate a different ranking when considering one or another model, which may lead to changes in the choice of bulls by the farmers depending on the model used (Sawalha et al., Reference Sawalha, Keown, Kachman and Van Vleck2005; Sun et al., Reference Sun, Madsen, Lund, Zhang, Nielsen and Su2010). In this sense, the evaluation of reproductive efficiency using AR model could create more opportunities for selecting the best bulls and, therefore, promoting greater genetic progress.
Rank correlation for Top 100 and Total population of bulls (a) and percentage of coincidence between the best 1, 10 and 50% of bulls (b) ranked according to solutions predicted by the Autoregressive repeatability model (AR) and the traditional repeatability model (REP) models for interval between calving to first service (ICF), days open (DO), calving interval (CI) and daughter pregnancy rate (DPR).
To quantify the EBV re-ranking between the two models, we also computed the percentage of coincidence in the top 1, 10 and 50% of the bull population for each trait (Fig. 1). The percentage of coincidence increased with the increase in the number of bulls compared (from 80 to 84.6% in the top 1% and from 90 to 93% in the top 50%). The high percentage of coincidence observed indicates that both models captured almost the same top animals, which emphasizes the importance of the reliabilities of EBV predictions as an auxiliary tool to increase the selection pressure.
The bulls' EBV reliabilities predicted by REP and AR models were significantly different (P < 0.05) in all traits evaluated, favoring the AR (Table 4). According to Chegini et al. (Reference Chegini, Hossein-Zadeh, Moghaddam and Shadparvar2019), when genetic values are predicted with higher accuracy and precision, the greater will be the selection differential and, consequently, a higher response to selection are expected. Thus, the results presented in our study indicate that greater genetic progress may occur with the use of EBVs predicted by the AR model in the genetic evaluations of the reproductive traits evaluated.
Means of EBV reliabilities (Rel) for the best 1%, 10% and 50% of bulls obtained from the autoregressive repeatability model (AR) and the traditional repeatability model (REP) for interval between calving to first service (ICF), days open (DO), calving interval (CI) and daughter pregnancy rate (DPR)
* Significant difference (P < 0.05) between RelAR and RelREP within reproductive trait.
Genetic parameters
The variance components and autocorrelations estimated with AR and REP models for ICF, DO, CI and DPR, are presented in Table 5. In general, the estimates of the variance components were similar between models, with heritabilities ranging from 0.08 to 0.11 for AR and from 0.06 to 0.09 for REP. In the present study, the HYS was the most auto-correlated effect (ranging from 0.7 to 0.9) probably because of the shorter time lag between events (between seasons) compared with the PE effect where the events are much more distant (between calvings). A similar pattern was reported in autoregressive analysis of production traits (e.g., milk yield) where the autocorrelations between lactation orders (long term effect) were negligible compared with the short term autocorrelations between test-days within lactations (Carvalheira et al., Reference Carvalheira, Pollak, Quaas and Blake2002; Costa et al., Reference Costa, Carvalheira, Cobuci, Freitas and Thompson2009; Silva et al., Reference Silva, Costa, Silva, Silva, Lopes, Santos, Thompson and Carvalheira2019b).
Variance components, autocorrelation coefficients and genetic parameters (±se) estimated by the Autoregressive repeatability model (AR) and the traditional repeatability model (REP) for interval between calving to first service (ICF), days open (DO), calving interval (CI) and daughter pregnancy rate (DPR)
${\rm \sigma }_{\rm a}^2 $, additive genetic variance; 
${\rm \sigma }_{{\rm HYS}}^2 $, herd-year-season variance; 
${\rm \sigma }_{{\rm pe}}^2 $, permanent environmental variance; 
${\rm \sigma }_{\rm e}^2 $, residual variance; ρHYS, autocorrelation coefficient for herd-year-season effect; ρpe, autocorrelation coefficient for permanent environmental effect; h 2, heritability; r, repeatability.
The heritabilities estimated in our study are in agreement with previous results reported in the literature, which presented values ranging from 0.02 to 0.05 for ICF, from 0.02 to 0.07 for DO and CI and from 0.03 to 0.07 for DPR (Kadarmideen et al., Reference Kadarmideen, Thompson, Coffey and Kossaibati2003; Ghiasi et al., Reference Ghiasi, Pakdel, Nejati-Javaremi, Mehrabani-Yeganeh, Honarvar, González-Recio, Carabaño and Alenda2011; Haile-Mariam et al., Reference Haile-Mariam, Bowman and Pryce2013; Frioni et al., Reference Frioni, Rovere, Aguilar and Urioste2017; Gernand and König, Reference Gernand and König2017).
Genetic trends for bulls and cows
The genetic trends of reproductive traits resulting from the REP and AR analysis are depicted in Fig. 2, for bulls and cows, respectively. It is possible to identify two distinct periods for each class of animal. The intervals between 1980 and 2000 (in bulls) and from 1980 to 2008 (in cows) are characterized by a drop in reproductive efficiency of the population. However, in the years following those periods, there was a recovery in these traits with a decrease in the EBV means of ICF, DO and CI and an increase for DPR in both, bulls and cows. As expected, the shape and magnitude of the curves are very similar between models, reflecting the high correlation between EBVs. The genetic gains expressed by regressing yearly means for EBV on year of birth were also close between the two models evaluated, as presented in the Table 6.
Genetic trends in the bull (a) and cows (b) populations estimated from EBVs obtained by the Autoregressive repeatability model (AR) and the traditional repeatability model (REP) for interval between calving to first service (ICF), days open (DO), calving interval (CI) and daughter pregnancy rate (DPR).
Genetic gains expressed by regression coefficients (b), with the respective standard deviations and coefficients of determination (R 2) estimated from analyses with the Autoregressive repeatability model (AR) and the traditional repeatability model (REP) for interval between calving to first service (ICF), days open (DO), calving interval (CI) and daughter pregnancy rate (DPR) considering two different year intervals for bulls and cows
* P < 0.05.
It is worth noting that until the present, the reproductive performance was never included in any selection index for the improvement of dairy cattle in Portugal. This suggests that the genetic trends observed along those years were a result of indirect selection on production traits. In fact, the trends here observed are almost the opposite of those of milk production for the same period (Silva et al., Reference Silva, Silva, Silva, Costa, Lopes, Caetano, Thompson and Carvalheira2019a). As in other countries, the Portuguese dairy breeding program has as its main goal, the increase of milk production and its components. Thus, the positive genetic gains between 1980 and 2000 for bulls, and between 1980 and 2008 for cows, are probably consequences of the genetic progress achieved with selection for improvement of productive traits, emphasizing the negative genetic correlation between these two classes of traits. Similar results were found in the dairy herds of Canada, United States, Denmark and Australia, where improved productive performance also led to a fall in the reproductive potential of herds (VaRaden et al., Reference VanRaden, Sanders, Tooker, Miller, Norman, Kuhn and Wiggans2004; Sewalem et al., Reference Sewalem, Kistemaker and Miglior2010; Sun et al., Reference Sun, Madsen, Lund, Zhang, Nielsen and Su2010; Haile-Mariam et al., Reference Haile-Mariam, Bowman and Pryce2013).
The negative genetic gains observed from the year of birth of 2000 for bulls (with practical effects 6 to 8 years after) and from the year 2008 for cows (with their first calving about 2010), may be indirectly associated with the economic scenario during that period. Starting in the year 2008–2010, Portugal underwent a strong economic recession, affecting all sectors of the economy, including the importation of genetics (frozen semen and heifers). Silva et al. (Reference Silva, Silva, Silva, Costa, Lopes, Caetano, Thompson and Carvalheira2019a) reported a drop in the genetic tendencies for productive traits for this period. Therefore, the response observed in the genetic trends for reproductive traits could be actually a correlated response from the selection on the productive traits.
Although the genetic improvement observed in the reproductive traits during the last years, the expectation is that this scenario will not remain given the tendency to prioritize the improvement of productive traits that will end up causing deterioration on the animals' reproductive performance. In general, reproductive traits tend to have low heritability estimates, indicating a strong influence of management and environmental effects on the respective phenotypic expressions. However, due to their impact on the profitability of dairy farmers, different breeding programs have considered the reproductive traits in the selection process, and genetic gains have been observed with the inclusion of these traits within selection indexes (Gernand e König, Reference Gernand and König2017). In this sense, the planning of strategies for inclusion of reproductive traits in the genetic evaluation process becomes of great importance for the Portuguese dairy farming, by making possible genetic gains for the productive traits while at the same time softening the negative effects on the reproductive performance.
In conclusion, results from the application of AR models outperformed the REP models in evaluating time reproductive traits for all statistical fitting criteria providing lower AIC estimates, AW close to unit and lower MSE values. AR model also provided higher EBV reliabilities, which for selection purposes, will permit a better discernment among the best bulls. Thus, we conclude that AR model represents a viable approach for genetic evaluations of longitudinal reproductive traits in dairy cattle, which correspond to more reliable results and consequently, faster genetic progress.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0022029919000931.
Acknowledgements
The authors acknowledge the Portuguese Dairy Cattle Breeders Association for providing data for this study. This study was funded by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES/FCT, no. 99999.008462/2014-03 and 88887.163433/2018-00), Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq 465377/2014-9 – PROGRAMA INCT).
