Genetic analysis of litter size and number of kids born alive across parities in Spanish goat breeds using a random regression model

Abstract Genetic parameters of litter size (LS) and number of kids born alive (NBA) were estimated across parities in Florida and Payoya goat breeds using a threshold random regression model (RRM). We analysed the reproductive records for the Florida and Payoya breeds separately, and a total of 130.849 and 67.478 reproductive records from the first to the seventh parities from Florida and Payoya, respectively, were included in the analysis. Random regressions on Legendre polynomials of standardised parity were included for permanent environmental and additive genetic effects. We based our estimation of all the covariance components on Bayesian inference under categorical distribution and considered heterogeneous residual variances in the model. The estimates of heritabilities ranged from null to 0.37 for LS and from 0.02 to 0.20 for NBA, while the repeatability estimates were between 0.008–0.60 and 0.03–0.44 for LS and NBA, respectively. Phenotypic correlations within-trait along parities ranged widely from −0.31 to 0.97, while the genetic correlations within-trait through parities showed a wider range and magnitude of values, from −0.98 to 0.99 for LS and from −0.72 to 0.99 for NBA. The estimates of genetic correlations between LS and NBA in the different parities were positive in the first two parities and negative in the later ones. Due to the noticeable changes in variance and the incomplete genetic correlations, the use of a random regression model to analyse LS and NBA is recommended in goat breeds. HIGHLIGHTS A bivariate threshold random regression model has been used for the first time to estimate genetic parameters of litter size and number of kids born alive across parities in Spanish goat breeds. The evolution of genetic parameters along parities is different between Florida and Payoya that have different production systems, being more intensive in the case of the Florida breed. Estimates of variance were not constant across parities for both traits and breeds. Genetic correlations within-trait across parities were different from unity. Thus, LS and NBA should be treated as different traits in successive parities. Given that the estimates of variance and covariances are not constant across parities for both traits and breeds, the use of RRM is highly recommended.


Introduction
Over the last years, genetic improvement in Spanish dairy goats has focused on milk production and milk composition traits (Men endez-Buxadera et al. 2008Molina et al. 2018). Nevertheless, there are other important economic traits in small ruminants, such as litter size and number of kids born alive, especially for breeds facing extinction like the Payoya breed. In the Florida breed, they already represent a significant source of income for breeders by producing a large average litter size. Florida and Payoya breeds are native Spanish dairy goats reared under different production systems varying from semi-extensive to semi-intensive for Florida and from semi-extensive to extensive for Payoya. Florida is spread all over the country (Andalusia, Castilla la Mancha, Extremadura, etc) and in different countries in the world (Italy, Portugal, Venezuela, etc) whereas the Payoya breed is present only in the Andalusian Spanish region. Several approaches have been used to carry out genetic evaluation of litter size in goat, the most common being repeatability and multiple-trait models. A repeatability animal model assumes that genetic correlations between litter sizes are equal to the unity and that the variance is constant in successive parities (Mrode and Thompson 2005). Nevertheless, it has been argued in several studies using multiple-trait model and random regression model (RRM) that phenotypic and additive genetic correlations between litter sizes in adjacent parities are often lower than unity (Alfonso et al. 1994;Irgang et al. 1994;Fern andez et al. 2008;Skorput et al. 2014), especially between the first and the later parities (Serenius et al. 2003;Oh et al. 2005;Skorput et al. 2014). Multiple-trait model and random regression model consider subsequent measures to be different traits. Nevertheless, a random regression model has some advantages in comparison with multiple-trait model: the possibility of estimating covariance structure and breeding values at any point of the trait trajectory, and the reduction on the number of parameters to be estimated (Kirkpatrick et al. 1990;Meyer and Hill 1997). Also, RRM offers more accurate modelling of covariance structure, giving more accurate predictions of breeding values (Huisman 2002).
In goat, RRM has been used for analysing growth and feed intake (Barazandeh et al. 2012 in Raini;Kheirabadi and Rashidi 2016 in Markhoz;Desire et al. 2018 in mixed-breed Saanen, Alpine, and Toggenburg dairy goats) and milk yield and composition traits (Men endez-Buxadera et al. 2008 in Payoya;2010 in Murciano-Granadina;Andonov et al. 2013 in Norwegian goats;Silva et al. 2013 in Alpine;Arnal et al. 2019 in Saanen); however, their application to genetic evaluation of litter size has not yet been implemented, and estimation of genetic parameters for litter size or number born alive using RRM has been reported only in pigs (Lukovic et al. 2003(Lukovic et al. , 2004(Lukovic et al. , 2006Fern andez et al. 2008;Chen et al. 2010;Skorput et al. 2014;Ogawa et al. 2019). Litter size is measured more than once in an animal's lifetime and, despite being different from longitudinal traits because the parity number is a discrete variable, the application of RRM to analyse this trait is feasible (Schaeffer 2004).
Therefore, the main aim of this study was to estimate for the first time in Florida and Payoya goat breeds genetic parameters and their changes across parities for litter size (LS) and number of kids born alive (NBA) treated as threshold traits using a bivariate random regression model (RRM).

Data description
The reproductive and genealogical records of the Florida and Payoya goats were provided by ACRIFLOR (National Association of Florida Goat Breeders) and ACAPA (Association of Payoya Goat Breeders), respectively. The original datasets comprised 143,002 and 72,514 records, respectively. Data from herds with less than 10 records or scarce genetic connections were removed. Genetic links were due to the use of artificial insemination sires. After data editing, a total of 130,849 records from 48,984 Florida females collected between 1986 and 2019 (an average of 2.67 kiddings per female) and 67,478 records from 14,942 Payoya females recorded between 2003 and 2019 (an average of 4.51 kiddings per female) were used in this analysis. Traits considered were litter size (LS), estimated as the total number of kids born per kidding and number of kids born alive (NBA), both recorded on the day of kidding from the first to the seventh parity in both breeds. For the Florida pedigree, only animals with approved parentage tests using DNA microsatellites were considered. The pedigree was traced back for as many generations as available in the breeds herd book, being 7.3 equivalent generations for Florida and 4 generations for Payoya. The total number of animals in the pedigree was 55,748 for Florida and 26,392 for Payoya.

Statistical analysis
The analyses were performed for the Florida and Payoya breeds separately, fitting the following bivariate random regression model (RRM) for LS and NBA: where y is a vector with n observations of LS and NBA, 'fixed' is the combined effect of herd-year-month of kidding (2629 levels for Florida and 803 for Payoya), bk are k coefficients of fixed regression of the LS and NBA traits across parities with orthogonal secondorder Legendre polynomials. Seven parity classes were considered as the unit of time t for both traits and breeds corresponding to parity numbers (t ¼ 1-6 and t ! 7). aik andlik are the random regression coefficients associated with the additive genetic function of animals and the permanent environmental function of the goat with LS and NBA records, respectively, and e is the residual modelled with heterogeneous variance according to parity with three measurements for Florida (parity 1, 2, and ! 3) and five measurements for the Payoya breed (parity 1, 2, 3 and 4, 5 and ! 6). X, Z, and W are incidence matrices relating y with the new parameters in b, a, and m, respectively.
To deal with the non-normal distribution of LS and NBA, these traits were treated as a threshold trait, divided into two categories in both breeds (one; more than one).
The (co)variance components, heritabilities (h 2 ), repeatabilities (t), and additive genetic correlations (r g ) throughout the scale of values of parity were estimated by the expression proposed by Jamrozik and Schaeffer (1997).
Estimation of all the (co)variance components was based on the Bayesian inference methodology using the Gibbs Sampling algorithm in the THRGIBBS3F90 software (Misztal et al. 2016). Flat priors were assumed for systematic effects of herd-year-month of kidding, as well as for additive genetic, permanent environmental, and residual effects. Chains of 100,000 samples were used, with a burn-in period of 20,000. One sample per 100 was saved to avoid high correlations between consecutive samples. A post-Gibbs analysis using the POSTGIBBSf90 program (Tsuruta and Misztal 2006) was performed in order to calculate posterior means and high posterior density intervals and to check convergence with the Geweke test.

Results and discussion
Up to now, there have been no studies on the estimation of genetic parameters for litter size with random regression model in goat: previous studies on the genetic analysis of litter size have normally used repeatability and multiple-trait approaches. In this work, we analysed the change in genetic parameters of LS and NBA across female parities with a threshold random regression model in Florida and Payoya goat breeds.
The distribution of number of records and average female age at different parities in Florida and Payoya breeds is provided in Table 1. The least-squares means for litter size (LS) and number of kids born alive (NBA) over the seven parities are summarised in Figures 1  and 2, respectively. The average values were 1.61-1.6 and 1.6-1.57 for LS and NBA, for Florida and Payoya breeds, respectively. The maximum values of LS and NBA were observed in the fourth parity in Florida and in the sixth parity in Payoya. The phenotypic trajectories of both traits (Figures 1 and 2) reflected the high reproductive performance of Florida and Payoya breeds and their high persistence. However, the initial increase and the final decline in both traits in the Florida breed were more remarkable and occurred earlier than in Payoya, possibly due to its high production level and the fact that the Payoya breed is a more rustic (less improved) breed and is, therefore, less precocious.
The results of variance components (additive genetic, permanent environmental, and residual), heritabilities (h 2 ), and repeatabilities (t) across the seven parities in Florida and Payoya breeds for litter size and number of kids born alive are presented in Tables 2 and 3, respectively. All estimates of variance components, heritabilities, repeatabilities, and additive genetic correlations in this study for LS and NBA in both breeds were presented on a liability scale. Estimates of variance components for LS in both breeds presented the same tendency, decreasing from the first parity to the second parity which presented the lowest values, and increasing up to the latest parities, reaching the highest values in the sixth and seventh parities. Heritability and repeatability estimates of LS in both breeds showed the same pattern as the variance components. Their estimated values oscillated between null to 0.22 and 0.008 to 0.42 in Florida and from 0.01 to 0.37 and 0.02 to 0.60 in Payoya for heritability and  repeatability, respectively (Table 2). Variance components, h 2 , and t for NBA in both breeds did not present a clear pattern from the first to the seventh parity. Nevertheless, in the Payoya breed, the first parity presented the lowest estimates. Besides, the highest values for variance components, h 2 , and t were observed in the latest parity (h 2 was 0.19 and 0.20 and t was 0.30 and 0.44 in Florida and Payoya, respectively, Table 3).
The estimated variance components and h 2 for LS and NBA changed over parities and differed slightly between Florida and Payoya breeds, showing low values, as expected for reproductive traits. In general, all HPD 95% for heritabilities for LS and NBA in both breeds were small.
To date, there are no works on the estimation of genetic parameters of reproductive traits in goat using threshold RRM. In previous studies conducted in other goat breeds, a wide range of h 2 values has been reported for LS (Kasap et al. 2013     there are no available studies on estimating its genetic parameters, and this work is the first to estimate it in goat using RRM. The high estimated repeatabilities observed in the last parities could be explained by the fact that at this stage, the least reproductively efficient animals have already been discarded. The differences found between Florida, Payoya and other breeds may be attributable to estimation errors due to differences in the used methodologies, the size of the datasets, the genetic structure of the populations, or management practices. Estimates of the phenotypic (r p ) and additive genetic (r g ) correlations within-trait over the seven parities in Florida and Payoya breeds are provided in Tables 4  and 5 for LS, and Tables 6 and 7 for NBA, respectively. Phenotypic correlations within-trait along parities ranged widely from À0.26 to 0.96 for LS and À0.24 to 0.96 for NBA in Florida (Tables 4 and 6, respectively) and from À0.31 to 0.97 for LS and À0.30 to 0.97 for NBA in the Payoya breed (Tables 5 and 7, respectively). The estimated genetic correlations for LS were generally similar in both breeds and covered a very wide range of values varying from À0.90 to 0.96 in Florida (Table 4) and from À0.98 to 0.99 in Payoya (Table 5). These estimates were high between the adjacent parities (in the diagonal) and close to unity, except for r g between the second and third parities in both breeds, and decreased as the interval between parities increased, except for r g between the seventh parity and the first and second parities, which showed a small increase. In both breeds, the decrease in r g from the first-second parities to the second-third parities was sudden, with clearly contrasting values. For NBA, estimates of r g varied from À0.66 to 0.99 in Florida (Table 6) and from À0.72 to 0.99 in Payoya (Table 7). Estimated r g among the first five parities, unlike LS, presented only positive values. In subsequent parities, the values were high, except for r g between the fifth and sixth parities (0.28 in Florida and 0.17 in Payoya, Table 7). Unlike LS, the variation in r g as the interval between parities increased did not present a clear pattern in both breeds.
A low genetic correlation within-trait in different parities indicates that LS and NBA records from different parities should be treated as different traits (Irgang et al. 1994;Roehe and Kennedy 1995;Skorput et al. 2014). Direct additive genetic correlations were highest between adjacent parities and decreased as the interval between parities increased. This decreasing tendency of r g with the increasing distance between parities was in agreement with the findings of Lukovic et al. (2006) and Skorput et al. (2014) in pigs. This could be explained by the smaller amount of data in later parities and the imprecise recording of parity order (Skorput et al. 2014).
Additive genetic and phenotypic correlations between LS and NBA in Florida and Payoya breeds across parities are presented in Table 8. The estimated r g between both traits covered a wide range and magnitude of values, being positive in the first two parities with the highest values in the first parity for Florida (0.81, Table 8) and the first and second ones in Payoya (0.93 and 0.66 for first and second parity, respectively, Table 8), and negative in the later parities. Therefore, a selection made at any parity will have a positive or negative impact on the other trait Table 5. Estimates of additive genetic correlations (above diagonal) and phenotypic correlations (below diagonal) between litter size (LS) across parities in Payoya breed using a threshold random regression model.  throughout the rest of the productive life of the goat, particularly after the third parity. On the other hand, the phenotypic correlations were much smoother compared to the genetic ones. The wide range and magnitude of r g estimates between LS and NBA at different parities implies that selection for increasing one of these traits will produce a genetic response (limited given the low magnitude of heritabilities) whose magnitude will depend on the moment it is selected. The maximum direct and correlated responses can be obtained by selecting during the first two parities in Payoya and the first parity in Florida. It is also important to recognise litter size as a complex trait (P erez-Enciso et al. 1994, Rauw et al. 1999Men endez-Buxadera et al. 2004). LS is the best natural index as it combines ovulation rate and survival rate of the embryos (directly dependent on the maternal uterine environment). Litter size is therefore considered as a trait genetically related only to the mother of the litter, although some authors (Men endez- Buxadera et al. 2004) have shown that there may be a small paternal component (probably due to the fact that there is an additive component of the embryo or embryos that make them more or less robust), consistent with Falconer's (1960) hypothesis.
Random regression model has become commonly used to analyse longitudinal data or repeated records of individuals over time. RRM is a suitable way of dealing with repeated traits since it describes performance over time for each animal to allow for individual variation in the course of the trajectory (Meyer 1998). To estimate the genetic parameters for LS and NBA in Florida and Payoya goat breeds, the use of RRM is justified since both traits show a very wide variation in variances and genetic correlations throughout parities. Even for lower values of h 2 , as in the case of reproductive traits, and for negative genetic correlations, the possibility of optimally selecting for increasing female prolificacy is suggested.
The wide range of estimates of genetic correlations for LS and NBA between the different parities confirms that these traits are partially controlled by different genes at the first two parities and at the later parities. Genetic selection of animals for reproductive traits could be implemented involving LS and NBA in the first and second parities in Florida and in the first parity in Payoya.
Finally, this study provides estimates of genetic parameters for LS and NBA across the female parities, using a random regression model that can be used in the future to support decisions about selection aims within the current selection programs of Florida and Payoya breeds. It also suggests that future research should evaluate the traits of litter size and number of kids born alive in conjunction with the current selection criteria (milk production and milk composition traits).

Conclusions
Using a threshold model in which a continuous variable underlying LS and NBA traits across parities by a bivariate RRM in Florida and Payoya goats, our results indicate that litter size and the number of kids born alive in subsequent parities are genetically not the same traits because estimated variances and covariances were not constant over parities and genetic correlations between repeated measures of the traits differed from unity. Hence, the use of a random regression model to genetically analyse litter size or the number of kids born alive over parities is recommended in both breeds. Although low heritability estimates were obtained, improvement of litter size and number of kids born alive by selection is feasible. Given the genetic correlations between both traits across parities, an early genetic response can be achieved if selection occurs in the two first parities for Payoya and the first one for Florida breed.

Ethical approval
This study did not require manipulation or modification of the usual handling of the animals, since we have worked directly with the routine records provided by the breeders' associations of the Florida and Payoya breeds. Table 8. Estimates of phenotypic (r p ) and additive genetic (r g ) correlations between litter size (LS) and number of kids born alive (NBA) across parities in Florida and Payoya breeds using a threshold random regression model. ACAPA (Asociaci on de Criadores de Raza Caprina Payoya -Association of Payoya Goat Breeders) for providing the data for this study.