Using High-Density SNP Array to Reveal Selection Signatures Related to Prolificacy in Chinese and Kazakhstan Sheep Breeds

Simple Summary Genetic improvement of litter size trait in domestic animals is an appealing way to improve production efficiency. In our study, the selection signatures between multiparous and uniparous sheep populations are identified, so that potential pathways and candidate genes related to litter size were screened out. Our findings help better understand the mechanisms of selection underlying the prolificacy trait in sheep and other mammals. Abstract Selection signature provides an efficient tool to explore genes related to traits of interest. In this study, 176 ewes from one Chinese uniparous breed and three Kazakhstan multiparous breeds are genotyped using Affymetrix 600K HD single nucleotide polymorphism (SNP) arrays, F-statistics (Fst), and a Cross Population Extend Haplotype Homozygosity Test (XPEHH). These are conducted to identify genomic regions that might be under selection in three population pairs comprised the one multiparous breed and the uniparous breed. A total of 177 and 3072 common selective signatures were identified by Fst and XPEHH test, respectively. Nearly half of the common signatures detected by Fst were also captured by XPEHH test. In addition, 1337 positive and 1735 common negative signatures were observed by XPEHH in three Kazakhstan multiparous breeds. In total, 242 and 798 genes were identified in selective regions and positive selective regions identified by Fst and XPEHH, respectively. These genes were further clustered in 50 gene ontology (GO) functional terms and 66 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in enrichment analysis. The GO terms and pathways were relevant with reproductive processes, e.g., oxytocin signaling pathway, thyroid hormone synthesis and GnRH signaling pathway, vascular smooth muscle contraction and lipid metabolism (alpha-Linolenic acid metabolism and Linoleic acid metabolism), etc. Based on the findings, six potential candidate genes ESR1, OXTR, MAPK1, RYR1, PDIA4, and CYP19A1, under positive selection related to characteristics of multiparous sheep breeds were revealed. Our results improve our understanding of the mechanisms of selection that underlies the prolificacy trait in sheep, and provide essential references for future sheep breeding.


Introduction
From the point of view of population genetics, when a novel mutation is subjected to the selection pressure over a long time, it will generate "selection signature", demonstrating some distinguished features on the genome, e.g., unusual linkage disequilibrium (LD) and changed population frequency [1]. Therefore, identifying the selection signatures underlying phenotypic difference can contribute to target causal variants for breeding, as well as explore the mechanisms of evolution. Furthermore, it can also help us to reveal the genetic basis of complex traits with phenotypic difference [2,3].   [4] detected strong signatures of selection in genes associated with the local adaptation of Tibetan sheep;   [5] genotyped 78 meat Lacaune and 103 milk Lacaune sheep to identify the selection signature related to ovine milk traits. Many methods have been proposed to detect pre-mentioned selection signatures, including Fst test based on population differentiation [6], the integrated Haplotype Homozygosity Score (iHS) [7] and the Cross Population Extend Haplotype Homozygosity Test (XPEHH) [8] based on linkage disequilibrium, etc. These methods have been widely applied in many studies to identify selection signatures. In this study, we employed the Fst and XPEHH test to detect the selection signatures between populations.
Litter size is one of the most important reproductive traits in sheep, as well as in other domestic animals, and has always been regarded as a critical index affecting the reproductive performance and productivity. Sheep presents variable litter size within and among breeds attribute to the natural and artificial selection for higher prolificacy during the long breeding history. For instance, Booroola Merino is well known as its reproductive characteristic with two offspring [9], and the Chinese local breeds, Hu sheep and small tail Han sheep are also characterized by high-prolificacy trait [10]. Many studies reported a series of genes or causal mutations associated with litter size through different genetic analyses in sheep [11]. Galloway et al. (2000) [12] detected FecXI and FecXH mutations in gene BMP15 associated with litter size in Romney sheep.   [9] reported a mutation FecBB in gene BMPR1B in Booroola Merino associated with the highly prolific phenotype. Vage et al. (2013) [13] found that FecGF mutation in gene GDF9 was related to larger litter size in Norwegian White sheep and Finn sheep through a genome-wide association study. Miao et al. (2016) [14] identified a set of differential expressed genes in different sheep breeds explaining the variation of fecundity by the integrated analysis of miRNAs and lncRNAs. Although the fecundity of sheep could be influenced by the interaction of environment (i.e., climate, nutrition, and stocking density) [15], previous researches indicated that genetic factor usually plays critical roles in the differential reproductive performance of sheep.
Four sheep populations (from one local Chinese breed and three local Kazakhstan breeds) were sampled, with the aim of exploring specific selection signatures for litter size in sheep. In addition, the geographic distribution of these sampled populations was shown in Figure A1. Of these, Aletai sheep, has a close relationship with Kazakstan sheep [16], mainly distributed in Xinjiang province, China, usually with one offspring [17]. Three Kazakhstan populations, Aldabas, Kurdyuchnyj and Karakul, with a high frequency of two offspring, were not only adapted to the local environment with harsh conditions but were also important genetic resources [18]. Our results will provide an essential reference for better understanding of the genetic mechanism of reproduction and further genetic improvement of litter size in sheep breeding.

Samples and Phenotype
According to the farms' reproduction records, ewes of Aletai sheep with single offspring and ewes of three Kazakhstan breeds with two offspring were collected. In total, 176 ewes were sampled in this study, including Aletai sheep (AL; n = 44), Aldabas (AD; n = 36), Kurdyuchnyj (KD; n = 48), and Karakul (KL; n = 48).

SNP Genotyping and Quality Control
DNA of each individual was isolated from blood, and all DNA samples were genotyped using 600 K Affymetrix Ovine HD genotyping arrays (Affymetrix, Santa Clara, CA, USA), including 633,619 SNPs across the entire sheep genome. The SNPs were re-mapped on the Ovis aries 3.1 genome assembly according to the position information provided by Affymetrix and SNPdb database. The whole procedure for collecting blood samples was carried out in strict accordance with the protocol approved by the Animal Welfare Committee of China Agricultural University (permit number DK996).
In order to improve the quality of genotyping data, genotype quality control was performed using PLINK v1.90 [19] for each group with the following criteria: (1) The individuals with >0.1 genotype missing rate were excluded; (2) SNPs with missing rate >0.1 were removed. After genotype quality control, three individuals (two from AL, one from KL) were removed, 591,189, 588,581, 570,692, and 560,856 SNPs remained for AD, KD, KL, and AL, respectively.

Principal Components Analysis (PCA), Population Admixture Analysis, and LD Decay
After filtering the individuals and SNPs, three methods were applied to analyze the genetic diversity among populations. Principal components analysis (PCA) was performed using the PLINK v1.90 software [19] to visualize patterns in relationships between individuals with filtered SNPs. Afterwards, population admixture analysis was performed by the ADMIXTURE program [20] to estimate the proportion of common ancestors among the four populations. Three scenarios of populations (K ranging from 2 to 4) were estimated using the cross-error estimation for genetic clustering, and the iteration times were set as 10. Additionally, levels of linkage disequilibrium (LD) for each sheep population were evaluated by genotypic correlation coefficient (r 2 ) using PopLDdecay software [21]. The options were set as: "-MaxDist 300 kb", and the visualization of LD decay among sheep populations across the whole genome was generated using self-written R script.

Identification of Selection Signatures
To detect potential selection signatures across the genome, Fst and XPEHH tests were employed to detect the selection signatures between breeds. The approaches were proved with high power in selection signatures with approximately fixed or fixed alleles. We firstly chose the uniparous AL as the common reference, and the other three multiparous populations as observed population, respectively. For each population pair compared of AL and one multiparous population, the common SNPs were used for calculating Fst and XPEHH values. The pairs of AD vs. AL, KD vs. AL, and KL vs. AL separately has 554,521, 555,037, and 547,884 common SNPs.
In this analysis, the single locus analysis method Fst proposed by   [6] was first employed to quantify the degree of population differentiation. The calculating of Fst was completed using VCFtools using non-window approach [22]. The Fst statistics ranges from 0 (identical population) to 1 (complete differentiation), and the top 1% SNPs were empirically considered as significant signatures in this study. Different from Fst, XPEHH is a haplotype-based method. In our study, haplotypes were firstly constructed in each breed by SHAPEIT [23], and XPEHH statistics based on the extended haplotype were calculated for each population pair using SELSCAN [24]. Since the XPEHH statistics approximately follow a normal distribution, the XPEHH values were normalized firstly. Then, the significance test of standard normal distribution (p < 0.05) was used to determine the variations caused by selection between populations, and the positive and negative XPEHH values represent the selection respectively occurred in the observed and reference population. Considering the selection regarding multiparous populations, the genes harbored in significant signatures with positive XPEHH values were used for further bioinformatics analysis.

Functional Annotation for Selection Signatures
According to the findings of selection signature, the common signatures of three population pairs were selected, and each core SNP of the common signatures was extended 200 kb towards upstream and downstream to be defined as selective regions. Candidate genes harbored in these regions were annotated based on the NCBI database (https://www.ncbi.nlm.nih.gov/). Because the annotation of the sheep genome is incomplete, corresponding human genomic information was regarded as a reference. The human orthologous genes were generated by the program of BioMart (http://www.biomart.org/). We further performed bioinformatics analyses to explore potential biological significance of genes harbored in these selective regions. A KOBAS 3.0 [25] (http://www.kobas.cbi.pku.edu.cn/kobas3) webserver was employed to perform enrichment analyses for biological processing GO (gene ontology) terms and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways. upstream and downstream to be defined as selective regions. Candidate genes harbored in these regions were annotated based on the NCBI database (https://www.ncbi.nlm.nih.gov/). Because the annotation of the sheep genome is incomplete, corresponding human genomic information was regarded as a reference. The human orthologous genes were generated by the program of BioMart (http://www.biomart.org/). We further performed bioinformatics analyses to explore potential biological significance of genes harbored in these selective regions. A KOBAS 3.0 [25] (http://www.kobas.cbi.pku.edu.cn/kobas3) webserver was employed to perform enrichment analyses for biological processing GO (gene ontology) terms and KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways.       Table 1   As shown in Figure 3b, the standardized XPEHH scores approximately followed normal distribution. Hence, outliers were identified through a normal test. The positive and negative XPEHH scores represent the selection happened in the observed and reference population, respectively. As shown in Table 1 and Figure 3a, for three population pairs, 16,935, 16,786, and 15,850 core SNPs were detected and correspondingly, 1399, 1409, and 1203 selective regions were identified in the observed populations AD, KD, and KL. Likewise, 13547, 13,332, and 12,789 core SNPs in 1416, 1471, and 1314 selection regions were detected in common reference population AL. Figure 3c,d show that 1337 common positive selection were detected in three observed populations, and 1735 common negative core SNPs were identified in AL. In total, 3111 common core SNPs were obtained in three population pairs (Figure 3c,d). Figure 3e presents the distribution of these core SNPs on chromosomes, reflecting positive selection signatures were mainly distributed on chromosome 2, 3, 4, 6, 14, and negative selection signatures were mainly observed on chromosome 1, 2, 3, 27. Particularly, chromosome 4 and 26 had positive selection signatures only, and chromosome 25 had only negative selection signatures. Additionally, nearly half of the common signatures detected by Fst were also captured by XPEHH test (Figure 4a). addition, these SNPs were harbored in 1705, 1716, and 1404 selective regions detected in corresponding population pair. The venn plot (Figure 2b) indicates a total of 177 common SNPs identified as selection signatures among three population pairs. Meanwhile, a bar histogram ( Figure  2c) illustrates the distribution of these common signatures on autosomes, indicating the selection mainly occurred on chromosome 1, 3, 7, and 14. As shown in Figure 3b, the standardized XPEHH scores approximately followed normal distribution. Hence, outliers were identified through a normal test. The positive and negative XPEHH scores represent the selection happened in the observed and reference population, respectively. As shown in Table 1

Functional Annotation
Through identification of selective signatures between one uniparous breed and other three multiparous breeds, 242 genes were identified in selective regions of Fst test, 798 and 1043 genes were detected in positive and negative selective regions of XPEHH test. Around 66% of genes identified by Fst were included in the gene sets by XPEHH. Of the genes identified by XPEHH, 30 genes were detected in both positive and negative selective regions (Figure 4b), it will be discussed later. Finally, 953 genes harbored in selective regions identified by Fst and XPEHH (positive selective regions) were used for further biological information analyses in the GO and KEGG databases using KOBAS 3.0 webserver.

GO Term Enrichment Analysis
According to the results of GO enrichment analyses, 953 genes under selection were further clustered in 50 GO functional terms ( Figure 5, Table S1), 35 out of 50 (70%) GO terms were classified as a biological process, such as cellular metabolic process, multi-organism process and developmental process, etc. Among the biological processes, the most abundant terms were biological regulation (GO: 0065007) with 92 genes, followed by the cellular metabolic process (GO: 0044237) with 88 genes. Membrane-bounded organelle (GO: 0043227) with 101 genes was the most abundant terms in the cellular component. Binding (GO: 0005488) with 114 genes was most dominant in molecular function. Among them, one important term closely related to reproduction, the multicellular organism reproduction was enriched, it contained 12 genes, including the previously reported ESR1, OXTR, and STAT5B (Table 2).

Functional Annotation
Through identification of selective signatures between one uniparous breed and other three multiparous breeds, 242 genes were identified in selective regions of Fst test, 798 and 1043 genes were detected in positive and negative selective regions of XPEHH test. Around 66% of genes identified by Fst were included in the gene sets by XPEHH. Of the genes identified by XPEHH, 30 genes were detected in both positive and negative selective regions (Figure 4b), it will be discussed later. Finally, 953 genes harbored in selective regions identified by Fst and XPEHH (positive selective regions) were used for further biological information analyses in the GO and KEGG databases using KOBAS 3.0 webserver.

GO Term Enrichment Analysis
According to the results of GO enrichment analyses, 953 genes under selection were further clustered in 50 GO functional terms ( Figure 5, Table S1), 35 out of 50 (70%) GO terms were classified as a biological process, such as cellular metabolic process, multi-organism process and developmental process, etc. Among the biological processes, the most abundant terms were biological regulation (GO: 0065007) with 92 genes, followed by the cellular metabolic process (GO: 0044237) with 88 genes. Membrane-bounded organelle (GO: 0043227) with 101 genes was the most abundant terms in the cellular component. Binding (GO: 0005488) with 114 genes was most dominant in molecular function. Among them, one important term closely related to reproduction, the multicellular organism reproduction was enriched, it contained 12 genes, including the previously reported ESR1, OXTR, and STAT5B (Table 2).

Pathway Enrichment Analysis
KEGG enrichment analyses detected a total of 66 pathways (p < 0.05) relevant to genes within selective regions (Table S2) and seven pathways exhibited the strong enrichment statistical signal (corrected p < 0.05). Figure 6 presents the top 20 significant pathways. Of these, for the multiparous sheep, the highly represented pathways were also associated with lipid metabolism (e.g., alpha-Linolenic acid metabolism, Linoleic acid metabolism, Arachidonic acid metabolism and ether lipid metabolism) and vascular smooth muscle contraction (e.g., the Vascular smooth muscle contraction and Calcium signaling pathway) ( Figure 6). Meanwhile, some known pathways related to reproduction were found to be significantly enriched. As shown in Table 2, five pathways were related to reproduction, and the genes involved in these pathways were also presented in Table 2, in all, 39 unique genes were identified relating to reproduction through these pathways.

Pathway Enrichment Analysis
KEGG enrichment analyses detected a total of 66 pathways (p < 0.05) relevant to genes within selective regions (Table S2) and seven pathways exhibited the strong enrichment statistical signal (corrected p < 0.05). Figure 6 presents the top 20 significant pathways. Of these, for the multiparous sheep, the highly represented pathways were also associated with lipid metabolism (e.g., alpha-Linolenic acid metabolism, Linoleic acid metabolism, Arachidonic acid metabolism and ether lipid metabolism) and vascular smooth muscle contraction (e.g., the Vascular smooth muscle contraction and Calcium signaling pathway) ( Figure 6). Meanwhile, some known pathways related to reproduction were found to be significantly enriched. As shown in Table 2, five pathways were related to reproduction, and the genes involved in these pathways were also presented in Table 2, in all, 39 unique genes were identified relating to reproduction through these pathways.

Candidate Genes Related to Litter Size
Through the functional annotation analyses, genes involved in reproduction-related processes or pathway were analyzed together to find candidate genes related to high prolificacy. Six candidate genes, presented in Table 3-OXTR and CYP19A1) were identified by both Fst and XPEHH, and another four genes were detected by XPEHH. All these genes were involved in GO terms or pathway. OXTR and ESR1 were presented in the reproduction-related GO term, multicellular organism reproduction. Meanwhile, OXTR, MAPK1, and RYR1 were involved in Oxytocin signaling pathway (hsa04921), ESR1, MAPK1, and PDIA1 in Thyroid hormone signaling pathway (hsa04919).

Candidate Genes Related to Litter Size
Through the functional annotation analyses, genes involved in reproduction-related processes or pathway were analyzed together to find candidate genes related to high prolificacy. Six candidate genes, presented in Table 3-OXTR and CYP19A1) were identified by both Fst and XPEHH, and another four genes were detected by XPEHH. All these genes were involved in GO terms or pathway. OXTR and ESR1 were presented in the reproduction-related GO term, multicellular organism reproduction. Meanwhile, OXTR, MAPK1, and RYR1 were involved in Oxytocin signaling pathway (hsa04921), ESR1, MAPK1, and PDIA1 in Thyroid hormone signaling pathway (hsa04919).

Discussion
In this study, genomic selection signatures related to litter size in sheep were detected using high-density SNP array in 176 sheep from one Chinese uniparous breed and three Kazakhstan multiparous breeds. The whole-genome single locus Fst statistics and haplotype-based XPEHH scores were calculated for each population pair comprising one multiparous breed and the uniparous breed. Compared to XPEHH, the selective SNPs and selective regions obtained by Fst were much less ( Table 1). As the haplotype-based method, XPEHH detected the core SNPs as the representative of the corresponding fixed regions based on haplotype, it can utilize more information of LD within one region, while the single locus Fst test only calculated the diversity of two loci between populations, even using the slide-window strategy, Fst could not make full use of the SNPs in the window together, and sometimes the division of window is not reasonable, due to the disperse of LD, resulting in the lower efficiency on the detection of selection signatures of Fst. In addition, different from using the normal test to determine outliers, SNPs located at the extreme 1% of the Fst values were considered as outliers empirically, and the threshold values for all three population pairs were near 0.15 (0.15, 0.13 and 0.15). This meant only the highly differentiated SNPs could be selected, and the moderate genetic difference (Fst values ranging from 0.05 to 0.15) [26] were ignored.
As the haplotype-based approach, XPEHH can identify positive and negative selection signatures in the observed and reference populations, respectively. In this study, no overlaps were found in positive and negative selection signatures, while 30 genes were detected both in positive and negative selective regions by XPEHH. It was mainly due to the overlaps between the positive and negative selective regions after extending the selected core SNPs, e.g., ENSOARG00000004311 gene (location: chr11, 26416218bp-26416934bp) was identified by the positive selective region (location: chr11, 26354799bp-27023145bp) and the negative selective region (location: chr11, 26107221bp-26565054bp). The corresponding core SNPs were located at chr11:26554799bp and chr11:26365054bp, the close distance of core SNPs lead to the overlap of positive and negative selective regions, further resulting in the same genes identified by positive and negative selection signatures.
The enrichment analysis found that the genes under selection were involved in 60 pathways, and five out of them were related to reproduction (Table 2). For example, oxytocin signaling pathway was known as the most well-established roles in stimulating uterine contractions during parturition and lactation containing 17 genes in all. Thyroid hormone synthesis was identified since it is essential for vertebrate embryogenesis and fetal maturation. Moreover, thyroid hormones triiodothyronine (T3) and thyroxine (T4) are critical for normal development, growth and metabolic homeostasis [27]. Besides, thyroid hormone deficiency during pregnancy was reported in rat, which caused a decrease of litter size [28]. GnRH signaling pathway, and Ovarian steroidogenesis are classical signaling pathways related to follicle development. GnRH signaling pathway has been shown to regulate gonadotropin-releasing hormone (GnRH), a precondition of the subsequent hormonal cascade which could induce the ovulation [29,30]. Ovarian steroidogenesis plays a role in normal uterine function, establishment and maintenance of pregnancy. The gene members of Ovarian steroidogenesis included the reported genes associated with sheep litter size, such as hormone regulate gene IGF1, and the famous main role gene of multiparous sheep, oocyte-derived factor BMP15 [31].
We also found that the genes under selection to be overrepresented in pathways related to lipids (e.g., alpha-Linolenic acid metabolism, Corrected p = 0.005; Linoleic acid metabolism, Corrected p = 0.02; Arachidonic acid metabolism, Corrected p = 0.05; ether lipid metabolism, p < 0.05), which were reported to be relevant to lipogenesis in Taihu pigs [32]. This finding indicated that the genes related to lipids traits have experienced intensive selection, and might be correlated with the fat tail trait of Aldabas and Kurdyuchnyj sheep in this study. In addition, vascular smooth muscle contraction pathway was found with the strongest enrichment statistical score (corrected p = 0.0027), and Calcium signaling pathway also was significantly enriched (corrected p < 0.05). Previous studies have reported that blood pressure is regulated by vascular smooth muscle contraction, which is triggered by an increase in intracellular free calcium concentration ([Ca 2+ ]) [33]. Although the data of blood pressure in multiparous sheep have not been reported during gestation, the higher blood pressure symptom is observed in human in the twin pregnancies than singleton pregnancies [34]. Moreover, vascular remodeling in the uterine and systemic circulation is important to meet the metabolic demands of the mother and developing fetus [35]. We inferred that these pathways overrepresented in multiparous sheep breeds may relate to the psychological changing in ewes in order to meet the higher metabolic demands during the twin pregnancies.
In this study, six potential candidate genes related to litter size experienced selection signatures in all three multiparous sheep populations (Table 3). All these genes were involved in the critical pathways or GO terms of the reproductive process. For example, the ESR1 (estrogen receptor-1) gene was found in thyroid hormone signaling pathway, which was a key gene affecting estrogen biosynthesis. Besides that, some studies reported that ESR1 plays a critical role in follicular growth and ovulation in ewes [36] which was also an important candidate gene of litter size in sheep [37]. As for MAPK1 (mitogen-activated protein kinase 1), it could mediates luteinizing hormone-induced breakdown of communication and oocyte maturation in rat ovarian follicles [38]. RYR1 (ryanodine receptor) was identified taking effect on regulating calcium release in oocytes [39], which also play roles in oocyte maturation [33,38]. The difference between the RYR1 genotypes were significant at the number of offspring in pigs [40]. Another candidate gene, PDIA4 (protein disulfide isomerase family A, member 4), one of the redox genes whose expression patterns are related to oocyte quality [41]. Also, PDIA4 gene was reported expressed in ovaries and associated with litter size in pigs [42]. According to our findings, OXTR (oxytocin receptor) and CYP19A1, which were identified by both Fst and XPEHH test could be the most potential genes affecting sheep litter size. Likewise, these two genes were repeatedly found in the results of enrichment analysis ( Table 2). In addition, the oxytocin receptor (OXTR) gene is known to be important during and after the ovulatory stimulus and expressing in many tissues, including brain, thymus, ovary, and testis [43]. During reproduction, OXTR could bind to OXT (oxytocin), while OXT plays a role in steroidogenesis, ovulation, luteinization, and luteal regression, in the mammalian ovary [44]. Recently research suggested that OXT is associated with larger litter sizes and signals of positive selection for OXTR forms were found in Cebidae sheep [43]. CYP19A1 gene in Ovarian steroidogenesis pathway encodes an estrogen-synthesizing enzyme aromatase, which is a mono-oxygenase and catalyzes many of the reactions associated with sterogenesis and the conversion of androgens to estrogen [45]. Previous studies reported that CYP19A1 played a critical role in gonadal development in sheep, and also involved in the development of follicular follicles and follicular atresia in bovine [46,47]. In Small Tail Han, one Chinese representative multiparous sheep breed, TIAN et al. (2019) [48] reported that CYP19A1 gene mainly expressed in ovary and hypothalamus, which suggested that CYP19A1 gene could promote the physiological function of ovary and negatively regulates hypothalamus during estrus or ovulation.
In order to explore the genetic mechanism of litter size on a molecular level, previous studies focused on genes relevant with ovulation rate, oocyte and follicle development [49,50]. These candidate genes (ESR1, OXTR, MAPK1, RYR1, PDIA4, and CYP19A1) may be the promising resource to explore the further mechanisms of high prolificacy in sheep, as well as other mammals.

Conclusions
With the help of whole genome selective sweep analysis, we detected selection signatures among the sheep genome using the Fst and XPEHH test. Gene enrichment analyses based on selection signatures suggested that six potential candidate genes related to litter size are worthy of further functional validation to reveal the underlying mechanisms of litter size in sheep. Our results will provide an essential reference for further genetic improvement of litter size in sheep breeding.