Inbreeding depression is a well-known biological phenomenon resulting from matings between genetically related individuals. The increase in homozygosity also increases the expression of deleterious recessive alleles while heterozygosity advantages are reduced [1]. In domestic animals, inbreeding depression has been shown to negatively affect fitness-related traits, including fertility, survival, body development, and functional performance [2,3]. These effects are of particular concern in closed populations subject to intensive selection pressure.

The Pura Raza Española (PRE) horse is a closed population, managed under a closed Studbook since 1912 [4]. Modern PRE horses are selected for both dressage aptitude and morphological traits aligned with breed standards, which has contributed to a strong fixation of type.

Several studies have reported evidence of inbreeding depression in morphological traits across horse breeds, including PRE [5,6]. In the PRE breed, the potential impact of inbreeding on dressage traits is still underexplored, but recent studies have highlighted the usefulness of modelling inbreeding depression load (IDL) to evaluate susceptibility in both morphological and functional traits [4,6]. Traditionally, Wright’s inbreeding coefficient (F) has been the standard for estimating individual inbreeding levels. However, it fails to account for the origin of inbreeding or the specific contribution of different ancestors. To overcome this limitation, García-Cortés et al [7] proposed the use of partial inbreeding coefficients (F_ij) by decomposing F into its ancestral sources. Later, the use of F_ij allowed the estimation of the individual IDL, successfully applied in PRE to reproductive and morphological traits, offering more informative tools to detect animals with high genetic value but reduced susceptibility to inbreeding depression [4,6,8].

Furthermore, complementary pedigree-based metrics such as recent inbreeding (F₆–over six generations) and ancestral inbreeding (F_k–Kalinowski inbreeding coefficient) provide additional information on the timing and persistence of inbreeding effects [9–11]. These indicators are particularly important in structured or subdivided populations like the PRE, where historical bottlenecks, line breeding, and founder imbalance have shaped the genetic architecture [12]. Evaluating the genetic IDL for both dressage traits and conformational defects in PRE horses using F_ij may provide useful insights for genetic evaluation and sustainable breeding strategies.

Therefore, the objectives were: (1) to determine the presence and magnitude of inbreeding depression in these traits using different inbreeding estimators; (2) to estimate genetic parameters for selected dressage traits and conformational defects in PRE horses; and (3) to estimate the individual IDL to support the selection of robust and genetically valuable individuals.

Due to their relevance in the functional and morphological evaluation of PRE horses, the variables studied were:

The dressage dataset comprised a total of 43,838 performance records from 4,546 PRE horses (4,221 males and 325 females), collected between 2004–2023 during dressage competitions. In these events, the scoring system utilized the average of three judges’ evaluations for the Walk movement, which was assessed on a continuous scale from 1 (poor quality) to 9 (excellent quality). Additionally, the PPR was recorded as a composite score summarizing the overall performance of the horse-rider combination throughout the entire dressage test. For this purpose, judges individually scored each movement of the test on a scale from 0 (not executed) to 10 (excellent execution). These scores were then weighted according to the relevance of each movement, and their weighted average was calculated. This final average was subsequently multiplied by 10 to obtain the final score, expressed on a continuous scale ranging from 10 (poor quality) to 100 (excellent quality). Although the score of each judge individually was not available, but rather the average score of the 3 judges, notably a key strength of the dataset lies in the interconnection of competitions through shared judges, as many of them officiated at multiple events. This overlap created a robust network of evaluations, enabling the modelling of judge effects across competitions and improving the consistency and reliability of the evaluations.

The defects traits dataset included conformation data from 57,949 horses (19,448 males and 38,501 females). The evaluations were performed by a group of 12 specially trained veterinarians who routinely conduct the basic aptitude tests for this breed. The phenotypic assessment of conformational defects was carried out while the horse stood on a hard, level surface, in a natural stance. Horses were positioned with their forelegs and hindlegs parallel and as perpendicular as possible, with hooves properly aligned. No sedatives were administered during the evaluation process. Hock conformational defects were assessed using a linear scoring scale from 1 to 3. Class 1 corresponds to the absence of the defect; class 2 corresponds to the slight presence of the defect and class 3 is the most pronounced degree of defect.

To ensure the consistency of the analysis, only horses with phenotype and both parents registered were included. Pedigree information was sourced from the Royal Association National of Purebred Spanish Horse Breeders (ANCCE) studbook. Based on the selected individuals, a pedigree comprising 398,866 horses (194,911 males and 203,955 females) was constructed with all the known generations.

Wright’s coefficient of inbreeding (F), also referred to as classical inbreeding, was calculated as the probability that an individual carries two alleles identical by descent (IBD) at a randomly selected locus. F and coefficient of inbreeding up to the sixth generation (F₆) were estimated following the methodology of Meuwissen and Luo [9], using the Endog v4.8 software [15]. The ancestral inbreeding coefficient proposed by Kalinowski (F_k), defined as the probability that an allele is currently autozygous and has been IBD in at least one previous generation [10], was computed using GRain 2.2 [11], a coefficient recently linked to inbreeding depression in this breed [4].

Moreover, F_ij, representing the joint probability that an individual i is autozygous for an allele inherited from a specific ancestor j [16,17], were determined for all individuals in the population, based on the methodology described by Casellas [18]. This strategy enabled the partitioning of total inbreeding into contributions derived from the coancestry between each individual’s parents, taking into account both the original founders and the Mendelian sampling variance among non-founders [7,19]. As a result, 193,214 F_ij coefficients exceeding 0.00001 were identified, corresponding to 3,994 common ancestors (including both founders and non-founders from paternal and maternal lines). Each of these ancestors contributed F_ij values to between 14 and 102 descendants. To assess the effect of inbreeding depression, simple linear regression coefficients were estimated between the inbreeding coefficients (F, F₆, and F_k) and the phenotypic values of the analysed traits. Two regression models were employed using R. The first, based on the lm() function, fits continuous response variables using ordinary least squares under the assumption of normally distributed errors (used for Walk and PPR). The second, utilizing the polr () function from the MASS package, is tailored for ordered categorical response variables by estimating the probability of the response falling at or below a given category in accordance with the proportional odds assumption (applied to Closed and Convergent hocks). The t-values resulting from these regression models serve as indicators of the statistical significance of each predictor, with higher absolute t-values denoting a stronger contribution of the corresponding variable to the model. All statistical analyses were conducted using the R software [20].

For the genetic parameters estimation, each trait (y) was analysed using a Bayesian linear model initially proposed by Casellas [18] and later reformulated by Varona et al [8]. This model incorporates both the standard breeding value and the IDL attributable to the individuals’ ancestors as random genetic effects, capturing the phenotypic variation while accounting for the additive nature of the individual IDL (i). The final model applied to dressage traits was expressed as:

The model used for defect traits followed a reduced form, represented as:

These models incorporate a vector representing the total inbreeding of the evaluated individual (f), a covariate associated with the overall inbreeding depression (c), systematic effects (b), infinitesimal additive genetic contributions (u), individual IDL effects (i), and residual terms (e). In addition, gait models included the permanent environment effect (h) and the rider random effect (r). The prior distributions assumed for the additive genetic effects, the IDL effects, and the residuals were as follows:

where G=(σu2σuiσuiσi2);

A is the numerator of relationship matrix; σu2,σi2, σ_h², σ_r² and σ_e² are the associated variance components for additive, IDL, permanent and rider effects, and σ_ui is the covariance between the additive genetic and the IDL effects, respectively. R, W, X, and Z are incidence matrices of rider and permanent effects, systematic effects and additive genetic contributions, respectively, and K = T(I−P). T is a lower triangular matrix in which each non-zero element is a F_ij which links the phenotype of an inbred individual to each ancestor causing inbreeding. For computational reasons, each F_ij was multiplied by 10 to obtain the IDLvariance for a F_ij of 10%. P is a projection matrix with a 0 diagonal and 0.5 in the elements that link an individual to its sire and dam, and I is the identity matrix.

Specifically, for dressage traits, the systematic effects vector (b) included the fixed effects of sex (2 levels), stud size (3 levels: small, fewer than 3 foals born per year; medium, between 3 and 9 foals per year; and large, more than 9 foals per year), location-date of competition (932 levels), and competition level (11 levels). Rider (1,395 levels) and permanent environmental effects (3,638 levels) were included as random effects, and evaluation age was incorporated as a covariate. For conformational defect traits, the model was simplified and only included sex and stud size as fixed effects, with evaluation age also included as a covariate. The qualifying effect was not included since it was not significant (data not shown), as had already happened with this type of variables.

The Markov chain Monte Carlo (MCMC) model was implemented using custom software developed in FORTRAN90 [8]. Convergence was assessed by visual inspection of trace and running mean diagnostic plots; no trend was observed after burn-in. A single chain consisting of 1,000,000 samples, following a burn-in period of 20,000, was used for each trait. This analysis further enabled the calculation of the IDL ratios, which can be interpreted as the proportion of the variation in each phenotypic unit attributable to the variation in inbreeding depression effects. This proportion is estimated in a theoretical population where each individual possesses a partial inbreeding coefficient derived from a single, specific ancestor [6,8]. Inference was based on posterior distributions and HPD95 credible intervals; therefore, classical frequentist residual diagnostics (e.g., QQ-plots, Shapiro–Wilk) are not required within this Bayesian framework. Finally, to assess the relationships among IDL among traits, Pearson correlation coefficients were calculated using R software [20] and the cor() function with the Pearson method.

To evaluate the potential impact of incorporating IDL into selection decisions, a series of selection index framework were implemented following Hazel and Lush [21], adapted for estimated breeding values (EBV) as described by Gutiérrez et al [22]. Three index scenarios were simulated: using the EBV selection criteria as the selection objective (index 1), using the IDL as selection criteria and the EBV as the selection objective (index 2), and combining the EBV and the IDL as selection and criteria objectives (index 3). For each index, expected genetic response (EGR) were computed, assuming equal additive genetic variance (unit scale), selection intensity of 1, and the same heritability-based accuracy across individuals [23]. Different economic weight vector (p′) were tested for index 3.

The descriptive statistics for traits analysed in the PRE population can be seen in Table 2. For the continuous gait traits, Walk has a mean of 6.67 with a very low standard deviation (0.005) and a coefficient of variation (CV) of 0.07%, indicating low dispersion relative to its mean. In contrast, PPR exhibits a mean of 65.84 with a standard deviation of 5.22 and a CV of 7.91%, showing a slightly lower relative variability compared to Walk. The range for Walk spans from a minimum of 1 to a maximum of 9, and for PPR from 11.85 to 98.10. Regarding the defect traits, which are assessed on an ordinal scale where lower scores indicate better conformation, both Closed and Convergent hocks are confined to a range of 1 to 3. For Closed hock, the mean is 1.32 with an SD of 0.60. The distribution is skewed toward the absence of the defect, with 74.76% of animals scored as 1 and 25.24% showing the presence of the defect (class 2+3). In contrast, Convergent hock displays a higher mean score of 1.93 with a standard deviation of 0.68. Its distribution is less favourable, with only 26.91% of animals scored as 1 and 73.10% showing the presence of the defect (class 2+3).

The regression coefficients (b) for dressage and defect traits with the different inbreeding coefficients in the PRE horse are displayed in Table 3. For the gait trait Walk, negative regression coefficients were observed across all inbreeding measures. Specifically, the coefficients for F, F₆, and F_k were −0.68, −0.74, and −0.75, respectively. These coefficients were highly significant (all p<0.001). Similarly, the PPR trait exhibited significant negative associations with inbreeding, with coefficients of −3.5 for F, −4.08 for F₆, and −3.10 for F_k. For the defect trait Closed, the ordinal logistic regression revealed positive coefficients for all inbreeding measures, evidencing inbreeding depression: 1.60 for F, 1.95 for F₆, and 2.18 for F_k. These values were statistically significant. Convergent hock trait presented a more variable pattern. The regression coefficient for F was not statistically significant, whereas F₆ yielded a positive coefficient of 0.43 and F_k produced a significant negative coefficient of −0.67.

Table 4 shows heritability estimates, IDL ratios (d²) and the variances attributed to the direct additive genetic effect, IDL (corresponding to an inbreeding value of 10%), the random effects (for dressage traits only) and the residual effect of the analysed traits. Analogously to heritability, d² were calculated as the relative magnitude of the IDL variances. The heritability estimates were 0.02 both for Walk and PPR, while higher values were obtained for Closed hock (0.15) and Convergent hock (0.22). d² calculated for an inbreeding value of 10%, were higher than heritability values for all traits (0.81, 0.89 and 0.21 for Walk, PPR and Closed hock) except for Convergent hock (0.18). The IDL variances were higher than the direct additive genetic variances as a direct consequence of the arbitrary scaling of F_ij by 10. Finally, as expected, repeatability (r) slightly exceeded heritability (h²) for the dressage traits with repeated records (Walk r: 0.05 vs h²: 0.02; PPR r: 0.04 vs h² = 0.02; Table 4).

Figure 1 displays histograms of the estimated IDL for each of the four traits within the evaluated PRE horse population. Across all traits, the distributions are largely unimodal and centered around zero, indicating that most individuals exhibit moderate IL. However, a subset of horses displays more extreme positive or negative values. The blue bars highlight IDL values surpassing (or falling below) the predicted regression slope threshold for an assumed 10% inbreeding level. Notably, these regions suggest that only a small fraction of the population experiences substantially higher or lower IDL which implies that the use of these animals as common ancestors of the future animals will contribute to an improvement of those traits even when they transmit a high inbreeding level. Specifically, the predicted regression slope threshold between IDL and the F for an assumed 10% inbreeding level were −0.01 (Walk, Closed hock and Convergent hock) and −0.02 (PPR), implying that the proportion of animals that can be used as a common ancestor with positive results are low, 2.94% (Walk), 0.77% (PPR), 1.30 (Closed hock) and 0.69% (Convergent hock).

Table 5 shows the Pearson correlation coefficients between IDL values and the number (and percentage) of individuals exhibiting simultaneous improvement in IDL values for each pairwise combination of traits. The correlation coefficients indicate a moderate positive relationship between Walk and PPR (0.45). The correlations between Walk and Convergent hock (0.07), and between Walk and Closed hock (0.07), as well as those between PPR and the defect traits, Closed hock (−0.02) and Convergent hock (0.06), indicate only a weak relationship.

Regarding the number and percentage of animals showing coincident favourable IDL across traits, the largest overlap was again between Walk and PPR, with 2030 animals (51%). In contrast, only 16% to 32% of the population showed overlapping improvement between one gait trait and one defect trait, and the lowest overlap was found between PPR and Convergent hock (16%).

Beyond these cross-trait overlaps, Figure 2 summarizes the temporal trajectories of F₆, EBV, and individual IDL by cohort. We focused on F₆ as the time-varying metric because it captures recent inbreeding and showed the strongest and most consistent associations with the phenotypes in our regressions (Table 3). F₆ increased from the early decades of the twentieth century, peaked around 1958, and then declined while remaining above early-cohort levels. Over the same period, mean EBV for dressage traits remained broadly stable, with slight improvements in recent decades, whereas IDL for PPR decreased (became more negative) and IDL for Walk remained essentially stable. For the hock defects, IDL increased over time for Convergent hock, while EBV for Closed hock stabilised with slight recent improvements and EBV for Convergent hock showed an overall decline.

The expected EBV responses for Walk and Closed hock, together with the relative change versus Index 1 (baseline = 100%), are shown in Table 6. Index 1 (EBV-only) delivered the largest response for both traits (0.14 for Walk; 0.39 for Closed hock). Index 2 (IDL-only) produced a marked reduction in EBV gain, especially for Walk (−589.46%). Hybrid indices that combine EBV and IDL (Index 3 variants) yielded progressively lower EBV responses as the IDL weight increased.

Descriptive results (Table 2) show a distribution typical of intensively selected populations, such as the PRE breed. The average Walk score (6.67) is like that of other dressage breeds, but its low variability (0.07%) suggests possible genetic fixation due to selection. This might also reflect a tendency among judges to favour intermediate scores, limiting the use of extreme values and contributing to score homogenization. In contrast, Warmblood horses show similar mean Walk scores but with variation coefficients above 15%, indicating more genetic variability [24,25]. The PPR had a mean of 65.84 and a moderate SD (5.22). Sánchez-Guerrero et al [26] reported a nearly identical average (65.33) and a CV of 6.94% in PRE horses, supporting trait consistency. While phenotypic variability in PPR appears stable across PRE studies, comparisons with other breeds are scarce. However, Solé et al [27] revealed that Warmblood horses, particularly those competing at the international level, tend to outperform both Lusitano and PRE horses in dressage disciplines.

For conformational defects, most animals were unaffected by for Closed hock (74.76%), but only 26.91% were unaffected for Convergent hock. These results match those by Ripollés-Lobo et al [14], who found 77.91% and 25.88%, respectively. This confirms the high frequency and functional relevance of Convergent hock defects in PREs. Compared to other breeds, the PRE shows intermediate prevalence. In Menorca Purebreds, Ripollés-Lobo et al [28] found fewer healthy animals for Closed hock (37.47%) and more for Convergent hock (61.17%).

Inbreeding depression has been studied in different horse breeds [5,29–31]. The magnitude of inbreeding depression is influenced not only by the overall inbreeding level but also by the timing and origin of inbreeding, which can be better captured through refined coefficients such as F₆ and Fk, rather than classical Wright inbreeding coefficient.

Comparatively, the magnitude of IDL and the negative regressions with recent inbreeding observed here align with inbreeding-depression evidence on performance in other sport breeds (notably Thoroughbreds) and with meta-analytic results across livestock [32,33]. However, the particularly strong response in PPR in PRE is consistent with a closed studbook and the historical concentration of elite sires [4]. In Lusitano horses, heritabilities for dressage and gait traits are typically higher [34], suggesting a lower expressed recessive load for these traits. Within PRE, morphological traits and fertility show detectable inbreeding effects but different sensitivities and responses to selection [6,31].

For the functional traits (Walk and PPR) all three inbreeding coefficients (F, F₆, F_k) showed consistent negative associations, suggesting that higher levels of inbreeding are linked to poorer performance in dressage-related abilities (Table 3). This agrees with prior studies on gait and locomotion traits across breeds like Andalusian horses [5]; Lusitan horses [35]; Thoroughbred horses [29,32]. The most pronounced effects occurred with F₆, supporting the idea that recent inbreeding increases the likelihood of phenotypic expression of deleterious recessives, as these have not yet been purged from the genome. This interpretation aligns with findings in Thoroughbred horses by Hill et al [32], who identified a significant association between a specific haplotype on chromosome ECA14 (THR14) and reduced probability of racing. In comparison, the weaker or even inverse regressions observed with F_k, in PPR and Convergent hock, suggest that ancestral inbreeding may be associated with partial purging of harmful alleles, especially in lines subjected to consistent selection pressure. This is supported by findings in other horse breeds as Lusitano horses [34], where ancestral inbreeding effects were less pronounced and often curvilinear.

Todd et al [29] reinforced this interpretation using the ancestral history coefficient (AHC) and F_k. While global inbreeding reduced performance, ancestral inbreeding improved it, and founder-specific effects highlighted the influence of inbreeding origin. Casellas [18], proposed a similar framework, noting that inbreeding depression depends on which genomic regions and ancestors are involved, not just on overall F.

Morphological defects examined, exhibited more heterogeneous responses. Closed hock showed a consistent positive regression with all inbreeding coefficients, suggesting the persistence of a stable additive load that has not been effectively purged over time. These results agree with Ripollés-Lobo et al [28] in Menorca Purebred. On the other hand, Convergent hock presented a distinct and biologically meaningful pattern, while the regression with F₆ was positive, indicating ongoing inbreeding depression, the association with F_k was negative, supporting a purging effect in animals with a history of ancestral inbreeding. This duality reinforces the idea that visible, easily penalized traits may be subjected to more effective purging dynamics, as suggested by Casellas [18]. This pattern has also been confirmed in the PRE population by Ripollés-Lobo et al [14], who found that animals with F>12.5% were significantly more likely to express limb defects, particularly Convergent hock, which also had a high prevalence (74.1%), yet many severely affected animals may be excluded from registration, masking the full impact of inbreeding.

The idea that not all conformation traits are equally sensitive to inbreeding is further supported by Vostrý [36] who analysed 22 linear conformation traits in Czech cold-blooded horses and found variable but detectable inbreeding effects, even at low mean F. Traits related to movement and overall body size were more affected by inbreeding than those associated with specific anatomical regions, such as localized limb structures. Moreover, they demonstrated that including inbreeding in the model altered heritability estimates and breeding value reliability, underscoring the relevance of modelling inbreeding explicitly in genetic evaluations.

The large-scale meta-analysis by Doekes et al [33] across seven livestock species showed that inbreeding depression is pervasive across all trait types, not just those related to fitness. Their study emphasized that refined metrics such as F_k or runs of homozygosity (ROH) improve the estimation of inbreeding effects and that even moderate inbreeding levels can significantly impact trait means and variances, especially in closed populations.

Heritability estimates for Walk and PPR were low (Table 4), contrasting with higher values reported by Sánchez-Guerrero et al [26] in the same breed (0.21 for Walk, 0.30 for PPR). These differences could be due to variations in data structure, model specification, or trait definition. In other breeds, Walk heritability ranges from 0.08 to 0.38, and PPR from 0.18 to 0.32 (Dutch Warmblood [37]; Hungarian sport horse [38]; Finnhorse and Standardbred foal [39]; Lusitano horse [35]). Heritability of hock defects varies by breed and method used. In PRE horses, Ripollés-Lobo et al [14] reported moderate to high values (Closed hock: 0.26; Convergent hock: 0.42) (Table 4), while in Menorca horses (0.12 and 0.24, respectively) [28]. The very low heritability estimates for Walk and PPR (0.02), together with the high d² ratios (0.81–0.89), are primarily a consequence of the modelling strategy adopted here. Unlike earlier PRE analyses, our Bayesian framework explicitly partitions the genetic variance into a direct additive component and an individual-specific IDL; accordingly, variance that in additive-only models would inflate heritability is reallocated to the IDL term, yielding lower apparent heritability values. This conceptual partitioning, rather than an abrupt biological change, likely explains most of the discrepancy with previous PRE reports [26,40]. Additional contributors may include differences in data structure and model specification (wider time span and stronger control of rider and permanent environmental effects), as well as judge-based subjectivity and environmental heterogeneity related to training level, all of which can further reduce the detectable additive genetic signal relative to studies in PRE or Warmblood populations [25,34,37].

The results highlight the relevance of IDL in the genetic architecture of performance and defect traits in PRE horses (Table 4). In Walk, PPR, and Closed hock; IDL variance exceeded additive genetic variance, with only Convergent hock showing higher additive values. This pattern aligns with findings in cattle [8] and equines [4,6], supporting that the transmitted recessive genetic load can explain a larger share of the phenotypic variance than additive effects alone, especially when partial inbreeding (F_ij) is considered.

The d², representing the proportion of phenotypic variance due to inbreeding load (assuming 10% F_ij) was particularly high for dressage traits: 0.81 for Walk and 0.89 for PPR, far exceeding their heritabilities (both 0.02). This indicates that most genetic variation in these traits originates from inbreeding effects rather than additive contributions. These high ratios are not only biologically meaningful but also potentially actionable, as they provide insight into how specific ancestral lineages might affect phenotypes through inbreeding. These figures reflect the underlying genetic load that possibly will exist in the population due to founder effects and pedigree bottlenecks, despite moderate additive control. Such results echo Casellas [18] who highlighted the importance of individual-specific IDL in identifying deleterious recessive variance often missed by the additive models.

Figure 1 illustrates the complexity of inbreeding depression through individual-level IDL distributions for each trait. All traits show roughly unimodal, symmetrical distributions cantered around zero, as also reported by Casellas [18], Poyato-Bonilla et al [6] and Perdomo-González et al [4]. This suggests most animals will have an intermediate IDL, although a subset displays extreme values, either positive or negative. Individuals with IDL values lower than the regression between IDL and the inbreeding are likely to transmit to their inbred offspring a potential overcome the negative effects of inbreeding on defects traits. Conversely, individuals with IDL values higher than the regression between IDL and the inbreeding may transmit to their inbred offspring a potential to offset the negative effects of inbreeding on dressage-related traits. This observation is crucial because it challenges the classical assumption that inbreeding effects are uniformly detrimental and instead supports the concept of individual heterogeneity in IDL, as discussed by Varona et al [8] and Perdomo-González et al [4]. Moreover, the Figure 1 highlights thresholds of IDL corresponding to an F value of 10%, providing a benchmark to identify tolerant or even favourable transmitters of IDL. For instance, only 0.77% of individuals had an IDL above −0.02 for PPR, suggesting that an individual with a F_ij of 10% or more, derived from one of those specific common ancestors, tends to perform better in dressage than its contemporaries. This implies that only a small elite group of animals can improve PPR performance in inbred offspring, despite the high susceptibility of this trait to inbreeding depression. Only 0.69%–2.94% of animals exceed the favourable IDL threshold under an assumed 10% F_ij, indicating that most lines still transmit a significant recessive load. This scarcity of resilient lineages (e.g., those predicted to perform well despite a 10% of inbreeding) is consistent with expectations from closed populations that have experienced historical bottlenecks and long-term overuse of a limited number of elite sires. Practically, this low resilience carries economic consequences, including reduced sporting performance and international competitiveness, and underscores the need to complement EBV based selection with IDL to preserve valuable lines while limiting the recessive load.

From a breeding perspective, these results provide a valuable tool: instead of discarding animals based only on inbreeding coefficients, breeders can select individuals with favourable IDL values (even at moderate inbreeding) preserving valuable lines while reducing harmful recessive effects [4,18].

Together, the results from Table 4 and Figure 1 provided strong evidence to include individual-specific IDL in PRE genetic evaluations. That IDL variance often exceeds additive variance, and that d² surpasses heritability in key traits (PPR and Walk), clearly indicated the concealed but substantial influence of recessive load. The distribution and thresholds of IDL further offer breeders a practical way to identify inbred individuals that avoid expected phenotypic decline and may even enhance performance or reduce defects. This suggest that managing inbreeding based solely on pedigree coefficients (F, F_k) may be insufficient. Selection based on IDL allows for more refined decisions, preserving valuable lines while limiting deleterious alleles, a scientifically grounded improvement to current PRE breeding strategies, particularly in traits with low heritability but high susceptibility to inbreeding depression.

Cross-trait correlations among IDL values were generally low to moderate and mostly positive (Table 5), consistent with some shared recessive burden across traits [33]. However, because the favourable IDL direction differs by trait (higher for dressage, lower for defects), a positive correlation does not necessarily imply a favourable correlated response. In practice, the proportion of animals showing simultaneous favourable IDL in both traits was limited, ranging from 16% to 32% across dressage–defect pairs and as low as 16% for PPR with Convergent hock (Table 5). These results support using IDL alongside EBV when ranking and mating, rather than relying on EBV alone.

As shown in Figure 2, despite the mid-century peak and later decline in F6, the recessive burden has not abated. In dressage, IDL for PPR became more negative, which is unfavourable because higher IDL denotes greater resilience to inbreeding, while IDL for Walk was essentially stable and EBV remained broadly stable with slight recent gains. For hock defects, IDL increased for Convergent, which is unfavourable since lower (more negative) IDL is desirable for defects, whereas EBV for Convergent hock declined, which is favourable, and EBV for Closed hock was largely stable with slight recent improvements; IDL for Closed hock showed no marked trend.

Two representative traits were selected based on the sign of their correlation between EBV and IDL for the selection index simulations: one with a positive correlation (e.g., Closed hock) and one with a negative correlation (e.g., Walk) (Table 6). For Closed hock, which shows a positive correlation between EBV and IDL, none of the designed indices appear to produce any improvement and index 1 will be recommended. In contrast, although a superficially similar pattern is observed for Walk, the negative correlation between EBV and IDL in this case results in a different outcome. The index 3 assigning 90% economic weight to EBV and 10% to IDL may be beneficial in the long term for Walk. Although it results in slightly lower genetic progress, it leads to an increase in IDL over time, which could be advantageous in inbred populations such as the PRE, by lowering the expected risk of inbreeding depression.

These patterns indicate trait-specific selection responses. Strong penalties and culling reduce the expression of Convergent hock (downward EBV), but the hidden recessive load persists (rising IDL), consistent with incomplete purging and with the negative F_k regression for Convergent hock. Conversely, the diffuse recessive burden in the polygenic and subjectively scored PPR proved difficult to remove. Practically, EBV only selection is insufficient; programmes should incorporate IDL aware decisions, for example EBV+IDL indices and IDL informed mate allocation, while diversifying lines and limiting the overuse of elite sires to curb the accumulation of recent inbreeding.

The findings of this study underscore the value of incorporating Individual IDL into the genetic evaluations of PRE horses. By going beyond traditional additive genetic values, IDL captures the specific impact of deleterious recessive alleles transmitted through shared ancestry, offering a more comprehensive explanation of phenotypic variability in both functional traits (Walk and PPR) and morphological defects (Closed and Convergent hocks). In many cases, the genetic load associated with inbreeding accounts for a larger portion of this variability than additive effects alone, which is particularly relevant in inbred or semi-closed populations such as the PRE. The identification of individuals with favourable IDL values opens the door to more refined selection strategies that maintain genetic diversity while enhancing performance and reducing defects.