Ethical approval was obtained at the ARC Animal Production Ethics Committee (APIEC23/07). The carcasses were classified as WCS-affected or unaffected based on abattoir grading and condemnation by the meat inspector, 24 h post slaughter at the abattoir. Also, both case and control animals were sampled from the same cohort to ensure consistency and comparability. Due to the sporadic nature of the condition, it is challenging to target specific breeds for sampling. Meat samples were collected from WCS-affected and unaffected carcasses at two commercial abattoirs: KLK Upington Abattoir (Northern Cape) and Cavalier Abattoir (Gauteng). Strict standard operating protocols were executed to ensure there was no cross contamination between tissue samples and between carcasses. The samples were taken from the hind leg (biceps femoris muscle) of the carcasses. Each sample was placed in an individual tube and stored at −80°C until further analysis.

The study included WCS affected animals (n = 164, cases) and WCS unaffected animals (n = 83, controls) that originated from farms in the Northern Cape, Free State, Mpumalanga, Namibia, Eastern Cape and Gauteng (Figure 2). Namibia is situated north west of South Africa, bordering the Northern Cape province of South Africa. However, animals from Namibia are transported to the Northern Cape for processing at the abattoir. The sampled population comprised of Dorper, unknown Crossbreds, Meatmaster, Van Rooyen, White Dorper and Dorper x Meatmaster, SA mutton Merino and Merino sheep. All animals were under 1 year old (classified as A-grade lambs) and were graded A1 or A2 for carcass quality immediately after slaughter. Additional demographic factors, such as weight and production system (extensive vs. intensive), were recorded when available. In addition to the animals genotyped in this study, genotypes for reference animals were obtained from the Grootfontein Agricultural Development Institute (GADI) and Dzomba et al. (2021). The reference population contained representatives of Afrino (n = 51), Blackhead Persian (BHP; n = 14), Dohne Merino (n = 50), Dorper (n = 22), Meatmaster (n = 48), Merino (n = 91), Namakwa-Afrikaner (n = 10), South African Merino (SAM; n = 10), and South African Mutton Merino (SAMM; n = 10). Sheep genome Oar_v4.0 was selected as the reference genome for this study.

For tissue lysis, 25 mg of meat sample was added to ATL buffer and Proteinase K, thereafter incubated at 56°C for 1 h. DNA was extracted from 200 μL of the resulting lysate using the Qiagen DNeasy Blood and Tissue extraction kit (Qiagen, Germany). Genomic DNA concentrations were measured using a Qubit^® 4.0 Fluorometer and standardized to 50 ng/μL (sample starting concentration between 10and 100 ng/μL). DNA samples were genotyped at the Agricultural Research Council - Biotechnology Platform, using the Ovine 50K SNP BeadChip (Illumina Inc., San Diego, CA, United States) that features 54,000 SNPs that are approximately equally distributed throughout the domestic sheep (Ovis aries) genome (Illumina, 2015; https://www.illumina.com). PLINK v1.9 (Purcell et al., 2007) was used for individual and SNP quality control (QC) for the genotyped and reference data. Quality control criteria included the following: SNP markers were excluded if minor allele frequency was less than 0.01 (--maf) and missing rate per SNP was greater than 0.05 threshold.

To investigate the genome-wide relationships and divergences among the sheep populations, population structure analysis was conducted. Nei’s genetic distances were used to construct a neighbour-joining tree to visualize relationship between the study population and reference groups and plotted using R package “ggplot2” v3.5.1 in R v4.0.2. Fixation index (Fst) values between each pair of breeds with 95% level of confidence were calculated in R using the package “StAMPP” v 1.6.3 (Pembleton et al., 2013). Principal component analysis (PCA) was performed to describe continuous ancestral heterogeneity within the study population (Greenacre et al., 2022). The Eigen values were generated using the “pca” parameter in PLINK v 1.9, and a PCA plot generated using the packages “tidyverse” v 2.0.0 (Wickham et al., 2019) and “ggplot2” v3.5.1 in R v4.0.3 (R Core Team, 2023). Thereafter, a cluster analysis using the k-means method was performed using “cluster” v2.1.8.1 in R 4.5.0, where both the Elbow method and Silhouette analysis were used to determine the optimal K for the number of clusters for this study population.

To assess the relationships between province, WCS status, breed, and clustering patterns, statistical analyses were performed using chi-square tests, Cramér’s V, Fisher’s Exact Test, and multinomial logistic regression. Contingency tables were constructed for each categorical variable pair, and chi-square tests were conducted to determine whether observed distributions deviated significantly from expected values. The strength of associations was quantified using Cramér’s V, with values interpreted as weak (<0.2), moderate (0.2–0.4), or strong (>0.4) (Kearney and Cramer’s, 2017). Fisher’s Exact Test was employed as an alternative for cases with low expected counts, using 10,000 Monte Carlo simulations to ensure robust p-value estimation. Additionally, multinomial logistic regression models were fitted to evaluate the predictive relationships between province, WCS status and breed with clustering. Model fit was assessed using residual deviance and Akaike Information Criterion (AIC), with lower AIC values indicating better model performance. All statistical analyses were performed in R v 4.3.3 using the “stats” v 4.1.1 and “nnet” v 7.3–20 packages (Ripley et al., 2016) for chi-square tests and multinomial logistic regression, respectively. The “vcd” v1.4-13 package (Meyer et al., 2024) and ‘“rcompanion” v 2.5.0 (Mangiafico, 2025) was used for Cramér’s V calculations, and Fisher’s Exact Test was conducted using the ‘fisher.test’ function with simulated p-values.

Ancestry matrix analysis was employed to determine whether certain genetic admixtures are more prevalent in WCS-affected vs. unaffected animals. Individual ancestry coefficients were calculated using the sparse non-negative matrix factorization (sNMF) with the R package, LEA v 3.14.0 (Frichot and François, 2015) with K (i.e., number of hypothetical ancestors) values ranging from 1 to 12, with five iterations each. The optimal number of ancestors was set as the lowest cross-entropy value. The resulting Q values (the ancestry coefficients) were plotted breed-wise to show proportion of ancestry from different breeds in each individual animal (Patel et al., 2022) using the ‘plotQMultiline’ function in “Pophelper” v2.3.1 package (Francis, 2017) in R v 4.1.1. The K value is the number of sub-populations that make up the total population, and is determined by the total number of reference populations that is included in the analysis (Alexander et al., 2015).

Following quality control, 34,528 SNPs and 540 animals remained. This comprised 164 affected animals and 376 unaffected animals which included the reference populations. The majority of affected were collected in the Northern Cape (47%), Free State (19%) and Namibia (21%), while the majority of unaffected were collected in the Northern Cape (43%) and Namibia (34%) (Figures 3, 4). The majority of affected animals of which the breed status was known, were Dorper (73%) and Merino (19%) sheep. The remaining breeds (8%), included Dorper x Meatmaster crossbred, Meatmaster, unknown Crossbred, White Dorper and Van Rooyen. Among the unaffected animals of which the breed status was known, included Merino (30%), followed by Dorper (15%), Afrino (14%), Dohne Merino (14%), Meatmaster (13%). The remaining breeds 29%, included SA Merino, Namakwa-Afrikaner, Crossbred, Blackhead Persian, as well as SA Mutton Merino sheep.

To understand the relationship between the WCS affected and unaffected populations in relation to reference populations, a neighbor-joining tree was constructed (Figure 5). The WCS affected and unaffected populations are grouped closely together. Furthermore, it shows the Dohne Merino and Dorper populations are close to the WCS-affected and unaffected populations. Figure 6 shows the results of PCA analyses that maximally explain the variance of all the variables across all SNP loci. The first (21.12%) and second (18.72%) principal components accounted for 39.84% of the variability. PCA analyses revealed 3 major clusters which was validated using k-means clustering analysis. Furthermore, the elbow method (Figure 7) helps determine the optimal number of clusters (K) by plotting within-cluster sum of squares (WCSS) against the number of clusters, while the average silhouette score (Figure 8) validates the quality of the clusters, with higher scores indicating better separation (Dalmaijer et al., 2022). The point where the curve starts to flatten out (the “elbow”) suggests the optimal K which is K = 3 in the study population (Figure 7). In addition, an average silhouette score of 0.7–1.0 indicates strong evidence for clustering and the location of the maximum is considered as the appropriate number of clusters, which is K = 3 in this study population (Figure 8). The PCA also showed the affected samples to reside in two of the clusters, supporting the interpretation that the origin of WCS is not breed-specific (Figure 9). With the exception of 2 samples, the affected and unaffected populations grouped in Clusters 1–3 (Figure 10). Since the affected animals were observed to reside in two major Clusters (2 and 3), suggested that additional factors beyond genetic background contribute to WCS occurrence. Even though the affected and unaffected populations grouped in Clusters 1, 2 and 3, the majority of affected animals were located in Cluster 2, while the majority of unaffected populations were located in Cluster 1. While the Merino were present in all 3 clusters, Dorper was only present in 2 clusters. It is interesting to note that the majority of affected Merino were observed in Cluster 1, while the unaffected Merino resided in Clusters 1 and 3. Majority of affected Dorper grouped in Cluster 3. Geographical location of the animals in each cluster reveal that each cluster has the largest proportion of animals from the Northern Cape province, with the highest present in cluster 1 (Figure 11). The second largest proportion for cluster 1 and 3 animals originating from Namibia while cluster 2 has animals from the Free State province. These clusters may reflect different production environments, regional variations, or management practices at the sampled abattoirs.

Ancestry-based clustering uncovered the diversity and breed relations of the populations under study. The cross-entropy criterion assists in selecting the number of ancestral populations or a best run for a fixed value of K. A smaller value of cross-entropy means a better run in terms of prediction capability (Frichot and François, 2015). The cross-entropy criterion that is based on the prediction of masked genotypes to evaluate the fit of a model with K population indicated that the sheep in the present study descended from 6 populations (Figure 12). The population structure revealed unique genomic signatures in the WCS affected and unaffected populations with increasing K values. The ancestry at the population level is quite diverse for most animals of a particular breed, like the merino-types, hence clear admixture patterns are difficult to distinguish in a single breed. One possible reason could be due the origin ancestry of the Merino breed, as intensive gene flow, founder effects and geographic isolation are the main factors that determined the genetic makeup of current Merino and Merino-derived breeds (Ciani et al., 2015). Therefore, the admixture analysis was run per individual and then sorted by ancestry probabilities, placing individual WCS affected and unaffected samples next to a reference breed. K = 6 was selected as it corresponds to the inflection point of the cross-entropy curve. The value of K = 6 corresponds closely to geographical distributions within the sample, that align with breeding history and population management practices in South African sheep. Thus, at the individual level of K = 6, the ancestry matrix shows the clustering of both affected and unaffected individuals with the individuals of the various reference breeds (Figure 13; Supplementary Table S1). The LEA program automatically assigned color codes to each individual based on K-value and ancestry. Most WCS individuals were clustered with the Merino-type and Dorper animals confirming the results observed in the neighbor joining unrooted tree (Figure 5). Also, the WCS affected displayed admixed signatures from all the Merino-type populations. The lack of a single ancestral population being significantly represented in an overwhelming majority of affected animals is viewed as evidence that WCS does not originate from one particular breed.

Population pairwise FST ranged from 0.008 (between the WCS affected sheep and Dohne Merino) to 0.036 (between WCS affected and Blackhead Persian sheep populations) (Table 1). This observation is also seen in Figure 4, where the WCS-affected population is genetically closer to the Dohne Merino than the Blackhead Persian population. The neighbour-joining tree (Figure 5) illustrates genetic differentiation among the populations that were studied. In contrast, estimates of FST account for genetic sampling and reflect differences in the allele frequency distributions. While Nei’s distance quantifies the overall genetic difference between populations and is often used to estimate divergence time. As a result, FST and Nei’s genetic distance measure different properties of the populations under study (Holsinger and Weir, 2009).

The statistical analysis revealed significant associations between province, status, and breed with clustering patterns (Tables 2 and 3). The chi-square tests indicated highly significant relationships across all comparisons (p < 0.0001), confirming that the observed distributions were not random (Table 2). The strength of these associations varied, with Cramér’s V values ranging from 0.35 to 0.36 based on province while, for breed it ranges from 0.25 to 0.60. Notably, breed exhibited the strongest association with clustering (Cramér’s V = 0.36), followed by WCS status (Cramér’s V = 0.25), suggesting that both genetic and management-related factors influence clustering (Table 3). Province also demonstrated a moderate effect (Cramér’s V = 0.32), indicating a geographic component in the clustering patterns (Table 2). The relationship between WCS status and breed was particularly pronounced, with a chi-square statistic of 166.92 (df = 11, p < 2.2e-16) and a strong effect size (Cramér’s V = 0.60), suggesting that breed composition varies significantly between the WCS affected group and the unaffected group (Table 2). Similarly, breed was significantly associated with clustering (X² = 126.79, df = 22, p < 2.2e-16), reinforcing the idea that genetic or phenotypic factors contribute to cluster formation. The Fisher’s Exact Test, performed with 10,000 simulations, consistently confirmed these associations, providing additional robustness to the statistical findings (Tables 2 and 3). Multinomial logistic regression models further supported these associations, with relatively low AIC values across comparisons, indicating a good model fit. The model for breed vs. cluster showed the lowest AIC (327.37), followed by WCS status vs. province (428.90), suggesting that these variables play a key role in defining clustering patterns. Collectively, these results highlight the significant influence of both geographic (Province) and biological (WCS status and breed) factors in shaping the observed clusters (Tables 2 and 3).