Genome-Wide Patterns of Linkage Disequilibrium Decay Across Seven Indigenous Goat Breeds of Punjab, Pakistan

1. Introduction

Linkage disequilibrium (LD) the non-random association of alleles at physically linked loci underlies both genome-wide association studies (GWAS) and genomic selection. The rate at which LD decays with physical distance determines the marker density required for a SNP panel to adequately tag causal variants: breeds or populations with slowly decaying LD can be characterized with sparser panels, whereas rapidly decaying LD demands denser genotyping to maintain equivalent prediction accuracy [1,2,3]. LD patterns are themselves shaped by effective population size (Ne), mating system, selection history and admixture, so characterizing LD decay also offers an indirect window onto a breed’s recent demographic history [4].

Pakistan maintains one of the world’s most diverse goat populations, with Punjab province alone hosting seven documented breeds Beetal, Daira Din Panah, Nachi (Bikaneri), Barbari, Teddi, Pahari (Kajli) and Pothwari that differ markedly in body size, coat type, production purpose and, presumably, effective population size [5]. These breeds were recently genotyped on the Illumina Caprine 50K SNP Bead Chip and used to identify genomic regions associated with body weight and conformation traits [6]. That study characterised trait-associated loci but did not examine the underlying LD architecture of the population information that is a prerequisite, rather than a byproduct, of any future genomic selection programme for these breeds.

No genome-wide LD decay analysis has previously been reported for any Pakistani goat breed. Comparable work in other populations including Italian and Moroccan native goats [7], Boer, Saanen, Alpine and Cashmere breeds [8], and Chinese indigenous breeds [9] consistently finds substantial inter-breed variation in LD decay rate, generally attributable to differences in Ne and breeding history rather than to marker density or genotyping platform. Given the contrasting demographic profiles of Punjab’s goat breeds Beetal is a large-scale, multi-strain commercial breed reared across several districts, whereas breeds such as Barbari, Daira Din Panah, Nachi and Local Pothwari are comparatively localised and numerically smaller we hypothesised that LD decay rate would differ substantially among them.

Beyond the conventional decay-curve and threshold-distance approaches used in prior goat LD studies, we additionally summarise overall LD magnitude for each breed using the area under the LD-decay curve (AUC) a single integrated statistic that is less sensitive than a single threshold crossing to noise in any one distance bin and quantify, at each physical distance, the deviation of each breed’s mean r² from the across-breed mean (Δr²), which highlights where in the genome breed differences are largest. We also explicitly test whether genome-wide marker density is sufficiently uniform across chromosomes that it could not, by itself, generate the observed inter-breed LD differences. Using the Illumina Caprine 50K genotypes from the Punjab goat population described by Moaeen-ud-Din, Muner & Khan (2022) [6], we therefore (i) quantify genome-wide and chromosome-wise LD decay for each of the seven breeds; (ii) estimate the physical distance at which mean r² falls below the commonly used genomic-selection thresholds of 0.5, 0.2, 0.1 and 0.05; (iii) rank breeds by overall LD magnitude (AUC) and quantify the pattern of between-breed divergence across the full 1-Mb window; (iv) confirm that this ranking is not an artefact of uneven SNP density; and (v) discuss the implications of breed-specific LD architecture for the design of future genomic selection and fine-mapping panels in Pakistani goats.

2. Results

2.1. Sample Composition and Quality Control

A total of 827 goats across seven Punjab breeds Beetal (n = 594), Teddi (n = 109), Pahari (n = 41), Barbari (n = 23), Daira Din Panah (n = 23), Nachi (n = 21) and Local Pothwari (n = 16) formed the starting dataset (Table 1). After breed-wise genotype quality control, 817 individuals and between 32,083 and 40,013 autosomal SNPs per breed were retained for LD estimation (Table 2). Beetal, as the numerically dominant breed, retained the largest sample (n = 588) and the largest SNP set after filtering (40,013), while the five minor breeds each retained between 16 and 23 animals.

2.2. Genome-Wide LD Decay Differs More than Ten-Fold Among Breeds

Mean pairwise r² across all SNP pairs within 1 Mb differed markedly among breeds (Table 3; Figure 1). Beetal showed by far the lowest genome-wide mean LD (r² = 0.0201, median = 0.0073), consistent with its large effective sample size and its status as a widely distributed, multi-strain commercial breed. Barbari showed the highest mean LD (r² = 0.2084, median = 0.1211) more than ten-fold greater than Beetal with Daira Din Panah (r² = 0.1310) and Nachi (r² = 0.1313) forming an intermediate-high tier, and Local Pothwari (r² = 0.1063), Pahari (r² = 0.0920) and Teddi (r² = 0.0801) forming an intermediate-low tier. The proportion of SNP pairs retaining strong LD (r² ≥ 0.5) mirrored this ranking, from 0.20% in Beetal to 12.46% in Barbari, while the proportion at r² ≥ 0.2 ranged from 0.89% (Beetal) to 37.03% (Barbari).

2.3. Overall Magnitude of LD (Area Under the Decay Curve) Confirms the Breed Ranking

To summarise overall LD magnitude with a single statistic that is less sensitive to noise in any one distance bin than a threshold-crossing distance, we computed the area under the LD-decay curve (AUC, r²×kb) for each breed across the full 0–1000 kb window (Table 4; Figure 2). This confirmed the ranking observed for mean r²: Barbari showed by far the largest AUC (212.1), roughly 8.6-fold greater than Beetal (24.5). Notably, Nachi (AUC = 134.2) and Daira Din Panah (AUC = 134.1) were almost indistinguishable from one another despite being independently sampled breeds, consistent with their overlapping home tracts along the Indus (Muzaffargarh–Layyah–Multan) [6] and a plausibly shared or recently diverged effective population. Local Pothwari (110.8), Pahari (96.2) and Teddi (83.5) formed a distinct lower tier, and Beetal (24.5) remained a clear outlier at the low end.

2.4. Distance-Class and Short-Range Decay Patterns

All seven breeds showed the expected sharp initial decline in LD within the first 10–25 kb (Supplementary Figure S1), but the rate and asymptote of decay diverged sharply thereafter. Mean r² in the 0–10 kb bin ranged narrowly from 0.499 (Beetal) to 0.634 (Daira Din Panah), showing that short-range haplotype structure at the shortest distances is broadly comparable across breeds. By the 25–50 kb bin, however, Barbari retained more than three times the LD of Beetal (0.278 vs. 0.080), and by the 500–1000 kb bin this ratio widened to more than fourteen-fold (0.193 vs. 0.014; Table 5).

2.5. LD at Specific Physical Distances Relevant to Marker-Panel Design

Table 6 reports mean r² for each breed at seven physical distances spanning the range typically relevant to SNP-panel design (10, 20, 50, 100, 200, 500 and 1000 kb), corresponding to the heatmap in Figure 3. At 10 kb, LD was substantial in every breed (r² = 0.25–0.45), showing that short-range haplotype blocks are preserved across the whole population regardless of breed-level demographic differences; divergence between breeds became visually and numerically pronounced from 20–50 kb onward. At 1 Mb, Barbari retained r² = 0.158, roughly nine times the residual LD retained by Beetal (r² = 0.018) at the same distance.

2.6. Distance at Which LD Falls Below Genomic-Selection Thresholds

The physical distance at which mean r² declined below the commonly used thresholds of 0.20, 0.10 and 0.05 differed by more than an order of magnitude between the two extreme breeds (Table 7). In Beetal, mean r² fell below 0.20 and below 0.10 within the first 20 kb, and below 0.05 by 70 kb. In Barbari, by contrast, mean r² did not fall below 0.20 until approximately 500 kb, and never fell below 0.10 or 0.05 within the 1-Mb window examined. Daira Din Panah and Nachi reached the 0.20 threshold by 30 kb but likewise did not reach the lower thresholds within 1 Mb. Pahari, Local Pothwari and Teddi were intermediate, reaching r² = 0.20 by 20 kb but requiring 90–400 kb to reach r² = 0.10, and only Pahari and Teddi reached r² = 0.05 within the 1-Mb window (at the outer edge, ~1000 kb).

2.7. Breed-Specific Decay Profiles Are Consistent Across the Full 1-Mb Window

Plotting each breed’s decay curve on an independent panel (Figure 4) confirmed that the qualitative decay shape — a steep initial drop followed by a long, shallow tail — was shared by all seven breeds, but that the height and slope of the shallow tail differed systematically. Barbari was unique in retaining a tail that hovered around, and for part of the range above, the r² = 0.20 reference line for the entire 1-Mb window; every other breed’s tail settled below 0.20 within the first 30–50 kb and thereafter declined only gradually.

2.8. Between-Breed Differences in LD Are Largest at Intermediate Distances and Persist to 1 Mb

To quantify where along the physical-distance axis breed differences were largest, we computed, at each 10-kb bin, the difference between the highest and lowest breed mean r² (Table 8; Supplementary Figure S6). This range was smallest at the shortest distance examined (0.135 at 0–10 kb, reflecting the relatively uniform short-range LD noted above) and rose sharply to approximately 0.196–0.209 by 20–300 kb, after which it declined only slightly, remaining above 0.16 out to 1 Mb. Beetal occupied the lowest-LD position at essentially every distance from 10 kb to 900 kb, while the highest-LD breed was Daira Din Panah at the shortest distance (0–10 kb) and Barbari at every distance from 10 kb outward; at the extreme edge of the window (1000 kb, based on very few SNP pairs) Local Pothwari and Pahari occupied the highest- and lowest-LD positions respectively, reflecting the reduced reliability of single-pair bins at this distance.

Figure 5. Deviation of each breed’s mean LD (Δr²) from the across-breed mean at each physical distance, illustrating that Barbari’s excess LD narrows only slightly with distance while Beetal’s deficit remains essentially constant across the full 1-Mb window.

Figure 5. Deviation of each breed’s mean LD (Δr²) from the across-breed mean at each physical distance, illustrating that Barbari’s excess LD narrows only slightly with distance while Beetal’s deficit remains essentially constant across the full 1-Mb window.

2.9. Chromosome-Wise Consistency of the Breed Ranking

Chromosome-wise mean LD (Figure 6) showed the same breed ranking observed genome-wide, on every one of the 30 goat autosomes examined. Barbari showed the highest per-chromosome mean r² throughout (range 0.164–0.247, overall mean 0.209), Beetal the lowest throughout (range 0.015–0.027, overall mean 0.020), and the remaining five breeds occupied the same relative order — Daira Din Panah and Nachi (both ~0.131 genome-wide) above Local Pothwari (0.107), which was above Pahari (0.093) and Teddi (0.080) — on essentially every chromosome (Table 9). No chromosome showed a qualitatively different pattern from the genome-wide average for any breed, indicating that the observed LD differences reflect population-wide processes (effective population size, mating structure) rather than localised structural or selective effects on any single chromosome.

2.10. Genome-Wide Marker Density Is Uniform Across Chromosomes and Does Not Confound the Breed Comparison

Because LD estimates from array data can in principle be influenced by uneven marker coverage, we examined SNP density across the genome using the full autosomal marker map (38,890 SNPs; Table 10; Supplementary Figures S4–S5). SNP density per chromosome ranged narrowly from 14.15 SNPs/Mb (chromosome 18) to 17.43 SNPs/Mb (chromosome 12), around a genome-wide mean of 16.20 SNPs/Mb — a difference of less than 25% between the sparsest and densest chromosomes, and with no systematic relationship to chromosome length. Density within 1-Mb windows along each chromosome was likewise reasonably even, without the large gaps or dense clusters that would be expected to bias LD estimation in specific regions (Supplementary Figure S4). This close-to-uniform marker coverage indicates that the more than ten-fold breed differences in LD magnitude reported above are very unlikely to be an artefact of uneven SNP placement on the array, and instead reflect genuine differences in breed-level haplotype structure.

3. Discussion

This study provides the first genome-wide characterization of LD decay in Pakistani goat breeds, and extends the conventional decay-curve and threshold-distance approach with two additional, complementary summaries: an integrated area-under-the-curve (AUC) statistic and a distance-resolved measure of between-breed divergence (Δr²). All three approaches — mean r², AUC and Δr² — converged on the same central finding: a more than ten-fold difference in overall LD magnitude between Beetal and Barbari, with a consistent five-tier ranking reproduced across every distance bin and every autosome. This convergence across independent summary statistics increases confidence that the ranking reflects a genuine population-genetic signal rather than an artefact of any single analytical choice.

The practical implication is direct. Under the classical relationship between marker density and effective population size (Ne ≈ 1 / (4c·Δr²) for recombination distance c) [13], populations with rapidly decaying LD, such as Beetal, require markers spaced far more densely than the ~60 kb average spacing of the 50K chip to fully capture haplotype blocks, whereas populations retaining LD to hundreds of kilobases, such as Barbari, could in principle be characterised adequately with a substantially sparser and cheaper panel. The AUC ranking (Table 4) makes this practical contrast especially concrete: Barbari’s AUC of 212.1 r²×kb is 8.6-fold greater than Beetal’s 24.5, a gap of similar order to the ten-fold difference in mean r², confirming that the result is not sensitive to whether LD magnitude is summarised by a simple mean or by the full area under the decay curve.

The rapid decay observed in Beetal is consistent with its status as Punjab’s largest and most widely distributed goat breed, comprising several documented strains (Faislabadi, Nuqri, Nagri, Gujrati and MakhiChini) reared across multiple districts [6], which together imply a comparatively large effective population size and reduced genetic drift; notably, the Δr² analysis (Figure 5) shows that Beetal’s deficit relative to the across-breed mean is essentially flat across the entire 1-Mb window (roughly −0.08 to −0.09 from 20 kb outward), indicating that its low-LD architecture is a stable, distance-independent property of the breed rather than a transient effect confined to any particular scale. Conversely, the slow decay and elevated LD observed in Barbari — which did not fall below r² = 0.20 until approximately 500 kb, and whose excess over the across-breed mean narrows only gradually from roughly +0.12 at 20 kb to +0.09 near 900 kb — is consistent with a smaller, more geographically restricted effective population subject to greater drift and possibly historical bottlenecks or closer within-flock mating.

A striking secondary observation is the near-identity of Nachi and Daira Din Panah on every summary statistic examined — genome-wide mean r² (0.1313 vs. 0.1310), AUC (134.2 vs. 134.1) and chromosome-wise mean (0.1308 vs. 0.1307). Both are hairy, meat-type breeds native to overlapping tracts along the Indus (Muzaffargarh–Layyah–Multan) [6], and this degree of convergence across independent breeds and independent samples is more plausibly explained by a shared or recently diverged effective population — consistent with reported gene flow or common ancestry between geographically adjacent Punjab breeds [5,11] — than by coincidence, and would be a natural target for a dedicated population-structure analysis (e.g., FST, admixture or identity-by-descent sharing) in future work.

The between-breed range analysis (Table 8; Supplementary Figure S6) showed that breed differences were smallest at the very shortest distance examined (0–10 kb, range = 0.135) and expanded rapidly to a broad plateau of approximately 0.18–0.21 from roughly 20 kb out to 1 Mb. This pattern indicates that very short-range haplotype structure (co-inherited blocks smaller than ~10–20 kb) is comparatively conserved across all seven breeds — plausibly reflecting shared ancestral recombination architecture common to the species — whereas breed-specific demographic history becomes the dominant determinant of LD from approximately 20 kb outward. This has a direct practical corollary: a marker panel dense enough to tag the shortest-range haplotype blocks (a few SNPs per 10 kb) would perform comparably across breeds, but a sparser panel relying on LD extending beyond 20–50 kb would perform very differently by breed, remaining effective in Barbari, Nachi and Daira Din Panah but failing in Beetal.

A necessary robustness check for any LD comparison based on SNP-array genotypes is that observed differences are not an artefact of uneven marker placement. Genome-wide SNP density in the present dataset was closely uniform across all 29 autosomes (14.15–17.43 SNPs/Mb, a less than 25% range around a genome mean of 16.20 SNPs/Mb; Table 10), with no systematic relationship to chromosome length or to the chromosome-wise LD pattern in Figure 6 (the highest-LD chromosome for Barbari, chromosome 6, and the lowest-density chromosome, chromosome 18, did not coincide with any single breed’s LD extremes). This makes it very unlikely that the ten-fold breed differences reported here reflect uneven SNP coverage rather than genuine population-level differences in haplotype structure.

A caveat that must be weighed against these interpretations is sample size. Five of the seven breeds (Barbari, Daira Din Panah, Nachi, Pahari and Local Pothwari) were represented by only 16–41 animals, and r² is a known upward-biased estimator of true population LD when sample size is small [14]. This bias likely inflates the absolute magnitude of LD in these breeds to some degree, and the wide confidence bands visible for the smaller breeds in Supplementary Figure S2 reflect this reduced precision directly; the single-SNP-pair bins at the extreme edge of the 1-Mb window (e.g., the 1000-kb points in Table 6 and Table 8) are particularly unreliable and should not be over-interpreted. However, several considerations argue that the qualitative ranking reported here is unlikely to be a pure sample-size artefact. First, Beetal — with the largest sample by far (n = 588) — showed a genome-wide LD level roughly consistent with what has been reported in other large, outbred commercial goat populations genotyped on the same 50K platform [7,8,9], arguing against a systematic inflation specific to this dataset. Second, the five small-sample breeds spanned a more than 60% range in mean r² (0.080–0.131) and nearly a three-fold range in AUC (83.5–134.2) despite comparable sample sizes (16–41), indicating that sample size alone does not explain the observed ordering among them. Third, the chromosome-wise analysis (Figure 6, Table 9) and the AUC ranking (Table 4) independently reproduced the identical breed ordering, which would be an unusual coincidence if the pattern were driven purely by estimation noise. Nonetheless, we recommend that future work prioritise expanding sample sizes for Barbari, Daira Din Panah, Nachi and Local Pothwari before these LD estimates are used to fix a genomic-selection panel density, and that r² estimates for these breeds be treated as upper bounds pending validation in larger cohorts.

These results complement the recent GWAS of growth and conformation traits conducted in the same Punjab goat population [6], which reported two SNPs (on chromosomes 8 and 16) associated with all six traits examined. The present LD analysis suggests that, in Beetal — the breed contributing 594 of the 827 animals in that GWAS — the effective resolution of the 50K chip for fine-mapping such pleiotropic loci is high, given that LD decays to background levels within tens of kilobases; associated SNPs in Beetal are therefore likely to lie close to, rather than merely linked with, the underlying causal variant. In the smaller breeds, however, the persistence of LD over hundreds of kilobases means that a significant SNP is markedly less likely to pinpoint the causal gene directly, and functional follow-up in those breeds should anticipate a substantially wider candidate interval.

This study has several limitations beyond sample size. First, LD was estimated using array-derived genotypes at a fixed marker density (32,083–40,013 SNPs per breed after filtering; 38,890 SNPs in the full pre-breed-QC marker map), which imposes a lower bound on the resolution at which very short-range LD (<10 kb) can be characterised; whole-genome resequencing would be required to resolve haplotype structure at base-pair resolution. Second, we report composite pairwise r² without formally partitioning the contribution of population stratification, admixture or recent bottlenecks to the observed patterns; a targeted analysis of runs of homozygosity and effective population size trajectories, ideally combined with pedigree records where available, would clarify the demographic mechanisms underlying the LD differences reported here and would be particularly informative for testing the shared-ancestry hypothesis raised above for Nachi and Daira Din Panah. Third, breed assignment relied on farmer-reported breed identity at the time of sampling [6], and any misclassification would tend to inflate apparent within-breed LD by mixing distinct genetic backgrounds.

4. Conclusions

Genome-wide LD decays more than ten-fold faster in Beetal than in Barbari, with Daira Din Panah, Nachi, Local Pothwari, Pahari and Teddi occupying a consistent intermediate ranking that is reproduced independently by mean r², by the area under the LD-decay curve, and on every autosome, and that is not attributable to uneven genome-wide marker density. Between-breed differences are smallest at the shortest distances examined (<10 kb) and largest from approximately 20 kb to 1 Mb, indicating that short-range haplotype structure is broadly conserved across breeds while longer-range LD is dominated by breed-specific demographic history. These breed-specific LD architectures should directly inform the design of future genomic selection and fine-mapping panels for Pakistani goats: Beetal will require considerably denser marker panels than the 50K chip to fully exploit LD-based genomic prediction, while the slower-decaying breeds may be tractable with sparser, lower-cost panels, subject to confirmation in larger sample sizes. This work provides the first such baseline for Punjab’s indigenous goat breeds and is intended to support subsequent genomic selection, effective population size and conservation-genetics studies in this population.

5. Methods

5.1. Study Population and Genotype Data

Genotype data were drawn from the Punjab goat population described in detail by Moaeen-ud-Din, Muner & Khan (2022) [6], comprising 827 unrelated animals from seven documented Punjab breeds — Beetal (including its five strains: Faislabadi, Nuqri, Nagri, Gujrati and MakhiChini; n = 594), Teddi (n = 109), Pahari/Kajli (n = 41), Barbari (n = 23), Daira Din Panah (n = 23), Nachi/Bikaneri (n = 21) and Local Pothwari (n = 16), sampled from districts including Fateh Jhang, Bhakkar, Layyah, D.G. Khan, Rajanpur, Faisalabad and Jhang. Genomic DNA was extracted from EDTA whole blood and genotyped on the Illumina Caprine 50K SNP Bead Chip (53,347 SNPs) by GeneSeek (USA), as originally reported [6]. All animal sampling and handling followed protocols approved by the institutional Ethics Committee of the collaborating university and complied with the ARRIVE guidelines.

5.2. Quality Control

Starting from the shared genotype dataset, quality control was applied independently within each breed, retaining SNPs and individuals meeting standard genotype call-rate, minor allele frequency and Hardy–Weinberg equilibrium criteria in PLINK v1.9 [15,16]. Breed-wise filtering retained between 32,083 (Nachi) and 40,013 (Beetal) autosomal SNPs and between 16 (Local Pothwari) and 588 (Beetal) individuals per breed (Table 2); genome positions were assigned according to the goat reference assembly ARS1.

5.3. Linkage Disequilibrium Estimation

Pairwise composite linkage disequilibrium (r²) between all SNP pairs within 1 Mb of physical distance was estimated separately for each breed in PLINK v1.9 (--r2). SNP pairs were assigned to 10-kb physical-distance bins from 0 to 1000 kb, and mean, median, standard deviation, standard error, interquartile range and approximate 95% confidence intervals of r² were computed within each bin. Bin-level estimates were further aggregated into eight broader distance classes (0–10, 10–25, 25–50, 50–100, 100–250, 250–500, 500–1000 and >1000 kb) for summary comparison across breeds (Table 5), and mean r² at seven selected physical distances (10, 20, 50, 100, 200, 500 and 1000 kb) was extracted for cross-breed comparison (Table 6). For each breed, the physical distance at which mean r² first declined below the thresholds of 0.50, 0.20, 0.10 and 0.05 was identified from the binned decay curve (Table 7); where a threshold was not reached within the 1-Mb window examined, this is reported as ‘not reached’.

5.4. Overall LD Magnitude and Between-Breed Divergence

To provide a single integrated summary of LD magnitude for each breed, the area under the LD-decay curve (AUC, in units of r²×kb) was computed by trapezoidal integration of mean r² across the 10-kb bins spanning 0–1000 kb (Table 4). To identify the physical distances at which breeds diverged most, the across-breed mean r² was computed at each 10-kb bin, and each breed’s deviation from this mean (Δr² = breed mean r² − across-breed mean r²) was calculated (Figure 5). Independently, at each 10-kb bin the highest and lowest breed mean r² and their difference (the between-breed range) were identified (Table 8; Supplementary Figure S6).

5.5. Chromosome-Wise LD and Marker Density

LD was additionally summarised separately for each of the 30 goat autosomes to assess whether genome-wide breed differences were consistent across chromosomes (Table 9). To assess whether the observed breed differences could instead reflect uneven marker coverage, SNP density was computed for each chromosome as the number of SNPs on the full pre-breed-QC autosomal marker map (38,890 SNPs) divided by chromosome length in megabases (Table 10), and additionally tabulated in non-overlapping 1-Mb physical windows across each chromosome (Supplementary Figure S4) to check for local gaps or clusters in marker coverage.

5.6. Data Visualization

LD decay curves, short-range decay plots, breed-faceted decay panels, confidence-interval bands, chromosome-wise LD plots, the LD heatmap, the AUC ranking plot, the Δr² plot and the SNP-density plots were generated in R (ggplot2). All summary statistics reported in the Results are computed directly from the accompanying data tables and are reproduced here to two or three significant figures for readability.