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

TL;DR

This study uses whole-genome sequencing to compare Royal White and White Dorper sheep, identifying genetic differences linked to traits like growth, immunity, and reproduction.

Contribution

The first whole-genome sequencing-based comparison of Royal White and White Dorper sheep, revealing breed-specific genetic markers and selection signatures.

Findings

Royal White sheep showed genetic differences in growth and health-related regions.

White Dorper sheep exhibited genetic variations in immunity and reproduction-related regions.

Over 21 million SNPs and nearly 2.9 million indels were identified in Royal White sheep.

Abstract

Royal White and White Dorper sheep are two breeds raised for meat in the U.S., yet little is known about their genomes. In this study, we examined the genome-wide DNA variation to identify unique regions that define each breed. We found that some of these genetic changes are likely linked to economically important traits such as parasite resistance, body growth, reproduction, and meat quality. For example, Royal White sheep showed genetic differences in regions related to growth and health, while White Dorper sheep showed differences in regions related to immunity and reproduction. These findings improve our understanding of the biology of these breeds and may help farmers refine breeding practices in the future.

Whole-genome sequencing (WGS) is a powerful tool for uncovering genome-wide variation, identifying selection signatures, and guiding genetic improvement in livestock. Royal White (RW) and White Dorper (WD) sheep are economically important meat-type hair breeds in the U.S., yet their genomic architecture remains poorly characterized. In this study, WGS was performed on 20 ewes (n = 11 RW, n = 9 WD) to identify and annotate SNPs and small insertions and deletions (indels). Functional annotation, gene enrichment, population structure, and selective sweep analysis were also performed. Selective sweep analysis was conducted by integrating the fixation index (FST), nucleotide diversity (π), and Tajima’s D to identify candidate regions under putative recent positive selection. A total of 21,957,139 SNPs and 2,866,600 indels were identified in RW sheep, whereas 18,641,789 SNPs and 2,397,368 indels were identified in WD sheep. In RW sheep, candidate genes under selection were associated with health and parasite resistance (NRXN1, HERC6, TGFB2) and growth traits (JADE2). In WD sheep, selective sweep regions included genes linked to immune response and parasite resistance (TRIM14), body weight (PLXDC2), and reproduction (STPG3). These findings were supported by sheep-specific quantitative trait loci (QTL) annotations and previously reported SNP–trait associations. This study provides the first WGS-based genomic comparison between RW and WD sheep, establishing a foundation for future genetic improvement, including targeted selection for enhanced immune function, disease resistance, and other economically important traits in these breeds.

Figures

Full text

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 [1]. In the U.S., the national flock totaled 5 million sheep and lambs as of 1 January 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. 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 suitability for arid production systems [5].

Given the commercial importance of RW and WD, 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 or 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 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 sequencing 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 genomic regions 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 procedures were performed in accordance with guidelines approved by 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 ewes, comprising Royal White (RW; $[eqn]$ ) and White Dorper (WD; $[eqn]$ ) breeds, were used in this study. All animals originated from the same privately owned flock in the Southern U.S., minimizing environmental variation. Because all animals were sampled from the same flock, genome-wide relatedness was evaluated to quantify potential kinship among individuals. Pairwise identity-by-descent (IBD) estimates were calculated using PLINK v1.9 (–genome option) [10] following linkage disequilibrium (LD) pruning, and relatedness was assessed using the PI_HAT statistic. Among the 20 individuals (190 pairwise comparisons), four within-breed pairs showed PI_HAT ≥ 0.0884 (approximately third-degree relatives or closer), including one second-degree pair (maximum PI_HAT = 0.1846). No cross-breed relatedness was detected. Because relatedness was confined within breeds and the primary objective of this study was between-breed genomic comparison, all samples were retained for downstream analyses. The average age was $[eqn]$ (mean ± SE) years for RW and $[eqn]$ (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 and processed for genomic DNA extraction. Peripheral blood mononuclear cells (PBMCs) were isolated and used as the source material for DNA extraction. For DNA extraction, up to $[eqn]$ PBMCs were centrifuged at 300 $[eqn]$ for 6 min 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, Hilden, Germany; lysis buffer), samples were incubated at 56 °C for 10 min, followed by the addition of 200 µL of 100% ethanol. The mixture was then transferred to a spin column and washed with wash buffers AW1 (Qiagen, Hilden, Germany) and AW2 (Qiagen, Hilden, Germany), and DNA was eluted in 100 µL buffer AE (Qiagen, Hilden, Germany). The concentration and purity of DNA were assessed using a NanoDrop ND-100 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA). Libraries were prepared using the TruSeq DNA PCR-Free Library Preparation Kit (Illumina, Inc., San Diego, CA, USA) following the manufacturer’s instructions. Libraries were sequenced on the Illumina NovaSeq 6000 platform (Illumina, Inc., San Diego, CA, USA) 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) [11] to remove adapter sequences, 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) [12] with default parameters. Aligned reads were first processed with GATK (version 4.6.1) [13] using 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. Breed-specific joint genotyping was performed to avoid allele frequency bias between breeds and to improve genotype accuracy for breed-specific variants. To obtain high-confidence variants, breed-specific VCF files were filtered using BCFtools (version 1.21) [14]. 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: 29 August 2024) 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) [15]. 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) [16] 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 Ovis aries gene lists were submitted to the Database for Annotation, Visualization and Integrated Discovery (DAVID, version 6.8) [17] for functional enrichment analysis. Gene identifiers were analyzed under the Ovis aries species setting, and enrichment was evaluated against the sheep-specific background provided within DAVID. Gene Ontology enrichment analysis focused on the biological process, cellular component, and molecular function categories, as well as KEGG pathways. Statistical significance was determined using false discovery rate (FDR)-adjusted p-values, with an FDR threshold of 0.05 considered significant. The top 10 enriched terms from each GO category and from KEGG pathways, ranked by FDR, were visualized in R using the ggplot2 package (version 3.5.1) based on their $[eqn]$ . Enriched terms were interpreted in the context of conserved biological processes relevant to ovine physiology and production traits.

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) [18]. The filtered and imputed variant datasets for both RW and WD breeds were converted to PLINK binary format using PLINK v2.0 [10]. Population structure was assessed using principal component analysis (PCA) based on the top 10 principal components. Principal component analysis was performed in PLINK v2.0 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 ( $[eqn]$ ), and Tajima’s D metrics. These metrics capture complementary signatures of selection: F_ST_ detects between-breed differentiation, $[eqn]$ identifies within-breed diversity reductions, and Tajima’s D reflects allele frequency spectrum skews consistent with recent positive selection. Integrating multiple statistics reduces false positives inherent to any single metric. 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 ( $[eqn]$ ) 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 $[eqn]$ values $[eqn]$ 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. Variant 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: 23 December 2024), 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 from 11 RW sheep and 625 million from 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. Results of 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 for SNPs and 0.998 for indels in RW, and 0.998 for SNPs and 0.992 for indels in WD. Annotation against the Ensembl Ovis aries variation database (release 113) revealed that 77.24% of SNPs and 78.17% of indels in RW were known, while 22.76% of SNPs and 21.83% of indels were novel. Similarly, 77.57% of SNPs and 78.13% of indels in WD were known, with 22.43% of SNPs and 21.87% of indels 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 $[eqn]$ bp (deletions) and $[eqn]$ 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 $[eqn]$ bp for insertions and $[eqn]$ bp for deletions. These results indicate that short indels are predominant in both populations.

3.3. Results of 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.47 million in RW; 23.39 million in WD), intergenic regions (11.97 million in RW; 10.14 million in WD), downstream gene regions (2.33 million in RW; 1.98 million in WD), and upstream gene regions (2.31 million in RW; 1.95 million 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.29 million in RW; 3.53 million in WD) and intergenic variants (1.74 million in RW; 1.42 million in WD), followed by downstream and upstream gene variants. High-impact functional classes such as frameshift variants (10,427 in RW; 9510 in WD), disruptive in-frame insertions/deletions, splice site variants, and stop-gained 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 or 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), complement and coagulation cascades (oas04610), cytoskeleton in muscle cells (oas04820), and the Fanconi anemia pathway (oas03460), indicating potential roles in transmembrane transport, cell–matrix signaling, immune regulation, structural integrity, and DNA damage repair. Enrichment in graft-versus-host disease (oas05332) and amoebiasis (oas05146) pathways likely reflects conserved immune signaling and host defense gene components shared across mammalian species. 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.

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, representing related but distinct ECM structural roles), 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, indicating potential roles in structural organization, cell membrane integrity, and extracellular scaffolding. Notably, RW showed unique enrichment for cytoskeleton (GO:0005856), suggesting potential breed-specific differences in cellular structural architecture. White Dorper showed additional enrichment for dynein complex (GO:0030286) and microtubule (GO:0005874), indicating potential differences in microtubule-based intracellular transport and cytoskeletal dynamics.

Biological process terms were also partially shared. Both breeds showed enrichment for homophilic cell adhesion (GO:0007156), indicating potential roles in intercellular adhesion and cellular communication. Royal White showed unique enrichment for regulation of cytokine production (GO:0001817), highlighting genes associated with modulation of immune and inflammatory signaling pathways. White Dorper showed unique enrichment for axon guidance (GO:0007411), which may reflect conserved cell signaling and migration pathway components, as well as complement activation (GO:0006958), indicating enrichment of genes related to innate immune pathways. Overall, functional enrichment analyses revealed shared impacts of HIGH and MODERATE impact variants on immune-related processes, cytoskeletal organization, and membrane-associated components, alongside breed-specific enrichment patterns that may reflect differences in immune regulation and cellular signaling pathways.

3.4. Results of 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 (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 (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 low but detectable genetic structuring between the two breeds (genome-wide mean F_ST_ = 0.029).

3.5. Results of Selective Sweep Analysis

Selective sweep regions were identified by integrating F_ST_, nucleotide diversity ( $[eqn]$ ), 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 were identified in genomic regions harboring TGFB2, TOX2, and HERC6, while in WD, they were detected in regions containing LAMC1, COLGALT2, TRIM14, and EPHA5. In addition to parasite resistance, RW showed putative selection signals encompassing genes previously reported to be associated with behavioral traits (GRM5, MAGI2), metabolic disease susceptibility (ALDH5A1), and growth- and quality-related loci (JADE2, PARP8, NIN, NRXN1) that have been linked in prior studies to body size, meat composition, milk production, and fiber-related characteristics (Table 5 and Table 6). White Dorper displayed additional putative selection signals in genomic regions previously associated with growth (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 signals in regions related to economically important traits and parasite resistance, RW appears to show a concentration of putative selection signals in regions containing genes previously linked to adaptive immunity and wool traits, whereas WD shows signals distributed across regions associated with 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. It is important to note that Ts/Tv ratios primarily reflect variant call quality and underlying mutational patterns rather than population structure or relatedness. Thus, while these values support the technical reliability of variant detection, they do not directly inform kinship or demographic relationships among individuals. 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, the variant metrics and quality indicators 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, elevated retinol-related gene expression has been shown to correlate 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]. Enrichment for protein binding (GO:0005515) and ATP binding (GO:0005524) in both breeds has also 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]. White Dorper sheep exhibited unique enrichment for microtubule (GO:0005874). This cellular component has been directly linked to reproductive tissue function in avian livestock. While direct ovine evidence is currently limited, microtubule-associated genes have been implicated in follicle growth and reproductive tissue function in avian livestock, providing a broad comparative context for a potential role in cellular transport and structural dynamics [49].

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 [50]. The enrichment of regulation of immune system process (GO:0002682) in RW sheep is consistent with findings from an ovine PBMC transcriptome study, in which GO:0002682 was enriched following adjuvant treatment, suggesting that genes within this category are commonly involved in immune activation processes [51]. In WD, unique enrichment was observed for neurodevelopmental processes, including axon guidance (GO:0007411) [52]. Axon guidance plays a central role in neuronal network formation and has also been reported in livestock selection studies, such as runs-of-homozygosity-based analyses in indigenous rabbit breeds [53], suggesting that neural development pathways may be subject to selection in domesticated populations. 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 [54].

4.3. Candidate Genes Under Selection

Selective sweep regions were identified by integrating population differentiation (F_ST_), nucleotide diversity ( $[eqn]$ ), and Tajima’s D metrics, a commonly used integrative approach for detecting recent positive selection [55,56,57]. We acknowledge that Tajima’s D is sensitive to sample size and demographic history. In our study, Tajima’s D was applied as a supporting filter following F_ST_ and $[eqn]$ screening, reducing reliance on any single metric. Because both breeds originated from the same flock, shared demographic effects are partially controlled. Identified regions represent candidates for future functional validation in larger cohorts.

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; this cross-species evidence is consistent with a conserved role in adaptive immunity, though ovine-specific regulatory mechanisms remain to be characterized.

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) [58] and fiber diameter coefficient of variation (associated SNP: rs429232758) [59]. 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 [60]. In addition, HERC6 has been associated with milk production, growth, and feed efficiency in various livestock populations [61,62]. 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 [63,64]. The SNP rs423531735, associated with immune regulation [65], 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 [66], providing supportive comparative evidence for a potentially conserved role in adaptive immune function, although regulatory contexts may differ across species. 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 [67]. 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 [68]. Together, these cross-species observations suggest a potentially conserved function in ruminant pathogen resistance, though direct experimental evidence in sheep is currently lacking.

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 [69]. 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 [70,71]; this association offers indirect support for a role in ovine behavioral regulation, with the caveat that functional effects may not translate directly across species. For MAGI2, while functional studies in sheep are lacking, its reported association with feed efficiency in cattle [72] warrants further investigation into MAGI2’s potential physiological relevance in sheep.

Milk production traits: 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 [73]. 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 traits: The gene JADE2 overlapped a QTL linked to 6-month body weight, supported by genome-wide association studies in Baluchi sheep [74]. In Djallonké sheep, JADE2 was located within a copy number variation region (CNVR) hotspot associated with lipid metabolism traits; this positional overlap offers additional circumstantial evidence for a potential role in ovine growth and energy metabolism [75].

Meat traits: SNP rs416975775, which influences the omega-6 to omega-3 fatty acid ratio in sheep meat, is located within the gene PARP8 [76]. 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 [77]; this family-level evidence raises the possibility that PARP8 may participate in analogous muscle-related processes, though direct evidence in sheep remains absent.

4.3.2. Candidate Genes in White Dorper Sheep

In WD sheep, selective sweep regions overlapped genomic regions containing genes previously associated with 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 [78], 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 [79]. Taken together, these findings suggest that TRIM14 may represent a candidate gene within a region under selection in sheep, though the extent to which these immune functions are conserved in ovine biology has yet to be established. 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) [80]. Additionally, COLGALT2 was among several glycosyltransferases identified via GWAS as candidate loci for milk oligosaccharide synthesis in Holstein and Jersey cattle [81], 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 [82]. Taken together, these cross-species findings suggest that selective pressure in this genomic region may involve COLGALT2, though whether these roles extend to ovine connective tissue biology requires species-specific investigation. The SNP rs596561468 (located in LAMC1) was identified as associated with resistance to Haemonchus contortus infection in sheep and goats [63]. 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 [83]. Similarly, SNP rs426828157 (located in the gene EPHA5) was identified as associated with low fecal egg count in sheep [84], 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 [85,86]. Additionally, in goats, EPHA5 has been associated with body length and implicated in insulin-mediated growth pathways [87]; these associations across livestock species offer contextual support, though regulatory differences between goats and sheep should be considered when extrapolating these findings.

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 [88,89]. Given the high energy demands of milk synthesis, these roles suggest a plausible functional relevance of BUD23 to lactation performance in ruminants, though empirical evidence for this role in ovine lactation has yet to be reported. 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 [90,91]. Given its involvement in cellular excitability and ion transport, SCN8A has been associated in other species with physiological processes relevant to energy metabolism and secretory activity; these associations offer a basis for hypothesizing relevance in ruminant secretory tissues, pending direct investigation in sheep. 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, providing positional and comparative evidence for a potential association with late-stage growth performance [92]. Similarly, SNP rs401963094, located in PLXDC2, was associated with body weight at 9 months in Lori-Bakhtiari sheep [93]. Additionally, PLXDC2 has been linked to reproductive traits in Holstein cattle [94], 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 [95]. 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 [96]. 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 [97]. These findings suggest that adducins, including ADD2, may be involved in carcass-related traits in livestock, though direct evidence linking ADD2 to such phenotypes in sheep has not yet been reported.

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 [98]. 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 [99,100]. These conserved cellular functions suggest that DYNC2H1 may represent a candidate gene within a region associated with sheep prolificacy, though how these functions translate to ovine reproductive regulation at the mechanistic level remains an open question. 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 [101], 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 [102]. 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 [103]. These findings suggest that STPG3, while not essential for male fecundity in mice, may be associated with conserved testis-specific pathways potentially relevant to sheep reproduction, though whether analogous pathways operate in ovine biology requires experimental confirmation.

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 [101], 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 [104]. Similarly, in humans and other mammals, it is enriched in the corpus region of the epididymis, suggesting a potentially conserved association with sperm development and epididymal function across species [105]. 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 [106]. 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 [107]. 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 [108]. 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 [109]. These findings suggest that NFKB1 may represent a candidate gene within regions associated with immune function, growth, and tissue development, though the extent to which these associations reflect sheep-specific biology remains to be determined.

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 moderate sample size (RW, n = 11; WD, n = 9) may limit representation of the full breed-level genetic diversity and reduce the power to detect rare variants or subtle selection signals. Although similar sample sizes are commonly reported in livestock whole-genome resequencing studies due to sequencing cost and depth considerations, larger and more geographically diverse cohorts would improve the generalizability of these findings. 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. Additionally, several functional interpretations were informed by studies conducted in other species; although many gene functions are evolutionarily conserved, regulatory contexts and phenotypic effects may differ in sheep, and therefore these interspecies comparisons should be interpreted cautiously pending species-specific validation. 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 delivers significant genomic resources through comparative whole-genome analysis of RW and WD sheep. We define distinct breed-specific genomic architectures that directly support applied breeding objectives: immune-associated regions in RW sheep provide targets for selecting enhanced parasite resilience, reducing dependence on chemical parasite treatments; growth- and reproduction-associated signatures in WD sheep enable selection for improved production efficiency. These findings facilitate concrete breeding strategies, including marker-assisted selection for trait-specific improvement and genomic prediction models tailored to hair sheep breeds. The comprehensive variant catalog generated here serves as an immediate reference resource for designing SNP panels and implementing genome-wide association studies. While validation across multiple herds will confirm broad applicability, this work establishes a significant, actionable genomic framework to implement targeted breeding programs and inform conservation priorities for emerging hair sheep breeds in the U.S.