Genome-Wide Association Analyses and Population Verification Highlight the Potential Genetic Basis of Horned Morphology during Polled Selection in Tibetan Sheep

Simple Summary The horned type of sheep, including the polled phenotype and scurs phenotype, emerged during the process of artificial domestication. Sheep serve as an optimal animal model for investigating the correlation between quality and quantity traits of horns, as well as the multi-gene regulatory mechanism underlying these traits. While significant progress has been made in understanding the effects of natural selection and sex selection on horn phenotypes, there is limited research on polymorphisms resulting from artificial selection. Through genome-wide association analysis, this study reveals the genetic mapping of horn phenotypes, polled genes, and superior haplotypes to provide valuable insights into studying the genetic mechanisms governing horn traits and other specific characteristics in ruminants while also laying a foundation for breeding. Abstract The types and morphology of sheep horns have been extensively researched, yet the genetic foundation underlying the emergence of diverse horn characteristics during the breeding of polled Tibetan sheep has remained elusive. Genome-wide association analysis (GWAS) was performed on 103 subtypes (normal large horn, scurs, and polled) differentiated from G2 (offspring (G2) of parent (G1) of polled) of the polled core herd. Six single nucleotide polymorphisms (SNPs) located on chromosome 10 of the relaxin family peptide receptor 2 (RXFP2) gene exhibited positive correlations with horn length, horn base circumference, and horn base interval. Furthermore, in genotyping 382 G2 individuals, significant variations were observed for each specific horn type. Three additional mutations were identified near the target SNP upstream of the amplification product. Finally, the RXFP2-specific haplotype associated with the horned trait effectively maintained horn length, horn base circumference, and horn base interval in Tibetan sheep, as confirmed by population validation of nine loci in a sample size of 1125 individuals. The present study offers novel insights into the genetic differentiation of the horned type during improvement breeding and evolution, thereby establishing a robust theoretical foundation for polled Tibetan sheep breeding and providing valuable guidance for practical production.


Introduction
In cases where strong directional selection is evident, genetic traits often fail to exhibit the expected evolutionary response [1].Therefore, elucidating the genes and genomic regions that contribute to trait variation provides an invaluable opportunity for comprehending the intricacies of genetic diversity processes within populations.The presence of adaptive variation at the genotype level may not be readily apparent when focusing on phenotypes for analysis, potentially due to genetic correlations between the phenotype of Animals 2024, 14, 2152 2 of 18 interest and other fitness-related traits [2].Genetic resources obtained through directional selection [3] were utilized in this study to identify the genetic polymorphisms associated with horn traits within the population.Additionally, the effects of genetic variation on specific traits were determined, thereby establishing a crucial foundation for understanding the relationship between genotype and phenotype traits following artificial intervention.
The cranial appendages of ruminant animals, commonly referred to as headgear, consist of a pneumatized bony core fused with the frontal bone of the skull and are covered by a continuously growing keratin sheath [4].Horn injuries to cattle caused by too-frequent and severe bumps are listed as one of the top 10 challenges to the U.S. beef supply [5].To reduce injury risk to handlers and other animals as well as economic losses, dehorning procedures are performed [6].Therefore, these appendages have become unpopular in modern captive breeding programs.The market-oriented sheep industry values its economic value because of animal welfare concerns [6][7][8].More and more artificial interventions have been employed in recent years to genetically breed sheep without horns; the occurrence of this has been documented in historical records dating back to ancient Egypt [9].However, there are still a few native breeds that naturally possess this trait.Prolonged human-mediated artificial selection has resulted in phenotypic variation within the genetic structure, which can be identified through patterns of genomic variation such as reduced nucleotide diversity and regions with very low heterozygosity.Importantly, the detection of numerous single nucleotide polymorphisms (SNPs) provides crucial information for identifying potential sites responsible for phenotypic variation.Therefore, the comprehensive identification and strategic deployment of specific SNPs via whole-genome sequencing (WGS) analysis promises substantial advancements in the domains of animal husbandry and molecular biology, particularly for the enhancement of production traits.A growing number of researchers are conducting extensive studies of various patterns of the presence and absence of horns, multi-horns, and deformed horns from around the world.With human domestication, domestic sheep breeds have also developed a hornless phenotype; breeders have been breeding polled sheep since the 16th century [10].Earlier studies proposed that the autosomal 10 locus Ho P , Ho + , and Ho hL alleles inheritance model controls the phenotype of horns, with the Ho P allele determining polledness [10][11][12].The horn allele has been mapped to OAR 10 such that there are at least two loci affecting the presence or absence of horns in Merino × Romney crosses [13].The horn locus has been fine mapped to a 50 kb region in the sheep genome and predictive haplotypes for the polled HP allele have been identified using a segregating (Merino × Romney) × Merino resource [13].RXFP2 is the principle candidate locus which supports a new model of horn-type inheritance in Soay sheep [14].The analysis of the sequence revealed a 1833-base-pair genomic insertion located in the 3 ′ -UTR region of the RXFP2 gene in polled animals [15,16].The main candidate gene associated with the type and size of sheep horns, RXFP2, was found to have several new mutation sites.Among these mutations, the synonymous mutation p.P375 (c.1125A>G) was identified as a potential indicator for determining the presence or absence of horns in Tan sheep [17].So far, the impact of two causal mutations on the process of horn formation remains unknown and, despite their influence on polled traits [18], no genetic variation has been identified that is causative of horn phenotype.Although the genetic loci associated with polledness have been mapped successfully via WGS, the observed variations in hornless sheep populations across different breeds are notable.Additionally, horn type in certain breeds is associated with variations in the RXFP2 gene region, indicating a highly complex genetic mechanism underlying the polled trait.These studies have only explained a small portion of the overall genetic variation; therefore, further experimentation is necessary to identify any causal sites related to the polled trait in ruminants.
Tibetan sheep are one of the prominent animal breeding resources unique to the Qinghai-Tibet Plateau.In order to elucidate the genetic characteristics of horn development, we established polled male and female groups (G1) based on phenotype in the adult phase and performed targeted mating, with the resulting offspring (G2) serving as our research subjects.Until now, there have been no reports on the genetic basis underlying horned differentiation enhanced breeding.
In this study, we conducted GWAS to analyze three horn phenotypes in G2 and subsequently identified specific haplotypes through population verification.By investigating candidate genes associated with polled traits at the genomic level, we were able to assess the molecular genetic characteristics of polled traits and provide fundamental data for localizing functional genes related to horn traits.This will serve as a crucial reference for the future molecular breeding of polled sheep.

Samples Collection and Whole-Genome Sequencing
The core group exhibited phenotypic variations in the second generation (G2), where polled rams and ewes (G1) were mated resulting in individuals with either normal large horns, scurs, or no horns.After continuously observing the horn status of a core population with a well-defined pedigree, we meticulously selected a total of 103 sheep for sequencing analysis, including 48 polled sheep, 29 sheep with normal large horns, and 26 sheep with scurs.The genotyping of 382 G2 individuals was conducted within a non-horned core group.The phenotypic values of the compiled horns were collected, and their whole blood was obtained and stored in EDTA vacuum tubes.The husbandry and management of these animals was consistently maintained at the same level.The QIAamp DNA Mini Kit (Qiagen, Hilden, Germany) was utilized to extract total genomic DNA from 200 µL of sheep blood, followed by assessment of DNA quality.Illumina's paired-end DNA Sample Prep kit standard library building procedure was employed to prepare sequencing libraries using a minimum of 3 µg genomic DNA.The Agilent 2100 bioanalyzer (Agilent, Santa Clara, CA, USA) was used for examination and real-time PCR quantification of the sequencing libraries prior to sequencing on the Illumina NovaSeq6000 NGS platform (Guangzhou, China).

Quality Control and Annotation
The Illumina data underwent filtration using FASTP (version 0.18.0) with the following criteria: (1) reads containing ≥10% unknown nucleotides (N) were removed; (2) reads with ≥50% of bases having Phred quality scores ≤ 20 were removed; and (3) reads with adapters were deleted.The resulting filtered clean reads were utilized for assembly analysis.The filtered reads were aligned to the reference genome using the Burrows-Wheeler Aligner (BWA) (version 0.7.12) alignment software (ncbi_GCF_002742125.1) with the settings 'mem 4 -k 32 -M', -k as the minimum seed length, and -M as an option used to mark shorter split alignment hits as secondary alignments [19].Variant calling was performed for all samples using GATK's Unified Genotyper (version 3.4-46) [20].SNPs and InDels were filtered using GATK's Variant Filtration with proper standards (-Window 4, -filter QD < 2.0 || FS > 60.0 || MQ < 40.0, -G_filter "GQ < 20"), and those exhibiting segregation distortion or sequencing errors were discarded.To determine the physical positions of each SNP, the variants were annotated using the ANNOVAR software (version 2) tool [21] to align and annotate SNPs.

Genome-Wide Association Analysis
The implementation of REGENIE [22] and SAIGE [23] was chosen due to the presence of rare variants and highly imbalanced case-control ratios.In an effort to enhance the power of the association test, we employed a gene-region-based approach (or other methods for generating variant sets) and utilized GMMAT (version 0.9.3) [24] software for conducting gene/set-based association analysis.We included three tests in the framework: the burden test, SKAT (sequence kernel association test) [25], and SKAT-O (optimal sequence kernel association test) [26].The SKAT method employed a multiple regression model to directly assess the relationship between the phenotypes and genetic variants in a variant set, as well as covariates.This approach allows for variations in both direction and magnitude of effects among different variants.In the burden test, SKATO estimated the weight and SKAT statistics that maximize power using the formula T SKAT−0 = (1 − ρ)T SKAT + ρT burden .Additionally, SKAT-O required a grid search over ρ to determine the minimum p-value.The Bonferroni correction threshold (p-value = 0.01/marker number or 0.05/marker number) was employed to identify statistically significant associations.Candidate genes (CAGs) located within a 50 kb region upstream or downstream of the significantly associated markers were identified.

Enrichment Analysis of Candidate Genes
The candidate genes were mapped to GO terms in the Gene Ontology database [27] (http://www.geneontology.org/,accessed on 20 April 2024).Gene numbers were calculated for each term, and significantly enriched GO terms were determined using a hypergeometric test.The p-value calculation formula is as follows: The number of all genes with GO annotation is represented by N, while n represents the number of genes/CAGs in N. Similarly, M denotes the number of all genes annotated to specific pathways, and m represents the number of genes/CAGs in M. The calculated p-value underwent false discovery rate (FDR) correction, with a threshold set at FDR ≤ 0.05.Pathways that meet this condition are defined as significantly enriched GO terms in genes/CAGs.

Group Validation of Significantly Correlated SNPs
The blood samples of 1125 sheep from a large group were collected in a random manner, and the horn phenotype data (Figure 1, Supplementary Table S11), including measurements of horn length, horn base circumference, and horn base interval, were recorded using the following methods: Animals 2024, 14, x FOR PEER REVIEW 5 of 19 reactions were performed in a volume of 5 µL including 0.5 µL SNaPshotMix, 3 µL Pooled PCR Products, 1 µL Pooled Primers, and 0.5 µL ddH2O with a program consisting of an initial denaturation at 95 °C for 30 s and 35 cycles comprising 95 °C for 5 s, 52 °C for 5 s, and 60 °C for 5 s.All SNaPshot ® -specific products were purified with Shrimp Alkaline Phosphatase-SAP (Thermo Fisher Scientific Inc., Waltham, MA USA) and incubated at 95 °C for 5 min.Finally, the Genemarker Software (v.1.90) was used to analyze all the SNaPshot ® -specific results.Allele frequency and genotype frequency were assessed using webbased software to determine their adherence to the Hardy-Weinberg equilibrium (http://scienceprimer.com/hardy-weinberg-equilibrium-calculator,accessed on 20 April 2024).
The association between various genotypes or haplotypes and horn traits was assessed using a sampling approach based on general linear models: Here: YijTl is the horned trait measured observation; ξ is the overall mean; αi is the genetic effect; βj is the individual age effect; γT is the year effect; δl is the permanent environment effect; θ is the random residual effect and is assumed to be independent, N (0, σ2) distribution.
The data were analyzed using the General Linear Model (GLM) in SPSS Statistics (version 29, IBM, Armonk, NY, USA), and statistical significance was determined at a pvalue of less than 0.05.Horn length: the arc length along the longitudinal edge of the back of the horn, from its base to its end; Horn base circumference: the circumference of the longitudinal base encircling the horn base; Horn base interval: the linear distance between the vertical bases on the posterior side of the horn base.
The nine SNPs detected in exons of genes identified by screening analysis were validated through Sanger sequencing and SNaPshot mini-sequencing, and primers were designed using Primer3 v0.4.0 (1) [28] while single-base extension primers were designed using the online tool Primer 3 plus (Supplemental Tables S1 and S2).
The amplification and screening of SNPs were conducted in 1125 sheep to analyze the correlation between candidate mutations and horn traits.Initially, the amplification of exons from genomic DNA was achieved using primers.PCR reactions were performed in 10 µL volumes containing 5 µL Taq PCR Mix, 1 µL Primer Mix, 1 µL DNA, and 3 µL ddH 2 O.A PCR program consisting of an initial denaturation at 94 • C for 5 min, followed by 35 cycles at 94 • C for 20 s, 60 • C for 30 s, and finally at 72 • C for 30 s was carried out using a Veriti™96-Well Thermal Cycler (Applied Bio-systems, Foster City, CA, USA).All PCR products were purified with ExoSAPIT ® (USB Corporation, Cleveland, OH, USA) at 37 • C for 40 min followed by incubation at 85 • C for 15 min.Subsequently, all SNaPshot ® reactions were performed in a volume of 5 µL including 0.5 µL SNaPshotMix, 3 µL Pooled PCR Products, 1 µL Pooled Primers, and 0.5 µL ddH 2 O with a program consisting of an initial denaturation at 95 • C for 30 s and 35 cycles comprising 95 • C for 5 s, 52 • C for 5 s, and 60 • C for 5 s.All SNaPshot ® -specific products were purified with Shrimp Alkaline Phosphatase-SAP (Thermo Fisher Scientific Inc., Waltham, MA USA) and incubated at 95 • C for 5 min.Finally, the Genemarker Software (v.1.90) was used to analyze all the SNaPshot ® -specific results.Allele frequency and genotype frequency were assessed using web-based software to determine their adherence to the Hardy-Weinberg equilibrium (http: //scienceprimer.com/hardy-weinberg-equilibrium-calculator,accessed on 20 April 2024).
The association between various genotypes or haplotypes and horn traits was assessed using a sampling approach based on general linear models: Here: Y ijTl is the horned trait measured observation; ξ is the overall mean; α i is the genetic effect; β j is the individual age effect; γ T is the year effect; δ l is the permanent environment effect; θ is the random residual effect and is assumed to be independent, N (0, σ2) distribution.The data were analyzed using the General Linear Model (GLM) in SPSS Statistics (version 29, IBM, Armonk, NY, USA), and statistical significance was determined at a p-value of less than 0.05.

Summary Statistics of Phenotype Data and Quality Control
We established a polled core population and observed polymorphisms in the horn phenotypes of G2 generation sheep through directional selection, including degenerated, deformed horns known as scurs, polled individuals, and individuals with normal large horns.The descriptive statistics of phenotypic values associated with horn-type traits in a population of 103 Tibetan sheep are provided in Supplementary Table S3.
After implementing quality control and management measures for 36,023,978 SNP loci, the following conditions were excluded: (1) 14,513,444 SNPs with a secondary gene frequency less than 0.05; (2) 200,718 SNPs with a deletion rate exceeding 0.5; and (3) 30,361 SNPs with a heterozygosity ratio greater than 0.8.Subsequently, a total of 21,279,455 independent SNPs remained for the subsequent analysis.The final count revealed a total of 33,065,497 single nucleotide polymorphisms (SNPs) and 2,958,481 insertions/deletions (indels).The majority of the high-quality single nucleotide polymorphisms (SNPs) (61.61%) were observed in intergenic regions, characterized by T/C and A/G substitutions, while only 0.70% were located within exonic regions.The remaining SNPs were situated upstream (0.61%) and downstream (0.63%) of the open reading frame, primarily found within intronic regions (33.97%) (Supplementary Table S4).

Population Structures and Genome-Wide Association Analysis
According to the population structure, admixture v1.3 was employed with a range of K values from 2 to 5 (Figure 2A).Specifically, a K value of 2 was utilized and the corresponding outcome is depicted in Figure 2B.Kinship estimation and principal component analysis (PCA) were conducted on all individuals, confirming the validity of the sampling (Figure 2C,D).The bivariables in the Kinship plot are perfectly aligned along a straight line, indicating a strong linear relationship.Moreover, the sample exhibits high repeatability, suggesting excellent data reproducibility.
According to the number of independent effective SNPs, the genome-wide significance and suggestive significance values were determined as a p-value of 4.06 × 10 −10 corresponding to a 1% significance level and a p-value of 2.34 × 10 −9 corresponding to a 5% significance level, respectively.SNPs with p-values lower than these thresholds were identified as being associated with Tibetan sheep horn traits.The p-value Manhattan map shows significant correlations between SNPs and horn length (Figure 3A), horn base circumference (Figure 3C), and horn base interval (Figure 3E), with statistically significant regions covering SNPs sites (Figure 3B,D,F).These aforementioned SNPs exhibit a strong positive correlation with the three horn-related traits (Table 1, Supplementary Tables S5-S7), surpassing the predetermined threshold.Further verification will be conducted subsequently.Concurrently, we conducted a comparison of the reported sites in the literature pertinent to RXFP2 and observed that only a single site (OAR10_29461968) aligned with our findings (Table 2).The genetic polymorphism of horn type was observed in both sexes, with three distinct morphological characteristics: large normal horns, scurs, and polled (polled accounted for 31.12% and 54.55% of all recorded cases in male and female sheep, respectively) (Table 3).Each horn type exhibited a significant number of quantitative variation sites.Differences in the dominance and expression of specific alleles were evident between the sexes, posing challenges to inferring an individual's genotype.The genetic polymorphism of horn type was observed in both sexes, with three distinct morphological characteristics: large normal horns, scurs, and polled (polled accounted for 31.12% and 54.55% of all recorded cases in male and female sheep, respectively) (Table 3).Each horn type exhibited a significant number of quantitative variation sites.Differences in the dominance and expression of specific alleles were evident between the sexes, posing challenges to inferring an individual's genotype.

Enrichment Analysis
To further elucidate candidate genes associated with significant SNPs, we conducted additional gene enrichment analysis.The pathways enriched in the RXFP2 gene were further subjected to additional recordings.The GO/KEGG pathway analysis revealed the presence of certain pathways associated with horn characteristics.Specifically, the relaxin signaling pathway and neuroactive ligand-receptor interaction were found to be significantly KEGG co-enriched in terms of horn length, horn base circumference, and horn base interval.Furthermore, these three traits collectively exhibited enrichment in 33 GO term pathways.The candidate genes associated with horn length (Figure 4A, Supplementary Table S8), horn base circumference (Figure 4B, Supplementary Table S9), and horn base interval (Figure 4C, Supplementary Table S10) exhibit co-enrichment during gonad development, oocyte maturation, sex differentiation, protein-hormone receptor activity, cell maturation and differentiation processes, as well as positive regulation of CAMP-mediated signaling pathways.These findings suggest that these three traits are likely regulated by shared biological pathways.The enrichment analysis revealed that the RXFP2 gene exhibited the highest level of enrichment among the Gene Ontology terms associated with biological processes.

Population Verification of SNPs Significantly Correlated with Tibetan Sheep Horn Traits
The SNPs identified by the GWAS mentioned above were subjected to significant testing between horned and polled Tibetan sheep, revealing a significant correlation between the RXFP2 gene and polled Tibetan sheep, which has been extensively reported in the literature.Furthermore, the descriptive statistical information of phenotypic values obtained from a large population through the random selection of 1125 ewes (Supplementary Table S11) provided population verification for these polymorphisms at the aforementioned RXFP2 sites.Genotyping, population genetic analysis, and association determination between SNPs and horn traits were conducted for all loci.
maturation and differentiation processes, as well as positive regulation of CAMP-mediated signaling pathways.These findings suggest that these three traits are likely regulated by shared biological pathways.The enrichment analysis revealed that the RXFP2 gene exhibited the highest level of enrichment among the Gene Ontology terms associated with biological processes.

Population Verification of SNPs Significantly Correlated with Tibetan Sheep Horn Traits
The SNPs identified by the GWAS mentioned above were subjected to significant testing between horned and polled Tibetan sheep, revealing a significant correlation between the RXFP2 gene and polled Tibetan sheep, which has been extensively reported in the literature.Furthermore, the descriptive statistical information of phenotypic values obtained from a large population through the random selection of 1125 ewes (Supplementary Table S11) provided population verification for these polymorphisms at the aforementioned RXFP2 sites.Genotyping, population genetic analysis, and association determination between SNPs and horn traits were conducted for all loci.GWAS predicted six individually detected SNPs, while three additional mutations were identified in close proximity to the target SNPs upstream of amplification products (Supplementary Table S12).The value of PIC at SNP001 was <0.25, which was indicative of low polymorphism.The other SNPs we identified exhibited moderate levels of polymorphism (0.25 < PIC < 0.5) (Supplementary Table S13).
Nine SNP verification results revealed significant correlations between three mutation sites and the horn length, horn base circumference, and horn base interval of Tibetan sheep (Table 4).The genotype of horns had a significant effect on horn length, horn base circumference, and horn base interval, while there might be a polygenic influence associated with the genotypic horn.

Linkage Disequilibrium and Haplotype Block Association Analysis
The analysis results of linkage disequilibrium (LD) indicate a strong correlation between the SNP sites of the two haplotype modules (D' > 0.85), with a coefficient of 0.81 observed for the correlation between the two haplotype modules (Figure 5A).Subsequently, a correlation analysis was conducted to examine the haplotypes of the two modules, revealing significant effects of different haplotypes on the horn type of Tibetan sheep (Tables 5 and 6).In the 19 kb LD block, twelve haplotypes containing three SNPs were identified, while in the 51 kb LD block, twenty-two haplotypes containing six SNPs were found (Figure 5B,C).In Module 1, three haplotypes, namely AAGGCT (n = 11, p-value < 0.01), GAGACT (n = 9, p-value < 0.01), and GAGGTT (n = 8, p-value < 0.01) (Table 5), exhibited significant correlations with polled traits.In Module 2, the haplotypes AACCAAGGCTCT (n = 2, p-value < 0.01), AATTGAAATTTT (n = 18, p-value < 0.01), AATTGGAATTTT (n = 3, p-value < 0.01) GACTGAGATTCT (n = 4, p-value < 0.01), and GATTGGAATTTT (n = 15, p-value < 0.01) (Table 6) were found to be significantly associated with polled traits.The maintenance mechanism of horn characteristic variation is determined by the co-breeding of haplotypes associated with these three traits.

Discussion
The vestigialization of horns occurs in most domestic lines due to the diminishing relevance of traits that ensure fitness in natural environments under artificial breeding conditions [30].Domestic animals have undergone domestication and improvement to produce similar phenotypic transformations [34], a phenomenon known as "domestication syndrome" [35][36][37].However, our findings indicate that artificial intervention has led to the differentiation of three distinct horn phenotypes (normal large horns, scur, and polled) in polled-orientated progeny.Furthermore, we have determined that all types of horns, including horn length, horn base circumference, and horn base interval, are regulated by the RXFP2 gene located on chromosome 10 [33,38].We can distinguish phenotypes that have the same horn but have different genotypes, suggesting that the genotype of the horn may be responsible not only for discrete variation in the diagonal type, but also for many quantitative genetic variations in the diagonal size [39], and that the quantified variation differences in the size of the horn appear to be caused by the genotype of the horn.This study establishes a foundation for investigating genotypic-level selection of horn types and morphology and enhances our understanding of the mechanisms underlying the maintenance and differentiation of horn polymorphisms in this population.
Interestingly, horn polymorphisms also occur in natural populations, such as Soay sheep (males exhibit both normal and scurred horns, while females can have normal, scurred, or polled horns) [39].The same phenomenon is observed in the natural population of Tibetan sheep.It remains unclear whether this phenotype is associated with genetic introgressions or a free combination of RXFP2 superior alleles linked to horn types within the population.This hypothesis can be tested in future studies by generating genotypic data for additional sheep breeds.
The development of horns in sheep is influenced not only by genetic and sexual factors, but also by hormonal and environmental influences.Seasonal horn growth typically occurs during spring, when prolactin secretion reaches its peak.There are notable variations in the timing and pace of horn development across different life stages [40].Sustained low levels of testosterone during the inactive period are essential for promoting male-type horn growth, while high levels of testosterone prior to estrus can actually impede growth-a mechanism that also regulates antler development in deer [41,42].Sheep exhibited cryptorchidism after the RXFP2 gene was knocked out using CRISPR/Cas9 technology [43].The growth rate of sheep horns is regulated by thyroid activity, with the peak concentration of thyroid hormone (T3, T4) coinciding with the rapid growth period [44].Our findings demonstrated a significant correlation between the concentration of RXFP2 and gonadal development and maturation, as well as the regulatory impact of reproductive hormones on sheep horn development.Further investigation into the hormonal regulatory mechanisms will be conducted.
The inheritance of horns in sheep is closely associated with sex, as indicated by Wood's findings [45].Specifically, the horned and hornless genes exhibit distinct effects on males and females.In this regard, it can be observed that large horns are dominant in males while recessive in females [45].A study conducted by Dominik et al. revealed the presence of a single nucleotide polymorphism located in a specific genomic region, which exhibited a prominent maternal imprinting effect in both males and females [12].QTL analysis of diagonal morphology revealed a significant interaction between horn shape traits and sex, indicating that the same chromosomal region accounted for genetic variation in both horn length and base circumference angle.Furthermore, it was determined that these two traits were subject to genotypic (synergistic) selection [38].The fixation of the ideal individual's characteristics can be rapidly achieved by mating a homozygous hornless male with a homozygous hornless female [45], resulting in all offspring being hornless.
Prediction within a single family will incorporate both linkage information and LD information, assuming that non-genetic factors have been accurately accounted for.However, in the case of cross-family prediction, reliance can only be placed on LD information [46].Haplotypes encompass a greater amount of information (heterozygosity) compared to a single SNP, thereby significantly enhancing accuracy and robustness in gene localization.Conversely, in small sample sizes, a single SNP fails to establish an association between a specific phenotype and gene, whereas haplotypes possess this capability [47].Horn phenotypes can be accurately predicted with an accuracy of approximately 0.7 when utilizing a SNP for prediction [33].Consequently, investigations into phenotypic correlation with genetic variation are more likely to yield successful outcomes if the SNPs employed have been demonstrated to exhibit linkage imbalance through methodologies such as haplotype analysis [48][49][50].The alleles located on OAR10 at position 29,461,968 were found to be genetically linked with the neighboring loci in the gene map [14,30].LD was observed to be lower between the 1.78 kb insert fragment of the RXFP2 gene and the SNP located at OAR10_29511510.1 (position 29,476,678) [16], whereas a higher LD was observed between the SNP inserted with OAR10_29546872.1)[33].In our study, the haplotype combinations of AAGGCT, GAGACT, and GAGGTT as well as AACCAAGGCTCT, AATTGAAATTTT, AATTGGAATTTT, GACTGAGATTCT, and GATTGGAATTTT are the predominant ones.The presence of fitness disparities may also occur among sheep exhibiting similar phenotypes but possessing different underlying genotypes due to their association with other alleles at adjacent sites in a gene-linkage imbalance.We have successfully identified the haplotype associated with the horn trait, indicating that the absence of horns in artificially bred sheep is influenced by the interaction of multiple SNPs.This finding also highlights the synergistic genetic basis underlying these three horn traits.It is important to note that single SNP loci alone cannot provide sufficient evidence for establishing a correlation between variations in the RXFP2 gene and the expression of the horn trait.
Polled rams exhibit relatively low breeding values and were predominantly heterozygous, the proportion of polled individuals has greater potential in the G2 population, where ewes outnumber rams.It is evident that there exists a common genetic basis for both horn type and length, with at least 76% of the genetic variation attributed to the "horn" genotype [14].Consequently, genome selection can significantly accelerate the rate of genetic improvement and positively impact both the breeding value and quantity of polled rams.Quantitative trait genes in animals were directly selected using molecular marker technology, and whole-genome association analysis based on the genome variation map was employed to enhance breed improvement more effectively.A breeding program thus may well produce polled sire with breeding values comparable to the top sire in the near future [51].Based on these findings, we anticipate that a majority of lambs born in generations G3-G4 will be polled, thus establishing a new prototype for a polled strain.

Conclusions
In summary, our GWAS identified six single nucleotide polymorphisms (SNPs) on the RXFP2 gene located on chromosome 10 that exhibited significant associations with horned traits.Each horn type displayed a substantial number of variation sites.Furthermore, three additional mutations were subsequently detected near the target SNP upstream of the amplification product.Population validation of the specific haplotype is considered indicative of distinguishing between horn and polled phenotypes while maintaining horn length, horn base circumference, and horn base interval.This study aims to provide novel insights into the molecular characteristics of polled sheep under human intervention and can offer practical guidance.

Figure 1 .
Figure 1.Measurement of horn length, horn base circumference, and horn base interval.The X axis represents the measured position, while the Y axis represents the measured value, which is depicted through a combination of a box plot and normal curve.

Figure 1 .
Figure 1.Measurement of horn length, horn base circumference, and horn base interval.The X axis represents the measured position, while the Y axis represents the measured value, which is depicted through a combination of a box plot and normal curve.

Figure 2 .
Figure 2. Population structure analysis in the nucleus herd of polled Tibetan sheep.(A) Population structure with K from 2 to 5; (B) cross-validation plot for determining the best K; (C) kinship plot of 103 Tibetan sheep; (D) principal component analysis.

Figure 2 .
Figure 2. Population structure analysis in the nucleus herd of polled Tibetan sheep.(A) Population structure with K from 2 to 5; (B) cross-validation plot for determining the best K; (C) kinship plot of 103 Tibetan sheep; (D) principal component analysis.

Figure 3 .
Figure 3. Manhattan plots and QQ plots of (A) horn length, (C) horn base circumference, and (E) horn base interval traits.(B) QQ plots of the horn length, (D) QQ plots of the horn base circumference, (F) QQ plots of the horn base interval.

Figure 3 .
Figure 3. Manhattan plots and QQ plots of (A) horn length, (C) horn base circumference, and (E) horn base interval traits.(B) QQ plots of the horn length, (D) QQ plots of the horn base circumference, (F) QQ plots of the horn base interval.

Figure 4 .
Figure 4. Gene Ontology (GO) and KEGG pathway analysis of the RXFP2 gene in terms of (A) horn length, (B) horn base circumference, and (C) horn base interval.

Figure 4 .
Figure 4. Gene Ontology (GO) and KEGG pathway analysis of the RXFP2 gene in terms of (A) horn length, (B) horn base circumference, and (C) horn base interval.

Figure 5 .
Figure 5. Linkage disequilibrium plot (r 2 ) and haplotype blocks for SNPs of the RXFP2 gene.(A) A 19-kb LD block and 51-kb LD block; (B) four haplotypes within the sampled sheep population; (C) five haplotypes within the sampled sheep population.

Figure 5 .
Figure 5. Linkage disequilibrium plot (r 2 ) and haplotype blocks for SNPs of the RXFP2 gene.(A) A 19-kb LD block and 51-kb LD block; (B) four haplotypes within the sampled sheep population; (C) five haplotypes within the sampled sheep population.

Table 1 .
Analysis of related SNPs associated with horn traits.

Table 2 .
Analysis of RXFP2 site characteristics.
a-this work.

Table 3 .
Phenotypic distribution and potential genotypes of horn type in genetically matched polled offspring.

Table 4 .
Association analyses of SNPs and genotypes in RXFP2 with horn traits.Within the same line, the mean values of different superscript lowercase letters were significantly different (p < 0.05).Within the same line, the mean values of different superscript uppercase letters differ significantly (p < 0.01).

Table 5 .
Association analysis of haplotypes in RXFP2 with horn traits in Module 1.Within the same line, the mean values of different superscript lowercase letters were significantly different (p < 0.05).Within the same line, the mean values of different superscript uppercase letters differ significantly (p < 0.01).

Table 6 .
Association analysis of haplotypes in RXFP2 with horn traits in Module 2.Within the same line, the mean values of different superscript lowercase letters were significantly different (p < 0.05).Within the same line, the mean values of different superscript uppercase letters differ significantly (p < 0.01).