Novel Heredity Basis of the Four-Horn Phenotype in Sheep Using Genome-Wide Sequence Data

Simple Summary We performed a genome-wide association study of 100 sheep genomes to investigate the hereditary basis of the four-horn phenotype (FHP) in sheep. The results revealed that a total of 3 InDels (e.g., CHR2: g.133,742,709delA and g.133,743,215insC) and 14 SNPs (e.g., CHR2: g.133,727,513C > T and CHR16: g.40,351,378G > A) in CHR2 and CHR16 were significantly associated with the four-horn trait. Furthermore, HOXD1 and ADAMTS12 annotated from both CHR sequences were the key candidate gene loci contributing to FHP in sheep. This study provided a novel understanding of the Mendelian genetic basis for four horns in sheep. Abstract Horns are an important breeding trait for sheep. However, no widely recognized viewpoint on the regulatory genes and mechanisms of horns is available, and the genetic basis of the four-horn phenotype (FHP) is unclear. This work conducted a genome-wide association study with 100 sheep genomes from multiple breeds to investigate the genetic basis of the FHP. The results revealed three significant associations (corrected as p < 1.64 × 10−8) of the InDels (CHR2: g.133,742,709delA, g.133,743,215insC, and g.133,743,940delT) for FHP in the intergenic sequence (IGS) between the MTX2 and the LOC105609047 of CHR2. Moreover, 14 significant associations (corrected as p < 1.42 × 10−9) of SNPs with the FHP phenotype were identified in CHR2 and CHR16, including five (e.g., CHR16: g.40,351,378G > A and g.40,352,577G > A) located in the intron of the ADAMTS12 gene, eight (e.g., CHR2: g.133,727,513C > T and g.133,732,145T > G) in the IGS between MTX2 and LOC105609047, and only one (CHR2: g.133,930,761A > G) in the IGS between HOXD1 and MTX2. Obvious divergence was also observed in genotype patterns between the FHP and others (two horns and hornless) in the HOXD1 and ADAMTS12 gene regions. An extremely significant linkage also occurred between Loci I and Loci II within 100 individuals (LD = −156.02186, p < 0.00001). In summary, our study indicated that the genomic sequences from CHR2 and CHR16 contributed to the FHP in sheep, specifically the key candidate genes HOXD1 and ADAMTS12. These results improved our understanding of the Mendelian genetic basis of the FHP in sheep.


Introduction
Sheep are important domestic animals that play an important role in regional economic development and meat supply security [1].To date, most studies on sheep have focused on the phylogenetic history of the population [2], the genetic diversity of indigenous sheep [3], and the genetic basis of important economic traits [4,5].
Horns are extensively considered to be an important Mendelian trait in many ruminants, and hornless herbivores accord better with the requirements of current intensive animal husbandry management [6].Numerous studies have increased our understanding Animals 2023, 13, 3166 2 of 8 of the formation basis of horn traits and the genetic basis of hornless phenotypes in ruminants [7,8].For example, RXFP2 expression is reportedly upregulated significantly in horn buds, indicating that RXFP2 may be a key gene in the regulation of horn development [9].Moreover, a 1.8 kb insertion in the 3 -UTR of RXFP2 contributes to the hornless phenotype in multiple domestic sheep breeds [10].
The four-horn phenotype (FHP) of sheep is a major research topic.The genetic locus of the FHP is located on chromosome (CHR) 2 [11,12].The genomic markers associated with this phenotype have also been narrowed down to the HOXD gene family member genome regions of CHR2 in Sishui fur sheep [13].However, previous works lack precise localization and overlook genetic background interference because of limited genetic markers (Illumina OvineSNP50 BeadChip) or a single population with four horns [11,13].
In the present study, an association analysis of multiple breeds and large samples was performed to identify FHP-associated genes with high-depth whole-genome sequencing.Our results can help further elucidate the Mendelian genetic basis of the FHP in sheep.

Materials and Methods
A total of 88 publicly available sheep genome datasets were downloaded from the NCBI SRA Database (PRJNA624020, Table S1), including 10 four-horn sheep (Sishui fur sheep), 26 two-horn sheep (e.g., small-tailed Han sheep), and 52 polled sheep (e.g., Hu sheep and large-tailed Han sheep).The blood samples of 12 Mongolian sheep with four horns were supported by the Institute of Animal Sciences of CAAS.Genomic DNA was extracted following a standard phenol-chloroform protocol (Sambrook and Russel, 2001).Sequencing libraries were constructed with Annoroad ® Universal DNA Library Prep Kit v2.0 (Illumina ® , 15,026,486 Rev. C, San Diego, CA, USA).Sequencing was performed using the BGI DNBSEQ-T7 platform (Beijing Genomics Institution, Shenzhen, China) with > 10× per individual.

Results and Discussion
An approximately 310.79 Gb dataset was generated from 12 sequenced individuals and uploaded to the Sequence Read Archives (PRJNA1003686).The sequencing depth of the 12 samples ranged from 10.43-to 18.04-fold (Table S2).First, a total of 3,042,046 InDels were obtained after filtering from all 100 individuals, and the GWAS results revealed that three InDels (CHR2: g.133,742,709delA, g.133,743,215insC, and g.133,743,940delT) exceeded the Bonferroni significance threshold (corrected as p < 1.64 × 10 −8 ; Figure 1B).Based on the physical location of the chromosomes, the three InDels were distributed in the intergenic sequence (IGS) between MTX2 and LOC105609047.In fact, MTX2 was confirmed to be related to bone development.For example, the splice site variant of the MTX2 gene (c.543+ 1G > T) can cause human mandibuloacral dysplasia, which is characterized by postnatal growth retardation, hypotonia, and dysmorphic facial features [17,18].
Obviously, the SNP (CHR2: g.133,930,761A > G) with the highest associated signal value was located in the region between HOXD1 and MTX2 on CHR2, consistent with previous findings [11,12].Near the threshold line, a series of high-signal SNPs was annotated to the introns and 3 -UTR of HOXD1.The genotype patterns of HOXD1 showed obvious differences between four-horn breeds and others (Figure 1F).HOXD1 haploinsufficiency leads to the division of horn sprout progeny [19].In particular, a number of studies have shown that the HOXD gene family of CHR2 may be related to the FHP [11,12,20].In fact, HOXD family genes have been identified as playing an important role in vertebrate limb development.For example, the targeted disruption of exon 2 in HOXD10 leads to changes in the hind limb, such as an anterior shift in the position of the patella.These changes contribute to defects in mouse locomotion [21].HOXD11 is highly expressed in the carpal, tarsal, and digital regions during limb development in mammals [22].HOXD11 has also been demonstrated to regulate chondrocyte differentiation at an early step in the zeugopods of mice [23].An SNP in HOXD12 (c.1359G > T) induces a change from alanine to serine and thus causes an abnormal limb phenotype in mice with macrodactyly [24].Severe brachydactyly with metacarpal-to-carpal transformation is also underlined by a missense mutation (c.938C > G) in HOXD13 [25].
Additionally, several SNPs located on CHR 16 with the highest associated signal values were annotated to the ADAMTS12 gene.Throughout the whole-genotype patterns, distinct differences existed in two regions of ADAMTS12 between four-horn breeds and other breeds (e.g., CHR16: 40,350,835-40,352,878 and 40,361,056-40,366,451; Figure 1G).The expression levels of ADAMTS12 and ADAMTS16 were reported to be related to reproductive capabilities in males, especially testicular development [26,27].Interestingly, plasma testosterone concentration was negatively correlated with horn growth, indicating that peripheral plasma testosterone levels may play a role in regulating horn development [28].ADAMTS12 also plays an important role in chondrogenesis regulation [29].The role of ADAMTS12 in chondrogenesis was further found to be related to the chondrogenic regulatory factor parathyroid hormone-associated peptide (PTHrP) [30].Other evidence shows that ADAMTS12 and PTHrP can inhibit chondrocyte differentiation by inhibiting chondrocyte differentiation marker genes, such as Col II and Sox9 [29].The proliferation and differentiation of chondrocytes are extremely important to antlers, similar to horns [31,32].
three InDels (CHR2: g.133,742,709delA, g.133,743,215insC, and g.133,743,940delT) exceeded the Bonferroni significance threshold (corrected as p < 1.64 × 10 −8 ; Figure 1B).Based on the physical location of the chromosomes, the three InDels were distributed in the intergenic sequence (IGS) between MTX2 and LOC105609047.In fact, MTX2 was confirmed to be related to bone development.For example, the splice site variant of the MTX2 gene (c.543+ 1G > T) can cause human mandibuloacral dysplasia, which is characterized by postnatal growth retardation, hypotonia, and dysmorphic facial features [17,18].Coincidentally, our analysis showed that ADAMTS12 was enriched in several GO terms (Table S4), such as the negative regulation of chondrocyte differentiation, extracellular matrix, extracellular matrix organization, cellular response to bone morphogenetic protein (BMP) stimulus, ossification, and the ossification involved in bone maturation.The extracellular matrix provides a living environment for chondrocytes to facilitate their adhesion, proliferation, and secretion of cartilage matrix, thereby affecting cartilage growth [33].BMP also plays an important role in the initiation of cartilage formation [34].For example, Animals 2023, 13, 3166 5 of 8 the loss of BMP2 and BMP4 leads to serious damage to bone formation [35].Therefore, the BMP signaling pathway is the main regulator of bone development [36].Endochondral ossification, a process of formation for many bones, is essential for the growth and development of bones [37].Interestingly, ossification was already reported to be an important process in the growth of horns and antlers [38,39].
Four haplotypes were obtained from the significant SNPs in CHR2 (Loci I, Table S5) and CHR16 (Loci II, Table S6) in this study.The Hap2 and Hap3 of Loci I occupied the highest frequency (88.64%) in the FHP animals (Figure 2A), and the heterozygous individuals (Hap2/Hap3) accounted for 72.73%.Moreover, Hap1 in Loci II accounted for the highest proportion of the FHP (Figure 2B), and the homozygous individuals (Hap1/Hap1) accounted for 81.82%.Extremely significant linkage was further observed between Loci I and Loci II within 100 individuals, although they originated from different chromosomes (LD = −156.02186,p < 0.00001).Eight haplotypes were reconstructed from Loci I and Loci II (Table S7), with Hap2 (43.18%) and Hap3 (45.45%) having the highest frequency in the FHP (Figure 2C).
growth [33].BMP also plays an important role in the initiation of cartilage formation [34].For example, the loss of BMP2 and BMP4 leads to serious damage to bone formation [35].Therefore, the BMP signaling pathway is the main regulator of bone development [36].Endochondral ossification, a process of formation for many bones, is essential for the growth and development of bones [37].Interestingly, ossification was already reported to be an important process in the growth of horns and antlers [38,39].
Four haplotypes were obtained from the significant SNPs in CHR2 (Loci I, Table S5) and CHR16 (Loci II, Table S6) in this study.The Hap2 and Hap3 of Loci I occupied the highest frequency (88.64%) in the FHP animals (Figure 2A), and the heterozygous individuals (Hap2/Hap3) accounted for 72.73%.Moreover, Hap1 in Loci II accounted for the highest proportion of the FHP (Figure 2B), and the homozygous individuals (Hap1/Hap1) accounted for 81.82%.Extremely significant linkage was further observed between Loci I and Loci II within 100 individuals, although they originated from different chromosomes (LD = −156.02186,p < 0.00001).Eight haplotypes were reconstructed from Loci I and Loci II (Table S7), with Hap2 (43.18%) and Hap3 (45.45%) having the highest frequency in the FHP (Figure 2C).
On the basis of the distributed frequency of eight haplotypes, the FST result confirmed negligible genetic differentiation between both FHP breeds, whereas significant genetic divergence was observed between FHP breeds and others (Figure 2D and Table S9).Previous studies have hypothesized that horn inheritance is controlled by two loci, with four horns occurring if the recessive allele is homozygous at either locus [40].The Loci I and Loci II from CHR2 and CHR16 identified in this study supported the findings of a previous one.Moreover, although the series of SNPs did not reach a significantly associated threshold, they still had high-correlation signals distributed throughout the intergenic region between ARNT2 and ABHD17C and the HOMER2 gene in CHR18.However, studies related to ARNT2 and ABHD17C have focused on the fields of cancer and immunity [41,42].Other researchers have shown that HOMER2 is expressed in osteoblasts [43,44].HOMER2 and HOMER3 can regulate bone metabolism [45].In particular, the FoxO signaling pathway, enriched with HOMER2, can promote chondrocyte maturation and extracellular matrix generation [46].Akt, also known as protein kinase B, controls the processes of On the basis of the distributed frequency of eight haplotypes, the F ST result confirmed negligible genetic differentiation between both FHP breeds, whereas significant genetic divergence was observed between FHP breeds and others (Figure 2D and Table S9).Previous studies have hypothesized that horn inheritance is controlled by two loci, with four horns occurring if the recessive allele is homozygous at either locus [40].The Loci I and Loci II from CHR2 and CHR16 identified in this study supported the findings of a previous one.
Moreover, although the series of SNPs did not reach a significantly associated threshold, they still had high-correlation signals distributed throughout the intergenic region between ARNT2 and ABHD17C and the HOMER2 gene in CHR18.However, studies related to ARNT2 and ABHD17C have focused on the fields of cancer and immunity [41,42].Other researchers have shown that HOMER2 is expressed in osteoblasts [43,44].HOMER2 and HOMER3 can regulate bone metabolism [45].In particular, the FoxO signaling pathway, enriched with HOMER2, can promote chondrocyte maturation and extracellular matrix generation [46].Akt, also known as protein kinase B, controls the processes of endochondral ossification and skeletal growth by regulating chondrocytes and cartilage matrix by tuning the activities of FoxOs [46,47].

Figure 2 .
Figure 2. Haplotype frequency of candidate genomic region (CGR) in CHR2 and CHR16 and pairwise difference in different horn phenotype populations with CGR haplotype data.(A) Haplotype frequency of Loci I, (B) haplotype frequency of Loci II, (C) haplotype frequency of Loci I and Loci II, (D) heatmap of FST for Loci I and Loci II among 12 breeds.

Figure 2 .
Figure 2. Haplotype frequency of candidate genomic region (CGR) in CHR2 and CHR16 and pairwise difference in different horn phenotype populations with CGR haplotype data.(A) Haplotype frequency of Loci I, (B) haplotype frequency of Loci II, (C) haplotype frequency of Loci I and Loci II, (D) heatmap of F ST for Loci I and Loci II among 12 breeds.