Longitudinal Reconstruction of Effective Population Size Reveals Severe Recent Genetic Erosion in Seven Indigenous Punjab Goat Breeds of Pakistan

1. Introduction

Effective population size (Ne) is one of the most informative parameters in conservation and quantitative genetics because it summarises, in a single value, the combined effects of genetic drift, inbreeding accumulation and the erosion of additive genetic variance within a population. Populations with small Ne lose genetic diversity faster, accumulate inbreeding more rapidly, and have reduced long-term capacity to respond to selection or to adapt to environmental and disease pressures. The Food and Agriculture Organization (FAO) and the Convention on Biological Diversity have therefore recommended Ne as a core indicator of genetic risk in livestock breeds, with Ne below 50 considered critical, Ne between 50 and 100 considered endangered, and Ne above 500 generally regarded as necessary to preserve long-term evolutionary potential [1,2].

Pakistan is the third-largest goat-producing country in the world, with an estimated national population of 78.2 million head, of which Punjab province accounts for approximately 37% [3,4]. Seven distinct, well-documented goat breeds are reared across Punjab: Beetal, Daira Din Panah, Nachi (Bikaneri), Barbari, Teddi, Pahari (Kajli) and Local Pothwari, each associated with a distinct home tract, morphology and production purpose ranging from large dual-purpose meat-and-milk types (Beetal) to small, highly prolific breeds favoured for religious sacrifice and smallholder systems (Teddi, Barbari) [5,6]. These breeds have recently been genotyped on the Illumina Caprine 50 K SNP Bead Chip and characterised for genetic diversity, population structure and genome-wide association with growth and conformation traits [7,8]. However, the historical demographic trajectory of these breeds—how their effective population size has changed from ancient founding events through to the present day—has not previously been reconstructed.

Genome-wide SNP data allow Ne to be estimated not only at the present generation but also at successive time points in the past, because the extent of linkage disequilibrium (LD) between markers separated by a given physical distance reflects the effective population size a corresponding number of generations ago: LD between closely linked markers persists over many generations and therefore reflects ancient Ne, whereas LD between more distant markers decays rapidly and reflects recent Ne. This LD-based approach, first formalised by Hayes et al. [9] and subsequently implemented in tools such as SNeP [10], has been widely applied to reconstruct the demographic history of cattle, sheep, pigs and, more recently, goat breeds, including a genome-wide survey of Spanish goat breeds [11]. Because the method requires no pedigree information, it is particularly well suited to smallholder-managed indigenous breeds such as those of Punjab, for which multi-generation pedigrees are rarely available.

The present study uses the same 50 K SNP genotype resource generated for the Punjab goat breeds to reconstruct breed-specific Ne trajectories spanning roughly 5,000 generations of history, quantify the magnitude of recent genetic erosion relative to ancestral population size, and rank breeds by their current genetic-risk status. These results are intended to complement previous diversity and association-mapping work on these breeds [7,8] and to provide an evidence base for breed-specific conservation and breeding-management recommendations.

2. Results

2.1. Sample Composition and Genotype Quality Control

A total of 827 goats belonging to seven Punjab breeds were genotyped on the Illumina Caprine 50 K SNP Bead Chip. Breed-specific quality control (SNP call rate ≥ 0.95, animal call rate ≥ 0.95, minor allele frequency ≥ 0.05) was applied independently within each breed prior to LD-based analysis, since per-breed allele frequency spectra and sample sizes differ substantially and pooling before filtering would bias rare-variant-driven LD estimates. Retained sample sizes ranged from 16 animals (Local Pothwari) to 588 animals (Beetal), and the number of SNPs passing quality control ranged from 32,083 (Nachi) to 40,013 (Beetal) (Table 1).

2.2. Historical Ne Trajectories

LD-based reconstruction across 18 physical-distance classes (10 kb to 4,500 kb midpoints, corresponding to approximately 11 to 4,992 generations ago) revealed a consistent demographic signature across six of the seven breeds: a moderate ancestral Ne of 3,300–4,500 individuals around 4,700–5,000 generations ago, a rise to a mid-history maximum of 3,800–8,300 individuals approximately 1,400–2,300 generations ago, and a subsequent, progressively steepening decline toward the present (Figure 1, Figure 2). Beetal showed the highest trajectory throughout the reconstructed history, peaking at Ne ≈ 8,308 around 1,729 generations ago before falling to Ne = 968 at the most recent distance class. Local Pothwari followed a similar rise-and-fall pattern with a lower peak (Ne ≈ 6,950). Daira Din Panah, Nachi, Pahari and Teddi showed a tightly clustered pattern with peaks of Ne ≈ 4,000–5,700 in the same 1,400–2,300-generation window.

Barbari was the only breed that did not display a mid-history expansion: its Ne declined almost monotonically from 2,748 at the earliest reconstructed time point to 63 at the most recent, indicating a long, sustained contraction rather than a recent-onset decline (Figure 3).

2.3. Recent Effective Population Size and Breed Ranking

At the most recent distance class analysed (mean physical distance ≈ 4,499 kb, corresponding to ~11 generations ago), estimated Ne differed by more than 15-fold among breeds (Table 2, Figure 4). Beetal retained by far the largest contemporary effective population size (Ne = 968), consistent with its status as the most numerous and widely distributed breed and its five recognised strains. Local Pothwari (Ne = 215) and Teddi (Ne = 139) occupied intermediate positions, while Nachi (119), Pahari (107), Daira Din Panah (91) and Barbari (63) all showed recent Ne below 150. Six of the seven breeds fall below the FAO Ne = 100 threshold associated with elevated short-term inbreeding risk, and Barbari and Daira Din Panah fall within the range flagged as approaching critical status.

2.4. Magnitude of Recent Genetic Erosion

Comparison of ancestral and recent Ne estimates (Table 2, Figure 5) showed that every breed has undergone a substantial contraction in effective population size, but the magnitude of decline was markedly uneven. Barbari (−97.7%) and Daira Din Panah (−97.2%) showed the most severe proportional erosion, followed closely by Nachi (−96.9%), Pahari (−96.8%) and Teddi (−96.8%). Local Pothwari, despite having a comparatively large ancestral Ne, also declined by 95.2%. Beetal showed the least severe proportional decline (−78.4%), consistent with its larger contemporary population and the maintenance of multiple recognised strains, but even this decline is substantial in absolute terms (a loss of 3,514 effective breeding individuals).

3. Discussion

This study provides the first genome-wide, LD-based reconstruction of historical effective population size for the documented goat breeds of Punjab, Pakistan. Two demographic signatures emerge from the analysis. First, six of the seven breeds share a common deep-history pattern of a moderate ancestral Ne, a mid-history expansion peaking roughly 1,400–2,300 generations ago, and a subsequent decline a trajectory broadly consistent with early post-domestication population growth followed by breed formation and, in the most recent period, contraction associated with modern livestock management. Second, and more strikingly, all seven breeds show a pronounced and accelerating loss of effective population size in the most recent generations, with proportional declines of 78–98% relative to ancestral levels. Because the LD-based method places this decline predominantly within the last few hundred generations well within the era of organised breeding, market-driven selection and, for several breeds, restriction to a narrow home tract the pattern points to recent breeding and management practices, rather than ancient founder effects, as the principal driver of contemporary genetic erosion.

The breed-level heterogeneity in current Ne has direct conservation implications. Beetal, the largest and most widely distributed breed with five recognised strains (Faisalabadi, Nuqri, Nagri, Gujrati and MakhiChini) and preferred for its size at Eid-ul-Azha, retains a comparatively robust recent Ne of 968, likely reflecting both its larger census population and greater effective breeding-male diversity across strains and sub-populations. In contrast, Barbari and Daira Din Panah both minor breeds with restricted home tracts in D.G. Khan, Rajanpur, Layyah and Muzaffargarh—show recent Ne below 100, placing them in the FAO endangered category and indicating an elevated risk of inbreeding depression, reduced fertility and diminished response to future selection if current management continues unchanged. Nachi, Pahari and Teddi occupy an intermediate but still concerning position, with recent Ne between 100 and 150.

The unique demographic profile of Barbari a long, continuous decline without the mid-history expansion phase seen in the other breeds may reflect its small and geographically restricted population throughout its history, its comparatively recent formal recognition as a distinct breed, and its selection primarily for milk and meat traits within closed, small-scale flocks resembling a deer in conformation. This continuous-decline signature contrasts with the expansion-then-contraction pattern of the other breeds and suggests that Barbari's conservation priority should not be interpreted solely from its percentage decline (which is comparable to several other breeds) but also from its consistently small absolute Ne across its reconstructed history.

These findings complement previous work on the same genotype resource that characterised genetic diversity, population structure and genome-wide associations with body weight and conformation traits in these breeds [7,8]. Taken together, the two analyses indicate that despite meaningful phenotypic and genomic differentiation among the seven breeds, several are simultaneously experiencing rapid contemporary loss of effective population size—a combination that increases the urgency of coordinated, breed-specific conservation action. Practical measures indicated by the present results include: (i) expanding the pool of registered breeding males and equalising male reproductive contribution, particularly in Barbari, Daira Din Panah and Nachi; (ii) establishing rotational or cooperative mating schemes across flocks within a breed's home tract to reduce co-ancestry; (iii) prioritising cryopreservation of semen and, where feasible, embryos from genetically distinct sires as an ex situ safeguard, especially for breeds with recent Ne below 100; and (iv) incorporating genomic Ne monitoring into future sampling rounds to track whether management interventions arrest the observed decline.

Several limitations should be considered when interpreting these results. LD-based Ne estimates for the most ancient distance classes are sensitive to sample size, and the smallest breed sample (Local Pothwari, n = 16) yields wider uncertainty around its historical trajectory than the larger Beetal and Teddi samples; the sample-size correction applied here (adjusted r² = mean r² − 1/n) mitigates but does not eliminate this effect. LD-based Ne also reflects the breeding population contributing to the genotyped sample and may not fully capture unsampled sub-populations or recent cross-breeding. Nonetheless, the consistency of the qualitative pattern—universal recent decline, breed-specific magnitude tracking known differences in population size and management—across breeds of very different sample size and history lends confidence to the overall conclusions.

4. Conclusions

Genome-wide LD-based reconstruction shows that all seven documented goat breeds of Punjab have experienced a severe contraction in effective population size in the recent past, with declines of 78–98% relative to ancestral levels and six of seven breeds currently below the FAO Ne = 100 threshold for elevated inbreeding risk. Beetal remains comparatively secure, but Barbari, Daira Din Panah, Nachi, Pahari and Teddi warrant urgent, breed-specific conservation attention, including expansion of the effective breeding-male pool, structured mating schemes, and ex situ genetic banking. These findings establish a quantitative baseline against which the effectiveness of future conservation interventions in Punjab's indigenous goat genetic resources can be measured.

5. Methods

5.1. Ethics Declaration

This study reused genotype and phenotype data collected from a survey of goat-keeping farmers together with blood sampling of their animals. Participants provided verbal informed consent for animal blood sampling and for related survey questions. Collection of blood samples and all animal-handling procedures were performed in accordance with relevant guidelines and regulations for animal welfare, and the study was carried out in compliance with the ARRIVE guidelines.

5.2. Breeds and Sampling

Blood samples were collected by jugular venepuncture into EDTA vacutainers from unrelated animals of all seven documented Punjab goat breeds: Beetal (including its five recognised strains Faisalabadi, Nuqri, Nagri, Gujrati and MakhiChini), Teddi, Pahari (Kajli), Local Pothwari, Daira Din Panah, Barbari and Nachi (Bikaneri). Sampling was carried out across the principal home tracts of each breed in Punjab, and animal identity was recorded to ensure that sampled individuals were unrelated and representative of each breed.

5.3. Genotyping and Quality Control

Genomic DNA was extracted from whole blood and genotyped using the Illumina Caprine 50 K SNP Bead Chip, which covers 53,347 SNPs distributed across the whole caprine genome (ARS1 assembly). After merging, 827 unrelated animals across the seven breeds passed initial quality control. For the present LD-based Ne analysis, quality control was subsequently applied independently within each breed in PLINK [12], removing SNPs with a genotyping call rate below 0.95 (--geno 0.05), animals with a call rate below 0.95 (--mind 0.05), and SNPs with a minor allele frequency below 0.05 (Table 3). Breed-specific filtering was used because pooling animals of very different sample size prior to MAF filtering would differentially bias the retained SNP set and downstream LD estimates for the smaller breeds.

5.4. LD-Based Estimation of Historical Effective Population Size

For each breed, pairwise LD (r²) was calculated between autosomal SNPs across 18 physical-distance classes spanning 0–5,000 kb (bin boundaries: 0–20, 20–30, 30–40, 40–50, 50–75, 75–100, 100–150, 150–200, 200–300, 300–400, 400–500, 500–750, 750–1,000, 1,000–1,500, 1,500–2,000, 2,000–3,000, 3,000–4,000 and 4,000–5,000 kb). Mean and median r² were computed per distance class per breed, and mean r² was corrected for finite sample size using adjusted r² = mean r² − 1/n, where n is the number of genotyped animals retained for that breed after quality control. For each distance class, the corresponding genetic distance (c, in Morgans) was derived from the mean physical distance assuming a constant recombination rate of 1 cM/Mb, and effective population size at the generation corresponding to that distance class was estimated as Ne = (1/r² − 1) / (4c), following the standard LD-based Ne framework of Hayes et al. [9], with the number of generations before present corresponding to a given distance class estimated as t = 1 / (2c). This approach yields, for each breed, a continuous trajectory of Ne estimates from approximately 11 generations ago (largest physical distance, 4,000–5,000 kb) to approximately 5,000 generations ago (smallest physical distance, 0–20 kb). ‘Recent Ne’ in the Results refers to the estimate at the most recent distance class available for each breed; ‘ancestral Ne’ refers to the estimate at the 0–20 kb (deepest) distance class.

5.5. Statistical Software

Quality control and LD calculations were performed in PLINK [12]. Historical Ne trajectories, breed-specific plots, and comparative ancestral-versus-recent summaries were generated in R.