Whole-Genome Sequencing Reveals Breed-Specific SNPs, Indels, and Signatures of Selection in Royal White and White Dorper Sheep

1. Introduction

The domestic sheep (Ovis aries) is a globally important livestock species, contributing to food, fiber, and income security through the production of meat, wool, milk, and hides. According to statistical summaries from the Food and Agriculture Organization (FAO) of the United Nations for 2024, the global sheep population reached approximately 1.5 billion head worldwide [1]. In the U.S., the national flock totaled 5 million sheep and lambs as of January 1, 2025, including 3.7 million breeding sheep and 1.4 million market lambs [2]. Among U.S. sheep breeds, Royal White (RW) and White Dorper (WD) are prominent meat-type hair sheep valued for efficient meat production and adaptability, and they are also increasing in popularity in the U.S. due to their climate adaptability and lack of shearing requirements. Royal White^® Sheep is a U.S.-developed composite breed created by Bill Hoag in the early 2000s, through the crossbreeding of Dorper and St. Croix sheep. The breed was developed to combine desirable traits such as carcass quality, parasite resistance, a clean-shedding hair coat, and adaptability to diverse production environments [3]. White Dorper is a South African meat-type hair sheep developed through strategic crossbreeding beginning in the 1930s. The breed originated from Dorset Horn rams imported from Australia and crossed with Blackhead Persian ewes, with later contributions from Van Rooy sheep. White Dorper shares identical breed standards with Blackhead Dorper, differing only in coat color and pigmentation [4]. Today, the breed is valued for its rapid growth, high fertility, adaptability, and broad use in arid production systems [5].

Given the RW and WD commercial importance, understanding the genomic architecture of these breeds is essential. Despite the growing economic relevance of U.S. meat production systems, genomic studies on RW and WD sheep remain limited. To date, no published studies have comprehensively characterized RW sheep and U.S. populations of WD sheep at the whole-genome level, leaving a critical gap in our understanding of their genetic architecture and selection history. While WD has been evaluated genomically in regions such as South Africa and Hungary using SNP chips [6,7], few genomic studies have focused on U.S. populations, limiting our understanding of how this breed adapts and performs under American production conditions. This gap is particularly relevant given that environmental pressures and selection objectives may differ across geographic regions, potentially shaping distinct genomic signatures. To address this gap, whole-genome sequencing (WGS) provides a powerful approach for capturing genome-wide variation, offering deeper insights into the genetic basis of breed adaptation, performance, and selection under specific production environments.

Whole-genome sequencing has substantially advanced livestock genomics by enabling comprehensive characterization of genetic variation across the entire genome. Unlike marker-based approaches, WGS captures both common and rare genetic variants at single-base resolution, including single nucleotide polymorphisms (SNPs) and small insertions and deletions (indels). These variants provide valuable insights into genetic diversity, population structure, breed evolution, and the genetic basis of traits under selection [8,9]. In livestock research, WGS facilitates the identification of putatively functional mutations and signatures of selection, offering a powerful framework for improving animal breeding, conservation, and adaptation strategies. Applying WGS to under-characterized breeds such as RW and U.S.-based WD sheep enables the discovery of breed-specific genomic features that may contribute to performance and environmental resilience. By applying WGS to both RW and WD sheep from the same flock, this study aims to uncover breed-specific genomic features and provide foundational insights for future breeding and conservation efforts.

The objectives of this study were to (1) identify and compare genome-wide genetic variants (including SNPs and indels) in RW and WD sheep using WGS, (2) annotate the functional impact of these variants and explore enriched biological pathways through Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses, and (3) detect regions of the genome putatively under recent selection in each breed and identify candidate genes associated with key traits.

2. Materials and Methods

2.1. Sample Collection, DNA Extraction, Library Preparation, and Sequencing

All animal samples were obtained and experiments conducted according to the Virginia Polytechnic Institute and State University Institutional Animal Care and Use Committee (IACUC) #17-233 and Institutional Biosafety Committee (IBC) #20-067. A total of 20 ewe sheep, RW (Royal White; $n=11$) and WD (White Dorper; $n=9$), were used in this study. All animals originated from the same privately owned flock in the Southern U.S., minimizing environmental variation. The average age was $2.7\pm 0.3$ (mean ± SE) years for RW and $2.89\pm 0.35$ (mean ± SE) years for WD. For this study, a single blood sample was collected from each of the 20 animals. Whole blood (approximately 4 mL) was collected via jugular venipuncture into Ethylenediaminetetraacetic acid (EDTA) tubes under approved IACUC protocols. Samples were shipped overnight on ice, processed the next day, and peripheral blood mononuclear cells (PBMC) were isolated using density gradient centrifugation. Isolated PBMCs were washed, cryopreserved in CryoStor CS10, stored at –80°C for 24 h, and then transferred to liquid nitrogen until DNA extraction. For DNA extraction, up to $5\times {10}^6$ thawed PBMCs were centrifuged at 300 $\times g$ for 6 minutes to obtain a cell pellet, which was resuspended in 200 µL phosphate-buffered saline (PBS) and mixed with 20 µL proteinase K. After the addition of 200 µL buffer AL (without ethanol; Qiagen; lysis buffer), samples were incubated at 56 °C for 10 minutes, followed by the addition of 200 µL of 100% ethanol. The mixture was then transferred to a spin column, washed with wash buffers AW1 (Qiagen, first wash buffer) and AW2 (Qiagen, second wash buffer), and DNA was eluted in 100 µL buffer AE (Qiagen, elution buffer). The concentration and purity of DNA were assessed using a NanoDrop ND-100 spectrophotometer. Libraries were prepared using the TruSeq DNA PCR-Free Library Preparation Kit (Illumina) following the manufacturer’s instructions. Libraries were sequenced on the Illumina NovaSeq 6000 platform with a paired-end 151 bp read length using the NovaSeq FASTQ workflow.

2.2. Alignment and Variant Calling

Raw reads were processed using Trimmomatic (version 0.39) [10], where adapter sequences were removed, low-quality bases (quality score < 30) were trimmed using a sliding window approach, and reads shorter than 50 bp were discarded. High-quality paired-end reads were aligned to the Ovis aries reference genome (ARS-UI_Ramb_v2.0) using BWA-MEM (version 0.7.18) [11] with default parameters. Aligned reads were first processed with GATK (version 4.6.1) [12] MarkDuplicates to identify and mark PCR duplicates. The resulting deduplicated BAM files were used for variant calling with GATK HaplotypeCaller in GVCF mode. Variants were initially called on a per-sample basis, and individual GVCF files were then combined by breed, RW and WD, using GATK CombineGVCFs. Joint genotyping was performed within each breed using GATK GenotypeGVCFs to produce breed-specific multi-sample VCF files. To obtain high-confidence variants, breed-specific VCF files were filtered using BCFtools (version 1.21) [13]. Variants were retained if they passed the following criteria: QUAL > 30, depth (DP) > 8, Fisher strand bias (FS) < 60.0, mapping quality (MQ) > 40.0, mapping quality rank sum (MQRankSum) >−12.5, quality by depth (QD) ≥ 2.0, read position rank sum (ReadPosRankSum) >−8.0, and strand odds ratio (SOR) ≤ 3.0. Filtered variants were compared against the Ensembl variation database (Ovis_aries_rambouillet; version 113; file date: 2024-08-29) to distinguish known and novel variants. The reference database contained a total of 89,897,984 variants, including 83,083,034 SNPs and 6,763,747 indels. Venn diagrams illustrating shared and unique variants between breeds were generated using the VennDiagram R package (version 1.7.3) [14]. The indel length distribution was visualized using a histogram generated with ggplot2 (version 3.5.1).

2.3. Functional Annotation and Gene Enrichment Analysis

Functional annotation of SNP and indel variant files was performed using SnpEff (version 5.2c) [15] with the Oar_ARS_UI_Ramb_v2_0 database. Annotated VCF files were generated for each breed and variant type. Following annotation, variants were filtered to retain those with predicted functional effects classified as HIGH or MODERATE impact by SnpEff. Genes affected by these variants were extracted from the annotated VCF files. To focus on genes with stronger predicted functional burden, only genes harboring at least five HIGH or MODERATE impact variants were retained for downstream analysis. The resulting gene lists were submitted to the Database for Annotation, Visualization and Integrated Discovery (DAVID, version 6.8) [16] for functional enrichment analysis. Gene ontology enrichment analysis focused on the biological process, cellular component, and molecular function categories, as well as KEGG pathways. The top 10 enriched terms from each GO category and from KEGG pathways, ranked by false discovery rate (FDR), were visualized in R using the ggplot2 package (version 3.5.1) based on their minus log10(FDR) values.

2.4. Population Structure Analysis

Prior to population structure analysis, variants were filtered to remove those with a call rate below 80% and a minor allele frequency (MAF) less than 0.05. Missing genotype values were imputed using Beagle (version 5.4) [17]. The filtered and imputed variant datasets for both RW and WD breeds were converted to PLINK binary format using PLINK2 (version 2.0) [18]. Population structure was assessed using principal component analysis (PCA) based on the top ten principal components. Principal component analysis was performed in PLINK2 using genotype data with standard allele-frequency–based scaling to account for differences in allele frequencies between breeds. The first two components were visualized using ggplot2 (version 3.5.1) in R (version 4.3.0) [19], with breed-specific coloring and ellipses to illustrate group separation. In addition to PCA, discriminant analysis of principal components (DAPC) [20] was conducted using the adegenet (version 2.1.11) package in R to further investigate breed-specific genetic differentiation. The first discriminant function was visualized as a density plot using ggplot2, with coloring consistent with PCA plots to enhance visual comparability.

2.5. Selective Sweep Analysis

Selective sweep regions were identified by integrating the population differentiation fixation index (F_ST), nucleotide diversity ($\pi$), and Tajima’s D metrics. Window-based F_ST values were calculated using VCFtools (version 0.1.16) [21] with a sliding window of 50 kb and a step of 10 kb. Breed-specific nucleotide diversity ($\pi$) and Tajima’s D were computed separately for RW and WD using the same window parameters. The top 10% of genomic windows based on F_ST were selected as candidate regions of differentiation. Within these, regions with reduced diversity in one breed relative to the other were identified by calculating the natural log ratio of $\pi$ values [$ln({\pi }_{\mathrm{RW}}/{\pi }_{\mathrm{WD}})$] and selecting the top and bottom 10% as WD- and RW-specific sweep candidates, respectively. To strengthen the evidence for putative recent selective sweeps, an additional filter was applied based on Tajima’s D. For each breed, only candidate regions with negative breed-specific Tajima’s D values were retained, indicating an excess of low-frequency alleles consistent with recent positive selection. Genomic positions in these regions were annotated with gene information. Genotypic identifiers were retrieved by cross-referencing positions with the Ensembl Ovis aries Rambouillet variation database (version 113). Trait associations were then incorporated by mapping SNP IDs to the sheep quantitative trait locus database (QTLdb; release 55; file date: 2024-12-23), enabling the identification of candidate genes linked to economically important traits.

3. Results

3.1. Read Quality, Mapping, and Depth Coverage

Whole-genome sequencing generated approximately 843 million raw paired-end reads across 11 RW sheep and 625 million raw paired-end reads across 9 WD sheep. After trimming, an average of 94.98% of RW reads and 92.95% of WD reads were retained. The average mapping rate for both breeds was consistently high at 99.93%, indicating excellent alignment quality. Royal White samples showed an average sequencing depth of 7.19× (ranging from 5.34× to 10.00×), while WD samples had an average depth of 6.33× (ranging from 5.01× to 7.41×) (Table 1).

3.2. Variant Calling

Variants were identified through within-breed joint genotyping using GATK and filtered to retain high-confidence SNPs and indels. A total of 26,555,583 SNPs and 3,703,099 indels were initially identified in RW, and 24,792,950 SNPs and 3,396,380 indels in WD. After quality filtering, 21,957,139 SNPs and 2,866,600 indels were retained in RW, and 18,641,789 SNPs and 2,397,368 indels in WD. The transition-to-transversion (Ts/Tv) ratio was 2.30 in RW and 2.16 in WD. The heterozygous-to-homozygous (Het/Hom) genotype ratio across called variants was 0.999 (SNPs) and 0.998 (indels) in RW and 0.998 (SNPs) and 0.992 (indels) in WD. Annotation against the Ensembl Ovis aries variation database (release 113) revealed that 77.24% of SNPs and 21.83% of indels in RW were known, while 22.76% and 78.17%, respectively, were novel. Similarly, 77.57% of SNPs and 21.87% of indels in WD were known, with 22.43% and 78.13% being novel (Table 2). Venn diagram analysis showed that RW and WD shared 13,498,534 SNPs and 1,350,346 indels, while 8,458,605 SNPs and 1,516,254 indels were unique to RW, and 5,143,255 SNPs and 1,047,022 indels were unique to WD (Figure 1). The length distribution of indels was examined to characterize insertion and deletion patterns in both breeds. As shown in Figure 2, the majority of indels in RW and WD sheep were short, with sizes concentrated between $-5$ bp (deletions) and $+5$ bp (insertions). One base-pair indels were the most frequent in both breeds, followed by two base-pair changes. The overall distribution was symmetric around 0 bp, with a peak at $+1$ bp for insertions and $-1$ bp for deletions. These results indicate that short indels are predominant in both populations.

3.2.1. Functional Annotation and Gene Enrichment Analysis

A comprehensive functional annotation of SNP and indel variant files was performed separately for RW and WD sheep using SnpEff to evaluate the potential biological effects of genomic variation within each breed. For SNPs, both breeds exhibited a high proportion of variants in non-coding regions, including introns (27.47M in RW; 23.39M in WD), intergenic regions (11.97M in RW; 10.14M in WD), downstream gene regions (2.33M in RW; 1.98M in WD), and upstream gene regions (2.31M in RW; 1.95M in WD). Functionally important categories such as missense variants were also prevalent (142,906 in RW; 124,264 in WD), suggesting protein-altering potential (Table 3). Additional variants were found in synonymous, splice region, and UTR regions. For indels, the most abundant categories included intron variants (4.29M in RW; 3.53M in WD) and intergenic variants (1.74M in RW; 1.42M in WD), followed by downstream and upstream gene variants. High-impact functional classes such as frameshift variants (10,427 in RW; 9,510 in WD), disruptive in-frame insertions/deletions, splice site variants, and stop-gainedd variants were also detected (Table 4). These findings indicate that both breeds harbor a substantial number of variants with the potential to affect gene function and regulation.

Gene ontology and KEGG enrichment analyses of genes harboring HIGH and MODERATE impact variants revealed broadly consistent functional categories and pathways in both RW and WD sheep (Figure 3 and Figure 4). In KEGG, both breeds were strongly enriched for ABC transporters (oas02010), ECM–receptor interaction (oas04512), graft-versus-host disease (oas05332), complement and coagulation cascades (oas04610), cytoskeleton in muscle cells (oas04820), amoebiasis (oas05146), and the Fanconi anemia pathway (oas03460), indicating potential roles in transmembrane transport, cell–matrix signaling, immune regulation, structural integrity, host–pathogen interactions, and DNA damage repair. Royal White-specific enrichment included retinol metabolism (oas00830) and linoleic acid metabolism (oas00591), suggesting potential breed-specific adaptations in vitamin A utilization and fatty acid processing that may influence growth and health. White Dorper-specific enrichment included Staphylococcus aureus infection (oas05150) and motor proteins (oas04814), pointing to putative differences in pathogen defense mechanisms, intracellular transport, and cytoskeletal regulation.

Table 4. Summary of indel functional annotation categories in Royal White and White Dorper breeds based on SnpEff.

Variant Type	Royal White Indels	White Dorper Indels
3’ UTR truncation	1	1
3’ UTR variant	67,521	55,407
5’ UTR truncation	1	3
5’ UTR variant	21,424	17,480
Bidirectional gene fusion	78	68
Conservative inframe deletion	820	617
Conservative inframe insertion	654	479
Disruptive inframe deletion	1,473	1,098
Disruptive inframe insertion	671	509
Downstream gene variant	386,904	318,194
Exon loss variant	5	7
Frameshift variant	10,427	9,510
Gene fusion	59	30
Intergenic region	1,744,317	1,418,580
Intragenic variant	920	793
Intron variant	4,292,930	3,529,191
Non-coding transcript exon variant	22,460	19,116
Non-coding transcript variant	304	235
Splice acceptor variant	1,339	1,106
Splice donor variant	977	759
Splice region variant	12,640	10,121
Start lost	86	54
Start retained variant	24	17
Stop gained	135	127
Stop lost	64	56
Stop retained variant	10	10
Transcript ablation	1	2
Upstream gene variant	380,196	309,740

Table 4. Summary of indel functional annotation categories in Royal White and White Dorper breeds based on SnpEff.

Variant Type	Royal White Indels	White Dorper Indels
3’ UTR truncation	1	1
3’ UTR variant	67,521	55,407
5’ UTR truncation	1	3
5’ UTR variant	21,424	17,480
Bidirectional gene fusion	78	68
Conservative inframe deletion	820	617
Conservative inframe insertion	654	479
Disruptive inframe deletion	1,473	1,098
Disruptive inframe insertion	671	509
Downstream gene variant	386,904	318,194
Exon loss variant	5	7
Frameshift variant	10,427	9,510
Gene fusion	59	30
Intergenic region	1,744,317	1,418,580
Intragenic variant	920	793
Intron variant	4,292,930	3,529,191
Non-coding transcript exon variant	22,460	19,116
Non-coding transcript variant	304	235
Splice acceptor variant	1,339	1,106
Splice donor variant	977	759
Splice region variant	12,640	10,121
Start lost	86	54
Start retained variant	24	17
Stop gained	135	127
Stop lost	64	56
Stop retained variant	10	10
Transcript ablation	1	2
Upstream gene variant	380,196	309,740

Figure 3. The top 10 enriched terms from each functional category, including Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, Gene Ontology (GO) molecular function, cellular component, and biological process, are shown, ranked by false discovery rate (FDR). Bars represent –log10(FDR) values. Enrichment analyses were based on genes carrying at least five SNPs or indels with HIGH or MODERATE impact.

Figure 3. The top 10 enriched terms from each functional category, including Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, Gene Ontology (GO) molecular function, cellular component, and biological process, are shown, ranked by false discovery rate (FDR). Bars represent –log10(FDR) values. Enrichment analyses were based on genes carrying at least five SNPs or indels with HIGH or MODERATE impact.

Figure 4. The top 10 enriched terms from each functional category, including Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, Gene Ontology (GO) molecular function, cellular component, and biological process, are shown, ranked by false discovery rate (FDR). Bars represent –log10(FDR) values. Enrichment analyses were based on genes carrying at least five SNPs or indels with HIGH or MODERATE impact.

Figure 4. The top 10 enriched terms from each functional category, including Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, Gene Ontology (GO) molecular function, cellular component, and biological process, are shown, ranked by false discovery rate (FDR). Bars represent –log10(FDR) values. Enrichment analyses were based on genes carrying at least five SNPs or indels with HIGH or MODERATE impact.

At the molecular function level, both breeds were significantly enriched for protein binding (GO:0005515), calcium ion binding (GO:0005509), ATP binding (GO:0005524), ATP hydrolysis activity (GO:0016887), ABC-type transporter activity (GO:0140359), lipid transporter activity (GO:0005319), microtubule binding (GO:0008017), and extracellular matrix structural constituent (GO:0005201 in WD; GO:0030020 in RW), suggesting potential common roles in energy metabolism, membrane transport, cytoskeletal interactions, and structural integrity. Royal White-specific terms included four-way junction helicase activity (GO:0009378), transmembrane signaling receptor activity (GO:0004888), and carbohydrate binding (GO:0030246), indicating possible breed-specific differences in DNA repair, signal transduction, and carbohydrate recognition. White Dorper-specific enrichment included cadherin binding (GO:0045296) and endopeptidase inhibitor activity (GO:0004866), pointing to potential variation in cell–cell adhesion and proteolysis regulation.

In the cellular component category, enrichment for collagen trimer (GO:0005581), myosin complex (GO:0016459), plasma membrane (GO:0005886), and extracellular matrix (GO:0031012) was observed in both breeds, reflecting structural and signaling roles commonly affected by coding variants. Notably, RW showed unique enrichment for the catenin complex (GO:0016342), microtubule organizing center (GO:0005815), and cytoskeleton (GO:0005856), suggesting enhanced roles in cytoskeletal organization and cell–cell adhesion. White Dorper showed additional terms like dynein complex (GO:0030286) and microtubule (GO:0005874), indicating potential differences in microtubule-based transport and structural dynamics.

Biological process terms were also partially shared. Both breeds showed enrichment for homophilic cell adhesion (GO:0007156), regulation of immune system process (GO:0002682), positive regulation of T cell activation (GO:0050870), and regulation of cytokine production (GO:0001817), highlighting immune-related and intercellular communication pathways. Royal White showed unique enrichment for immunoglobulin production (GO:0002381) and peptide antigen assembly (GO:0002503), suggesting possible stronger signals in adaptive immunity. In contrast, WD showed unique enrichment for neurodevelopmental and physiological processes, including axon guidance (GO:0007411), membrane depolarization during action potential (GO:0086010), and tissue development (GO:0009888), along with complement activation (GO:0006958), the latter indicating a potential innate immune component. Overall, functional enrichment analyses revealed shared impacts of HIGH and MODERATE impact variants in RW and WD sheep on immune function, cytoskeletal organization, and membrane-associated transport, alongside breed-specific putative signatures highlighting metabolic, immune-regulatory, and neurophysiological differences.

3.3. Population Structure Analysis

Population structure was assessed using both PCA and DAPC to explore genetic differentiation between breeds. As shown in Figure 5, the PCA plot (Panel A) displays the first two principal components, which explained 14.83% and 13.69% of the total genetic variation, respectively. While overlap between individuals from the two breeds is observed, breed-level clustering is evident, with RW and WD forming partially distinct groups along PC1. The 95% confidence ellipses further illustrate the separation pattern between breeds. To complement PCA, DAPC was performed to maximize between-breed variation. The density plot of the first discriminant function (Panel B) reveals a clearer distinction between RW and WD individuals, indicating that DAPC was able to enhance the separation observed in PCA. These results collectively suggest detectable but moderate genetic structuring between the two breeds.

3.4. Selective Sweep Analysis

Selective sweep regions were identified by integrating F_ST, nucleotide diversity ($\pi$), and Tajima’s D metrics. Genes within these regions were annotated, and existing SNPs overlapping with sheep QTL records were mapped to known trait associations. Both RW and WD exhibited putative signatures of selection related to parasite resistance, but the underlying genes differed between breeds. In RW, parasite resistance signals involved TGFB2, TOX2, and HERC6, while in WD, they were associated with LAMC1, COLGALT2, TRIM14, and EPHA5. In addition to parasite resistance, RW showed putative selection on behavioral genes (GRM5, MAGI2), metabolic disease susceptibility (ALDH5A1), and growth- and quality-related loci (JADE2, PARP8, NIN, NRXN1) linked to body size, meat composition, milk production, and fiber-related characteristics (Table 5, Table 6). White Dorper displayed additional putative selection on growth loci (PLXDC2, HYDIN), milk composition (TENM2, BUD23, SCN8A), reproduction (STPG3, DYNC2H1), and morphology (LCN8, NFKB1) (Table 6). Collectively, these patterns indicate that while both breeds show putative selection for economically important traits and parasite resistance, RW may have stronger signals related to adaptive immunity and wool traits, whereas WD shows a broader balance across growth, reproduction, and parasite resistance.

4. Discussion

This study aimed to characterize genomic variations, assess population structure, and identify putative signatures of selection in two sheep breeds, RW and WD, which have distinct origins but share similar meat production purposes. Leveraging WGS data, we comprehensively analyzed SNP and indel variants, their functional impacts, and genomic regions under putative selection. Our findings provide valuable insights into the genetic architecture differentiating these breeds, with particular attention to loci related to parasite resistance, growth, and reproductive traits that may contribute to subtle phenotypic differences despite their overall production similarity.

4.1. Genomic Variant Characteristics

In the current study, WGS of 11 RW and 9 WD sheep identified approximately 21.96 million and 18.64 million SNPs and 2.87 million and 2.40 million indels, respectively. The sequencing coverage depth in this dataset ranged from 5.01× to 10×. These variant counts fall within the expected range when compared to other sheep WGS efforts. For example, a multi-species study involving domestic and wild sheep populations, including 18 domestic sheep and multiple wild relatives, reported 125.98 million SNPs and 13.04 million indels in total. Per-breed SNP counts ranged from approximately 13 million in European mouflon (n = 3) to 53 million in Asiatic mouflon (n = 16), with indel counts ranging from about 3 million to 7 million per breed. These samples were sequenced at coverage depths between 12.2× and 36.9× [22]. Furthermore, the Ts/Tv ratios were calculated as 2.30 for RW and 2.16 for WD in the current study. These values are consistent with typical mammalian genomes, which often fall within the range of 2.0 to 2.5, and specifically align with ratios reported in other sheep population genomic studies [23,24]. The consistency of these ratios further supports the high quality and accuracy of our SNP calls. Differences in variant counts across studies may be influenced by several factors. One key factor is sequencing depth, which affects the sensitivity of variant detection. Another important factor is sample size, as larger cohorts are more likely to capture rare and population-specific variants. On a broader scale, a larger study of 297 Duolang sheep identified 43.97 million SNPs and 6.50 million indels at approximately 13.35× coverage [25]. This higher variant yield is expected, given the substantially larger sample size and deeper sequencing depth, which together increase the power to detect both common and rare variants. The comparison underscores how sequencing depth and cohort size can influence variant discovery, and supports the interpretation that the variant counts observed in RW and WD are appropriate for the study design and technical parameters. Importantly, these results confirm that the sequencing and variant calling pipeline captured sufficient polymorphic sites for downstream analyses, including population structure, functional annotation, and sweep detection. Such comprehensive variant catalogs are essential for understanding the genetic basis of trait differentiation and provide a foundation for breed-specific genomic selection strategies.

4.2. Functional Annotation and Enrichment

Functional annotation using SnpEff revealed that the vast majority of SNPs and indels were located in non-coding regions, such as introns and intergenic regions, consistent with findings in other complex genomes, including those of sheep. This pattern has been well-documented in livestock genomics, where over 95% of detected variants typically lie outside coding regions due to the large proportion of non-exonic DNA in mammalian genomes [26,27]. Although the majority of SNPs and indels were located in non-coding regions (e.g., introns and intergenic areas), a smaller but biologically relevant subset occurred in coding regions or at exon-intron boundaries. These included missense, frameshift, splice site, and stop gain/loss variants, all of which are predicted to affect protein structure or gene regulation and may contribute to phenotypic variation [28,29]. Variants with predicted functional consequences are especially important in livestock genomics because they often underlie key traits like body weight and health. For example, body weight has been associated with specific SNPs and QTL regions in Merino sheep [30], while milk production traits have been linked to high-impact variants in crossbred dairy sheep [31]. Similarly, a genome-wide association study (GWAS) in meat sheep revealed associations of production traits such as birth weight, weaning weight, scan weight, and fat and muscle depth, alongside health traits including footrot and mastitis, demonstrating the polygenic and multifaceted nature of livestock traits [32]. The identification of predicted HIGH and MODERATE impact variants in both RW and WD sheep highlights changes primarily affecting protein-coding regions through amino acid substitutions, premature stop codons, or splice site disruptions. These findings suggest that selective processes, whether natural or artificial, continue to shape breed-specific genomic landscapes. Similar observations have been made across livestock species, where selection often acts on coding or regulatory variants to promote adaptation and improve performance traits [28].

Functional enrichment analyses based on KEGG pathways and GO terms revealed both shared and breed-specific biological processes in RW and WD sheep, providing insights into the molecular mechanisms underlying immunity, metabolism, and other key physiological traits. Specifically, KEGG pathway analysis showed that in both RW and WD breeds, enrichment of ABC transporters (oas02010) likely reflects roles in substrate transport and immune function, consistent with studies showing ABC transporter involvement in antigen processing in cattle [33]. ECM–receptor interaction (oas04512) is central to tissue remodeling and cellular communication; it is notably enriched in the ovine mammary gland during lactation, where it regulates epithelial cell adhesion and remodeling [34]. The complement and coagulation cascades (oas04610) pathway is a central component of innate immunity and has been shown to mediate early defense responses in sheep against Haemonchus contortus [35]. Retinol metabolism (oas00830), uniquely enriched in RW, has been associated with parasite resistance in sheep. Specifically, it elevated retinol, related gene expression correlates with resistance to Echinococcus granulosus infection [36]. The WD-specific enrichment for Staphylococcus aureus infection (oas05150) may reflect genetic adaptations linked to immune defense against bacterial pathogens. Similar KEGG enrichment was reported in bovine mammary gland transcriptome analyses, where the Staphylococcus aureus infection pathway (bta05150) was significantly enriched among genes differentially expressed in cows with subclinical Staphylococcus aureus mastitis [37]. Complementary to the KEGG results, GO term analysis further revealed shared and unique functional categories in RW and WD sheep, providing an additional layer of insight into the biological significance of the identified variants.

At the molecular function level, both RW and WD showed strong and highly significant enrichment for ATP binding (GO:0005524), protein binding (GO:0005515), calcium ion binding (GO:0005509), and ABC-type transporter activity (GO:0140359). The term ATP binding (GO:0005524) has also been detected in genome-wide selection scans in Dorper and Hu sheep, implicating it in growth and metabolic regulation [38]. Both breeds’ enrichment for protein binding (GO:0005515) and ATP binding (GO:0005524) have been reported in Chinese indigenous sheep adapted to warm climates, where these terms are functionally linked to metabolic regulation and heat loss mechanisms [39]. Calcium ion binding (GO:0005509), another shared term, has been observed in goat muscle development, with Leizhou goat fetal muscle studies showing differentially expressed genes enriched for this function [40]. ABC-type transporter activity (GO:0140359) underscores roles in membrane transport and detoxification; ABC transporters are well-characterized in veterinary pharmacology and pathogen defense, such as in drug absorption and xenobiotic handling [41]. Royal White-specific enrichment included transmembrane signaling receptor activity (GO:0004888), which involves membrane-bound receptors that detect extracellular cues. In wild cervids, this GO term was significantly enriched in antler-related genomic regions, suggesting a role in regulating tissue growth and regeneration [42]. Another RW-unique term, carbohydrate binding (GO:0030246), has direct ovine evidence from structural studies of the secretory glycoprotein SPS-40, which demonstrated specific carbohydrate-binding properties and conformational switching upon binding chitin-like oligosaccharides [43]. In WD, cadherin binding (GO:0045296), a function crucial for cell–cell adhesion, was significantly enriched, and has been identified in epigenomic studies of tissue remodeling in other mammals, such as cattle rumen during weaning [44]. White Dorper also showed unique enrichment for endopeptidase inhibitor activity (GO:0004866), which may imply regulation of proteolysis, a function critical in tissue remodeling and inflammation across vertebrates. This term has been reported in marine male fish Cyprinodon variegatus, where its expression was significantly altered following environmental chemical exposure, suggesting its sensitivity to physiological and environmental stressors [45].

In the cellular component category, enrichment was observed for structural elements such as collagen trimer (GO:0005581), myosin complex (GO:0016459), plasma membrane (GO:0005886), and extracellular matrix (GO:0031012) in both RW and WD sheep. Similar GO terms have been reported in transcriptomic studies of goats and sheep, where structural components are commonly enriched in tissues under selective pressure, including muscle and skin [46,47]. The enrichment of these cellular structures reflects a shared influence on genes involved in tissue organization, cellular architecture, and interactions between cells and their environment. Functional studies in livestock have shown that collagen and myosin-related components are essential for muscle fiber formation, extracellular support, and mechanotransduction processes that contribute to animal growth and performance [48]. Royal White sheep exhibited unique enrichment for the catenin complex (GO:0016342), microtubule organizing center (GO:0005815), and cytoskeleton (GO:0005856). The catenin complex, which is integral to adherens junctions and cell–cell adhesion, plays a key role in tissue integrity and differentiation. In transgenic mice expressing ovine $\mathsf{\beta }$-catenin, increased hair follicle density was observed, underscoring its potential impact on sheep morphological traits [49]. White Dorper sheep exhibited unique enrichment for the microtubule (GO:0005874). This cellular component has been directly linked to reproductive tissue function in avian livestock. In Sichuan white geese, microtubule (GO:0005874) enrichment was identified in ovarian tissue through integrated DNA methylation and transcriptomic analyses, implicating microtubule-associated genes such as EML6 in follicle growth and development [50].

Biological process terms highlighted both common and breed-specific functional enrichments in RW and WD sheep. Among the shared terms, homophilic cell adhesion (GO:0007156) was enriched in both breeds and has been reported among the top biological process GO terms in protoscoleces from sheep liver cystic echinococcosis cysts, accounting for 18% of annotated genes and involving plasma membrane adhesion molecules critical for cell–cell recognition during host–parasite interactions [51]. The enrichment of regulation of cytokine production (GO:0001817) and regulation of immune system process (GO:0002682) in both breeds is consistent with an ovine PBMC transcriptome study, where GO:0001817 was enriched post-vaccination and GO:0002682 was enriched post-adjuvant treatment, underscoring their central roles in immune activation [52]. Royal White-specific enrichment for immunoglobulin production (GO:0002381) and peptide antigen assembly (GO:0002503) points to emphasis on adaptive immunity. This is supported by deep sequencing of ovine immunoglobulin repertoires, which revealed extensive CDR3 diversity and active somatic hypermutation, indicating the capacity for enhanced and efficient humoral immune responses against infections [53]. In WD, unique enrichment was observed for primarily neurodevelopmental and excitability pathways, including axon guidance (GO:0007411) [54] and membrane depolarization during action potential (GO:0086010), processes also identified in livestock selection studies such as runs of homozygosity-based analyses in indigenous rabbit breeds [55]. While these terms are not classically immune-related, WD also exhibited unique enrichment for complement activation (GO:0006958), a core innate immune process. Complement activation has been reported as a key component of resistance to Haemonchus contortus infection in parasite-resistant sheep breeds (e.g., Canaria Hair Breed), linking this WD-specific term to protective immune functions [56].

4.3. Candidate Genes Under Selection

Selective sweep regions were identified by integrating population differentiation (F_ST), nucleotide diversity ($\pi$), and Tajima’s D metrics, a commonly used integrative approach for detecting recent positive selection [57,58,59].

4.3.1. Candidate Genes in Royal White Sheep

In RW sheep, selective sweep analysis revealed putative candidate genes associated with a range of traits, including health, behavior, growth, meat quality, and milk production. Many of these genes overlapped with known QTLs and have functional support from studies in sheep or other species, suggesting potential roles in breed-specific adaptation and productivity.

Health traits: Genes in this category include NRXN1, HERC6, TGFB2, TOX2, and ALDH5A1. NRXN1 was identified in two differentiated regions, associated respectively with red blood cell distribution width (associated SNP: rs409057468) [60] and fiber diameter coefficient of variation (associated SNP: rs429232758) [61]. These traits are consistent with findings from QTL studies in Alpine Merino and fine-wool sheep breeds, suggesting pleiotropic effects on health and fleece quality. HERC6, linked to fecal egg count, has been implicated in the host response to parasitic infection in an Australian sheep population [62]. In addition, HERC6 has been associated with milk production, growth, and feed efficiency in various livestock populations [63,64]. TGFB2, located in a region containing two different SNPs (rs162057314 and rs160759291) on chromosome 12, was associated with resistance to gastrointestinal nematodes (Haemonchus contortus), supported by multiple studies of resistance loci in sheep and goats [65,66]. The SNP rs423531735, associated with immune regulation [67], maps to the gene TOX2. While its specific function in sheep immunity has yet to be explored, studies in mice and humans indicate that TOX2 is integral to germinal center T follicular helper (GC TFH) cell formation and memory responses [68], suggesting a cross-species role in adaptive immune function. The gene ALDH5A1 harbored a SNP (rs421181203) that was identified in sheep as significantly associated with susceptibility to Mycobacterium avium subsp. paratuberculosis and antibody titer levels, suggesting a potential role in immune defense [69]. Supporting evidence from dairy cattle indicates that ALDH5A1 expression is linked to antibody-mediated immune responses, as individuals with higher expression showed traits consistent with enhanced immunity [70]. Together, these findings suggest that ALDH5A1 may play a conserved role in pathogen resistance across ruminant species.

Behavior traits: Genes in this category include MAGI2 and GRM5. Notably, two SNPs (rs424244818 and rs424837012) located within the GRM5 gene and one SNP (rs429561404) within the MAGI2 gene overlapped with QTL associated with vocalization and locomotion responses, as identified in studies of social and handling reactivity in sheep [71]. Although experimental evidence for these two genes in sheep is limited, this positional overlap highlights their candidacy as behavior-related loci. In particular, GRM5, which encodes a glutamate metabotropic receptor, has been associated with movement patterns and grazing behavior in beef cattle [72,73], suggesting a conserved role in behavioral regulation across ruminants. For MAGI2, while functional studies in sheep are lacking, its reported association with feed efficiency in cattle implies potential relevance to broader physiological or behavioral traits [74].

Milk production trait: The gene NIN (Ninein) harbored SNP rs410734119, which overlapped a QTL associated with milk yield in sheep. Although this positional evidence suggests potential involvement in lactation traits, current functional studies across species have not linked NIN to milk production or mammary gland biology. NIN encodes a centrosomal protein involved in microtubule anchoring and epithelial cell organization, with well-characterized roles in neural development and cytoskeletal dynamics in humans and mice [75]. However, no direct evidence currently supports its role in lactation, either through gene expression profiling or functional assays. Further studies are needed to determine whether the observed association reflects a causal relationship, a regulatory linkage to nearby lactation-relevant genes, or an indirect positional effect.

Growth trait: The gene JADE2 overlapped a QTL linked to 6-month body weight, supported by genome-wide association studies in Baluchi sheep [76]. In Djallonké sheep, JADE2 was located within a copy number variation region (CNVR) hotspot associated with lipid metabolism traits, further suggesting its potential role in growth and energy regulation [77].

Meat trait: SNP rs416975775, which influences the omega-6 to omega-3 fatty acid ratio in sheep meat, is located within the gene PARP8 [78]. Although the direct involvement of PARP8 in meat traits in sheep remains unconfirmed, members of the same poly(ADP-ribose) polymerase gene family, such as PARP1, have been implicated in post-mortem muscle tenderization mechanisms [79], supporting a potential role for PARP8 in muscle-related phenotypes.

4.3.2. Candidate Genes in White Dorper Sheep

In WD sheep, selective sweep regions uncovered functionally relevant genes linked to growth, immunity, reproduction, and milk production traits.

Health traits: The genes TRIM14, COLGALT2, LAMC1, and EPHA5 were identified in sweep regions linked to immune-related traits. The gene TRIM14, containing SNP rs422296454, was identified within a selective sweep region in WD sheep and is associated with increased hematocrit levels during gastrointestinal nematode infection [80], which found the same SNP and the same gene in the current study. In humans, TRIM14 has been described as a regulator of innate immune signaling and a putative tumor suppressor, modulating interferon pathways in non-small cell lung cancer [81]. These findings indicate that TRIM14 is under selective pressure in sheep due to its immunological and potentially pleiotropic functions. The SNP rs402132699, located in COLGALT2, was identified in a Brazilian Morada Nova sheep study as being associated with hematocrit (packed cell volume after nematode challenge) and fecal egg count (fecal egg count after nematode challenge) [82]. Additionally, COLGALT2 was among several glycosyltransferases identified via GWAS as candidate loci for milk oligosaccharide synthesis in Holstein and Jersey cattle [83], suggesting its potential influence on the nutritional quality and functional properties of milk. Moreover, COLGALT2 has been shown to be overexpressed in human ovarian cancer, where it interacts with PLOD3, suggesting a role in collagen glycosylation and extracellular matrix organization [84]. Taken together, these cross-species findings indicate that selective pressure on COLGALT2 in sheep may reflect its conserved function in glycan metabolism, growth regulation, and tissue adaptation. The SNP rs430289425 (located in LAMC1) was identified as associated with resistance to Haemonchus contortus infection in sheep and goats [65]. In dairy cattle, a novel QTL was discovered near the LAMC1/2 locus (BTA16:63823597), which was associated with variation in teat width, suggesting a potential role in tissue organization and mammary gland morphology [85]. Similarly, SNP rs426828157 (located in the gene EPHA5) was identified as associated with low fecal egg count in sheep [86], indicating a potential role in parasite resistance. Although direct functional validation in immunity is limited from this work, previous studies have highlighted EPHA5 as a candidate gene for wool traits in Chinese Merino and Kirghiz sheep populations [87,88]. Additionally, in goats, EPHA5 has been associated with body length and implicated in insulin-mediated growth pathways [89], suggesting broader physiological functions across livestock species.

Milk production traits: Genes under selection included TENM2, BUD23, and SCN8A. SNP rs409487914 in TENM2 overlapped a QTL for milk fat yield in sheep. SNP rs430795622, associated with 180-day milk fat yield [31], was located in BUD23, a gene identified in the current study. Although no livestock-specific studies have directly linked BUD23 to milk traits, it encodes an 18S rRNA methyltransferase known to regulate mitochondrial function and lipid metabolism in mice and humans [90,91]. Given the high energy demands of milk synthesis, these roles suggest a plausible functional relevance of BUD23 to lactation performance in ruminants. The gene SCN8A, harboring SNP rs419496265, intersects with selective sweep signals and QTLs for milk fat percentage. Although its role in sheep lactation remains unconfirmed, livestock studies have demonstrated that SCN8A is expressed in spermatozoa, specifically localized to the flagellum and neck of mammalian sperm cells, and is associated with sperm motility traits in pigs and horses [92,93]. Given its involvement in cellular excitability and ion transport, SCN8A may contribute to broader physiological processes relevant to energy metabolism and secretory activity in ruminant tissues. SNP rs409487914 (located in TENM2) was associated with milk fat yield in sheep, although its specific role in ovine lactation remains unverified in any species, indicating a need for further research.

Growth traits: Selective sweep regions identified PLXDC2 and HYDIN as candidates linked to body weight. The SNP rs410323459, located in the gene HYDIN, was associated with 8-month body weight in Iranian sheep populations, suggesting its role in late-stage growth performance [94]. Similarly, SNP rs401963094, located in PLXDC2, was associated with body weight at 9 months in Lori-Bakhtiari sheep [95]. Additionally, PLXDC2 has been linked to reproductive traits in Holstein cattle [96], highlighting its broader developmental importance.

Meat traits: The gene ADD2, harboring SNP rs417859328, was located within a selective sweep region in WD sheep and was associated with dressing percentage [97]. While direct evidence linking ADD2 to meat traits in sheep is limited, its paralog ADD1 offers valuable insights. In beef cattle, multiple SNPs within ADD1 were significantly associated with growth traits and are potentially useful for marker-assisted selection in breeding programs [98]. Similarly, in pigs, polymorphisms in ADD1 were linked to meat quality differences between Meishan and other commercial breeds, suggesting that the adducin gene family plays a role in adipose and muscle development [99]. These findings support the hypothesis that adducins, including ADD2, may influence carcass-related traits in livestock through cytoskeletal regulation and tissue organization pathways.

Reproductive traits: Genes were represented by loci such as DYNC2H1 and STPG3. SNP rs413723884 (located in DYNC2H1) was identified in association with offspring number across four parities [100]. While direct functional studies of DYNC2H1 in livestock are unavailable, the gene is known in humans and mice to encode a cytoplasmic dynein heavy chain essential for retrograde intraflagellar transport in primary cilia, which is critical for Hedgehog and Wnt signaling pathways that regulate ovarian follicle development and reproductive tissue morphogenesis [101,102]. These conserved cellular functions support its potential involvement in sheep prolificacy, consistent with its positional mapping in reproductive trait QTL. The gene STPG3, harboring SNP rs430682724, was situated within a selective sweep region linked to litter size in sheep, overlapping QTL evidence for offspring number in global breeds [103], suggesting a role in prolificacy. It is also known to be abundantly expressed in the testes of both humans and mice, as identified in a CRISPR-based screening study targeting testis-enriched genes for contraceptive development [104]. Although knockout of STPG3 did not impair male fertility in mice, its high expression in reproductive tissues supports a potential role in gametogenesis or sperm function. Additionally, a related gene, STPG2, has been implicated in male infertility, specifically azoospermia, in a Taiwanese cohort study investigating microtubule-associated gene clusters [105]. These findings suggest that STPG3, while not essential for male fecundity in mice, may participate in conserved testis-specific pathways that are potentially relevant to sheep reproduction.

Morphological traits: Genes under selection included LCN8 and NFKB1. The LCN8 gene, identified in a selective sweep region and associated with horn number through SNP rs415039972 [103], presents an interesting case of potential pleiotropy or positional linkage. Although its QTL association relates to horn phenotype, functional studies primarily describe LCN8 as a member of the lipocalin family involved in male reproduction. In sheep, LCN8 is highly expressed in the caput epididymis, where spermatozoa begin maturation [106]. Similarly, in humans and other mammals, it is enriched in the corpus region of the epididymis, suggesting a conserved role in sperm development and epididymal function [107]. This apparent functional divergence may indicate a pleiotropic influence or a neighboring regulatory element that affects both horn development and reproductive traits, which warrants further functional validation. Morphological selection was further supported by NFKB1, identified in a selective sweep region containing SNP rs416625889, which has been associated with bone area in QTL mapping studies of Scottish Blackface lambs [108]. Although primarily known for its role in immune regulation, NFKB1 is a key transcription factor in the nuclear factor kappa B signaling pathway, which mediates inflammatory responses and cellular stress signaling. In sheep, NFKB1 is actively expressed in maternal inguinal lymph nodes during early pregnancy, indicating its involvement in immune modulation at the maternal–fetal interface [109]. Furthermore, a retrospective SNP analysis of host resistance and susceptibility to ovine Johne’s disease, caused by Mycobacterium avium subsp. paratuberculosis, identified significant variants near genes involved in immune-related pathways, including the nuclear factor kappa B and mitogen-activated protein kinase signaling pathways, underscoring their role in host defense mechanisms against infection [110]. Broader livestock studies also implicate NFKB1 as a key transcription factor in the regulation of immune and inflammatory responses, playing a major role in mastitis susceptibility in beef cattle [111]. These findings suggest that NFKB1 may influence multiple traits under selection, including immune function, growth, and tissue development, either directly or through pleiotropic effects.

4.4. Limitations and Future Directions

This study provides foundational information on the genomic architecture of RW and WD sheep; however, several limitations should be acknowledged. The relatively small sample size may have limited the detection of rare variants and reduced the power to identify subtle signals of selection. Candidate genes were prioritized based on the presence or overlap of SNPs, recorded in the Sheep QTL Database, that were located within identified selective sweep regions; however, functional validation of these genes is still lacking. Future research should include larger and more diverse populations across multiple breeds, integrate gene expression and functional assays, and broaden variant discovery to encompass structural variants and epigenetic modifications, thereby enabling a more comprehensive understanding of breed-specific adaptations.

5. Conclusions

In conclusion, this study presents the first WGS-based comparative analysis between RW and WD sheep, identifying over 40 million SNPs and 5 million indels across the two breeds. Through variant annotation, population structure analyses, and selective sweep detection, we revealed breed-specific genomic regions and candidate genes associated with traits such as health, reproduction, growth, and production. These findings enhance our understanding of genetic differentiation between hair sheep breeds and provide a valuable foundation for future research and genomic selection in sheep breeding programs. Furthermore, the variants identified in this work can be investigated in additional populations and incorporated into genome-wide association studies to further elucidate the genetic architecture of economically important traits, ultimately supporting selective breeding strategies aimed at improving production and health in these breeds.