Non-genetic effects affecting fertility traits in local Reggiana cattle

This work is a preliminary study to develop a selection index for fertility traits in the local Reggiana cattle breed. The study aimed to investigate non-genetic sources of variation for fertility and to identify the best model in terms of variance explained. Moreover, the variation of the target fertility traits in the different parities and months of first service have been considered. The fertility traits under investigation were the interval between calving and first insemination, the interval between calving and conception, the number of inseminations per conception, and the calving interval. The dataset included 22,731 records of 10,502 cows, collected between 1986 and 2019. Four different models were tested: Model 1 included the fixed effects of herd-year-month of first service and parity; Model 2 separately accounted for herd-year and month of first service, in addition to parity. Additionally, Model 1a and Model 2a presented the same effects of Model 1 and Model 2 with the addition of the age at first insemination as linear covariate. The best fertility performances were observed in March and April, whereas for parity effect the best performance was in third lactation. Regarding the Model 1, the coefficient of determination was on average 0.40 for the four traits, and this value increased to 0.48 in Model 1a. Moving to Models 2 and 2a, the coefficient of determination dropped to the average values of 0.21 and 0.27, respectively. Therefore, Model 1a had the best fitting performances for all studied traits, although data editing for age at first insemination was very strict. On the other hand, in Models 2 and 2a the variance absorbed by fixed effects was low, and this could potentially bias the final estimates. Model 1 was therefore the best one in terms of trade-off between data loss and predictivity.


Introduction
Reggiana is a local Italian cattle breed mainly reared in Emilia Romagna region (North of Italy). During the 17 th century, Reggiana was the most common breed in this area. After the Second World War the number of animals decreased because of the substitution with cosmopolitan Holstein and Brown Swiss breeds, and it reached 985 heads in 1981 (Forabosco et al., 2011). In the last decades the number of cows slowly increased (Forabosco et al., 2011), as a result of the commitment of some passionate breeders, who adopted strategies to increase the intrinsic value of Reggiana milk within the PDO cheese named "Parmigiano Reggiano delle Vacche Rosse". Currently, more than 2,700 animals are reported in Italy (including 1,200 lactating cows), mainly in the province of Reggio Emilia (https://www.regionalcattlebreeds.eu/breeds/Reggiana.html).
The Reggiana breed is characterized by red color coat with dark or pale changing, toned down in the inner parts and in the lower limbs, at the eye's edges, around the face and the tail. The average 305days milk yield, protein content and fat content are 5,700 kg, 3.41% and 3.70% (http://bollettino.aia.it/). The milk of Reggiana breed is particularly suitable to produce Parmigiano Reggiano cheese, due to its high casein percentage.
Due to the limited number of animals and the recent interest for this breed, the development of selection strategies for Reggiana breed has started in 1996 (https://www.razzareggiana.it/en/). However, selection for fertility traits has not been implemented yet. The application of selection plans for fertility could imply lower costs for animal management due to lower replacement rate (Groen et al., 1997).
The aim of this study was to investigate non-genetic fixed factors affecting fertility traits in Reggiana local cattle breed. The identified factors could be then included in genetic models to estimate genetic parameters and breeding values for fertility of Reggiana cattle.

Data and editing
Fertility traits were analyzed combing information from two data sources, the insemination dataset and the test-day dataset, both provided by the National Reggiana Cattle Breeders Association (ANABoRaRe, Mancasale Reggio Emilia, Italy). The insemination dataset contained repeated events (n = 59,904) of 20,537 cows collected from 1986 to 2019, and the test-day dataset contained 29,291 records of 10,693 cows collected in the same period during the routine controls of milk production performed by the Italian Breeders Association (Rome, Italy). In addition to milk traits, the test-day dataset included information on herd, date of calving and parity of the cows. The whole herd book of Reggiana breed, dated back to 1943, was provided by ANABoRaRe.
Animals with an inconsistent date of birth were removed. Also, when multiple dates of calving for the same parity were present, only the one that guaranteed a calving interval of 287 ± 5 days was retained. In case the cow changed herd during a lactation, the herd in which she spent more time was considered as the herd effect for that lactation. After these checks, the test-day and the insemination datasets were merged according to the following steps.
From the test-day dataset, one row per each parity of a cow was created. In this row, the date of calving of the target lactation was considered as the starting date, while the date of subsequent calving was considered as the ending date. From the insemination dataset, the events that occurred between these two dates were joined. If the date of a subsequent calving was absent, all insemination events that occurred after the date of calving were joined. Insemination events that occurred before the first date of calving were also joined and considered as data for heifers. Subsequently, the insemination events showing a distance greater than 305 days with the subsequent date of calving were removed. First insemination events per each lactation that occurred outside the interval 20 to 180 days after calving were also removed.
Then, the following fertility traits were calculated on the edited data: i) interval between calving and first insemination (P_1), which represents the renewal of the estrus cycle; ii) interval between calving and conception (days open, DO); iii) number of inseminations to achieve pregnancy (NI); and iv) calving interval (CI), i.e. the interval between two consecutive calvings. Cows of parity higher than 6 were removed, as well as records without herd information. Following González-Recio and Alenda (2005), animals with age at first calving outside the interval 18 to 40 months were also removed. Calving intervals shorter than 300 days or longer than 600 days were discarded from the final dataset.
After the previous data editing, the merged dataset consisted of 22,731 lactations of 10,502 cows. A final editing was performed separately for each phenotype (four traits, but two following the same data editing process) and for each model (four models). Each phenotype had indeed a different number of data, apart DO and P_1 (Table 1). This was due to the different way phenotypes were calculated. The CI had the lowest number of data because it was necessary to retain the information of two consecutive dates of calving for each record, while for DO and P_1 this was unnecessary. The NI exhibited the highest number of data because it was possible to calculate this phenotype also for heifers. In addition to this editing by phenotype, an editing for each of the four models considered in the study was done separately for each trait (apart DO and P_1, included in the same dataset as explained above). This was done to have data with at least two records per each categorical effect included in the models described in the next chapter.

Statistical analysis
Descriptive statistics of fertility traits were calculated, and the normality of the phenotypic distribution was checked. As previously mentioned, four different models were tested for each trait using a linear model analysis in R environment (R Core Team, 2017) including different combinations of categorical fixed effects. Parity was retained as fixed effect in all 4 models, whereas the other effects differed, as reported below. Age at first insemination was considered as linear covariate. The distribution of data at different ages at first insemination is depicted in Figure 1. The least squares means of the cross-classified effects included in the models were also calculated. The four models considered were: Model 1: y ijk = μ + NL i + HYM j + e ijk , Model 1a: y ijkl = μ + NL i + HYM j + AGE_1 k + e ijkl , The number of levels for each cross-classified effect included in the models is reported in Table 1.

Descriptive statistics
Descriptive statistics of the fertility traits by parity are reported in Table 2. Overall, the third parity had a bit lower average values for P_1, DO, and CI, which means better reproductive performance. Traits P_1, DO, and CI were not highly skewed, therefore they can be considered normally shaped. However, NI was highly skewed (values close to 3), as more than 66% of cows had one service per conception. Skewness and kurtosis were also analyzed within each parity ( were slightly lower in the present study because our dataset included more recent data and animals showed negative (i.e., favorable) phenotypic trends in terms of fertility performance, despite the absence of an actual selection for fertility. The value of age at first insemination within lactation was consistent with the study of Gandini et al. (2007) performed on the same breed.

Linear model analysis
Results of the analysis of variance for the four models are summarized in Table 3 (Table 3 -last page of this article), whereas coefficients of determination (R 2 ) for each model are reported in Table 4. The inclusion of age at first insemination among the effects already included in Model 1 and in Model 2 increased the R 2 of Model 1a and Model 2a, but the high correlations between age at first insemination and all the 4 traits suggest that the use of this effect is questionable. Looking at the single phenotypes, CI was the trait with the highest R 2 in all 4 models (40%, on average). The other traits had an average R 2 of 0.32 considering all models.
Least squares means for the parity and month of first service effects are plotted in Figure 2 and Figure  3, respectively. The best results for DO and CI were observed for the second and third parity. On the contrary, in cosmopolitan breeds such as Holstein and Brown Swiss, the best performances were obtained in first-parity animals due to high metabolic demand of milk production in higher parities (Tiezzi et al., 2012). This does not happen in Reggiana since it is a rustic dairy breed, not showing high levels of daily milk production. Figure 3 suggests that in March and April reproduction performances were at their best (97.1 days, 1.4, 374 days, and 78.3 days for DO, NI, CI, and P_1, respectively). This is in accordance with Cavestany et al. (1985), who reported that the month affects reproductive performance for reasons related to photoperiod and temperature.

Fertility in Reggiana and future perspectives
This study aimed to carry out a preliminary analysis on fertility traits in Reggiana breed through the evaluation of the quality of data editing and making some considerations about fixed effects to be accounted for in future genetic analysis, phenotypic distribution and predictability of different models.
As expected, this breed exhibited good performances in terms of fertility compared with the specialized Italian dairy breeds. Gandini et al. (2007) also reported this discrepancy considering Reggiana and Italian Holstein. They found that NI averaged 1.52 ± 0.06 in Reggiana and 1.70 ± 0.07 in Holstein. Reggiana in the study of Gandini et al. (2007), presented values of 85.70, 103.66, and 381.13 days for P_1, DO, and CI, respectively. In Holstein, the same traits averaged 96, 122.52, and 414.43 days. As mentioned before, it is interesting to note that even if no selection to improve fertility was applied during the years in Reggiana cattle, reproductive performances have slightly improved over time. However, emphasis on fertility traits in the aggregate selection index is essential from an economic point of view in addition to cow well-being. The setting up of a good model for fertility traits in the Reggiana breed is more complicated than in other breeds, due to the limited data available. The number of animals in the population is much lower than in other breeds and, furthermore, the collection of phenotypes regarding fertility has started more recently than in other cosmopolitan or local breeds.
Removing 40% of animals could be a problem in terms of variance components estimation, since important information can be lost, and biases introduced because of the reduced sample size. Even if the majority of studies treated the herd-hear-month or herd-year-season as fixed effects (e.g., Van Bebber et al., 1997, Tiezzi et al., 2011, 2012, they could be considered as random. In general, fitting this effect as random may be advantageous for the estimation of breeding values, however, fixed time effects like this one should be included in the model to adjust for time trends (Schaeffer, 2018).

Conclusions
In conclusion, all fixed effects, especially the herd-year-month of first service, explained a significant amount of variance. Better reproductive performances were observed in second-and third-parity cows and in March and April. Overall, this study demonstrated the difficulty to implement a good model for complex phenotypes as fertility traits in a local breed as Reggiana. However, a good compromise between a rigorous data editing and a good R 2 is necessary. For this reason, Model 1 can be considered as the best model among the different ones analyzed in the study. Model 1a presented a higher proportion of variance explained by the factors included, but the data editing reduced too much the sample size. On the contrary, in Model 2a the editing was not too restrictive, but the amount of variance explained resulted low, although a low R 2 is common for functional traits like fertility.