Hotan sheep, sourced from Hotan Huimin Animal Husbandry Technology Co., Ltd., were chosen as the research subjects, and their reproductive phenotypes, including the number of lambs per ewe, were recorded and statistically analyzed. Based on lambing records, multiparous ewes were categorized into single and twin lamb groups. Blood samples were collected from 180 healthy ewes aged 3 to 4 years (Supplement 1) and preserved at −80°C until further analysis. All experimental animal protocols were performed in accordance with the “Guide to Animal Experimentation” and were approved by the Use Committee under the norms of the Ethics Committee of Tarim University of Science and Technology (SYXK 2020-009). Thirty ewes with single lambs and 30 ewes with twin lambs, along with their corresponding lambing trait records (Supplements 2,3), were used to examine the relationship between single nucleotide polymorphism (SNP) loci and lambing. Genomic DNA was extracted using the standard phenol/chloroform extraction method. DNA integrity was verified using 0.75% agarose gel, and the purity of DNA (OD 260/280 ratio) was determined using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The DNA concentration was measured using a Qubit3.0 (Life Technologies, Carlsbad, CA, USA).

After the DNA fragments passed the test, they were interrupted using a Covaris ultrasonic crusher (Thermo Fisher Scientific) to concentrate the DNA fragments to approximately 200 to 400 bp, and a DNA library was prepared. Subsequently, the construction of the library is completed, the nucleic acid concentration is preliminarily quantified by using Qubit2.0, the library is diluted, and then the inserted fragments of the library are detected using Agilent2100 (Agilent, Beijing, China). When the DNA fragments met expectations, the effective concentration of the DNA library was accurately quantified using quantitative polymerase chain reaction to ensure quality. Arrange computer sequencing for qualified libraries (DNBSEQ): Single-stranded circular DNA molecules replicate by rolling rings, and these DNBs were then sequenced by combined probe-anchored polymerization.

Fastp software was used for quality control and data filtering. BWA v0.7.17 software package [9] was used to align the filtered clean reads with the reference genome. Samtools v1.11 [10] was used to change the comparison results from SAM files to sorted BAM files, Samtools was used to sort and remove duplicates, and Python v3.9.2 scripts were used to count the comparison rate and coverage. Mutation sites were detected using GATK v4.2.5.0 [11]. All mutation sites were annotated using ANNOVAR (last accessed July 5, 2024) [12]. GATK v4.2.5.0, was used to filter the SNP and insert deletion mutation (INDEL) sites. Subsequently, vcftools v0.1.16 [13] were used to exclude loci with a minor allele frequency (maf) of less than 3% and a genotype deletion rate greater than 20%.

ADMIXTURE v1.3.0 [14] was used to infer the population structure and number of genetic clusters (K) ranging from 2 to 5. Using PLINK v1.9 [15], we performed quality control of the genotype dataset using multiple filter options (−geno: 0.1; −maf: 0.01). Subsequently, we used PHYLIP v3.695 to construct a phylogenetic tree and MEGA v11.0.13 to visualize the kinship data. Principal component analysis (PCA) was conducted using PLINK v1.9 [15] to derive the principal component values for each sample. A scatter plot of principal components was generated using the scatter plot 3d function from the R package to further explore the genetic structure of the population, and the decay of linkage disequilibrium (LD) analysis was performed using PopLDdecay v1.5 software package [16].

To identify the selection signal, we used vcftools v0.1.16 to calculate the whole-genome FST and nucleotide diversity (π) using the sliding window method. The parameters used were an FST-window size of 100,000 and FST-window step of 50,000, with each interval between the two populations being 100 kb as the window size and 10 kb as the step size. To minimize false positives, we focused on regions with substantial FST values and extremely low or high θ π ratios (5% left tail and 95% right tail) as candidate genomic regions. These regions were visualized using the R package cmplot.

Furthermore, we annotated gene functions within the selected regions by comparing them with the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases. Based on the literature, genes related to the lambing trait of sheep were screened. To explore protein interactions among the screened genes, protein interaction analysis was performed using STRING. Following these steps, we aimed to identify the genomic regions associated with selection signals, annotate gene functions, and explore protein interactions related to lambing trait in sheep.

To ensure that 60 DNA samples extracted from Hotan sheep could be used for subsequent database construction, the purity and integrity of DNA samples were detected by 1% agarose gel electrophoresis and ultraviolet spectrophotometry. As shown in Figure 1 (gel electrophoresis results of some DNA samples), the DNA extracted for constructing the Hotan sheep resequencing DNA library has a single bright band, which shows that the DNA integrity is good and can meet the needs of sequencing. The DNA extracted from serum had clear bands, no obvious tailing, and high purity and concentration, and all DNA samples met the requirements of sequencing and database construction.

Table 1 presents the phenotypic differences in reproductive characteristics among the various lambing trait of Hotan sheep. Upon analysis of the phenotypic data from Hotan sheep used for blood sample collection, it was observed that the rate of twin lamb (TLE) was notably low, at only 21.8%. Although no substantial differences in pregnancy rates were observed among single lamb (SLE), TLE, and Hotan sheep, variations were noted in estrus onset and cycle duration. TLE Hotan sheep exhibited earlier estrus than SLE Hotan sheep, and the duration of the estrous cycle was longer in TLE Hotan sheep than in SLE and Hotan sheep (Table 1).

The reference sheep genomes used in this study are shown in the following table (Table 2). In the present study, 60 Hotan sheep were re-sequenced (Table 3). The data range of the raw data was 90,554,556 to 186,357,182, and the data range of clean reads was 90,551,892 to 186,347,082 after software quality control filtration, with a filtration rate of 99.98% to 99.99%. The sequencing quality was high (Q20≥97.35%, Q30≥ 91.50%). The average cytosine-guanine content was 41.95%, with an average of 1.15 genes in each interval of each chromosome, and an average of 1617 SNP loci was distributed in each interval (Figure 2).

Hotan sheep data were compared with the sheep reference genome using the BWA software Men program, and the comparison rate and average coverage were used as indicators to evaluate data reliability and accuracy. A general comparison of the samples is presented in Table 4. After alignment, the number of reads was between 20.72 and 47.44 GB, the comparison rate of samples was between 99.81% and 99.98%, and the average sequencing depth was 6.77× (4.66×–9.90×), indicating that the database of samples is well established, with an average coverage of 97.70% (94.93% to 98.96%), and indicating a good comparison that could be used for the next step of population genetic structure analysis.

According to the distribution in the chromosome, we found that the distribution trend of INDEL was consistent with that of SNPs, and both were mostly distributed in intergenic regions. To better understand the distribution of SNPs and InDels in sheep chromosomes better, 100Kb windows were set up, and the density maps of SNPs and InDels in each window were counted. It could be seen that SNPs and InDels had good homogeneity (Figures 3A, 3B).

As can be seen from Table 5, according to the referenced annotation information of the sheep genome, the number of intergenic regions in which variation is annotated to coding genes is the largest, with 31,778,335. Second, the intron region of the coding gene included 16,102,157, with 264,715 mutations in the upstream 1kb region of the coding gene, 255,853 mutations in the downstream 1 kb region of the coding gene, and 5,261 mutations in both the upstream and downstream 1 kb regions of the coding gene. There were 5 mutations in the 5′ untranslated region of the coding gene were 5, and 253,259 mutations in the exon region of the coding gene. There were 1,442 regional variations in the shear sites of the coding genes. There were four mutations that simultaneously existed in the exons and cleavage sites of coding genes.

For the functional annotation of mutations located in exons, the number of synonymous mutation regions was the highest at 130,992. Followed by nonsynonymous mutation regions, including 109,724; the number of frameshift mutation regions was 5,922; The number of frameless mutation regions was 3,241; The number of unknown mutation regions was 2814; The number of termination mutation regions was 2,432; The number of non-terminating mutation regions is 252.

PCA analysis based on the SNPs in autosomes (Figure 4) showed that most SLE and TLE groups were not obviously separated, and only a few groups were scattered in the opposite direction. Considering the existence of outliers, a few subgroups may not represent the whole sample, and these samples may not have an obvious structure. The results of the adjacent phylogenetic trees are shown in Figure 5. Geographically, the individual compositions of the two Hotan sheep populations are not clear.

To detect LD attenuation in Hotan sheep, the r2 values of the SLE and TLE groups were calculated. The results of the LD attenuation in this experiment are shown in Figure 6. These findings indicate that the SNP distance between the SLE and TLE groups in Hetian sheep decreased rapidly in the range of 0 to 50 kb, with r2 reducing from 0.4 to 0.1, and the LD attenuation rate between the two groups was nearly identical. The decay rates of Hetian sheep with different litter sizes were very similar, which suggests that differentiation was not considerable. This may be attributed to individuals originating from large populations with similar origins, and not being subjected to varying degrees of (artificial or natural) selection.

Structural analysis showed that the Cross-Validation error was the smallest when K = 2, and SLE and TLE began to separate when K = 3. As the k value decreased, the separation gradually became obvious, and the SLE and TLE were distinctly separated when K = 2 (Figure 7).

Based on selective elimination analysis (Figure 8), we analyzed the common areas of FST and Pi rates in SLE and TLE, wherein the regions above the threshold line (top 0.05%) were considered. Subsequently, 332 genes were annotated within these areas (Supplement 4). The GO and KEGG pathway analysis results are shown in Figure 9. Upon review of relevant references for 332 genes, it was found that 13 genes (BMPR2, ID2, SMAD7, PLCZ, PAK1, ELOC, RAD50, H2A, H2B, THBS1, RBX1, H4, and ABCC4) were associated with litter size. These genes were enriched in the TGF-β signaling pathway, renal cell carcinoma, inositol phosphate metabolism, non-homologous terminal junction, ATP-binding cassette (ABC) transporters, and NET pathway. Furthermore, protein interaction analysis revealed interactions among BMPR2, ID2, SMAD7, PAK1, ELOC, RAD50, THBS1, and RBX1 (Figure 10), providing an in-depth understanding of the interrelationship of these genes in the context of the regulation of lambing trait in Hotan sheep.

Population Pi can accurately analyze the ancestral group and individual lineage composition of livestock and poultry populations. Techniques such as PCA, genetic map analysis, and population genetic structure analysis can provide comprehensive insights, and the results can be verified. Such analysis is crucial for the effective utilization and protection of population genetic resources. In this study, the population structure of Hotan sheep with different numbers of lambs and Tian sheep was studied. It was found that the population structures of the two groups were not completely separated, and a small number of sheep individuals were genetically mixed, which may be because some individuals in the two groups had the same genome ancestor. For example, the phenotype of the individual “105” is the low fertility group, however, the results of population genetic structure show that it is a high fertility group, and the experimental results are similar to those of Wang et al [17]. Considering that the ewes sampled had the same feeding environment and management mode, it may be due to individual or environmental factors that ewes with a high reproductive genetic basis have fewer lambs. The characteristics of natural selection events and reproductive behavior experienced by different species or groups can be reflected through a linkage imbalance. Domestication and selection will lead to a decline in population Pi and strengthening of linkage between loci. The higher the degree of domestication, the greater the selection intensity and the lower the LD decay rate of the population. However, the attenuation curves of the SLE and TLE groups were similar, the attenuation speed was fast, and the attenuation distance was almost the same, which indicates that the degree of domestication of the Hetian sheep population was low, and the selection intensity was small. The above results show that there is a certain degree of population genetic differentiation between the SLE and TLE groups, but the degree of differentiation is small, which may be from the same population in the same sheep farm or the population has not been strongly artificially selected, which is in line with the actual situation.

Selection elimination analysis is primarily used to identify the genes selected during domestication and those selected during the adaptive evolution of species. The lambing trait, a quantitative trait regulated by multiple genes and factors, has been subjected to positive selection to a certain extent during the evolution and development of sheep. In this study, Hotan sheep of the same population with different numbers of lambs were selected for selection elimination analysis to screen for candidate genes related to the lambing trait. The FST values of the high- and low fertility groups were calculated, the Fst values were corrected to ZFst values, and the Pi values were calculated using the sliding window method with a window size of 100 kb and step size of 10 kb. Subsequently, the θ π ratio (pi.SLE.windowed.pi/pi.TLE.windowed.pi) was calculated, and the common window of the joint analysis of ZFst and θ π ratio was selected as the candidate window. Setting the top 5% as the threshold, windows with both the top 5% ZFst and θ π values were designated as candidate regions, ultimately yielding 290 shared genomic regions. From these 290 regions, 332 candidate genes were identified by annotation. These genes are primarily concentrated in TGF-β, renal cell carcinoma, phosphoinositide metabolism, non-homologous terminal junction, ABC transporter, and the NETs pathway. The association of PAK1 with the Hotan sheep lambing trait is further supported by a study on Qira Black sheep [18], which also suggested a possible influence of PAK1 on precocious puberty onset.

Among the pathways detected in this study, transforming growth factor-β (TGF-β) was the most closely related to the lambing trait. TGF-β is a pleiotropic cytokine that regulates many cellular processes, such as cell growth, differentiation, apoptosis, migration, cell adhesion, and immune response [19]. TGF-β family ligands are intricately related to the control of ovulation and fertilization and the establishment and maintenance of pregnancy. In mammals, even the earliest stages of reproductive development, including the specification of male and female species, are controlled by TGF-β-related proteins [20]. SMAD7 has been proven to play a role as an intracellular antagonist of TGF-β family signaling, which is a key regulator of TGF-β bone morphogenetic protein (BMP) signaling through a negative feedback loop, and inhibits TGF-β signaling by competing with receptor-regulated SMADs, which are involved in oocyte-somatic cell interaction and regulation of granulosa cell function [21]. Gao et al [22] regulated the expression level of SMAD7 in ovarian granulosa cells by overexpression and small interfering RNA knockdown, thus affecting follicular development in mice. SMAD7 signal transduction can reduce the dependence of TGF-β on ID2 [23]. As the inhibitor of the DNA binding (ID) gene is the main target of BMP/SMAD signal transduction, downregulation of ID2 leads to the inhibition of TGF-β. da Silveira et al [24] and others screened exosomes isolated from mares’ follicles in the middle estrus and preovulatory periods and identified microRNA (miRNA)s that are expected to regulate TGF-β/BMP signal members. ID2, a predicted target of exosome miRNAs, was found to exist in granulosa cells and exosomes in the follicular fluid of the middle estrus and preovulatory periods. Simultaneously, BMPR2, the receptor for BMP4, regulates germ cell activity and activates ID2 [25]. BMPR2 is a type II receptor of the TGF-β family that exhibits serine/threonine protein kinase activity. BMPR2 can be detected at all stages of follicular development and plays an important role in promoting sex hormone secretion, ovarian development, follicular growth, and granulosa cell proliferation. Studies have shown that BMPR2 plays an important role in the growth and function of mammalian and chicken follicles [26]. THBS1 is an important activator of TGF-β, which can change the conformation of the TGF-β protein, expose its binding site with cell receptors, and activate the TGF-β signaling pathway [27]. Overexpression of THBS1 substantially activates the TGF-β pathway, and the TGF-β signaling pathway can be inhibited by down-regulating THBS1. At the same time, the expression of THBS1 in fibroblasts increases, promoting the growth and migration of fibroblasts and TGF-β expression [28].

TGF-β inhibits the expression of RAD50 and stimulates TGF-β expression in cultured fetal ovarian fibroblasts [29]. RAD50 not only plays a role in the repair of double-stranded breaks and homologous synapses during female meiosis but also affects the elimination of oocytes that have not been repaired after birth. This suggests that similar to spermatocytes, it may play a role in monitoring meiotic recombination during the prophase of meiosis. Knockdown of RAD50 leads to a decrease in the expression level of two proteins in spermatocytes (GC-2spd), and downregulation of RAD50 inhibits the proliferation of GC-2spd cells and promotes apoptosis, which is consistent with the simultaneous decrease in RAD50 in the testes of spermatogenic failure patients [30]. Changes in the expression of RAD50 in the testis may lead to an imbalance in male germ cell proliferation and apoptosis, ultimately damaging spermatogenesis. This also indicated that RAD50 may play a key role in the repair, proliferation, differentiation, and apoptosis of male germ cells. P21-activated kinase 1(PAK1) is an effector of the RhoGTP enzymes RAC1 and CDC42, which are involved in regulating mitotic events such as chromatin condensation and subsequent chromosome capture, movement, and separation. High expression of PAK1 inhibits TGF-β function [31]. Previous reports have shown that PAK1 is a key regulator of various neuronal processes and binds to the BMPI receptor, ALK2 [32]. Sperm-specific phospholipase Czeta (PLCz) is a candidate sperm-derived oocyte-activating factor that can trigger a series of characteristic physiological stimulation fertilization processes through cytoplasmic Ca2+ [33]. PLCz is expressed in the epididymis, acrosome of sperm, equatorial segment, head-middle junction, and the main part of the flagella. Some evidence suggests that sperm deliver PLCz when they fuse with oocytes. It has been reported that PLCz can penetrate the oocyte in the sperm acrosome, which is helpful for successful fertilization [34], and is the key protein in the process of sperm egg fusion in mammals. ABCC4 is mainly localized in the uterine cavity epithelial and glandular epithelial cells, and is expressed in the endometrium in a pregnancy and stage-dependent manner. ABCC4 plays a key role in supporting the establishment and maintenance of pregnancy by regulating prostaglandin E2 transport in the mother and fetus. Studies have shown that ABCC4 is expressed in fertilized eggs of pigs during early pregnancy and that ABCC4 may act on the endometrium and fertilize eggs in an autocrine or paracrine manner [35], thus promoting the establishment and maintenance of pig pregnancy. RBX1 plays a role in E3 ligase complexes, which regulate various biological processes including DNA repair and mitotic fidelity. During nucleotide excision and repair, histone H2A ubiquitination can repair damage regulated by RBX1 [36]. RBX1 is also an important part of the SCF complex, and its subcellular localization indicates that RBX1 plays an important role in the meiosis and maturation of mouse oocytes [37]. An imbalance in RBX1 leads to late stagnation and maintains cells in the germinal vesicle (GV) stage. RBX1 can bypass oocytes, the GV stage, and metaphase MII (II) stage.

Core histones are composed of four proteins (H2A, H2B, H3, and H4), and the study of histone modification modes is important to understand their functions in biological processes. H2A and H4 are conserved in the nucleus of mature mouse sperm, indicating that a certain amount of H2A and H4 conserved in the nucleus may represent a more accurate marker for mature mouse sperm [38]. In the prophase of meiosis, all H4 modifications are very high, and elongated sperm cells exhibit acetylation, which increases histone H4. Modification of the terminal tail of histone H4N changes considerably in male germ cells during spermatogenesis, suggesting that modification of histone H4 plays an important role during this stage of spermatogenesis. H2A-H2B dimer is located on both sides of the nucleosome, stabilizing the entire structure of the nucleosome and regulating its dynamics such as DNA unpacking and sliding. This epigenetic regulation is achieved through post-translational modifications and the incorporation of histone variants [39]. The bending dynamics of H2A-H2B dimers are regulated by the interaction of the H3-H4 tetramer, distortion defects in nucleosome DNA, and amino acid sequences of histones. H2A-H2B modified fast short-term memory is helpful for H3-H4 chromatin to return to a stable state. We propose that the short-term guidance provided by H2A-H2B in the process of DNA replication is based on the long-term maintenance of H3-H4 and that the exchange rate of H3A-H4 before and after H2A-H2B replication is higher [40].