Age-related linear and nonlinear modelling of semen quality parameters in Miranda donkeys

Abstract

In endangered donkey populations, reproductive performance is one the most important areas to study, provided its connection with conservation opportunities. Among reproductive traits, the knowledge on the evolution of seminal qualitative and quantitative parameters may achieve a particular relevance given its implication with improved or more efficient fertility and fecundation chances. The aim of the present study was to model the gonadosomatic index and the quantitative and qualitative seminal parameters along age in Miranda donkeys. Forty semen collections from 8 Miranda donkeys from 2.8 and 21.6 years of age, were fitted to eleven linear and nonlinear models. Cubic function modelling reported the best-fitting properties for almost all variables, except for sperm concentration (exponential model). Cubic model presented a slightly higher ability to capture interindividual variability than less parametrically complex models and more computationally complex models, such as exponential or logarithmic. Senility may course with and larger numbers of abnormal morphology spermatozoa, and also reduced sperm concentration and motility.

Highlights

• Testis reach mature size rapidly which may be ascribed to a breed factor.

• Age promotes ejaculate volume increase and concentration reduction.

• Reduced motility may ascribe to increased abnormal morphology spermatozoa.

• Major age-related changes are linked to testicular functions degeneration.

Keywords:

• Equus asinus
• flexibility criteria
• seminal evaluation
• testicular ultrasound
• equine reproduction

Introduction

Understanding and controlling the progression of seminal parameters in endangered populations becomes especially important in species that present a relatively long reproductive cycle (Fielding 1988; Tosi et al. 2013; Quaresma 2015). The results derived from related research may be of direct applicability and transference to veterinarians, equine practitioners, breeders or owners derived from the optimisation of assisted reproduction programmes developed under the scope of in situ conservation strategies (Quaresma et al. 2014; Navas et al. 2017).

The Burro de Miranda breed is one of the two officially recognised donkey breeds in Portugal, with its origin in the most northeastern region of the Country. The breed is endangered of extinction, with just around 780 females and 60 males registered. Although the breed is under a conservation plan and the number of animals has slightly increased in the last decade, the current reproductive rates are not high enough for changing the endangered situation and the population is ageing. The animals are distributed by around 450 owners. The average age of the owners is 61.8 ± 18.2 years old (SPREGA, 2020).

Knowledge of the normal progression of events before and after puberty is necessary if the maximisation of the reproductive handling in the species is sought (Miragaya et al. 2018). Furthermore, the determination and ability to predict sperm quality and its relationship with fertilising ability has a great economic importance to breeding herds when artificial insemination is used, as it may enable to perform the early selection of those donkeys not only presenting a better reproductive performance, but also to them having interesting phenotypical or biological features such as increased levels of diversity or reduced levels of inbreeding (Quaresma et al. 2015). However, these analyses are expensive and time‐consuming, and cannot be applied under field conditions, at least as efficiently as to be routinely performed at a commercial level.

Classical methods of semen evaluation generally measure seminal parameters as sperm concentration and volume, motility, percentage of morphologically normal and abnormal cells, among others (Miró et al. 2005; Canisso et al. 2011; Dorado, Acha, Ortiz, et al. 2013). However, these studies or evaluations are time-limited rather focus on explaining the trends or the findings of data at the moment that the individuals are sampled, while, efforts to predict for future seminal quality evolution related outcomes in males remains scarce (Quartuccio et al. 2011).

Mathematical regression applications can help to approach and solve this problem or at least can shed light on critical stages of seminal evolution along the life of particular males (Blanchard et al. 2013). This may not only be important from a conservation point of view but also from a clinical or biological perspective, which may enable providing an improved attention to reproductive management in donkeys, to choose the most appropriate moment for reproductive handling or to detect anomalies which may compromise the reproductive efficiency of individuals (Quartuccio et al. 2011).

The donkey is a reproductively efficient species for which reproductive artificial management has resulted to be more scientifically and biologically challenging than for other species (González et al. 2015). Under field conditions, the use of simple, non‐expensive and accurate statistical evaluations is required. This goal has been sought after over decades but still research has not been able to report a definitive solution.

Up to the date, no systematic study has approached the changes occurring in seminal parameters in the donkey from puberty until reproductive senectitude occurs. The present study aims to evaluate the trends described by seminal parameters from prepuberal to post puberal stages in Miranda donkeys. Software-automatised Linear and nonlinear regression analytical strategies were developed to determine the possible relationship between sperm parameters and age.

Materials and methods

Animals

Ultrasound assessment was performed on 23 clinically healthy Miranda jackstocks, with ages ranging from 7 months to 21.5 years and weighting from 120 to 400 Kg. Seminal parameters were evaluated on a random sample subset comprising 8 Miranda mature donkeys between 2.8 and 21.6 years old, with a live weight variation of 214–400 kg. Examinations were performed during the autumn-winter in 2018–19 and spring-summer in 2019. The donkeys considered in the study sample were housed at the Veterinary Teaching Hospital of UTAD (VTH-UTAD, Vila Real, Portugal). It was not possible to perform semen extraction in all the animals on which ultrasound was performed due to many reasons linked to the handling of the animals. The main limitation was the fact that the availability of the animals was subject to the fact that some of the owner required their animals back to their farms and that some of them were not trained to be collected. Hence, extraction could not be performed several times as to eliminate extragonadal reserves of sperm for an ejaculate to be considered representative. Furthermore, seasonality was not included because our aim was to determine the evolution of semen parameters in time, and our sample did not allow to infer the effect of seasons due to it being unequally distributed. Formal consent and collaboration of the Association for Study and Protection of the Donkey Breed Burro de Miranda (AEPGA) was received. The study was carried out under the framework of a scientific agreement of cooperation signed between both institutions. Clinical procedures were conducted following national laws for animal welfare and experimentation as with the European Directive 2010/63/EU for animal experiments and the Directive Hospital Committee (Approval Ref. 408/VTH-UTAD).

Clinical and genital examinations were performed previous to semen collection. Body weight (BW) (Kg) was recorded using an equine digital floor weighting scale (Salter Brecknell® PS-3000HD Floor – Veterinary Scale, USA). Only animals with both testis symmetrical and located at scrotal position, presenting a regular size respecting to their age, normally consistent and for which no echogenic changes in the testicular parenchyma, epididymis and spermatic cord were detected, were included in the present study.

Gonadosomatic index calculation

Ultrasonography examination of testicular measurements were performed using a Philips® CX30 Portable Ultrasound (Philips^®, Amsterdam, Holland) with a sectorial 3.0–7.0 MHz transducer, following the premises developed by Love (2013) . Longitudinal and transversal plans were performed on each testicle, excluding epididymis measurements. Three scan rounds were performed placing electronic cursors at the limit of tunica albuginea. Right and left testicle length (L), height (H) and width (W) (cm) were computed. Right and left testicular volumes (TV) were then calculated using the Lambert formula, TV = L×W × H × 0.5233 (Love, 2013). Total testicular volume (TTV) or the sum of the right and left TV was also individually calculated. In the calculation of the gonadosomatic ratio (GSI) (%), i.e. testicular weight/BW, the TTV (cm³) was transformed into grams, considering testicular density in mammals is very close to one (Johnson and Neaves 1981).

Semen collection and evaluation

A total of 40 ejaculates (5 ejaculates per donkey) were collected. All individuals accounted with a successful record of previous natural mating. Prior to experiment, collections were performed during three consecutive days to minimise the collection of sperm from extra-gonadal reserves (Quartuccio et al. 2011). During the experiment (autumn-winter in 2018–19 and spring-summer in 2019), semen collections were performed at 2 days intervals. and using an artificial vagina (AV) (Hannover model – Minitub Iberica S.L., Tarragona, Spain) lubricated with non-spermicidal gel (ReproJelly – Minitub Iberica S.L., Tarragona, Spain), using a jenny in heat as a mount. AV inner temperature was maintained at 40–45 °C temperature using a warm water filling, so as to enable sterile semen collection, a new bottle was used at each collection. A nylon filter was used to remove gel fraction by filtering the whole ejaculate (Minitub Iberica S.L., Tarragona, Spain).

Gel-free ejaculate was immediately evaluated for volume (mL), motility (%), concentration (×10⁶/mL), and percentage of morphologically normal (%). A scaled semen collection bottle was used to determine gel-free volume. Afterwards, sperm motility and concentration were evaluated for each ejaculate. For sperm motility evaluation, an aliquot of gel-free ejaculate was immediately extended 1:1 (vol/vol) with INRA 96 (IMV Technologies, L’Aigle, France) extender at 37 °C. Sperm motility was blind and subjectively estimated by the same trained operator after the evaluation of spermatozoa motility (%) in 5 different fields under light microscopy (× 200), placing a semen droplet in a prewarmed (37 °C) slide covered by a cover slip. Concentration was determined by using an improved Neubauer haemocytometer (Marienfeld, Lauda-Königshofen, Germany). Total sperm number (TSN, ×10⁹) was obtained by the product between the volume of gel-free ejaculates and sperm concentration, whilst total motile sperm count (TMS, ×10⁹) was calculated computing the product between motility and TSN. Sperm morphology anomalies (head, intermediary piece, tail) were evaluated in eosin-nigrosine stained smears using light microscopy in oil immersion objective lens (× 1000), counting 200 sperm cells (Kenney 1983).

Statistical analyses

Approach determination: pre-evaluation of parametric assumptions

Normality assumption was tested on our study sample. Shapiro–Wilk’s test was used to test for normality (for n < 50 samples) in the sperm parameters datasets using the Explore procedure of SPSS Statistics for Windows statistical program, Version 25.0 (IBM Corp 2017).

The Durbin-Watson test (Durbin 1970) was conducted on the residuals of each model to test for possible first-order autocorrelations among residuals (using the Linear regression test of regression procedure in SPSS version 25.0). The Durbin-Watson statistic ranges from 0 to 4. The Durbin-Watson (DW) test is reliable for sample sizes larger than 15 (Greenberg et al. 2020). Durbin Watson statistic is only suitable for ordered time or spatial series (Chen 2016) such as the evolution of testicular morphometry and seminal parameters in time. A value of DW equal to 2 is indicative of the lack of autocorrelation (critical value). If DW statistic is below this critical value, there is statistical evidence that the data is positively autocorrelated. If DW is over this upper critical value, there is no statistical evidence that the data is positively correlated. However, if DW is in between the lower and upper critical values, the test is leads to no conclusion (Gujarati and Porter 2003). Another assumption of ordinary least squares regression is that the variance of the residuals is homogeneous across levels of the predicted values, also known as homoscedasticity. Heteroscedasticity appears when the lower values of the standardised predicted values tend to have lower variance around zero. Durbin Watson and Scatter plots were performed using SPSS Statistics for Windows statistical program, Version 25.0 to test homoscedasticity. Residual values are computed after the result of the difference between observed and predicted values.

Considerations prior to regression analyses

When outliers are present, predictive power weakens and the violation of the assumption of normally distributed data. The identify outliers procedure of the Analyze/Built-in analysis of the column analyses package of GraphPad Prism version 8.3.0 was used to detect the likelihood of outliers that could distort fitting properties. To this aim the ROUT method was implemented as it combines Robust regression and Outlier removal, and can be used to fit a curve that is not influenced by outliers. The ROUT method is based on the False Discovery Rate (FDR). To this aim a level of Q coefficient must be predefined. This level is the maximum desired FDR.

The value of coefficient Q determines how aggressively the ROUT method defines outliers. The mathematical details are explained in reference (Motulsky and Brown, 2006). This value is set in the Weights tab of the Nonlinear regression dialog. Unless you have a strong reason to choose otherwise, the default value of 1% is recommended. If there really are outliers present in the data, Prism will detect them with a False Discovery Rate less than 1%. A Q coefficient of 1% was used (Motulsky and Brown 2006).

Mathematical models and curve shape parameters

One linear model and ten non-linear models were used to describe the evolution of calliper and ultrasound testicular measurements of the 23 donkeys and seminal parameters of 8 of these donkeys considered in the study sample. Table 1 shows the equations for the 11 models used, the abbreviation used to identify each model and the model syntax equations used. SPSS model syntax equations were designed to facilitate the automatised application of the models in this study. The model syntax presented is ready to be copied and pasted in the non-linear regression task from the Regression procedure of SPSS version 25.0 (IBM Corp 2017).

Table 1. Models and model syntax.

Model	Initials	Model syntax/equation
Linear	LIN	Y = b0 + (b1 * t)
Logarithmic	LOG	Y = b0 + (b1 * ln(t)
Inverse	INV	Y = b0 + (b1/t)
Quadratic	QUA	Y = b0 + (b1 * t) + (b2 * t**2)
Cubic	CUB	Y = b0 + (b1 * t) + (b2 * t**2) + (b3 * t**3)
Power	POW	Y = b0 * (t**b1) or ln(Y) = ln(b0) + (b1 * ln(t))
Compound	COM	Y = b0 * (b1**t) or ln(Y) = ln(b0) + (ln(b1) * t)
S-curve	SCRV	Y = e**(b0 + (b1/t)) or ln(Y) = b0 + (b1/t)
Logistic	LGI	Y = 1/(1/u + (b0 * (b1**t))) or ln(1/y-1/u) = ln (b0) + (ln(b1) * t)
Growth	GRO	Y = e**(b0 + (b1 * t)) or ln(Y) = b0 + (b1 * t)
Exponential	EXP	Y = b0 * (e**(b1 * t)) or ln(Y) = ln(b0) + (b1 * t)

t represents time expressed in year units. Accessed from IBM Corp (2015).

The curve estimation task from the Regression procedure of SPSS version 25.0 (IBM Corp 2017) was used to iteratively specify the parameter bounds of each model (b0, b1, b2, b3, and b4 parameters) using the Levenberg-Marquardt method of iteration. The iterative process considered as many rounds as it was necessary until a tolerance convergence criterion (error sum of squares between successive iterations (Dematawewa et al. 2007)) of 10⁻⁸ was reached as suggested by other authors, as stricter criteria such as 10⁻⁶ or 10⁻⁸ have been suggested to report the same outcomes out of a slightly higher number of iterations (Feistauer et al. 2004; Arora 2017). After convergence was reached, initial parameters were predefined and considered to run the mechanised protocols for model fitting.

Model selection criteria

The selection criteria used to determine the best explicative and predictive models included Residual Sum of Squares (RSS), Mean Squared Prediction Error (MSPE), Adjusted R Squared (Adj. R²) and its standard deviation across does, Akaike (AIC), corrected Akaike (AICc), and Bayesian information criteria (BIC). The residual sum of squares (RSS) is a statistical technique used to measure the amount of variance in a data set that is not explained by a regression model. RSS computes the explanatory ability of the model. Cross-validated Minimum Mean-Square Residual or error (MMSE) (Asherson et al. 2008), was chosen to determine error variation as an alternative to cross-validated Mean Squared Error (MSE) which has been suggested to be influenced by the number of parameters (Val-Arreola et al. 2004) if sample sizes are limited.

In comparison to R², Adjusted R squared or modified R squared (Adj. R²) is a measure of the models’ ability to predict responses for new observations while it compensates (penalises) for the overfitting event occurring after the inclusion of high numbers of predictors.

Following the premises of information theory, several methods have been described (AIC and AICc) and (BIC) of the model designed for the data being modelled. Akaike information criterion (AIC) and Corrected Akaike information criterion (AICc) were computed to compare models in regards their ability to explain or capture the variability observed in the data set being studied and Bayesian information criterion (BIC) was calculated to determine the predictive potential of each model as suggested in Karangeli et al. (2011).

In cases of a limited number of data points (observations) (N) or relatively complex models (high number of parameters), corrected AICc may be more accurate. Still, AIC and AICc become similar to a higher number of observations is studied. AICc should be used when N/K < 40 (Burnham and Anderson 2004). As Adj. R², Bayesian information criterion (BIC) is a model order selection criterion and penalises more complicated models for the inclusion of additional parameters and was computed after Leonard and Bayesian (2001).

Curve shape parameters computation for the best fitting model

Peaks were computed as described by the papers referenced in Table S1. If the computation of a peak was not possible, change in variable units per event was computed as suggested in Table S1. For the cases in which no specific manner to compute curve, shape parameters were found in literature, Symbolab® Mathematical calculation tool for education (Eqsquest 2020) was used to determine relative maxima (peak).

Results

Pre-evaluation of parametric assumptions and descriptive statistics and prior to regression considerations

Table S2 suggests all variables resulted to be non-normally distributed and residuals heteroscedastic (p < .01) except for gel-free volume (mL) and sperm concentration (× 10⁶/mL). This, together with the relatively limited sample indicated the necessity to test nonlinear regression models and Bayesian inference statistics to compare their outputs. Table S3 presents the descriptive statistics for seminal parameters. No outliers were detected for any of the seminal parameters studied.

Table 3. Summary of the measures for model fitness and flexibility criteria of the best fitting linear or nonlinear regression models for each particular trait.

Traits	SS	df	MS	F	Sig.	RSS	MSE	MMSE	AIC	AICc	BIC
Gel-free ejaculate volume (mL)	37553.61	3.00	12517.87	15.96	0.00	28232.33	784.23	3.61	264.37	264.38	10498.64
Ejaculate concentration (× 10^6/mL)	10.73	1.00	10.73	73.64	0.00	5.54	0.15	0.00	−77.08	−77.07	−3159.33
Total sperm number (TSN) (× 10^9)	2617107713.84	3.00	872369237.95	9.71	0.00	3235844535.54	89884570.43	414214.61	730.35	730.35	29137.57
Motility (%)	5886.41	3.00	1962.14	15.17	0.00	4657.97	129.39	0.60	192.30	192.30	7615.62
Normal morphology (%)	2276.62	3.00	758.87	17.11	0.00	1596.48	44.35	0.20	149.47	149.47	5902.37
Abnormal morphology (%)	1966.86	3.00	655.62	14.42	0.00	1636.92	45.47	0.21	150.47	150.47	5942.40
GSI in 1:1	4.69	3.00	1.57	29.97	0.00	1.88	0.05	0.00	−120.30	−120.30	−4888.48
TMS (× 10^9)	1790755678.46	3.00	596918559.49	13.23	0.00	1624416224.92	45122672.91	207938.58	702.78	702.78	28034.94

GSI: gonadosomatic index; TSN: total sperm number; TMS: total motile sperm.

General outputs derived from regression analyses

Table 2 shows a summary of the results for determination coefficients, curve shape parameters, and peak for best fitting models for all the variables in the study. Table S4 shows the outputs obtained after fitting the eleven models tested in the study.

Table 2. Summary of the results for determination coefficients, curve shape parameters, and peak for best fitting models for all the variables in the study.

Traits	Model	R square	F	df1	df2	Sig.	Constant	b1	b2	b3	R/Peak
Gel-free volume (mL)	Cubic	0.571	15.96	3	36	0.000	−29.439	2198.00	−0.013	0.000	0.9365
Sperm concentration (× 10^6/mL)	Exponential	0.660	73.64	1	38	0.000	485.332	−0.01	4560.000	59360.000	0.9232
TSN (× 10^9)	Cubic	0.447	9.71	3	36	0.000	−36037.895	1636.67	−12782.000	0.028	27779.1210
Motility (%)	Cubic	0.558	15.17	3	36	0.000	74.765	0.11	−0.001	0.000	78.0657
Morphologically normal (%)	Cubic	0.588	17.11	3	36	0.000	100.077	−0.35	0.004	0.000	95.6015
Morphologically abnormal (%)	Cubic	0.546	14.42	3	36	0.000	−0.006	0.35	−0.004	0.000	10.4817
GSI	Cubic	0.714	29.97	3	36	0.000	1.817	−0.03	0.000	0.000	1.3009
TMS (× 10^9)	Cubic	0.524	13.23	3	36	0.000	−26411.888	1207.80	−9.346	0.020	21023.5586

GSI: gonadosomatic index; TSN: total sperm number; TMS: total motile sperm.

Seminal parameters modelling

Determination coefficient for best fitting model was 0.571 for the cubic model (CUB) and progressively decreased in quadratic (QUA), Logarithmic (LOG) and Linear (LIN) models (Figure 1). The model suggested a progressively increasing trend was followed by gel-free ejaculate volume with age, with an inflection point occurring between 12.50 and 16.67 years. The rate of increase in gel-free ejaculate volume per event unit is 0.9365.

Figure 1. Gel-free ejaculate volume (mL) modelling using a cubic function.

Ejaculate concentration modelling

Compound (COM), Growth (GRO), Exponential (EXP) and Logistic (LGI) models presented the highest determination coefficients (best fitting model) of 0.660, followed by the cubic model (CUB) (Figure 2). A progressively decreasing trend was followed by ejaculate concentration with age. No local maxima were found; hence ejaculate concentration can be assumed to decrease constantly from the moment that the animal 2.5 years on. The rate of increase in ejaculate concentration per event unit was 0.9232. Only cubic (CUB) and quadratic (QUA) models reported a moderately acceptable determination coefficient for total number of spermatozoa (Figure 3). A progressive reduction in the number of total spermatozoa occurs from around 7.5 years on when the peak (27779.1210) takes place. The peak in TSN occurs at around 7.5 years and is 27779.1210. Cubic (CUB) and quadratic (QUA) model presented best-fitting properties followed by Growth (GRO), Exponential (EXP) and Logistic (LGI) (Figure 4). The peak in motility is 78.0657 and progressively decreases from 3.33 years on. Peak in normal sperm morphology 95.6015 occurs at around 12.5 years months while the peak in abnormal morphology 10.4817 occurs at around 20.83 years months. The progression of normal sperm morphology is somehow exponential as shown provided the best fitting values (highest determination coefficients) were reported for cubic (CUB), quadratic (QUA) and compound (CM) models but followed by Growth (GRO), Exponential (EXP) and Logistic (LGI) (Figure 5). However, for abnormal sperm morphology, cubic (CUB) and quadratic (QUA) model were followed by linear (LIN) model with the improved performance of linear model being potentially attributed to the later stages of the curve from 16.67 years on when abnormal sperm morphology steadily increases and normal morphology decreased (Figure 6). Cubic (CUB) and Quadratic models were best-fitting models for US GSI with the peak (2.2750) for this parameter occurring around 14.16 years (Figure 7). Total motile sperm were best-fitted by cubic (CUB) and quadratic (QUA) models. Peak was attained around 6.25 years and presented a value of 21023.5586 (Figure 8).

Figure 2. Ejaculate concentration(×10⁶/mL) modelling using a compound function.

Figure 3. Total Spermatozoa number (TSN) in × 10⁹ modelling using a cubic function.

Figure 4. Sperm motility (%) modelling using a cubic function.

Figure 5. Normal sperm morphology (%) modelling using a cubic function.

Figure 6. Abnormal sperm morphology (%) modelling using a cubic function.

Figure 7. Gonadosomatic Index modelling using a cubic function.

Figure 8. Total motile sperm (× 10⁹) modelling using a cubic function.

Evaluation of model fitness and flexibility criteria

Table 3 shows a summary of the measures for model fitness and flexibility criteria through mean square residual or error and minimum square residual or error (MSE and MMSE), variability explanation power through Akaike Information Criterion (AIC) and Corrected Akaike Information Criterion (AICc) and Predictive power through Bayesian Information Criterion of the best fitting linear or nonlinear regression models for each particular trait.

Discussion

GSI modelling

GSI was reported to increase with age, following a pattern similar to body weight (Coulter and Foote 1977). These findings suggested that testis tended to reach mature size more rapidly, as indicated by the curvilinear relationship between scrotal circumference and body weight. Coulter and Foote (1977) would suggest a substantial correlation of 0.81 was found between scrotal circumference and body weight in bulls.

The best fitting models found in our study were also determined to be the best fitting ones in rams in the study by Yakubu and Musa-Azara (2013). For instance, these authors reported estimation of body weight from testicle length and diameter appeared to be better under the cubic model based on higher R² values. Among the single traits in the linear, quadratic and cubic models, scrotal circumference was the most important trait in predicting body weight, which appeared to be best done under the cubic model (R² and adjusted R²= 0.890 and 0.887, respectively), which may be comparable to our results in Miranda Donkeys. This was also supported by the high correlations in the range of over 0.8 found for body weight and testicular dimensions in Martins-Bessa et al. (2020) for Miranda donkeys.

The findings by Aissanou and Ayad (2020) revealed testicular and epididymal biometry increases as donkeys age. Furthermore, the results support the trends in Miranda donkeys for gonadosomatic index presented in Figure 7. As suggested by Aissanou and Ayad (2020) for Algerian donkeys, the peak of gonadosomatic index may be of 1.6 g/kg and may be attained at 12 years. The higher peak reported by our results may be ascribed to the highest size and body weight of Miranda donkeys in comparison to Algerian ones (from 140 to 250 kg), which was also consistent with the increase in testicular measurements. In these regards, the values of gonadosomatic index for these donkeys below 6 years are half the values of the Miranda foals of around 4.16 years to 8.33 years. Gonadosomatic index was higher (1.27 ± 0.41) than that of the Algerian local breed donkeys (0.82 ± 0.07) and much higher than in Pêga donkeys (0.15) (Neves et al. 2018). The results of studies in Algerian donkeys (Ayad et al. 2019; Aissanou and Ayad 2020) suggest a faster growth in Miranda donkeys.

Seminal parameters modelling

Dorado, Acha, Gálvez, et al. (2013) reported the correlations of the age with the measures of semen quality (sperm concentration, total motility and abnormal sperm morphology) were low and not significant (p > .05) in Andalusian donkeys, but this may be ascribed to the range of males used in the study (4–15 years), which may coincide with the section of the curve during which the slope is rather constant (from 4 to 15 years). For instance, the study by Kidd et al. (2001), suggested that increased male age is associated with a decline in semen volume, sperm motility, and sperm normal morphology but not with sperm concentration in humans.

This contrasts the finding in horses by Wilson et al. (2019), who suggested that younger stallions aged 2–4 years had 24% lower total and 21% gel-free ejaculate volume than stallions aged 5–9 years, while 10–14 years stallions presented 27% lower total and gel-free ejaculate volume. However, in young stallions (2 to 4 years) sperm concentration was 30% higher when compared to stallions aged 10–14 years. These results are in line with our findings in Miranda donkeys, whose ejaculate’s volume increased from puberty, stabilising around the age of 12 years, to increase again as age progresses from 12 years on. Dowsett and Knott (1996), suggested stallions below the age of 3 years and over the age of 11 reported poorest values for gel-free ejaculate volume, sperm concentration, total sperm numbers and sperm morphological abnormalities, with the breed of the animals conditioning the values reported. This has also been reported for boars (Savić et al. 2013). As suggested by Savić et al. (2013), there is an increase in the volume of ejaculate with age in boars. The main reason for this increase in ejaculate is probably due to the testicular parenchyma growth associated with animal age.

Several mechanisms have been proposed to explain how ageing in males may cause changes in semen parameters. These changes can be related to seminal vesicle inadequacy, which reduces semen volume or prostate atrophy such as reduction in water and protein content which might affect sperm motility and ejaculate volume (Kidd et al. 2001). However, the same authors reported a significant reduction in semen volume and sperm motility of only 3 and of 4% in morphology, which may explain the lowest value for R² of 0.5708 in our study and the modelled trend which was increasing until the age of around 21.67 years. Also, major changes associated to age must be related to testicular degeneration (Nipken and Wrobel 1997), which impacts spermatogenesis and may justify deterioration of seminal quality, as supported by decreased motility (Miró et al. 2005).

The opposite trends followed by gel-free ejaculate volume and sperm concentration in our study may be indicative of the fact that as age increases volume does too, which contributes to a reduction on ejaculate concentration. As model computational complexity increases, variability capturing properties increase too as it has been suggested in literature (Pizarro Inostroza et al. 2020a, 2020b). Our results suggest that as R² increases, the decreasing trend of seminal parameters becomes more patent. Although some authors (Aleisa 2013) have reported the factor age as not to have/having a direct negative effect on sperm concentration, viscosity, or morphology, an inverse effect of age was observed on sperm motility and volume in humans. Hence, it may be the changes in volume (increasing trend as suggested by Figure 1), which, besides age-related idiopathic testicular degeneration, may contribute to the decrease in/of concentration (Figures 2 and 4).

For horses (Wilson et al. 2019), age was found to significantly influence the quality of the seminal characteristics investigated, with a marked improvement in semen quality (sperm progressive motility, sperm concentration, sperm count and total progressively motile sperm) associated to stallions aged 20 years and above. These results differ from previous research within commercial breeding stallions, for which declines in semen quality by the age of 10 years have been consistently reported. This decline has been attributed to ageing stallions being more susceptible to substandard spermatogenesis and testicular degeneration. However, although previous research had reported stallion fertility to be optimal at 5–9 years, Wilson et al. (2019) observed an increase of 23% in the total progressively motile sperm of stallions below the age of 20 years old, when compared to stallions aged 5–9 years.

In the equine species, puberty occurs at an age which ranges from 56 to 97 weeks, with an average of 83 weeks (Naden et al. 1990) and sexual maturity, is achieved 2–4 years later, since testicular size and daily sperm output increase for several years after puberty (Johnson et al. 1991). However, a patent breed dimorphism may occur. For instance, the studies by Gamboa et al. (2009) in Lusitano horses evaluated in a retrospective study suggested that despite the lack of a clear age effect on the evolution of seminal parameters, a significant positive relationship between sperm motility and sperm vitality as well as between sperm motility and morphology was reported as described by other authors (Voss et al. 1982). This may be supported by the study by Aurich et al. (2003) who reported total sperm count in Noriker stallions may be higher than values for other breeds, however, a reduction in motility and percentage of normal sperm was also patent. Hence, an increased sperm number may probably be an attempt for testicle to compensate the latter two parameters (Aurich et al. 2003).

The study by Dorado, Acha, Ortiz, et al. (2013) used Spearman correlation analysis to evaluate the relationship between sperm quality parameters in the Andalusian donkey and age. Dorado, Acha, Ortiz, et al. (2013) suggested that the correlations found between motility and age were very low in absolute value, either for total motility or progressive motility. These findings (Dorado, Acha, Ortiz, et al. 2013) contrast our results as a progressive decrease was found when total motility was fitted using a quadratic or cubic function in our study, but also may contrast results reported for other species such as rats (Lucio et al. 2013) and humans (Henkel et al. 2005) for which a deceasing trend was reported alongside with age.

Dorado, Acha, Ortiz, et al. (2013) also suggested mean gel-free volume was similar in Andalusian donkeys at all ages apart from 11 to 13 years of age, where it was significantly lower (p < .05). In males of other species such as bulls (Brito et al. 2002), ejaculate volume, total number of spermatozoa and number of viable spermatozoa has been suggested to increase with age and has been associated with increased testicular volume which is in line with our results. These authors suggested a progressive degenerative process of the testicle germinal tissue may occur after 8 years, still, certain recovery of this tissue may occur as a result of compensatory action of one testicle over the other. This has also been suggested for donkeys in the study by Omar et al. (2013) who reported that unilateral orchiectomy had compensatory effects on the weight and volume of remaining testis. However, the same authors explained that adverse effects on sperm motility and viability may still affect the fertility of the animal.

Additionally, azoospermia stages have been described in donkeys due to plugged ampullae causing total ampullae obstruction (Segabinazzi et al. 2018) or after immunocastration procedures (Rocha et al. 2018), but also may be physiologically mediated by undifferentiated embryonic cell transcription factor 1 (UTF 1) and deleted in azoospermia‐like (DAZL) expression in donkey testis (Lee et al. 2017; Choi et al. 2020) and has been suggested to be an effect derived from ageing (Liu et al. 2021) which was supported by our results as shown in Figure 3. Figure 8 also describe a decrease in total spermatozoa motility which is concomitant with the azoospermia stage detected in Figure 3.

Aurich et al. (2003), suggested total sperm count and motility were significantly affected by age. These authors suggest that the percentage of stallions with semen abnormalities increased with age too. No abnormalities were found in stallions around the age of 2–3 years old. Motility reduced in 14.1% and the percentage of morphologically normal spermatozoa in 18.8% of the animals from 4 to 9 years old. 45% stallions around an age equal or over 10 years old had an abnormally low semen motility and 56% an increase in the percentage of morphologically abnormal spermatozoa. Both seminal parameters were reported to be strongly correlated, while the same was described by progressive motility of spermatozoa and percentage of morphologically normal spermatozoa, whose correlation was more patent as age increased from 4 to 9 years old. This supports our findings as suggested by Figures 4–6, which describe the same aforementioned trends in Miranda donkeys.

Aissanou and Ayad (2020) found a significant increase of the gonadosomatic index and testicular dimensions with age and body weight. Additionally, a seasonal variation in gonadosomatic index and testicular dimensions was reported with these being significantly higher in the autumn period compared to other seasons. Contrastingly, these authors did not report significant differences for body weight as age increased in donkeys. These findings were supported by Ayad et al. (2019) who reported that the testis develop until donkeys reach aged over 11 years old. However, as it has been noticed through our review of existing literature, the differences between the average of many biometric measurements can be explained by breed difference, geographical locations, nutritional level, agroecological areas, husbandry practices or even the function for which the animals are used.

Conclusions

Testis tended to reach mature size rapidly, although a faster growth in Miranda donkeys may be ascribed to breed particularities as suggested by the values of GSI. Age-related changes in donkeys may be concomitant with an increase in ejaculate gel-free volume, which may also reduce its concentration. Simultaneously, there is a reduction in motility which may rather be ascribed to the increase in abnormal morphology spermatozoa. Major changes associated to age must be related to degeneration of the testicular functions, which in turn conditions spermatogenesis and translates into a worse seminal quality which may be characterised by reduced motility patterns. Consequently, based on the results obtained herein, donkey males are more likely to be at their best reproductive capacity between the ages of 4 and 15 years old, with important individual variations.

Ethical approval

The study follows the national guidelines and premises described in the Declaration of Helsinki. Protocols applicated were permitted by the regulations of the European Union (2010/63/EU) in their transposition to the Royal Decree-Law 53/2013 and its credited entity the Ethics Committee of Animal Experimentation from the University of Córdoba.

Acknowledgements

Authors thank the Veterinary Degree’ students for their willingness to participate in the present study. Also, authors would like to thank the support of AEPGA – Association for the Study and Protection of Donkeys (Miranda do Douro, Portugal) and the Veterinary Teaching Hospital of UTAD (Vila Real, Portugal) to carry out the present work.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

Data will be made available from the corresponding authors F.J.N.G. and A.M.B. upon reasonable request.

Additional information

Funding

This work was supported by the project UIDB/CVT/00772/2020 funded by the supported by the Portuguese Science and Technology Foundation (FCT).