Selection Signature Analysis for Important Economic Traits in Sheep (Ovis aries)

During the long process of domestication, sheep have shown a wide range of behavioral habits, physical characteristics and production performance, mainly in terms of meat, wool and milk production, as people have domesticated them for different purposes (Zeder, 2008; Bernardo et al., 2009). Due to differences in geography, environment, culture, customs and other natural and human factors, sheep for meat have in turn begun to differentiate into different breeds. Various studies demonstrate that genomic selection became a preferable choice for selection of appropriate candidate for animal breeding research owe to potential benefits and advantages over traditional breeding methods (Lopez et al., 2019; Budhlakoti et al., 2021). The targeted selection of different sheep breeds for target genes during anthropogenic selection has led to an increase in the frequency of dominant alleles in their populations. According to the neutral genetic hypothesis, closely linked neutral loci in their vicinity exhibit the same pattern of variation due to factors such as genetic drift, a phenomenon known as selection signature or selective sweep (Nielsen et al., 2007). Selection signals are the imprinting of selection on the genome, in the form of reduced polymorphism at certain loci or DNA fragments, increased linkage disequilibrium values or pure genes.

Nucleotide diversity is often used to measure the degree of genetic variation in a population. In its original form proposed by Nei, nucleotide diversity is defined by the average nucleotide differences per site between any two sequences (Nei et al., 1979). A sequence is taken at random from multiple samples of DNA obtained from a population and the differences in bases at the same locus in these sequences are averaged to obtain a value for nucleotide diversity. This value is denoted as “π” or Pi. The Pi method is still commonly used for selective sweep analysis in individual populations.

If regions subject to selection between different populations or breeds are to be detected, commonly used methods include the population differentiation index (Fst) method. The fixation index (Fst) is one of the most important indicators for detecting genetic differentiation between populations and estimating Fst values for individual loci on a genome-wide scale allows analysis of the degree of differentiation at each locus, ultimately enabling selection signal detection. Detection of selection signals between populations based on Fst values was originally employed in human evolutionary genetics (Akey et al., 2002) and is now widely used in a variety of livestock species (Liu et al., 2012; Gao et al., 2020; Tang et al., 2020; Ahmed 2021).

With the introduction of SNP microarrays, it has become possible to examine population genetic differences at the whole genome level. Selective sweep studies are now widely used across livestock species with the aim of explaining the molecular mechanisms of inheritance hidden beneath complex economic traits. Many scholars have done extensive sequencing and selection signals in different sheep and goat breeds and have identified many candidate genes associated with the trait (Guan et al., 2016; Yuan et al., 2019; Guo et al., 2018; Zhu et al., 2021). It is clear from the above description that the vast majority of studies have been based on a genome-wide selective sweep of extreme traits, including altitude, reproductive traits, tail type, in indigenous sheep breeds. However, little research has been reported on the use of selective sweep approaches to elucidate the molecular mechanisms of differentiation between introduced and indigenous sheep breeds.

The introduced breed, including Dorper and German Merino and indigenous breed Sunit were all considered to be excellent meat sheep breeds. However the growth development, including body height, chest circumference and reproductive performance of the introduced breeds of Merino and Dorper sheep are found to be better than that of the indigenous breeds of Sunit sheep. During breed formation, selection leaves a selection signature on the genome as a result of selection. It is hypothesized that their good meat performance and reproductive performance may be linked to marks left on the genome over time by evolutionary processes. Therefore, this paper takes as a starting point the differentiation mechanisms between introduced and indigenous breeds in meat performance and reproductive performance and attempts to investigate the underlying genetic molecular mechanisms by means of a selective sweep strategy, with the expectation that the selected candidate genes will provide some ideas and insights for future sheep germplasm resource mining and foreign introductions.

Experimental materials

In this study, Illumina OvineSNP50 BeadChip containing 54241 SNPs from three breeds, consisting of 20 Sunit sheep, 20 German Mutton Merino sheep and 20 Dorper sheep, were downloaded from public databases (www.sheephapmap.org).  The three sheep breeds were named GMM, DOR and SUN, for subsequent bioinformatics analysis. The experiments were conducted at the Cloud Server Centre of Qinghai University and spaned the period from June 2021 to March 2022.

Analysis of the genetic structure of the population

Quality control was performed using plink software (9.2) with the criteria of (1) minimum allele frequency (MAF) > 0.05, (2) individual detection rate (call rate) > 0.9, (3) SNP detection rate > 0.9 and (4) Hardy Weinberg equilibrium P-value > 10-6. In order to better annotate the existing sheep genome (Ovis aries v3.1), the SNP positions were readjusted to Ovis aries v3.1 using plink software, resulting in a total of 43,299 SNPs for subsequent analysis. After the quality control was completed, PCA analysis was carried out using the plink -pca command and the glPca function in R. The ancestral components of the population were analyzed using admixture software. The distance tree was constructed using the brewer.pal function and the minimum spanning tree was analyzed using the bitwise.dist and poppr.msn functions.

Comparative analysis of nucleotide diversity in different populations

In this experiment, nucleotide diversity p was calculated for each of the three populations using bcftools and vcftools software, with a sliding window (-window-pi) set to 100000. All data were visualized using ggplot2.

Genome scans for selection using Pi within the group

The three breeds are considered as a population and the Pi method is used to find regions within the population that are subject to selection. The principle is to calculate the population nucleotide diversity p. The population nucleotide diversity p was calculated using vcftools software, with a sliding window (-window-pi) set to 100000. All data were visualized using ggplot2 and qqman, with the sugestiveline value in Manhattan set to -log10 (2e-06).

Genome scans for selection using Fst between species

In this study, the Fst method was used to detect selection signals between breeds. The allele frequency difference between the natural and free-mating populations for each SNP was measured using Pair-wise FST, as described by B.S. Whir (Weir et al., 1984).

Annotation and functional enrichment analysis of candidate genes

In order to accurately identify the genes selected in the breeds, we further intersected the candidate genes annotated using the Fst with those annotated using the Pi to create a Venn diagram of the selected candidate genes within and between breeds. This study is based on the Gene ontology (GO) database and KEGG annotation data to annotate candidate gene for gene function (Huang et al., 2009). In this study, R and Linux were used for data processing and analysis, statistical mapping and text information mining.

Population genetic structure

A total of 47015 SNPs and 60 individuals were selected for further analysis. After quality control, 43299 SNPs were used for PCA, distance tree and ADMIXTURE analysis. The PCA results showed that all animals can be divided into three groups (Fig 1a). The results from the distance-tree analysis showed that these samples could be appropriately divided into three subgroups (Fig 1b). The results from the ADMIXTURE analysis showed that the least amount of cross-validation error occurred when K=3  Fig 1c indicating that K=3 was the optimal modeling choice. Therefore, these samples could be appropriately divided into three subgroups (Fig 1), which was consistent with PCA and distance-tree analysis.

Fig 1: Population structure of the three sheep populations. (a) principal-component analysis for population stratification in three sheep breeds; (b) distree-analysis for three sheep breeds; (c) when K=3, the least amount of cross-validation error occurred; (d) makes it fairly clear that K=3 was the optimal modeling choice. The green background represents the Sunit sheep group (SUN); the red background represents the German Mutton sheep group (GMM); the blue background represents the Dorper sheep group (DOR).

Genome scans for selection using Pi within the group

The Pi calculations showed that the nucleotide diversity of the indigenous breed of Sunit sheep was higher than that of the introduced breed represented by the Dorper and German Merino sheep (Fig 2a), suggesting that the nucleotide diversity of the introduced breed was fully exploited compared to the native, while for the indigenous breed of Sunit, its higher nucleotide diversity implies a greater potential for future selection of this breed. Here we speculate that genomic nucleotide diversity may be the potential reason for differences in growth development and reproductive performance between introduced and indigenous breeds. In order to explore the molecular mechanisms underlying the differences, it is an urgent task to adopt selection signals to explore candidate loci across the genome. The Pi method was used to treat the three breeds as a population and then calculate the population nucleotide diversity p. The results showed that there were large differences in the population SNP density distribution on different chromosomes, mainly concentrated on chromosomes 1 to 6 (Fig 2b). Smaller values of p indicate greater susceptibility to selection. Manhattan results using the Pi strategy revealed the presence of 167 selected QTL regions within the population (threshold of 2e-06) and found 293 candidate peak SNP loci. The genomic location, size, peak SNP and peak genes in the selection QTL regions identified using Oar_v3.1 were summarized in Table 1. Eventually, 248 candidate genes were identified in 293 peak SNP loci (Fig 2c).

Fig 2: Population nucleotide diversity analysis and selective sweep analysis using Pi strategies. (a) Nucleotide diversity comparison among three sheep population; (b) SNP density distribution on chromosomes in sheep population; (c) PI-based selective sweep analysis in sheep population.

Genome scans for selection using Fst between species

The Fst strategy is to consider each of the three species as a group and then use the vcftools software to calculate the Fst value, the larger the Fst value, the more likely it is to be selected. For the DOR_VS_GMM comparison group, 256 selected QTL regions (threshold 0.4) were successfully obtained and found 1804 peak SNP loci. Eventually, 961 candidate genes were identified from 1804 peak SNP loci (Fig 3a). For DOR_VS_SUN comparison group, 78 selected QTL regions were successfully obtained and found 469 peak SNP loci. Eventually, 252 candidate genes were identified from 469 peak SNP loci (Fig 3b). For GMM_VS_SUN comparison group, 84 selected QTL regions were successfully obtained and found 487 peak SNP loci. Eventually, 217 candidate genes were identified from 487peak SNP loci (Fig 3c).

​

Fig 3: Analysis of selective sweep based on the Fst strategy. (a) DOR_VS_GMM; (b) DOR_VS_SUN; (c) GMM_VS_SUN.

Functional enrichment analysis of candidate genes

In this experiment, two strategies (Fst, Pi) were used for selection sweep analysis in order to screen for more precise target genes. The genes identified by the two strategies were then intersected, which greatly reduced the QTL range and thus the number of false positives. As depicted in Fig 4 and Table 1, we identified a total of 32 genes in the DOR_ VS_GMM group, including MAD2L1, PDE5A, ARHGAP42, ZBTB42, DOCK3, UNCX, CHRM3, CHORDC1, SV2A, EHBP1,  RUNDC3B, SNORA70, PRMT6 and ALCAM. Of these genes, the MAD2L1 gene was reported to be significantly associated with stemness characteristics of uterine corpus (Zhang et al., 2021). Previous studies have shown that the PDE5A gene is strongly associated with sexual function, regulation of vascular tone, Leydig cell secretory function and sperm motility (Made et al., 2010; Gebska et al., 2011; Dimitriadis et al., 2009). It is hypothesized that PDE5A and MAD2L1 are responsible for the mechanism of difference in growth development, reproductive performance between Dorper and Merino sheep.

Fig 4: Venn diagram based on Fst and Pi method for annotation of candidate genes. PI_population represents candidate genes obtained through the Pi strategy, DOR_VS_GMM, DOR_VS_SUN and GMM_VS_SUN represent candidate genes obtained through Fst strategy.

Table 1: List of candidate genes obtained in the sheep population based on selective sweep.

A total of 26 genes were identified in the DOR_VS_SUN group, such as SV2A, EHBP1, RUNDC3B, SNORA70, PRMT6, ALCAM, EPHA5, CSRP1, IPO9, TNNI1, SHISA4 and MBNL1. Of these genes, EH domain binding protein 1-like 1 (EHBP1), a protein with calponin homology domain, was reported to be highly expressed in the testis tissue of rat (Venditti et al., 2020). SNV and haplotype analysis reveals new CSRP1 variants associated with growth and carcass traits in cattle (He et al., 2013). Slow skeletal muscle troponinI 1(TNNI1) located in slow-twitch myofibers that provides a calcium-sensitive switch for the contraction of striated muscle and plays an important role in muscle growth (Ji et al., 2018). Some studies found that MBNL1 isoforms exhibit differential influence on enhancing brown adipogenesis and found mice lacking MBNL1 exhibit sudden cardiac death and molecular signatures recapitulating myotonic dystrophy (Hung et al., 2019; Lee et al., 2022; Huin et al., 2013). It is hypothesized that EHBP1, CSRP1, TNNI1, MBNL1 may be responsible for the differences in growth development between Dorper and Sunit sheep.

A total of 25 genes were identified in the GMM_VS_SUN group, including MOV10L1, HDAC10, SNCA, KAZN, AADACL3, WFDC8, TNNT1, RDH13, SLC6A3, TRIM67, TUBGCP6, TNNI3, SNORA19, BMP2, PPP1CC and MBNL1. Of these genes, the MOV10L1 gene was found to be involved in spermatogenesis and its gene polymorphisms were significantly associated with male sterility (Sarkardeh et al., 2014; Zhu et al., 2015). Similar to the results of this analysis, 460 fine-wool sheep were resequenced and performed GWAS study of body weight traits and identified a number of candidate genes, including AADACL3, suggesting that the AADACL3 gene may be an important candidate molecule for controlling body weight traits (Lu et al., 2020). It was found that the Wfdc8 gene was highly expressed in the male reproductive tract of the rat (Rajesh et al., 2011). Similar to the results of this analysis, 16 bulls with significant differences in daily weight gain and feed intake were performed transcriptome analysis and identified 1831 differential genes, including the RDH13 gene (Lindholm-Perry  et al., 2017), suggesting that the RDH13 gene may be an important candidate molecule for controlling daily weight gain and feed intake traits in sheep. Related studies have shown that showed the TUBGCP6 gene responsible for microtubule organization (Hasan et al., 2019). Troponin (TNNI3) is a member of a multigene family of calcium-binding proteins that play a key role in muscle contraction and influence meat quality traits in livestock (Lipska et al., 2009; Yang et al., 2009). The TNNI3 gene was also obtained in the results of this analysis, suggesting that this gene may be a candidate molecule for differences in meat-producing traits between introduced and indigenous breeds of sheep. Some literature showed that BMP2 played a major role in the muscle development and bone formation, also found it can be influence the proliferation of sheep myoblast (Ma et al., 2017; Lu et al., 2018). PPP1CC2 levels are highly regulated in testis to ensure normal spermato genesis and male fertility (Dasgupta et al., 2012; Sinha et al., 2012; Sinha et al., 2013). It is hypothesized that AADACL3, RDH13, TNNI3, BMP2, Wfdc8, PPP1CC2 and MOV10L1 may have left deep imprints on the genome during long-term domestication and evolution, which may be the underlying cause of the differences in growth development and reproductive performance between Merino and Sunit sheep.

In order to further extract known biological meanings from these candidate genes, GO and KEGG enrichment analyses were implemented using the OmicShare tools. GO functional enrichment indicates that these genes are mainly enriched in muscle development related terms, such as skeletal muscle contraction, smooth muscle contraction, multicellular organismal movement, musculoskeletal movement, dopamine biosynthetic process, cardiac muscle tissue morphogenesis (Fig 5). KEGG enrichment analysis shows that these genes are mainly enriched in pathways related to muscle development and reproductive regulation, such as regulation of actin cytoskeleton, ECM-receptor interaction, platelet activation, oocyte meiosis, adrenergic signaling in cardiomyocytes, hippo signaling pathway, cAMP signaling pathway, calcium signaling pathway (Fig 6). The genes identified above are only bioinformatic predictions and their functions are not known. Future gene function validation tests are needed to determine the authenticity of the candidate molecules identified. These genes can be used as candidate genes for growth development and reproduction in sheep and provide useful information for the production and genomic selection of meat sheep.

Fig 5: GO enrichment of candidate genes.

Fig 6: KEGG enrichment of candidate genes.

Selective sweep based on Pi and Fst strategy initially elucidated the genetic-molecular mechanisms underlying the differentiation of growth development and reproductive traits between introduced and indigenous breeds. In the DOR_VS_GMM group, 32 genes were identified. In the DOR_VS_SUN group, 26 genes were identified. In the GMM_VS_SUN group, 25 genes were identified. GO, KEGG enrichment analysis showed that these genes are mainly involved in the pathways of muscle development and reproductive regulation. We speculated that EHBP1, CSRP1, TNNI1, MBNL1, AADACL3, RDH13, TNNI3, BMP2, Wfdc8, PPP1CC2 and MOV10L1 may have left deep marks on the genome during long-term domestication and evolution and may be responsible for the productive and reproductive performance in the Merino and Sunit sheep. These genes have played an important role in the evolution of selection in the three breeds. It is expected that the selected candidate genes will provide some insights for future sheep germplasm mining and foreign introductions.

This study was supported by the Youth Fund of Qinghai Provincial Science and Technology Department (2021-ZJ-978Q), the Key Project of Youth Foundation of Qinghai University (No.2020-QNY-1) and the “Thousand Talents Plan for High-end Innovative Talents” of Qinghai Province.

The authors declare no conflicts of interest.