Modelling for the inheritance of endangered equid 2 multiple births and fertility : determining risk factors 3 and genetic parameters in donkeys ( Equus asinus ) 4

Multiple births or twinning in equids are dangerous, undesirable situations that 37 compromise the life of the dam and resulting offspring. However, embryo vitrification and freezing 38 techniques take advantage of individuals whose multiple ovulations allow flushing more fertilised 39 embryos from the oviduct to be collected, increasing the productivity and profitability of such 40 techniques. Embryo preservation is especially important in highly endangered populations such as 41 certain donkey (Equus asinus) breeds; for which conventional reproductive techniques have 42 previously failed. For instance, becoming an effective alternative to artificial insemination with 43 frozen semen to preserve the individuals’ genetic material. The objective of this study was to 44 examine the historical foaling records of Andalusian donkeys to estimate genetic parameters for 45 multiple births, assessing the historical foal number born per animal, maximum foal number per 46 birth and multiple birth number per animal. We designed an Animal Model with single records 47 considering the fixed effects of birthyear, birth season, sex, farm, and husbandry system, and age as 48 Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 5 November 2018 doi:10.20944/preprints201811.0084.v1 © 2018 by the author(s). Distributed under a Creative Commons CC BY license. a linear and quadratic covariate. Restricted maximum likelihood reported heritability estimates 49 ranging from 0.18±0.01 to 0.24±0.01. Genetic and phenotypic correlations ranged from 0.01±0.01 to 50 0.83±0.01 and 0.12±0.01 and 0.53±0.01, respectively. These estimates enable the potential for selection 51 against/for these traits, offering a new perspective for donkey breeding and conservation. 52


Introduction
The occurrence of multiple births has been addressed one of the principal causes of fetal and neonatal loss in equids [1].The majority of twin pregnancies in horses (73%) terminates in abortion or stillbirth of both twins from eight months to term.The likelihood of one or both twin foals to be born alive or to survive after birth complications reduces to around 11%.The foals are usually born stunted or emaciated, what does not allow them to survive further from 2 weeks of age [1].
In the case of the donkey species, Quaresma et al. [2] addressed the overall neonatal mortality for the first month of life to be of near 9%.These authors would also report that the percentage of twin foaling at full term was only around 3%, with a neonatal foal mortality rate of 40%.Hence, the selection of individuals that may be less prone to present multiple ovulation could be a preventive alternative to decrease the risks attached.
By contrast, the donkey is a species for which the most of its breed populations have been classed as endangered [3] and that has been reported to be highly reproductively compromised as it happens to many other endangered populations [4] what may be attributed to the deleterious effects of inbreeding in such populations [5].The long gestation cycle (a norm of 12 months to give birth in the 13 th month [6]), a fertility that steadily decreases over generations [2] and the highly inbred status of donkey breed populations [5,7] only contribute to worsening the endangerment risk situation that donkey breeds face worldwide.Furthermore, highly standardised reproduction techniques in horses and other equids [8] such as artificial insemination with frozen semen and embryo transfer still represent a challenge in donkeys [9][10][11].
Under this context, embryo vitrification and freezing arise as new possibilities that may enable the preservation of the genetic material of donkeys belonging to populations for which the numbers rarely exceed 1000 individuals.This is supported as the pregnancy rates of 50% and 36% after the transfer of fresh and vitrified embryos, respectively [12], overcome the best currently reported results for pregnancy rate (28%) obtained for uterine horn insemination using frozen-thawed semen [13].
The efficiency of such reproductive techniques could be improved relying on the higher ability of certain animals to develop multiple ovulations, even more, when those animals may be genetically prone to develop them more often.
Studies of the genetic background of multiple pregnancies are anecdotal as fertility, in general, has a very low heritability as common to other reproductive traits.These studies are even more limited when we focus on studying equids such as horses [14,15] or donkeys, for which no study has been reported.This was a retrospective study over a period of 38 years (the birth year of the oldest animal registered in the studbook was 1980) aimed to investigate the occurrence of multiple pregnancies in the historical population of Andalusian donkeys and the influence that non-genetic factors such as farm, husbandry system, year of birth, birth season, birth month or age may have had on the prevalence of these multiple gestations.Second, we estimated the genetic parameters of fertility and multiple births through the analysis of historical foal number born per animal, maximum foal number per birth and multiple birth number per animal using REML methods.Last, we estimated predicted breeding values for all the traits as a way to assess the potential implementation of a bidirectional breeding scheme.This scheme may simultaneously consist of animals selected against multiple births because of the gestation complications that they involve, while other individuals may promote the occurrence of multiple births, seeking higher conservation profitability based upon an increased number of embryos to collect while implementing assisted reproduction plans.

Sample size and background
We studied the foaling recordings of 765 individuals registered in the historical pedigree record of the Andalusian donkey breed (181 jacks and 584 jennies).As age range was not normally distributed (P≤0.01 for both Kolmogorov-Smirnov and Shapiro-Wilk's tests for normality), we used minimum, Q1, median, Q3 and maximum to describe the age range in our sample.Minimum age in the range was six months, Q1 age was six years, median age was ten years, Q3 age was 14 years, and the maximum age was 29 years.Such wide age range was considered, given the fact that we assess an endangered breed from which the information belonging to each individual is indispensable, using Restricted Maximum Likelihood REML methods and after including age of birth in our model to correct for the cases in which the animals were too young to have given birth to their first foal/s.
The youngest age at which both jacks and jennies gave birth for the first time was three and four years old, respectively [5].It is often a decision of owners in particular not to breed the animals until they have been recognised as apt for reproduction and included in the studbook of the breed what takes place when the animals turn 3 years old.
The donkeys in the sample were the progeny of 93 jackstocks and 253 jennies.All the donkeys were registered in the breed's Spanish studbook.The relationships in the pedigree of the breed are routinely genetically tested through microsatellite-assisted genotyping and parentage tests for the resulting offspring of each mating.Parentage tests for each offspring had previously been performed with microsatellite molecular markers to ensure the reliability of the information in the pedigree as a way to counteract the small size of the sample tested.The DNA used for parentage tests was obtained from hair samples that are routinely taken when the inscription of each new animals takes place and from the historical bank of samples of the breed kept at the laboratory of applied molecular genetics of the University of Córdoba.All tests were carried out using a pedigree file provided by the Union of Andalusian Donkey Breeders (UGRA).The pedigree file included 1017 animals (272 males and 745 females) born between January 1980 and July 2015 from which only 914 donkeys, 246 males, and 668 females, were alive during the development of the study.
The pedigree of the donkeys in our sample was traced back six generations providing indirect information from 930 connected ancestors (91% of the historical population registered) and reporting an average inbreeding of 0.7%.

Birth related traits
First, we studied multiple birth occurrence assessing three different traits for each jack or jenny.
To obtain this information, we contrasted the registries of the pedigree file with interviews with the 145 owners whose animals participated in the study.We decided to interview the owners as due to the early abortion of multiple gestations it is very likely for them not to be included in the pedigree.
Second, from this initial sample of owners, we only considered the ones who affirmatively responded to the question in block 2 for the estimation of genetic parameters (90 out of 145 owners interviewed) as a veterinarian or theriogenologist had issued an official gestation diagnosis (simple or multiple).
This excluding criterion was applied as a way to consider those cases when abortions had occurred.
Many twin (and triplet) pregnancies in equids are already lost at very early stages and the aborted material stays mostly undetected by the owners, what could have distorted the true number of pregnancies with multiple conceptuses.
First, we summarized the historical foal number born per animal.That is for the total of 765 individuals, the number of offspring foaled (resulting from natural mating or artificial insemination) by each of 584 jennies or born to each of 181 jacks, either over their reproductive lifetime or up to July 2015.Second, the maximum foal number per birth, or the maximum number resulting at any of all the deliveries through the life of each jenny, considering which jack was used to breed.That is, for the same 765 animals, the maximum number of offspring born in a single foaling event in which the individual (male or female) was part of either over its reproductive lifetime or up to July 2015.Third, multiple birth number per animal, that is for the same 765 animals, the sum of all mating events resulting in multiple gestations either over the reproductive lifetime of the individual (male or female) or up to July 2015.
The units of study considered for descriptive statistics and populational data were the births occurring in the Andalusian donkey population and their characteristics.Hence, births were only computed once per couple of breeding animals to avoid redundancies, instead of computing them twice (one per jack and one per jenny involved).However, for genetic analyses, the unit of analysis that we considered was the lifetime parentship record of each animal separately.That is to say; we summed every molecularly confirmed jack and jenny's birth registries separately as a single value for each so that for the data considered reliable and irredundant.Given the BLUP methodology was applied [16], data can either belong directly from field observations and registries or indirect, where individuals are directly genealogically linked to common individuals.

Interview description
A telephone survey was carried out to 145 different owners whose farms were located in Andalusia (southern Spain).The survey took place in June 2017.We interviewed owners regarding the specific foaling registry of all the animals historically present at their farms since the 1980s until 2017 and registered in the stud-book of the breed at the moment that the survey took place.The oldest donkey from which there was information available had been born in 1984.All the interviews comprised a battery of 18 questions that were asked by the same interlocutor and each interview lasted for a mean time of 10 minutes.Despite the lack of prevalence of multiple births or gestation in their farms stated by the owners, all the questions were asked indistinctly.A description of the questions and options asked to the owners is shown in Table S1.Table S2, defines the categories (levels) of the husbandry system level.There were open questions (regarding the location of the farms, the age of the animals or the number of animals present in the farms at the moment that the interviews took place) and closed questions (regarding the sex, the husbandry system under which the animals were handled, and the prevalence of this condition from the Past up to the date when the interview was performed).All the information provided by the owners was contrasted with the information provided by UGRA and the information present in the official stud-book of the breed.

Records description and scales
We organized the questions into three blocks (Table S1).The first block aimed at describing the farms of the owners' interviewed in order to statistically assess the possible effects that may condition the prevalence of multiple gestations or births.We included the questions asked to the owners to classify the husbandry system under which their farms were managed in Table S2.These questions based on the extension of territory to which the donkeys had access, whether the donkeys were reproductively handled and whether the owner held daily contact with them or they were handled just for minimum punctual health inspection and stud book inclusion.The second block comprised a single question related to whether the diagnosis by a veterinarian or theriogenologist had been requested.The second block comprised an excluding question as only the owners affirmatively responding to it were included in the statistical and genetic analyses.The third block consisted of questions regarding the assertive diagnosis of the multiple births, and the care and preventive measures taken in each case.When the animals had never given birth, had suffered from an undetected early embryonic loss nor had carried any embryo, we gave them a score of 0.

Previous statistical analysis (screening)
The average number of foals born per year was 28.19, reaching the highest number (71) in 2003.
The mean prevalence of multiple births per hundred births in the Andalusian donkey population was 9.73%.The 11.18% of the population had not given birth to any foal when the registries were studied.The proportion of single, twins and triplets' pregnancies detected (all triple pregnancies were interrupted) was 58.28%, 30.46%, and 0.08%, respectively.The percentage of females with progeny selected for breeding was 10.76% and 25% for males in the historical population.Historically breeding jacks were 2.98 years older than breeding jennies on average.The average age of parents when their offspring was born was 8.08 years (8.03 for jennies and 8.16 for jacks).The average generation interval or the average age of the parents at the birth of their offspring that in their turn will produce the next generation of breeding animals for the historical population was 7.40 years [5].
A Shapiro-Wilk test was applied to the data to check the fitness degree of the variables in the model to a normal distribution.Second, the highly statistical significance of all the elements in the model (P<0.001),revealed that the data significantly deviated from a normal distribution.Kurtosis values supported these results (Table S3).Thus, we carried out a cross-sectional study employing Chi-square analysis to determine whether the categorical independent effects of birth year, birth season, birth month, sex, farm/owner, and husbandry system and the covariate of the age may randomly influence the dependent variables of historical foal number born per animal, maximum foal number per birth and multiple birth number per animal.We performed a Kruskal-Wallis H test to study the potentially existing differences between levels of the same factor except for age, as it is measured on a continuous scale (Table 1).We present Kruskal Wallis H Ranks for all the levels of the factors affecting historical foal number born per animal, maximum foal number per birth and multiple birth number per animal in Table S4.Simultaneously, we studied the pairwise comparisons for any dependent variables for which the Kruskal-Wallis test was significant, aiming at assessing whether there were differences between groups of the same factor affecting them.We used the Mann-Whitney U Test for sex, as it only has two levels, jack and jenny, and Dunn's test for the rest of the factors.
If we test multiple comparisons (hypotheses), the likelihood of incorrectly rejecting a null hypothesis increases, that is rejecting the existence statistically of significant differences between two or more groups (for instance, making a Type I error).The Bonferroni correction compensates for that increase by testing each hypothesis (comparison) at a significance level of α/m, where α is the desired overall alpha level, and m is the number of hypotheses.
Once we test for the differences in the distribution of the levels for each category, an independent-sample median test was carried out to assess the differences in the median between levels within the same factor.
After conducting a Kruskal-Wallis H with three or more groups (k), we computed the strength effect of the factors on the variables tested.Kruskal-Wallis H produces chi² values with k-1 degrees of freedom.Then, we can transform chi² into an F value with k-1 numerator degrees of freedom (dfn) and N-k denominator degrees of freedom (dfd) using the expression F(dfn,dfd)=chi²/(k-1), modified from Murphy et al. [17].From F(dfn,dfd) we calculate partial eta squared [18] from the following equation ηp²=(F•dfn)/(F•dfn + dfd).
Eta was computed to measure the strength of association between each categorical independent factor from the first set with the numerical dependent variables of historical foal number born per animal, maximum foal number per birth and multiple birth number per animal using the Crosstabs procedure from SPSS Statistics for Windows, Version 24.0,IBM Corp. (2016) (Table 1

) [19] [20] [20] [20] [20] [20] [20] [20] [20] [20]
[20] [20].Similarly, for the case of age, Spearman's rho was computed to measure the strength of association between it and the numerical dependent variables of historical foal number born per animal, maximum foal number per birth and multiple birth number per animal using the Bivariate procedure from SPSS Statistics for Windows, Version 24.0,IBM Corp. ( (Table 1).The Spearman's correlation, denoted by the symbol ρ (rho), is a nonparametric measure of the strength and direction of association that exists between two variables measured on at least an ordinal scale.The test is used for either ordinal variables or for continuous data that has failed the assumptions necessary for conducting the Pearson's product-moment correlation, as age in our study (Table 1
Categorical regression (CATREG) was used to describe how the variables in our study depended on the factors considered (Table 2 and 3

Genetic model, phenotypic and genetic parameters
As we only considered one measure per animal, the model used in the analysis of birth related variables was a simple Animal Model with single records.The fixed effects that were submitted to the above described statistical procedures and comprised the mixed model consisted of the birth year (from 1984 to 2017); birth season (summer, spring, autumn and winter); birth month (from January to December); sex (jack or jenny); the farm (90 farms/owners) and husbandry system (intensive, semiintensive, semi-extensive and extensive).
At a previous stage of the study, we computed the double interaction between herd and year of birth (Herd*Birthyear) and the triple interaction between the herd, the year of birth and season of birth (Herd*Birthyear*Birth season) as these were the most regularly included in literature for the same kind of studies in other species such goats or sheep.Adjusted R-squared was computed both including and without including the interaction, and we include a summary of the results in Table S4.We computed expected prediction error of regression with 0.632 Bootstrap ("leave-one-out bootstrap") from 200 bootstrap samples [20,21].Adjusted R-squared tells better than R square (percentage of variance explained by the model) about whether we should add a new variable or not as it tries to correct for the fact that the improvement in R square could be attributed to chance alone, considering the ratio (N-1)/(N-k-1) where N = number of observations and k = number of variables (predictors).The triple interaction was statistically nonsignificant (P>0.05)so that it was not included in the model.Although, the Herd*Birth year double interaction was statistically significant P<0.01, its inclusion within the model distorted the results in the following way so that we decided not to include such interaction.The model for historical foal number born per animal explained a higher percentage of the variance in the sample when we included the interaction.However, the estimation of the genetic parameters reported almost three times the standard error of the same model without including the interaction as stated below, that may have its basis on the high amount of possible levels of the interaction matched to a proportionally small sample.For maximum foal number per birth, there was a reduction in Adjusted R squared from 0.421 to 0.406 and the expected prediction error increased from 0.113 to 0.198 when we included the Herd*Birth year interaction.For multiple birth number per animal, one or more levels for the categories were supplementary (only used by cases that do not occur in the sample).These results suggested that the inclusion of this interaction in the model may result in potentially distorting effects which were highlighted at the statistical level as expected prediction error could not be computed.The results of the genetic and phenotypic parameters estimated by a preliminary model including Herd*Birth year iteration supported such distorting effects, as there was an increase in the standard errors from the current model (without the interaction) 0.081 to 0.128 to 0.154 to 0.643 (including the interaction).As the previous statistical analysis had reported, the basis for such distorting effects may be the fact that the number of categories considered for herd*year interaction was 441, while the whole sample size was 765.This data may generate an statistical imbalance that may result in an overestimation of the effect of the interaction as it has been reported by literature [22], making it impossible to test for its effects properly, due to the lack of enough animals in the pedigree between whom to compare.
The multi-trait animal model used for the analyses is as follows: where Yijklmnop is the separate record of ith trait for jth donkey (historical foal number born per animal (1 in matrix below), maximum foal number per birth (2 in matrix below) and multiple birth number per animal for a given donkey (3 in matrix below)); μ is the overall mean for the trait; aij is the additive genetic effect of the jth donkey for ith trait, Yeak is the fixed effect of the kh birth year (k=1984-2017); Monl is the effect of the lth birth month (January, February, March, April, May, June, July, August, September, October, November, December); Seam is the fixed effect of the mth birth season (m=summer, spring, autumn, winter); Sexn is the fixed effect of the nth sex (n=jack, jenny); Faro is the fixed effect of the oth farm/owner (o=1-90); Sysp is the fixed effect of the pth husbandry system (p=intensive, semi-intensive, semi-extensive, extensive); b1 and b2 are the linear and quadratic regression coefficients on age when the tests took place (Aq and A 2 q) and eijklmnop is the random residual effect associated with each record.No maternal effect was computed because of the low completeness level found in the pedigree, as 53.36% of the dams in the study were unknown [5].Such a lack of information could have represented a problem when performing the genetic analyses.
However, as our sample provides direct or indirect information from 91% of the animals included in the pedigree, we could save the possible drawback meant by the missing information.Then, the quality of the predicted genetic values estimated was quantified by reporting their reliability.
We included the age of the animals expressed in years as a linear and quadratic covariate to correct the variables measured according to the lifetime of each animal and specifically the cases in which the animals were too young to have given birth to their first foal/s.We included the effect of sex on our model to save the imbalance between sexes, even more, when we consider the vast differences between the offspring of males and females given the long duration of the gestation of the species.
In matrix notation, the multi-trait model used was: where y1 to y3 represent the phenotypical observation for each trait and animal.The vectors of fixed effect for the three different traits considered (β1 to β3,) include all the effect related in the model described above and the vectors α1 to α2 and ε1 to ε2, are random additive genetics and residual effects for each trait, respectively.The incidence matrices X1 to X3 and Z1 to Z3 associate elements of β1 to β3 and α1 to α2 with the records in y1 to y2.
If A is the matrix of additive genetic relationships among individuals, the mixed model equation (MME) used is as follows: Proxies of prolificacy (i.e.number of offspring produced in a single parturition) are calculated as sums over random time periods eventually censored by nature, and/or the will of the owner, and/or the timeframe of the study (each donkeys' lifetime, especially in animals that are too young to have given birth).Hence the importance of including, assessing and controlling factors such as owner and age of birth as reported above.

Genetic assessment software
The first stage of our analyses was to obtain estimates of genetic parameters and correlations (phenotypic and genetic) for historical foal number born per animal, maximum foal number per birth and multiple birth number per animal for a given donkey in Andalusian donkeys.Estimations were performed using the MTDFREML Multiple Trait Derivative Free Restricted Maximum Likelihood software [23], a package employing the simplex procedure to locate the maximum of the log likelihood (logL).We considered having reached convergence when the variance of function values (-2 logL) in the simplex was less than 10 -12 .We run the model several times restarting the program with initial values from the previous round (run) to ensure convergence of logL to a global maximum.
When the -2 logL and estimates were similar in successive analyses, we concluded that convergence reached the global maximum.After we reached convergence, the second stage aimed at predicting breeding values (BLUP, best linear unbiased predictors for animal random effects), their accuracies and reliabilities for each animal sorted by sex.The analyses were carried out including the relationship matrix of animals with direct records related through at least one known ancestor, considering the 1017 donkeys in the historical pedigree.Covariate and fixed effects posterior means (BLUES, best linear unbiased estimators for fixed effects and covariable) were computed with MTDFREML, as well.The standard errors of genetic correlations among fertility and multiple birth traits were obtained directly from the MTDFREML analyses.

Interview results
Out of the 145 owners interviewed, we considered the information from 90 farms/owners.These owners had affirmatively responded to the question in the second block as they were they only who had requested information concerning diagnosis by their veterinarians or theriogenologists and therefore, were the only ones providing reliable information.Due to the particularities of the species and the breeding routines carried by the owners, the artificial insemination with fresh semen of the animals registered in the studbook was infrequent, and almost all the matings were performed naturally.No productive artificial insemination using frozen semen was registered.The matings of only 66 animals out of the 765 donkeys from which there was information (8.63% of the total sample) had resulted in multiple gestations.Out of this percentage, 1.04% of the animals developed multiple gestations in more than one occasion through their lives and only one of the animals was responsible for 0.13% of multiple gestations in the population (five multiple births out of 40 births through his life).
Shapiro-Wilk Test (P<0.001) and higher or lower kurtosis values than three on all the fixed effects, the covariate and interaction showed that they highly significantly did not fit a normal distribution.The variability observed for the two traits analyzed was from moderate to high, with a coefficient of variation of 21.3% for the husbandry system effect and 82.2% for the effect of the farm/owner.

Statistical analyses
The results of Chi-Square, Eta (for each independent categorical-dependent numerical pair of variables) and Spearman's rho correlation coefficient (for the effect of age on the numerical dependent variables studied), testing for the existence of linear correlation are shown in Table 1.Eta effectively and statistically significantly measured the strength of collinearity that the sex and farm factors have on continuous variables of historical foal number born per animal, maximum foal number per birth and multiple birth number per animal for a given donkey.Husbandry system reported highly statistically significant collinearity with the historical foal number born per animal (Table 1).
Kruskal-Wallis H test and Chi square reported the effects birth year, birth season and birth month to be statistically nonsignificant (P>0.05) for the three dependent variables considered.The same test reported the rest of independent variables (sex, owner/farm and husbandry system) to be statistically significant (P<0.05) for all dependent variables except for husbandry system on maximum foal number per birth and multiple birth number per animal for a given donkey (P>0.05)(Table 1).
From the results of the Mann-Whitney U Test, we can conclude that historical foal number born per animal and maximum foal number per birth in jacks was statistically significantly higher than in jennies (U=46363.500,P<0.001 and U=50364.000,P<005).However, the opposite trend was described by multiple birth number per animal for a given donkey, which was statistically significantly higher in jennies than in jacks (U= 47730.000,P<0.05) (Table S5).
The results of the Dunn test in our study reported the fact that there were highly statistically significant differences for 14.89% of pairwise comparisons of farms/owners (Table S5).The same test reported statistically significant differences between extensive, semiextensive and semiintensive husbandry systems (P<0.05) for maximum foal number per birth.
CATREG was performed to the 7 qualitative independent variables (birth year, birth month, birth season, sex, farm, husbandry system) and age as a covariable with the three birth-related continuous variables (historical foal number born per animal, maximum foal number per birth and multiple birth number per animal for a given donkey) as dependent variables.Then stepwise linear regression to the data with the resulted quantifications was applied, and we present the summary results with the significant variables in Tables 2 and 3.The standardized coefficients (β) are listed in Table 3. CATREG reported all of the independent variables except for birth year and sex to be significant for historical foal number born per animal.Sex was nonsignificant for maximum foal number per birth.Birth season and husbandry system were nonsignificant for Multiple birth number per animal.
There was a small to moderate monotonic (whether linear or not) significant (P<0.05)correlation between age and the three variables tested (Table 1).For the month of birth, Chi square values were non-significant (P>0.05).Thus, there was not any statistical difference between the values of the dependent variables for each of the twelve levels of For sex, Chi square values were all significant (P<0.05),thus there were statistical differences between the values of the dependent variables for each of the two levels of sex.Eta values ranged from 0.090 to 0.240 what reported a low to moderately high association between sex and the dependent variables of maximum foal number per birth and multiple birth number per animal for a given donkey.CATREG standardized coefficient for sex and historical foal number born per animal and maximum foal number per birth were non-significant.However, CATREG standardized coefficients for multiple birth number per animal (0.445) reported a high increase of the standard deviation of sex was needed to increase a unit of standard deviation in multiple birth number per animal.
Owner/Farm, Chi square values, were all significant (P<0.001),thus there were highly significant statistical differences between the values of the dependent variables for each of the 90 levels of the farm/owner factor.Eta values ranged from 0.346 to 0.634 what reported a high association between owner/farm and the dependent variables of historical foal number born per animal and multiple birth number per animal for a given donkey.CATREG standardized coefficient for owner/farm were all highly statistically significant (P<0.001).CATREG standardized coefficients ranging from 0.598 to 0.873 reported a high increase of the standard deviation of owner/farm was needed to increase a unit of standard deviation in all three dependent variables measured.
Husbandry system Chi square value was only significant (P<0.001) for historical foal number born per animal, thus there were highly significant statistical differences between the values of that dependent variable for each of the four levels of the husbandry system factor.Eta value for this dependent variable was 0.187 what reported a moderately low association between husbandry system and the independent variables of historical foal number born per animal and multiple birth number per animal for a given donkey.CATREG standardized coefficient for owner/farm was statistically significant (P<0.05) for historical foal number born per animal and maximum foal number per birth.CATREG standardized coefficients ranging from 0.532 to 0.553 reported a moderately high increase of the standard deviation of the husbandry system was needed to increase a unit of standard deviation in all three dependent variables measured.
We show the factors affecting the three birth related variables in order of importance according to the CATREG standardized coefficients (β) in Table 4. Since we used the stepwise method, there was no multicollinearity problem.The standardized solution for the regression equations can be found in Table 4 as well.Non-significant effects for each variable were not included (P>0.05)

Model predictive power
CATREG R squared coefficient obtained ranged from 0.0.994 to 0.431 for multiple birth number per animal and maximum foal number per birth, respectively (Table 2)

Genetic model, variance components, genetic and phenotypic correlations, predicted Breeding Values and prediction accuracy
We show the estimates for heritability, genetic and phenotypic variance estimated with REML in Table 5.Table 6 shows the genetic and phenotypic correlation chart.The results for the estimates of predicted breeding values (PBV) for both jacks and jennies are shown in Table 7.

Covariate and fixed effects posterior means
We show the results for the best linear unbiased estimators (BLUEs) obtained from the REML quantitative genetic analysis, including age as a linear and quadratic covariate, the fixed effects of birth year, birth month, birth season, sex, farm/owner and husbandry system in Table S8.

Discussion
According to literature, donkeys have a 13% higher fertility than horses [24], reaching an incidence for multiple ovulations of 61% in Mammoth and standard jennies.This higher incidence of multiple ovulation in donkeys translates in twinning occurring more frequently, especially in mammoth and standard donkeys.Although researchers have found the incidence of twins to be as high as 40% (via ultrasound at day 21 in standard donkeys), for endangered donkey breeds such as Asinina de Miranda, the percentage of twin foaling at full term was 2.85% [2].According to our results, only 8.75% of the Andalusian donkey breed population has produced multiple offspring, and only around 10% of total births considered in the study (controlled) were multiple.The rate of multiple ovulations in the donkey species varies with the reports from literature, ranging from 5.3% to almost 70% [7] so that our results fall within the range reported for other donkey breeds.
The reproductive trends of this polygynous species have been reported to highly depend on the owner tastes for certain morphological or coat characteristics and local availability of the animals.
Navas et al. [5] suggested the typical excessive contribution of few ancestors to the gene pool of small critically endangered donkey populations may lead to narrow bottlenecks shortly whose hidden effects can only be controlled by tracking the populations.Among such hidden effects, the compromises exerted on the reproductive and immunity system of the animals have been addressed to be some of the determinants of the difficulties experimented to conceive by individuals [25].
Such reproductive compromises have been suggested to be a direct cause of inbreeding depression in donkeys.However, the lack of completeness of the pedigree of endangered donkey populations and the irregular distribution through great extensions of territory makes the estimation of this parameter little reliable [5].Quaresma et al. [2] reported the numbers obtained in 40 captive mammalian populations indicated an average value of 3.14 of lethal equivalents with 50% due to recessive lethal alleles.
Taberner et al. [26] stated that multiple ovulations tend to repeat in several estrous cycles, what may support the existence of animals that present a certain cyclical predisposition towards multiple births.The relative frequencies for multiple pregnancies of certain donkeys were higher than for others, what suggested a genetic background behind multiple births, as it had previously been reported by Ginther [27].Similarly, Quaresma et al. [2] suggested an indirect selection of certain family lines may have been carried out in the Mammoth donkey, what may have resulted in the higher incidence of multiple ovulations reported by Blanchard et al. [28].
Specific studies have assessed the possible repercussion of certain environmental factors on the fertility of donkeys.For example, in our study, the Chi square values for the birth season were nonsignificant (P>0.05).Thus there was not any statistical difference between the values of the dependent variables for each of the four levels of birth season.The findings by Contri et al. [29] support our results.These authors would report estrous cycle can be detected during the whole year in jennies, with no differences in the estrous cycle length among seasons.Paralelly, the pattern of the plasma concentration of certain hormones such as E2 and P4 during the estrous cycle did not report any difference among seasons, although a larger diameter of the ovulating follicle was reported for spring and summer.
Breeding season and month significantly affected gestation and estrous cycle length in donkeys [30].However, these authors did not study whether the effect of the month may condition the occurrence of multiple births and fertility.Quaresma and Payan-Carreira [31] reported the incidence of single, double, and triple ovulations to be 57.58%,36.36%, and 6.06%, respectively.The same authors stated, multiple ovulations affected neither the length of the interovulatory interval nor the individual cycle stages (P > 0.05) but lengthened the interval from the beginning of estrus to the last ovulation (P = 0.01), what may support the results found by our study and those found by Galisteo and Perez-Marin [30] as well.
No paper has reported the higher prevalence of multiple births or a higher likelihood of presenting a higher maximum number of foals depending on the husbandry techniques carried in the farms.The results found in our study for Dunn's and independent samples median tests suggested donkeys located at semiextensive farms presented a higher likelihood of presenting higher maximum foal numbers per birth, followed by semiintensive farms and extensive farms, respectively.
The criteria used to classify the husbandry systems of the farms in the study (Table 2) may suggest that the access to more extensive territories, when owners provide regular reproductive care to the animals and the daily contact with the owners may have an increasing importance in the occurrence of a higher number of foals per birth.The higher strength effect of the farm factor on all the variables tested ranging from 0.598 to 0.873, for multiple birth number per animal and historical foal number born per animal, respectively supported the finding.
A higher relevance was attributed to jennies in having a historically higher number of foals, a higher number of multiple offspring and a higher maximum number per birth.These values balanced (providing an equal relevance to jacks and jennies) as the number of foals and multiple births increased, as we can observe in the charts in Table S5.However, still there seem to be a very slight effect of specific jacks on promoting the obtention of a higher historical number of foals.This could be attributed to the reproductive characteristics of the jenny and breeding strategies of donkey owners, as it has already been suggested by Bresińska et al. [14] and is addressed by the results of the Mann-Whitney U test of our study (Table S5).According to our results, the fact that foal number born per animal and maximum foal number per birth in jacks was statistically significantly higher than in jennies could be attributed to the fact that jacks can act as the sire for several jennies at the same time, while jennies are going to be reproductively blocked for a whole year when they have become pregnant.The same test suggested that although jacks were likely to significantly reach a higher number of foals on a certain gestation through their lives when compared to jennies, jennies were statistically significantly more prone to develop multiple gestations through theirs.This could be supported by the greater chance of jacks to mate and the fact that multiple ovulations are a female trait, usually associated with endocrine changes that originate a sort of independence from the falling FSH values, that allow two (or more) dominant follicles to ovulate.
Using Restricted Maximum Likelihood REML methods, as we consider the relationship among the individuals present in the pedigree, we can estimate genetic information for the animals from which we have direct observations, and predict such information for animals assessing the additive indirect observations obtained from their ancestors.Hence, we can get the information for a particular trait of an individual when it is naturally impossible or potentially difficult to obtain it.For instance, prolificacy in foals that are too young to give birth, milk production from a male or when fertility rates are unbalanced between sexes (i.e., the number of offspring that a male can produce compared to the number of offspring a female can give birth to) [16].
Estimates of additive genetic variance for maximum foal number per birth and multiple birth number per animal for a given donkey were around the lowest margin of the values reported for twinning and fertility in horses.However, the models used in such circumstances differed from ours.
By contrast, the estimate of additive genetic variance for historical foal number born per donkey was around the highest margin reported for fertility in horses (Table 3), what resulted in higher heritabilities [32].The low genetic component of variance did not affect heritability estimates which were moderate and ranged from 0.18 to 0.24, for multiple birth number per animal for a given donkey and historical foal number born per donkey, respectively.Furthermore, these heritability values were very accurate as suggested by the negligible estimation error found (Table 5).Moioli et al. [33] found similar SE for the same parameters and traits in the Maremmana local cattle breed whose sample size was similar to the one in our study.Among the common factors to the two studies, microsatelliteassisted genotyping of the pedigree relationships may have played an essential role in the estimation of such reliable genetic parameters.
Several authors have suggested Bayesian inference Threshold models to be more suitable to analyse non-normally distributed functional traits in small samples [15,[34][35][36].However, when we analyse the results, the REML estimates tend to be included within the credible interval of the estimates obtained using Bayesian methods [32].This finding suggests more accurate estimates when the prediction error using REML methods is low as it shortens the credible interval.Hence, a previous analysis should be carried out testing for different models and using both threshold and REML methodologies to assess what is the one that fits better the statistical properties and biological characteristics of the situation that we want to model.
Our estimates for phenotypic and residual variance are almost 4 to 6 times higher than genetic variance estimates.As it has been reported in horses [32], the current analysis assumes that fertility and multiple births are determined by an infinite number of loci that contribute each with a minimal effect in what is called infinitesimal mode of inheritance.Hence, we can suppose, fertility may complexly depend on many physiological processes each of which is controlled by specific biochemical pathways.
We could have expected the high genetic, and the moderate phenotypic correlation between maximum foal number per birth and multiple birth number per animal for a given donkey as the fact that an animal is more prone to have multiple births may make it more prone to have a higher maximum number of foals per birth.We found low genetic and phenotypic correlations between maximum foal number per birth and historical foal number born per donkey.This finding may mean a weak relationship between animals having a high historical number of offspring through their lives and the same animals having a high maximum number per birth, what may suggest a lower reproductive life for those animals producing multiple offspring.Genetic and phenotypic correlations between the number of multiple per animal and historical foal number born per donkey were moderately high, what suggests the higher the number of total offspring through the life of a given donkey is (that is the more fertile), the more likely these animals are to produce multiple births.
These correlations have been described as well in humans [37,38].For instance, all the findings by Mbarek et al. [39] point to spontaneous twinning being a heritable trait and suggest the potential for polygenic inheritance.The genetic correlations found by our analyses support it as well.Mbarek et al. [39] reported that consistent with its effects on higher circulating FSH levels; the rs11031006-G allele also associates with higher total lifetime number of children.Boomsma et al. [40], reported an increased frequency of the S allele in fathers of dizygotic twins.However, this may be a secondary effect of assortative mating for family size.Assortative mating is a form of sexual selection in which individuals with similar phenotypes mate with one another more frequently than would be expected at random.The Andalusian donkey is a highly standardized breed for which assortative mating may have played an indirect role when seeking for obtaining specific phenotypical characteristics what may account for the low genetic variance for maximum foal number per birth and multiple birth number per animal for a given donkey.
Despite its demographic bottlenecks, the Andalusian donkey still maintains considerable levels of genetic variability for fertility and multiple birth traits [5].Given the favourable existing genetic relationships between the traits involved, these traits can play an essential role in a selection program aimed at improving the breeding efficiency of the animals.The potential opportunities arising from the incorporation of genomic information in the selection program should be investigated and implemented carefully in the future.Their contribution to reducing generation intervals and enhancing selection accuracy could result in extraordinary benefits for genetic progress, avoiding to detrimentally increase the inbreeding problems and endangerment risk from which the species suffers [41].PBVs for multiple births and fertility show considerable variability, indicating a possibly effective selection based on genetic merit objective estimates.The moderate heritability values balance the high existing phenotypic variability, resulting in a moderately wide PBV distribution (Table 7).Implementing a systematic genetic evaluation procedure through the genetic information available, allowing the early selection of breeding animals becomes then one of the main aims of the study.However, the reduction of generation intervals, enhancing selection accuracy through multivariate animal models for functional traits, and thus, the reduction in the number of breeding jackstocks to compatible levels with an increased selection response, must consider the detrimental problems that are likely to appear because of an increase in inbreeding in breeds with such a low effective population number.In these breeds, the protection of genetic variability and minimizing inbreeding are primary concerns as they may prevent population bottlenecks form occurring.The incorporation of genetic markers in the functional selection against or for donkeys for multiple births or fertility is a still a developing possibility.Hence, the exceptional importance of the implementation of these validated assessment tools and new methods and the perspective to develop routinely studies assessing the same animals over several years.

Conclusions
High values for genetic parameters enable the potential inclusion of these traits within breeding programs seeking the genetic progress of donkey breeds.Positive and moderate genetic correlations enable the combined selection for maximum foal number per birth and historical foal number born per donkey, with low detrimental effect for either one.Selection for multiple births or fertility in donkeys may have traditionally been carried out indirectly.Thus, the routine application of the assessment including a higher number of animals is required to standardize the valuation methodology implemented.However, this is a difficult task to achieve, considering the current extinction risk of donkey breed endangered populations.Functional traits related to fertility and prolificacy can play an essential role in a selection program aimed at improving the suitability of donkeys for their inclusion in embryo vitrification, or freezing assisted reproduction programs.The present results enable a bidirectional selection strategy.On one hand, the specific nature and the magnitude of the existing genetic relationships may make interesting to consider the possibility of developing and maintaining specialized lines relying on the ability of particular donkeys to develop multiple births within the Andalusian donkey breeding program, hence, increasing the productivity of assisted reproduction techniques.On the other hand, when embryo collection is not the purpose aimed at, selection could focus on the obtention of those individuals that may be less prone to develop multiple births, thus, avoiding the risks of multiple gestations, what in the end translates in the improvement of the reproductive welfare of the individuals.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1.Table S1.Questions asked to the owner interviewed during the survey carried out regarding the fertility and prevalence of multiple births in donkeys; Table S2.Description of the levels included in the husbandry system fixed effect; Table S3.
Descriptive statistics for fixed effects (yellow), interaction (green), covariates (red) and birth related variables (blue) in Andalusian donkeys (N=765); Table S4.Kruskal Wallis H Ranks for all the levels of the factors affecting historical foal number born per animal, maximum foal number per birth and multiple birth number per animal; Table S5.Summary of the results of the Independent-Samples Mann-Whitney U Test and Dunn's test and Bonferroni's significance correction for the effects for whose levels there were statistically significant differences on the historical foal number born per animal, maximum foal number per birth and multiple birth number per animal; Table S6.Summary of the results of the independent sample median test for the effects for whose levels there were statistically significant differences on the historical foal number born per animal, maximum foal number per birth and multiple birth number per animal compared to the significances for the results of the Dunn's test for the same parameters; Table S7.Comparison of the model summary of stepwise linear regression with transformed variables including and without included the interaction of herd*birthyear; Table S8.Estimates of non-genetic effects obtained from the REML quantitative genetic analysis, including age as a linear and quadratic covariate, the fixed effects of birth year, birth month, birth season, sex, farm/owner and husbandry system.
).The resulting regression equation could be used to predict historical foal number born per animal, maximum foal number per birth and multiple birth number per animal for any combination of the seven independent factors.Categorical Regression was carried out using the Optimal Scaling procedure from the Regression task from SPSS Statistics for Windows, Version 24.0,IBM Corp. (2016).

Table 1 .
Summary of the results for the Kruskal-Wallis H test for fixed effects and the covariate included in the model to test for birth related traits in Andalusian donkeys.

Table 2 .
Model summary of stepwise linear regression with transformed variables.

Table 3 .
Standardized Coefficients and significance of CATREG model.

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 5 November 2018 doi:10.20944/preprints201811.0084.v1 birth
month.Eta values ranged from 0.117 to 0.146 what reported a moderate association between month of birth and the dependent variables of multiple birth number per animal for a given donkey and historical foal number born per animal.The low to moderate CATREG standardized coefficients reported a low to moderate increase of the standard deviation needed to increase a unit of standard deviation of each of the three variables on this factor (from 0.092 to 0.170, for historical foal number born per animal and multiple birth number per animal, respectively).
For the birth season, Chi square values were non-significant (P>0.05).Thus there was not any statistical difference between the values of the dependent variables for each of the four levels of birth season.Eta values ranged from 0.064 to 0.095 what reported a low association between birth season and the dependent variables of maximum foal number per birth and multiple birth number per animal for a given donkey.CATREG standardized coefficient for the birth season and multiple birth number per animal was non-significant.However, CATREG standardized coefficients for maximum foal number per birth (0.100) and historical foal number born per animal (0.071) reported a low increase of the standard deviation of the birth year was needed to increase a unit of standard deviation in both dependent variables.

Table 4 .
Regression equations for maximum foal number per birth, multiple birth number per animal for a given donkey and historical foal number born per donkey.

Table 5 .
Estimated components of variance, heritability (h 2) and standard error (SE) for maximum foal number per birth, multiple birth number per animal for a given donkey and historical foal number born per donkey obtained from multivariate analyses using REML methods in Andalusian donkeys.

Table 6 .
Estimated phenotypic (rP) (above diagonal) and genetic (rG) (below diagonal) correlations for maximum foal number per birth, multiple birth number per animal for a given donkey and historical foal number born per donkey obtained in bivariate analyses using REML methods in

Table 7 .
Descriptive statistics of predicted breeding values (PBV), standard error of prediction (SEP), accuracy (rTi) and reliability (RAP) for maximum foal number per birth, multiple birth number per animal for a given donkey and historical foal number born per donkey for all the donkeys included in the pedigree sorted by sex.
a Standard error for Skewness statistic was 0.148 and 0.090 and standard error for Kurtosis statistic was 0.294 and 0.179 for jacks and jennies, respectively for all factors assessed.