SNPs with High Linkage Disequilibrium Increase the Explained Genetic Variance and the Reliability of Genomic Predictions

TL;DR

Using SNPs with high linkage disequilibrium improves genetic variance explanation and prediction reliability in Holstein cattle traits.

Contribution

This study demonstrates that SNPs with high linkage disequilibrium enhance genomic prediction reliability and genetic variance explanation in cattle.

Findings

Using SNPs with high linkage disequilibrium increased genomic prediction reliability by 0.05 to 0.11 on average.

SNPs with high linkage disequilibrium doubled the explained genetic variance across all traits.

Haplotype-based SNPs showed larger individual marker effect estimates compared to general SNPs.

Abstract

This study aimed to evaluate and compare the proportion of explained genetic variance and the reliability of genomic breeding value predictions for six traits in Holstein cattle: milk yield, fat yield, protein yield, fat percentage, protein percentage, and somatic cell score. Three types of genomic information were tested. The first approach included 88,911 single nucleotide polymorphisms from 8290 animals. The second approach used haplotypes defined by strong linkage disequilibrium (r2 ≥ 0.80), encoded as pseudo-SNPs, with a total of 35,552 pseudo-SNPs from 8331 animals. The third approach analyzed only the SNPs forming haplotypes, resulting in 33,010 SNPs from 8192 animals. All analyses were performed using the single-step genome-wide association study method implemented in BLUPF90. The findings revealed that using single nucleotide polymorphisms with high linkage disequilibrium improved the reliability of genomic breeding value predictions compared with the use of single nucleotide polymorphisms in general, with average increases ranging from 0.05 to 0.11. Furthermore, analysis using single nucleotide polymorphisms with high linkage disequilibrium doubled the explained genetic variance across all traits, likely due to larger estimates of individual marker effects. Overall, the study highlights the advantages of haplotype-based information for improving prediction reliability and explaining genetic variance in Holstein cattle.

The objective of this study was to compare the proportion of explained genetic variance (EXGV) and the reliability of genomic breeding values (GBVs) predictions for milk yield (MY), fat yield (FY), protein yield (PY) fat percentage (FP), protein percentage (PP), and somatic cell score (SCS) in Holstein cattle. Three types of genomic information were evaluated. (a) SNP-ALL: this analysis included 88,911 single nucleotide polymorphisms (SNP) from 8290 animals. (b) HAP-PSEUDOSNP: haplotypes, defined based on high linkage disequilibrium (LD, r2 ≥ 0.80) between SNPs, which were encoded as pseudo-SNPs, with a total of 35,552 pseudo-SNPs and 8331 animals included. (c) SNP-HAP: analysis using only individual SNPs included in the haplotypes (without recoding); for this analysis, 33,010 SNPs and 8192 individuals were retained. All analyses were conducted using the single-step genome-wide association study method implemented in the BLUPF90 software package. The results showed that the inclusion of SNPs with high LD (SNP-HAP) increases the reliability of GBVs’ predictions compared to the SNP-ALL analysis; average reliability increased between 0.05 and 0.11. Moreover, the SNP-HAP analysis resulted in a twofold increase in the EXGV for all traits, likely due to increased estimates of individual marker effects compared to the SNP-ALL analysis.

Figures

Full text

1. Introduction

The use of SNPs has enabled the inference of haplotype structures, thereby enhancing our understanding of gene flow and population structure in both cattle and humans [1,2]. A haplotype refers to a set of contiguous alleles within a genomic region that shares a common genealogical origin. Additionally, haplotypes can be defined as groups of SNPs that are co-inherited due to high LD or based on predefined genomic windows of arbitrary length either in mega bases or in number of SNPs [2].

The use of haplotypes in genome-wide association studies (GWAS), aimed at identifying genomic regions associated with productive, reproductive, and health-related traits typically governed by polygenic effects, has demonstrated advantages in enhancing the detection power of these methods [3,4,5]. The relevance of incorporating haplotypes into GWAS lies primarily in the nature of genetic selection in livestock, which alters the frequencies of genomic blocks rather than individual loci [2,6]. Haplotypes-based approaches sometimes can increase the likelihood of detecting quantitative trait loci (QTL) of varying sizes, especially across genomic regions with heterogeneous allele frequencies. This is largely due to the high levels of LD among SNPs, which results in low or negligible recombination rates within these blocks [4,7]. When constructing haplotypes from SNP array data for GWAS, it is crucial to consider both the density of SNP arrays and the number of SNPs per haplotype. Low-density SNP arrays (<50 K) often exhibit reduced power in GWAS due to greater average physical distances between adjacent markers (>70 kb), leading to a loss of LD information [7]. Furthermore, the number of SNPs included in each haplotype must be carefully selected, as LD between markers diminishes with increasing physical distance. For instance, Salem et al. [8] reported a decline in LD measured by D’, from 0.815 at ~30 kb to 0.578 at ~500 kb, and a decrease in r^2^ from 0.283 to 0.090 over the same range.

The use of haplotypes defined as sets of SNPs in LD rather than individual SNPs, for the prediction of GBVs using statistical methods such as best linear unbiased prediction (BLUP), has shown potential to enhance the accuracy and reliability of genomic predictions in certain scenarios [9]. This improvement is also influenced by multiple factors that determine the predictive power of the model. Some of these factors are the parameters used in haplotype construction, including the number of SNPs per haplotype and the threshold LD level (commonly measured by r^2^) between adjacent SNPs [9,10]. Additionally, trait-specific factors such as heritability, the number and effect size of underlying genes or QTLs, and the polygenic nature of the trait also influence prediction accuracy [11]. Population-related factors further affect genomic prediction. These include genetic diversity, population stratification, and historical selection pressures [12].

In various livestock populations, the use of haplotypes instead of individual SNPs has led to small changes in the accuracy of GBV predictions. For instance, in Korean Hanu cattle, the evaluations of meat quality traits reported increase in prediction accuracy, ranging from 1.3 to 4.6 percentage points, when using haplotypes [9]. Similarly, studies in Simmental cattle observed improvements in GBV accuracy of up to 11.29% for carcass traits [10]. In contrast, studies in dairy cattle have reported smaller effects. For example, genomic evaluations for MY, FY, and PY showed an increase from 1 to 2% in reliability when using haplotypes of fixed length comprising five SNPs. Moreover, when haplotype length was increased to 20 SNPs, reliability slightly decreased by approximately 2% [13].

There is still no consensus on the optimal strategy for haplotype construction; however, some studies suggest that defining haplotypes based on LD levels may be more effective than using a fixed number of SNPs. The fixed-length approach can result in spurious haplotypes formed merely from random allelic combinations [14], without considering the degree of LD between SNPs. This method overlooks the key principle of haplotype that they represent blocks of alleles with strong LD, inherited together across some generations [2,11,15]. Constructing haplotypes based on LD thresholds may better capture the underlying genomic structure and more accurately reflect biologically meaningful inheritance patterns.

On the other hand, the effects of SNPs depend on the LD between the markers and the QTLs affecting the phenotype. Levels of LD between markers vary among populations, as they are influenced by factors such as breeding objectives, selection intensity, animal adaptability, migration, mating types, and breed history [2,4].

Consequently, evaluating the use of SNPs with high LD in genomic analyses of specialized cattle breeds involves prioritizing allelic variants that occur at high frequencies within the population. These SNPs are more likely to be associated with relevant genes because of long-term artificial selection. As a result, their inclusion may increase estimated marker effects and enhance the contribution of SNPs to the EXGV of the traits of interest [2,11].

Previous studies have reported modest improvements in the reliability or accuracy of genomic predictions when using haplotypes defined by fixed SNP windows or by LD, both in dairy cattle and other species [12,13,15]. However, these approaches generally rely on recoding haplotypes as pseudo-SNPs, which increases computational complexity and implementation time [15].

In contrast, the present study advances existing research by evaluating an alternative strategy that exploits LD information without recoding haplotypes. Specifically, we assessed that the use of individual SNPs characterized by high LD (r^2^ ≥ 0.80) could increase the proportion of EXGV and the reliability of GBVs in Holstein cattle, while remaining computationally simpler and more practical for routine genomic evaluations.

The objective of this study was to evaluate three different genomic sources: (a) SNP-ALL analysis using all available individual SNPs, (b) HAP-PSEUDOSNP analysis using haplotypes encoded as pseudo-SNPs based on LD structure, and (c) SNP-HAP analysis using only individual SNPs with high LD identified during haplotypes construction, and measure their impact on the proportion of EXGV by markers, the magnitude of marker effects on GWAS, and the reliability of GBVs for MY, FY, PY, FP, PP, and SCS in a Holstein cattle population.

2. Materials and Methods

2.1. Phenotypic Information

The analysis included 640,746 records of MY corresponding to the 1st (45.26%), 2nd (33.72%), and 3rd (21.02%) lactations, previously adjusted to 305 days and mature equivalent yields from 358,857 Holstein cows born from 1979 to 2019. These records included data on 204,873 records for FY, FP, PY, and PP and 171,890 records for SCS. The SCS was expressed on a linear scale ranging from 0 to 9 (where 0 = 12.5 × 1000/mL, increasing twofold for each point, until 9 = 6400 × 1000/mL) [16]. The cows were distributed across 353 herds located in 18 states of the Mexican Republic. The pedigree file comprised 470,695 animals, including 17,220 sires and 161,757 dams, with an average pedigree depth of five generations. The data were provided by the Holstein Association of Mexico, A.C. (AHM).

2.2. Genomic Information

The SNPs and pseudo-SNPs used in the three analyses corresponded to 8331 genotyped animals from 7809 cows and 522 sires, with information from ~116 K SNP markers previously imputed with FindHap V2 [17]. The imputation process incorporated genotyping arrays of different densities, including BovineLD v2.0 9K (4.15%), GGP Bovine LD v3.0 26K (5.91%), GGP Super LD v4.0 26K (1.38%), BovineSNP v3 50K (9.16%), GGP LD 77K (15.42%), GGP Bovine 100K (0.03%), GGP HD 150K (59.77%), Genome-Wide BOS 1 Bovine Array 640K (2.47%), and GGP HD 777K (1.71%), ILLUMINA, San Diego, CA, USA. NEOGEN, Lansing, MI, USA and Axiom, Waltham, MA, USA [18,19,20]. The genomic data was provided by the National Center for Disciplinary Research in Animal Physiology and Improvement of the National Institute of Forestry, Agricultural and Livestock Research (CENIDFyMA-INIFAP).

In the first analysis (SNP-ALL), 8290 genotypes with 88,911 individual SNPs were included. Prior to analysis, genomic quality control was applied to exclude animals with a call rate < 0.95, or with parent–progeny conflicts; the SNPs’ filtering included the parameters of minor allele frequency (MAF) < 0.05, call rate < 0.95, or a value of p for the Hardy–Weinberg equilibrium test < 0.15, or monomorphic [21].

For the second analysis (HAP-PSEUDOSNP), the construction of haplotypes as pseudo-SNPs was conducted in three sequential steps. In the first one, the SNPs from the imputation process were recoded into A/B format, where 0 = BB, 1 = AB, 2 = AA, and 3, 4 and 5 were coded as 00, because this number is for missing alleles (3 = B_, 4 = A_, 5 = __). The resulting dataset was then converted to binary BED format using PLINK 1.09 software [22,23] with the --make-bed option. In the second step, haplotypes were defined using PLINK v1.07 [1,22], applying the --blocks option for block detection as described by Taliun et al. [24]. Haplo-blocks consisted of at least two SNPs and were required to meet a minimum LD threshold of r^2^ ≥ 0.80 [11]; genomic quality control parameters included the following: --maf 0.01 --mind 0.05 --geno 0.02 and --allow-no-sex. In the third step, each haplotype was recoded into pseudo-SNPs in 2, 1, and 0, corresponding to the presence of two copies, one copy, or the absence of paternal and maternal alleles, respectively, allowing only haplotypes present in more than 1% of the population [15]. This step was also performed with PLINK v1.07, resulting in a dataset of 9095 genotypes and 43,026 pseudo-SNPs. Due to the recoding of haplotypes as pseudo-SNPs, multiple haplotypes with different allele combinations can be identified within the same genomic region at the population. However, because haplotypes represent phased allele combinations across multiple SNPs, different haplotypes within the same region cannot occur simultaneously within the same chromosome segment of an individual.

Genomic quality control was applied to the pseudo-SNPs dataset, to exclude animals with a call rate < 0.95 or with parent–progeny conflicts < 1%. Additionally, pseudo-SNPs with a minor allele frequency (MAF) < 0.05, a call rate < 0.95, or a Hardy–Weinberg equilibrium p-value < 0.15 were removed. After filtering, the final dataset included 8331 genotypes and 35,552 pseudo-SNPs, representing the constructed haplotype. This quality control was performed using the PreGSf90 1.22 program [21].

For the third analysis (SNP-HAP), a total of the 33,010 individual SNPs contained in the haplotypes (with an LD value of r^2^ ≥ 0.80) from 8182 animals were included. Both SNP and animal genotypes were subjected to the same quality control procedures as previously described to explain what the levels of EXGV would be for the evaluated traits if only high-LD SNPs were used. Notably, these SNPs were the same as those used in the construction of haplotypes in the second analysis.

After the construction of the genomic datasets for each analysis, the estimation of marker effects (SNPs or haplotypes) and their associations with the traits of interest was performed using the single-step genome-wide association study (ssGWAS) methodology. This approach is based on the prediction of genomic breeding values (GBVs) using the single-step genomic best linear unbiased prediction (ssGBLUP) method [25,26].

2.3. Prediction of Genomic Breeding Values with Haplotypes or SNPs

The GBV of the six traits in the three analyses were predicted with the single-step genomic BLUP method (ssGBLUP [27,28]), using BLUPF90 2.57 software [26] for each trait and dataset. The model to obtain the GBV was

$$[eqn]$$

where $[eqn]$ is each one of the six traits (MY, FY, PY, FP, PP, and SCS). The model’s fixed effects are as follows: $[eqn]$ is the general mean for each trait, $[eqn]$ is the i-th level of herd year season (6831 levels; 40 years and 2 seasons: January to June and July to December), and $[eqn]$ is the j-th age level at calving in months in each lactation classified from 1 to 9 according to its distribution: for the first lactation: 1 ≤ 23.4, 2 ≥ 23.5 and ≤25.5 and 3 ≥ 25.6; for the second one: 4 ≤ 35.9, 5 ≥ 36 and ≤39.5 and 6 ≥ 39.6; and for the third one: 7 ≤ 48.7, 8 ≥ 48.8 and ≤53.3 and 9 ≥ 53.4. The model’s random effects included $[eqn]$ is the k-th level of the permanent environment; $[eqn]$ is the l-th level of the sire–herd interaction (40,742 levels), $[eqn]$ is the animal effect, and $[eqn]$ is the residual effect.

2.4. Association Analyses for Haplotypes and SNPs

After the prediction of GBVs, the GWAS were carry out with postGSf90 1.70 program to estimate the haplotypes or SNP effects for each trait and genomic data source [21,27] and to compare whether haplotypes capture the same regions of significant SNPs.

The effects of haplotypes or SNPs were estimated as;

$$[eqn]$$

where $[eqn]$ is the estimated effect to a haplotype or SNP, $[eqn]$ , a scalar where $[eqn]$ is the A allele frequency for the i-th marker (SNP or haplotype), $[eqn]$ is the B allele frequency for the i-th marker, $[eqn]$ is a weight diagonal matrix which represents the variance or weight of each haplotype or SNP as proposed by VanRaden [29], $[eqn]$ is the transposed matrix of haplotype or SNPs adjusted by allele frequencies, $[eqn]$ is the inverse genomic relationship matrix, and $[eqn]$ is the GBV for each trait [28].

The p-values for the haplotype or SNP effects were calculated according to Aguilar et al. [30]:

$$[eqn]$$

where $[eqn]$ is the p-value for haplotype or SNP effects of the evaluated traits, $[eqn]$ is the cumulative standard normal function, $[eqn]$ is the estimated effect of the i-th marker (SNP or haplotype), and $[eqn]$ is the standard deviation of the effect for the i-th marker. The significant association threshold was defined with the Bonferroni adjustment [31] as follows: −log10(0.01/88,911 SNP) for SNP-ALL analysis, −log10(0.01/35,552 haplotypes) for HAP-PSEUDOSNP analysis, and −log10(0.01/33,010 SNP) for the SNP-HAP analysis.

2.5. Explained Genetic Variance by Haplotypes in HAP-PSEUDOSNP, SNPs in SNP-HAP and SNPs in SNP-ALL

The EXGV by one SNP or haplotype ( $[eqn]$ ) was defined as proposed Legarra et al. [32] and Abdel-Shafy et al. [33]:

$$[eqn]$$

where $[eqn]$ is the A allele frequency for the i-th marker (SNP or haplotype), $[eqn]$ is the B allele frequency for the i-th marker, and $[eqn]$ is the estimated effect of the i-th marker.

The EXGV for all SNP or haplotype was calculated, as done by Legarra et al. [32] and Lourenco et al. [28], with the following formula (components previously described in Formula (4)):

$$[eqn]$$

The EXGV calculated using this method assesses the contribution of SNPs or haplotypes to the additive genetic variance. This assessment is based on the estimated marker effect and the allele frequency, under the assumptions of the Hardy–Weinberg equilibrium and an additive polygenic model [33,34].

The calculation of the reliability for the GBV of each trait in the three analyses was carried out with the BLUPF90 methodology [35] and an ANOVA was performed to detect possible differences between the different genomic source information with SAS 9.4 [36]:

$$[eqn]$$

where

$[eqn]$ = Reliability of the GBV of the i-th animal for each trait.

$[eqn]$ = Prediction error variance for the GBV of the ith animal.

$[eqn]$ = Additive genetic variance.

$[eqn]$ = Inbreeding coefficient of the i-th animal.

3. Results

The SNP-ALL analysis was previously presented by Cortes-Hernández et al. [37], and this research addresses results to compare the differences with the other two analyses: HAP-PSEUDOSNP analysis, using haplotypes, and SNP-HAP analysis using only SNPs with high LD.

3.1. Descriptive Statistics of Haplotypes

A total of 11,788 haplotypes were identified, with an average of 3.48 ± 2.44 SNPs per haplotype, ranging from a minimum of 2 to a maximum of 59 SNPs. These haplotypes resulted in the generation of 35,552 pseudo-SNPs, a largest number than the original count of haplotypes due to the presence of multiple allele combinations within the population. The average physical length of the haplotypes was 41.06 Kb, with lengths ranging from 0.003 Kb to a maximum of 199 Kb. Longer Bos taurus autosomes (BTAs) 1, 2, 3, 4, and 5 harbored a higher number of haplotypes compared to the smaller BTAs, such as BTA 26, 27, 28, and 29. Notably, BTAs 1, 2, and 6 exhibited the highest number of haplotypes. Additionally, BTAs 14, 20, and 24 showed an elevated number of haplotypes relative to their neighboring BTA (Figure 1), these results are consistent with the distribution of the pseudo-SNP by BTA.

3.2. Trait-Associated Haplotypes and SNPs Identified by Three Genomic Analyses

In the GWAS analyses, no haplotypes or SNPs were found to be significantly associated with SCS (Figure 2, Figure A1, Figure A2 and Figure A3). In the HAP-PSEUDOSNP analysis, 53 haplotypes were significantly associated with milk production traits, based on the significance threshold (p > −log10(0.01/35,552)). These associations were distributed across BTA 3, 5, 6, 14, and 20 (Figure 2; Table A1). In comparison, the SNP-HAP analysis identified 69 SNPs significantly associated with milk production traits (Figure A1), while the SNP-ALL analysis revealed 162 significant SNPs associated with the same traits (Figure A2 and Figure A3).

From the associated haplotypes, 23 were only for PP and 16 for FP (Figure 3). Only two haplotypes (22111 and 11112, Table A1) in BTA 14 were associated with all traits (Figure 2), which were in the region of 1.51 to 1.69 Mb and 1.80 to 1.92 Mb.

In BTA 5, one haplotype was determined for two SNPs with two combinations highlighted, associated with FY and FP; in BTA 6, another haplotype of three SNPs with two combinations was associated with PP; and in BTA 20, one haplotype with nine SNPs and a single combination was associated with PP (Table A1).

In the SNP-HAP analysis, the 69 SNPs were associated with milk production traits (Table A2) distributed in BTAs 3, 5, 6, 12, 14, and 20, with only one SNP shared among the five milk production traits. The traits that presented unique SNP associated were FP and PP, whereas FY and PY did not show unique SNP associated with (Figure 4).

3.3. Explained Genetic Variance and Sum of Squared Effects by SNPs and Haplotypes in the Three Analyses

A higher estimation of EXGV across all traits was observed in the SNP-HAP analysis compared to the SNP-ALL analysis (Table 1). For instance, in the case of MY, the EXGV increased from 12,636.84 ± 0.233 in the SNP-ALL analysis to 16,012.56 ± 0.904 in the HAP-PSEUDOSNP analysis and further increased to 25,798.48 ± 1.266 in the SNP-HAP analysis even though the number of markers used in the SNP-HAP analysis was lower. This suggests that the inclusion of SNPs in high LD, as identified during haplotype construction, may enhance the estimation of EXGV by better capturing the joint effects of linked loci.

Similar to the increasing rates of the EXGV for all the traits with the SNP-HAP analysis, the sum of the squared effects of the markers to MY, FY, and PY was greater in the analyses of SNP-HAP compared with the SNP-ALL analyses, but not so different between HAP-PSEUDOSNP and SNP-HAP to FP, PP, and SCS (Table 2).

3.4. Reliability of the GBV for the Evaluated Traits in the Three Analyses

The average reliability of GBVs increased by up to 0.11 points in the SNP-HAP analysis compared to the SNP-ALL analysis, as observed for PY. This improvement is likely due to the high LD among SNPs used in the SNP-HAP analysis, which may better capture the underlying genetic architecture of the traits. Overall, the highest reliability levels were observed for percentage traits PP and FP followed by PY and FY. In contrast, SCS exhibited the lowest reliability, ranging from 0.61 ± 0.001 to 0.64 ± 0.001 across the analyses. The reliability estimates differed significantly among the three analytical approaches (p < 0.0001; Table 3).

4. Discussion

4.1. Associated Haplotypes and SNPs with the Evaluated Traits

The number of haplotypes identified per BTA in this study was consistent with findings reported in other Holstein cattle populations. For example, a study conducted on 450 Holstein sires in Spain found that the longer autosomes (BTAs 1 to 5) harbored a greater number of haplotypes compared to the shorter autosomes (BTAs 26 to 29) [38]. Similarly, the average number of SNPs per haplotype observed in our study (3.48 ± 2.44) was comparable to that reported in Portuguese Holstein cattle, where the mean was 4.21 across the population. In that study, the highest number of haplotypes was also observed on BTA 1 and BTA 27 [8].

The haplotypes identified on BTAs 6, 14, and 20 that were associated with the evaluated traits correspond to previously reported quantitative trait nucleotides (QTNs). One of the most notable is the DGAT1 gene (Diacylglycerol O-acyltransferase 1), located on BTA 14, which has been extensively studied and linked to MY, FY, PY, FP, and PP [39,40]. In addition, the GHR (Growth Hormone Receptor) gene on BTA 20, within the genomic region of 31.9 to 32.1 Mb, and the ABCG2 gene, located between 37.63 and 38.41 Mb on BTA 6, were regions where associated haplotypes and SNPs were detected. These genes have previously been reported to be linked to PP and productive lifespan [41].

A study conducted on U.S. Holstein cattle under approximately 40 years of artificial selection reported a high frequency of extended haplotype homozygosity on BTA 20, specifically within the 21 to 49 Mb region, because of intense genetic selection [42]. These findings align with those observed in the present study, likely due to the historical importation of U.S. genetic material into the Mexican Holstein population. Within this same genomic region (21–49 Mb), a haplotype spanning nine SNPs and located between 31.91 and 32.10 Mb was identified in the current population and found to be significantly associated with PP (Table A1). In contrast, Holstein cattle populations from China, which do not have a direct genetic link to the animals evaluated in this study, also exhibited associations between milk production traits and nine haplotypes located in the same region of BTA 20 [43]. This convergence of evidence suggests that this genomic region harbors haplotypes that are consistently associated with productive traits, likely due to the presence of key candidate genes involved in trait expression, such as SLC1A3, GHR, CCNB1, and NIPBL [40,43].

On BTA 3, at the 15.36 Mb region, two haplotype combinations (212222222 and 222212222) were identified, composed of a set of nine SNPs and associated with PP (Table A1). This genomic region has previously been reported to harbor over 20 QTLs associated with MY, PP, FY, and FP [44]. On BTA 5, within the 93.5 to 95.6 Mb region, haplotypes were associated with FY and FP. This region includes the ITPR2 gene, which has been linked to both MY and FY in earlier studies [45].

Although a high frequency of haplotypes was detected on BTA 24, none were significantly associated with the evaluated traits in this study. However, in Dutch Holstein cattle, loci on BTA 24 have been linked to natural antibodies in milk, which are associated with mastitis susceptibility, productive lifespan, and postpartum uterine health [46]. Additionally, a study conducted in Florida, USA, on Holstein cows reported loci in this same region related to pregnancy loss at 42 days in first- and second-lactation animals [47]. These findings may help explain the high haplotype frequency observed on BTA 24 despite the absence of trait associations in the current study.

Although no significant SNP or haplotype associations were detected for somatic cell score (SCS) across the evaluated analyses, this outcome is consistent with the highly polygenic architecture and relatively low heritability of SCS, which reduce the statistical power to detect individual marker or haplotype effects [33]. Therefore, the absence of significant associations for SCS should not be interpreted as a lack of genetic control, but rather as a limitation of marker-based association approaches for this trait.

Some haplotype blocks included a large number of SNPs (up to 59), with 12 haplotype combinations and population frequencies ranging from 1.3% to 21.8%. While such extended haplotypes may reflect genomic regions under strong historical selection, they may also increase the risk of overfitting by capturing multiple functional variants or genes with heterogeneous effects [48]. Consequently, the effects of long haplotypes should be interpreted with caution, particularly when their population frequency is low.

4.2. Explained Genetic Variance by SNPs, Haplotypes, and SNPs with High LD

The higher level of EXGV observed in the SNP-HAP analysis is likely attributable to the inclusion of SNPs with high LD (r^2^ > 0.80), which can enhance the detection of local additive effects among adjacent SNP [13]. These additive effects between closely linked markers have been shown to influence the overall EXGV estimation [49], supporting the findings obtained in the HAP-PSEUDOSNP analysis, where haplotypes associated with traits showed substantial explanatory power. For instance, a haplotype consisting of five SNPs located on BTA 14 (ARS-BFGL-NGS-4939, BovineHD1400000243, BovineHD1400000246, BovineHD1400000249, and Hapmap52798-ss46526455; Table A1) produced three combinations (22221, 11112, and 21112) with a population frequency greater than 1%. Among these, only the 11112 combination was significantly associated with all five evaluated traits. Notably, the individual SNPs comprising this haplotype were also independently associated with the same traits in both the SNP-ALL and SNP-HAP analyses (Table A2 and Table A3), reinforcing the importance of local LD and the additive interaction of tightly linked markers in explaining trait variation.

Unlike studies using genomic windows defined by a fixed number of SNPs [12,13], in which markers grouped within the same genomic interval may not necessarily be in LD or jointly influenced by the same underlying QTL, the use of LD-defined haplotypes may influence estimates of EXGV. Specifically, SNPs within high-LD blocks share correlated effects and are more likely to tag the same QTL [45], thereby increasing the contribution of selected markers to EXGV. As a result, EXGV values derived from LD-based groupings may be higher than those obtained using fixed-window approaches that include SNPs with weak or no LD relationships.

The estimation of EXGV was higher when using SNPs with high LD (SNP-HAP) and haplotypes coded as pseudo-SNPs (HAP-PSEUDOSNP) compared to the SNP-ALL analysis across most traits (Table 1). These results contrast with those reported by Abdel-Shafy et al. [33], who conducted a GWAS for SCS using daughter yield deviations (DYD) from German Holstein sires. In that study, the EXGV per significant SNP was approximately 2% higher than that obtained with significant haplotypes. The discrepancy with the German Holstein population may be attributed to differences in sample size (2354 genotypes), the lower number of SNPs used (44,576), and the nature of the phenotype (DYD), which may have limited the power to detect haplotype effects. Conversely, in Nelore cattle, it has been reported that haplotypes can capture epistatic interactions between variants within a haplotype locus, increasing additive genetic variance around 2% in some regions than individual SNPs for certain traits, such as meat tenderness [50]. However, the same study highlighted that longer haplotypes do not necessarily capture more genetic variance than shorter ones. This may be due to the use of fixed window sizes during haplotype construction, which does not account for regions of high LD. As a result, important marker interactions could be missed, leading to a reduced ability to detect genomic effects.

The magnitude of the EXGV observed for each trait was proportional to the magnitude of the phenotypic values. The gain in EXGV in the HAP-PSEUDOSNP analysis can be attributed to the fact that haplotype-based approaches profit larger marker effect estimates compared to individual SNPs (Table 2). This is primarily because haplotypes consist of multiple allelic combinations at a locus, increasing the likelihood that at least one of them is in linkage with a QTL. Moreover, haplotype-based models reduce the number of estimated effects while retaining relevant genomic information, effectively excluding rare alleles or those not associated with any QTL [4,51]. The greatest gains in EXGV were observed for MY, FY, and PY when using the SNP-HAP analysis. These improvements are likely due to the higher estimated effects of SNPs included in the SNP-HAP approach, which selectively incorporates markers in high LD and is more likely to capture the underlying genetic signal (Table 2).

4.3. Reliability of GBV in the Three Analyses and the Inclusion of SNPs with High LD in Genomic Prediction

The increase in the GBV reliability for FY, PY, and SCS from around 2% with the use of LD-defined haplotypes (HAP-PSEUDOSNP) compared to SNP-ALL analyses was found to be similar to that reported by other studies; for example, Cuyabano et al. [52] showed increases up to 3.1% for the reliability of GBV of PY when haplotypes were calculated with an r^2^ value for LD ≥ 0.75 in Nordic Holstein cattle, although they mention that increases in reliability of up to 3% can be achieved when a limit of r^2^ for LD ≥ 0.45 in low heritability traits such as fertility and mastitis. A factor to highlight in the Nordic population is the number of SNPs used for the construction of haplotypes, as they considered 492,057 SNPs and the methodology for estimated reliability. The greatest increases in the reliability of GBV were obtained with the SNP-HAP analyses, possibly since only SNPs with a high LD with r^2^ > 0.80 were included.

In contrast to the gains in reliability of GBV for the production traits with the use of haplotypes in Holstein cattle, in other cattle breeds, such as the Nelore, decreases were observed in the reliability of predictions in the same magnitude (1 to 2%) with the use of haplotypes for fatty acid traits of meat [3]. The authors attribute these results to the fact that the density of markers used was not enough to capture all LD between SNPs or haplotypes with some QTL region, because the total number of SNPs they used was 469,981 out of 893 genotypes.

The inclusion of SNPs with high LD in genomic prediction to increase the GBV accuracy and reliability of productive traits has been analyzed in other ways, for example, Mathew et al. [53] proposed the inclusion of a genomic relationship matrix adjusted by LD in the genomic prediction of different populations, showing positive results in increasing the reliability, and detecting higher EXGV, and therefore an increase of the heritability of the evaluated traits, mentioning that there are computational difficulties for the implementation of the process. Other studies suggest that the use of models with SNP selection, according to their stratification considering the LD level and the genetic structure of the trait, heritability level, and genes involved in the expression, would increase the percentage of EXGV [54], although these authors mention that these methods are not recommendable with medium-density (50 k) SNP arrays. In contrast to the results of this study about the increase in the reliability of GBV and EXGV for the traits evaluated in the SNP-HAP analyses, there were no computational problems and the number of markers used did not cause discrepant results between the evaluated traits.

On average, the use of SNPs with high LD (SNP-HAP) increases the EXGV twofold and the reliability of GBV up to 0.11, compared to the use of SNP-ALL. These increases were not found in analyses using haplotypes (HAP-PSEUDOSNP), but the findings are like those reported by other studies [51,52].

The highest increases in GBV reliability in the present study were observed in the SNP-HAP analysis, which selectively included only SNPs in strong LD (r^2^ > 0.80). This targeted selection likely enhanced predictive accuracy by focusing on SNPs more likely to tag relevant QTL regions.

The inclusion of high-LD SNPs in genomic prediction has also been explored through alternative approaches. Mathew et al. [53], for example, proposed incorporating an LD-adjusted genomic relationship matrix, which improved reliability and EXGV across diverse populations. While this approach enhanced trait heritability estimates, it was associated with notable computational challenges. Other research [54] has suggested that selecting SNPs based on LD stratification and the genetic aspects of traits (e.g., heritability, gene architecture) may improve EXGV. However, the authors cautioned against using such models with medium-density arrays (e.g., 50 K SNP chips), due to limitations in resolution and informativeness.

In contrast to those limitations, the present study demonstrated clear gains in both GBV reliability and EXGV using the SNP-HAP approach, without encountering computational difficulties or inconsistent results across traits even when using medium-density genotyping data, although the number of animal genotypes is not large.

Overall, the use of high-LD SNPs (SNP-HAP) resulted in approximately a twofold increase in EXGV and up to a 11% improvement in GBV reliability compared to the SNP-ALL analysis. While such gains were not fully replicated in the HAP-PSEUDOSNP analysis, the trends observed were consistent with findings reported by Cuyabano et al. [52] and Jónás et al. [51], reinforcing the potential of LD-informed marker selection strategies to enhance genomic prediction performance.

Although the observed 11% increase in GBV reliability may appear substantial, it is important to note that the set of markers used in the SNP-HAP analysis is not independent of the full SNPs dataset (SNP-ALL). Therefore, the observed improvement is likely driven by the strong LD among the selected SNPs.

Additionally, the higher EXGV estimated in the SNP-HAP analysis may be attributed to the larger estimated marker effects, which result from both the reduced number of markers and the strong LD structure among them. These factors directly influence the calculation of EXGV as defined in Formulas (4) and (5). Consequently, higher EXGV values obtained with this approach may not fully reflect the true genetic architecture of the evaluated traits, as the inclusion of a smaller subset of SNPs reduces genome coverage and may exclude relevant genomic regions contributing to trait expression.

5. Limitations

While the use of high-LD SNPs (SNP-HAP) can enhance genomic evaluations, this should be considered with caution because the set of markers in SNP-HAP represents a subset of LD-filtered SNPs, rather than an independent source of information; the improvements observed in GBV reliability are based on internal comparisons rather than external or cross-validation results, and the increase in EXGV could be reflecting inflation effects due to LD among SNPs. Additionally, several important limitations should be considered.

First, selection based on haplotypes may have antagonistic effects in certain genomic regions. For instance, a single haplotype that spans multiple genes or even different variants within the same gene may harbor alleles with opposing effects on the same trait. This could potentially reduce the overall genetic merit or increase the likelihood of expressing deleterious alleles [37,42,55]. As a result, caution is necessary when interpreting haplotype effects, especially in regions with complex gene architectures or pleiotropic effects.

Second, the generation of pseudo-SNPs from haplotypes is a labor-intensive process. It involves first defining haplotypes either by a fixed number of SNPs or by specific LD thresholds, and then reconfiguring both haplotype maps and chromosomal maps to incorporate these new markers into the genomic evaluation process. These additional steps increase the computational burden and may present challenges in routine applications.

Third, the use of SNPs in high LD for genomic evaluations requires constant updating of the selected SNP set. The addition of new markers or the inclusion of genotypes from newly analyzed animals can alter the LD structure, potentially impacting the effectiveness of previously selected SNPs. This ongoing need for recalibration makes the approach more complex and time-consuming [56]. And this also makes this process less predictable in populations with different levels of LD.

Fourth, the haplotype construction and pseudo-SNP encoding procedures were implemented using earlier versions of PLINK to ensure compatibility with established block-detection algorithms and downstream analytical pipelines. Although this choice does not affect the validity of the results, the use of more recent software versions (e.g., BEAGLE 5.5 [57]) for haplotype construction approaches may improve computational efficiency and reproducibility.

Lastly, although the use of haplotypes (like pseudoSNP) produced modest gains in EXGV and reliability of genomic predictions, the improvement over traditional SNP-based analyses was not always substantial. From a practical standpoint, the use of individual SNPs remains more straightforward, computationally efficient, and easier to implement in large-scale genetic evaluation programs.

6. Conclusions

Significant haplotypes in GWAS were found in the same regions as significant SNPs, and the presence of a greater number of haplotypes on BTAs 6, 14, and 20, compared to their neighboring BTAs, coincided with the presence of haplotypes associated with the evaluated traits; furthermore, associated haplotypes were also found in BTA 5.

The use of SNPs in high LD resulted in improved estimates of EXGV and marker effects across all six evaluated traits compared to using all individual SNPs (SNP-ALL). Moreover, this approach led to higher GBV reliability, particularly for FY and PY with reliability gains of up to 11%, suggesting that this method may offer a more efficient alternative for enhancing genomic predictions in Holstein cattle. However, it is essential to consider that LD patterns may change over generations, which could impact the consistency of this strategy in long-term selection schemes.

Although the use of haplotypes (like pseudoSNP) provided modest improvements in GBV reliability compared with the SNP-ALL analysis, these gains were limited and accompanied by increased computational and analytical complexity. Consequently, the HAP-PSEUDOSNP strategy should be regarded primarily as a methodological exploration rather than a practical alternative for routine genomic selection. In contrast, the SNP-HAP approach offers a more efficient and operationally feasible framework, as it leverages LD structure without requiring haplotype recoding, making it better suited for large-scale and routine genomic evaluation programs. However, further studies incorporating cross-validation or independent validation across different populations and traits are required to confirm the applicability of this methodology.