Genetic analysis of reproductive efficiency in Spanish goat breeds using a random regression model as a strategy for improving female fertility

Abstract This study aimed to estimate genetic parameters of reproductive efficiency over a wide age range of females using a random regression model (RRM) in Spanish goat breeds. A total of 138,139 and 64,638 reproductive records from the first to the sixth parities from Florida and Payoya, respectively, were included in the analysis. Random regressions on Legendre polynomials of standardised age were included for permanent environmental and additive genetic effects. Estimation of the covariance components was based on the Bayesian inference using the GIBBS3F90 software. Differences among genetic variance components for reproductive efficiency were observed over the animal’s lifetime. The estimates of heritabilities were moderate, ranging from 0.21 to 0.32 for Florida and from 0.25 to 0.35 for Payoya, while the fractions of phenotypic variance explained by the permanent environmental effect were high, varying between 0.45–0.68 and 0.58–0.71 for Florida and Payoya, respectively. The correlations for permanent environmental effect over age ranged from 0.37 to 0.99, while the genetic correlations between the different ages varied from 0.36 to 0.98 for Florida and from 0.80 to 0.99 for Payoya. The results from this study support the validity of using an RRM to genetically analyse reproductive efficiency in Spanish dairy goats following the changes in variances and the genetic correlations different from unity over the animal’s lifetime. Moreover, reproductive efficiency is a highly heritable trait that is expressed early in a female’s life, and it could be used as a precocious selection criterion to improve female fertility in Spanish dairy goats. Highlights Reproductive efficiency (RE) has been proposed as a trait to improve female fertility in dairy goats. RE was genetically analysed using an RRM over a range of ages in females from the Florida and Payoya breeds raised under different production systems. Variance components for RE were not constant over age, while the genetic correlations between RE in the different ages were different from unity. The use of RRM for RE is therefore justified in both goat breeds.


Introduction
Female reproductive efficiency is an economically important trait for the dairy goat industry and is considered the most important factor to ensure high productivity for certain environmental conditions (Hossain 2004). Until recently, improving the production of milk, fat, and protein per doe, together with morphology, were the main target breeding traits in genetic improvement programs of Spanish dairy goats (Menendez-Buxadera et al. 2008;Men endez-Buxadera et al. 2010;Molina et al. 2018). Placing too much emphasis on production traits, while ignoring other traits, may lead to undesirable consequences in the animals' health and fertility, which decreases longevity (Oltenacu and Broom 2010). The classic antagonistic genetic-environmental relationship between milk production traits and female fertility has been reviewed in various studies and species (Andersen-Ranberg et al. 2005 in cattle; David et al. 2008 in sheep; Montaldo et al. 2010 in goat). The importance of including fertility traits in dairy goats breeding programs has therefore increased, as a way of correcting negative effects in this parameter due to the selection of highly productive females.
Fertility in does is a complex trait that can be defined as the doe's ability to resume its ovarian function after kidding, show detectable oestrus, become pregnant, maintain gestation and succeed at kidding; this results in a variety of measures being used in genetic evaluations across the world (Diskin and Morris 2008). To date, despite the economic importance of reproductive traits in dairy goats, there are very few studies to estimate their genetic parameters and none have carried out a longitudinal analysis with a random regression model. The most commonly used traits in the international genetic evaluation of female fertility are (1) age at first kidding (Garc ıa-Peniche et al. 2012;Castañeda-Bustos et al. 2014;Ziadi et al. 2021a) and (2) different kidding intervals (Bagnicka et al. 2007;Garc ıa-Peniche et al. 2012;Castañeda-Bustos et al. 2014;Atoui et al. 2018;Ziadi et al. 2021a). Nevertheless, their very low heritability estimates make it difficult to determine which traits must be used as a selection criterion to improve female fertility. A reason for the low heritability of reproductive measures considered as fitness-related traits is natural selection. The fundamental theorem of natural selection states: 'The rate of increase in fitness of any organism at any time is equal to its genetic variance in fitness at that time' (Fisher 1930). It has been construed that characters with the lowest heritabilities are those most closely connected with fitness, while characters with the highest heritabilities are those that might be judged on biological grounds to be the least important as determinants of natural selection (Falconer 1982).
Genetic parameters of fertility in Spanish dairy goats were estimated only in one study for the age at first kidding and the different intervals between kidding in the Florida and Payoya populations (Ziadi et al. 2021a) and, as expected for this kind of traits, the magnitude of estimates was relatively low. However, this study used a repeatability model and did not deal with fertility over a time trajectory.
In Spanish dairy goats, random regression models have primarily been used to analyse milk production traits (Men endez-Buxadera et al. 2010) and, more recently, prolificacy traits (Ziadi et al. 2021b). The basic idea underlying all of these models consists of modelling the additive genetic values (or other random effects in the model) as a function of an observed dependent variable (i.e. time or weight) through a set of random coefficients.
The purpose of this study was to analyse genetically for the first time reproductive efficiency as a new fertility trait in Florida and Payoya goat populations using a random regression model over a range of ages in females. These breeds are prototypes of the Spanish breeds exploited in different production systems varying from semi-extensive to semi-intensive for Florida and from semi-extensive to extensive for Payoya.

Phenotypic data and pedigree
The data consisted of reproductive records of Florida and Payoya females and genealogical information provided by the National Association of Florida Goat Breeders and the Association of Payoya Goat Breeders, respectively. The original datasets comprised 143,002 and 72,514 reproductive records, collected between 1986-2019 and 2000-2019, respectively. Animals from the Payoya breed are reared under semi-extensive to extensive production systems, whereas the Florida breed presents a higher level of intensification. The data editing consisted of keeping herd-year contemporary groups with at least 10 records and removing herds with scarce genetic connections. Genetic links between herds were ensured by the use of artificial insemination and the sale of bucks. After data editing, the final datasets consisted of 138,139 records from 50,992 Florida females (average 2.7 records per female) and 64,638 records from 21,715 Payoya females (average 2.97 records per female). Female fertility was defined as the reproductive efficiency (RE), calculated as the deviation between optimal and real parity number of females at each age. Only records of females with up to six parities were considered in this analysis. The pedigree was traced back for as many generations as available in the breed herd book, which was 7.3 equivalent generations for Florida and four generations for Payoya. For the Florida pedigree, only animals with approved parentage tests using DNA microsatellites were considered. Finally, the total number of animals in the additive genetic relationship matrix was 55,748 for Florida and 26,392 for Payoya.

Statistical analysis
The significance of the non-genetic effects for the herd, year of parity, and month of parity and their inclusion in the model were determined using the 'GLM2' R package. All of them had a significant effect on reproductive efficiency at a 0.05 significance level. Two models with second and third-order Legendre Polynomials were tested and compared, and for both breeds, the Deviance Information Criterion (DIC) indicated that both models were equivalent.
We analysed the RE records to fit the following second-order Legendre polynomials, random regression model, separately for both breeds: where y is a vector with n observations of RE, fixed is the contemporary group effect of herd-year-season of kidding (1260 and 584 levels for Florida and Payoya, respectively), in which the season of kidding was coded as 1 if a goat gave birth in the period June-September and as 2 if a goat gave birth in the period October-May, bk are k coefficients of fixed regression of the RE trajectory in females over age with orthogonal Legendre polynomials of second order. Six age ranges were considered as the unit of time t (t ¼ 1-5 and t ! 6) and their effects were considered as a continuous variable expressed in a standardised form between À1 and þ1 according to the expression from Kirkpatrick et al. (1990). The terms aik and lik are sets of random regression coefficients of second-order associated with the animals' additive genetic function and the permanent environmental function, respectively, and e is the residual. X, Z, and W are incidence matrices relating observations with the new parameters in b, a, and m, respectively. In this RRM model, the expected (co)variance components are assumed as: in which G is an additive genetic (co)variance matrix among all animals; P and R are respective (co) variance matrices among permanent environments and residuals, A is the relationship matrix between all animals in the pedigree; I p and I n are identity matrices of respective order number of animals with own record (p) and number of records (n), and is the Kronecker product. The Bayesian inference approach was applied to estimate variance components and genetic parameters using the GIBBS3F90 software of the BLUPF90 family programs (Misztal et al. 2016). Flat priors were used for non-genetic effects and variance components. Chains of 100,000 samples were used, with a burn-in period of 20,000. One sample in 100 was saved to avoid high correlations between consecutive samples. The saved samples were used to obtain the posterior mean and posterior standard deviation of estimates and to check convergence with the Geweke test, using the POSTGIBBSF90 program (Tsuruta and Misztal 2006).
The (co)variance components, heritabilities (h 2 ), the fraction of phenotypic variance explained by permanent environmental variance (c 2 ), permanent environmental (r pe ), and additive genetic (r g ) correlations throughout the range of ages in females were estimated using the expression proposed by Jamrozik and Schaeffer (1997).

Results and discussion
In this work, we analysed reproductive efficiency in females as a fertility trait in the Florida and Payoya dairy goat breeds, estimating its genetic parameters across a range of ages using an RRM.

Phenotypic parameters of reproductive efficiency
A summary of the basic statistics of reproductive efficiency in the different ages in the Florida and Payoya populations studied is shown in Table 1. A RE value <100% indicates that the female has registered a delay in having her real parity compared to the optimal parity of her age. In the first age group, the average RE value was lower in Payoya than in Florida, indicating that most Payoya females start their reproductive life late since they are reared in a more extensive production system. Afterward, the average RE value increases with age in both breeds, reaching maximum reproductive efficiency in the fourth and fifth age groups because the less fertile females are removed from the herds at each stage. The decrease in RE in the last age group can be explained by the growing number of females with reproductive problems at these final stages. In general, the Florida breed showed a higher RE than the Payoya due to its intensive production system. The preliminary GLM analysis showed that all the effects included in the model were statistically significant to greater or lesser degrees. The least-squares means for reproductive efficiency for the herd, parity year, and parity month effects in both breeds are illustrated in Figures 1-3, respectively. Figures 1(a,b) showed a wide variation in the reproductive efficiency of the different herds of these breeds, with the average RE ranging from 68.42 to 129.75% and from 60 to 118.74% in Florida and Payoya, respectively. This illustrated how increased herd fertility can be achieved through environmental improvement. Also, RE registered a decline at the end of the last century and the first decade of this century, followed by a notable improvement over the last few years when it exceeded 100% (Figure 2(a)). RE has also improved over recent years for the Payoya breed, although its magnitude is still below 100% (Figure 2(b)). Furthermore, Figure 3 indicates that there is a seasonal pattern marking the RE. The minimum RE values were recorded in the hot seasons, from June to September.

Genetic parameters of reproductive efficiency
Tables 2 and 3 show a summary of the changes in the posterior distributions of variance components (additive genetic, permanent environmental, and residual), heritabilities (h 2 ), and the fraction of phenotypic variance (c 2 ) explained by permanent environmental effect variance for RE over the lifetime of females for Florida and Payoya, respectively. The RRM allows us to study the change in (co)variance with age. Estimates of additive genetic and permanent environmental variances tended to increase continuously with age at parity in both breeds and their magnitudes were Reproductive efficiency (  clearly more important in the case of Florida. The permanent environmental variances were larger than the additive genetic variances for all ages, indicating that the permanent environmental effects had a higher influence on the variation of reproductive efficiency among does of both breeds than the additive genetic effects. Permanent environmental effects contribute permanently to all phenotypes, are constant across repeated measures on the same individual (Kruuk and Hadfield 2007), and are cumulative over time (Schaeffer 2011). This result implies that improvements in reproductive management (especially in the organisation of the first service) and health care should lead to a greater increase in fertility in both breeds than by making a genetic selection alone. The heritability estimates were moderate for both breeds: for the Florida breed, h 2 decreased slightly from 0.23 in the first age group to 0.21 in the second age group and then increased steadily until the last age group with 0.32, making an overall average of 0.26 (Table 2). In the case of Payoya, the estimates of h 2 showed a decrease from 0.35 in the first age group to 0.25 at the third and fourth intermediate ages, followed by a slight increase until the last age group, making an overall average of 0.28 (Table 3). In all cases, HPD 95% values for heritabilities for RE in both breeds were small. The variance components for additive genetic and permanent environmental were not constant across the range of ages in females, indicating that from a genetic point of view, RE should not be considered as the same trait throughout the life of the animal.
The fraction of phenotypic variance (c 2 ) explained by permanent environmental variance varied from 0.45 to 0.68 in Florida (an average of 0.62) and from 0.58 to 0.71 in Payoya (an average of 0.67). Estimates of c 2 were not constant across the trajectory of age at parity in both breeds, implying that non-additive genetic effects or permanent environmental effects may not affect RE similarly at all ages.
The reasons for the apparent observed discrepancy between the Florida and Payoya breeds may be attributable to the differences in management practices and production systems, with Florida being raised under more intensive management, in addition to the genetic structure of the populations. For instance, for the process of culling females, 24.7% of Florida females are removed between the first and second  parities compared with 14.8% in the Payoya breed, and then 52.7% of Florida females have more than six parities (78.6% in Payoya) (Ziadi et al. 2021a).
The posterior means of the additive genetic (r g ), as well as the permanent environmental (r pe ) correlations over a range of ages in females in the Florida and Payoya breeds, are shown in Tables 4 and 5, respectively. The estimates of genetic correlations between RE for the different ages were all positive, varying from moderate (0.36 between the first and sixth age groups, Table 4) to high (0.98 between the fourth age group and the third and fifth, Table 4) for Florida. In the Payoya breed, they were high, ranging between 0.80 for the first and sixth age groups and 0.99 for the different r g of the diagonal, except for the r g between the first and second age groups (Table 5). In all cases, r g between any two age groups decreased, with the greater distance between RE in young does and RE in older does being more noticeable in the Florida breed. This result suggests that, if a doe fails to produce a kid, her chance of having a consecutive kidding decreases with an increase in the interval between successive parities. In addition, the permanent environmental correlations were not constant over age groups in both breeds, and varied from 0.37 for the first-sixth age groups to 0.99 for the second-third age groups in Florida, with a similar pattern in the case of Payoya. These differences in r pe between different ages would indicate that environmental effects and or non-additive genetic effects (dominance and epistatic effects) are not constant with age at kidding. Hohenboken (1985) reported that some environmental effects affecting the repeated measures of the trait in the same animal seem to be semi-permanent rather than permanent. Until now, the estimation of nonadditive genetic effects for fertility traits has not been reported in dairy goats.     Random regression models (RRM) are often more appropriate for the genetic analysis of the traits of repeated records on the same animal since they describe the changes in the additive genetic effects and permanent environmental effects over the time scale adjusting for covariance between the different points of the trajectory (Meyer 1998).
To our knowledge, in the literature, there are no estimates of genetic parameters for reproductive efficiency, as defined in this study, or for the application of RRM for fertility traits, and this work is the first work to present them.
Comparisons with other studies are difficult because of differences in the definition of fertility, in addition to the notable absence of studies of reproductive parameters related to fertility using a random regression model. The h 2 estimates for RE in both breeds were clearly higher than those for the classical traits used worldwide for the genetic evaluation of female fertility, such as the age at first kidding and the different intervals between kidding. For example, Jembere et al. (2017) in their meta-analysis of several studies on dairy goat breeds reported h 2 values of 0.17 ± 0.012 for the age at first kidding, 0.002 ± 0.018 for the first kidding interval, and 0.09 ± 0.01 for the interval between all kidding. Similarly, in our previous study on estimating genetic parameters for the classical fertility traits in the Florida and Payoya breeds, we obtained very low h 2 estimates, varying from 0.01 for the interval between second, third, and remaining kiddings in Florida to 0.11 for the age at first kidding in Payoya (Ziadi et al. 2021a). Besides, the higher expected genetic responses calculated using selection indices occurred when the interval between all kiddings was combined with the age at first kidding. However, the interval between all kiddings is expressed later in a female's life and direct selection for this trait will widen the generation interval, and then reduce the annual genetic response. Nevertheless, reproductive efficiency is measured from the first kidding and could be applied as a precocious fertility selection criterion. Given its high h 2 and genetic correlations with next parities, a faster genetic improvement of fertility would be achieved in both breeds.

Conclusions
According to the results obtained from this study, the trait of RE is highly heritable and should be included in the breeding goal to improve dairy goat female fertility in Spanish goat breeds. However, reproductive efficiency in subsequent ages should not be considered genetically as the same trait, because estimated variances and covariances were not constant in the different ages and the genetic correlations between repeated measures of the trait differed from unity. In genetic evaluations for the fertility trait, RE should therefore be treated as a different trait in females of different ages. In the same way, the use of a random regression model to genetically analyse this trait in different age groups is recommended in both breeds. Moreover, the magnitude of h 2 in early age groups and the genetic correlations between RE in subsequent age groups are high enough to apply a precocious selection for RE to improve dairy goat female fertility, which can be used to overcome the limitations of selection using classical fertility traits.

Ethical approval
All information used in this study was obtained from existing datasets provided by the National Association of Florida Goat Breeders and the Association of Payoya Goat Breeders. Therefore, no Animal Care and Use Committee approval was necessary.

Disclosure statement
The authors declare that there is no conflict of interest associated with the paper.

Data availability statement
The data that support the findings of this study are available on request from the corresponding author with the Table 5. Estimates of genetic a (above diagonal) and permanent environmental (below diagonal) correlations for reproductive efficiency across a range of ages in females in Payoya breed using a random regression model. permission of the ACRIFLOR (National Association of Florida Goat Breeders).