Skip to main content
Download PDF

Genome-wide association study provided insights into the polled phenotype and polled intersex syndrome (PIS) in goats

BMC Genomics volume 25, Article number: 661 (2024) Cite this article

Abstract

Background

Breeding polled goats is a welfare-friendly approach for horn removal in comparison to invasive methods. To gain a comprehensive understanding of the genetic basis underlying polledness in goats, we conducted whole-genome sequencing of 106 Xinong Saanen dairy goats, including 33 horned individuals, 70 polled individuals, and 3 polled intersexuality syndrome (PIS) individuals.

Methods

The present study employed a genome-wide association study (GWAS) and linkage disequilibrium (LD) analysis to precisely map the genetic locus underlying the polled phenotype in goats.

Results

The analysis conducted in our study revealed a total of 320 genome-wide significant single nucleotide polymorphisms (SNPs) associated with the horned/polled phenotype in goats. These SNPs exhibited two distinct peaks on chromosome 1, spanning from 128,817,052 to 133,005,441 bp and from 150,336,143 to 150,808,639 bp. The present study identified three genome-wide significant SNPs, namely Chr1:129789816, Chr1:129791507, and Chr1:129791577, as potential markers of PIS-affected goats. The results of our LD analysis suggested a potential association between MRPS22 and infertile intersex individuals, as well as a potential association between ERG and the polled trait in goats.

Conclusion

We have successfully identified three marker SNPs closely linked to PIS, as well as several candidate genes associated with the polled trait in goats. These results may contribute to the development of SNP chips for early prediction of PIS in goats, thereby facilitating breeding programs aimed at producing fertile herds with polled traits.

Peer Review reports

Introduction

Horns are cranial appendages unique to ruminant animals [1]. The presence of horns is crucial for self-defense [2, 3] and serves as a prominent sexual trait in intra-male competition [4]. However, the utilization of horns in inter-animal fights poses risks to the herd’s welfare and complicates management and transportation processes, consequently leading to increased economic burdens [5]. Therefore, exfoliation is increasingly gaining popularity in modern animal husbandry. Undoubtedly, breeding polled animals offers two primary advantages. Firstly, it enhances animal welfare by minimizing injuries resulting from inter-animal fighting and avoiding the significant distress caused by conventional dehorning methods during the early stages [6]. Secondly, it protects animal caretakers, particularly farmers and all people in contact with animals [7]. Although the polled trait follows a dominant inheritance pattern, selective breeding for polledness poses challenges due to the risk of inbreeding and loss of genetic variance, as polled animals do not have the genetic level achieved as horned populations [8]. With the development of science and technology, the application of modern animal breeding techniques, such as high-density genotyping and genome resequencing, will facilitate the resolution of this biological issue.

In recent years, extensive investigations have been conducted on the genetic basis underlying the horned and polled phenotypes in cows, sheep, and goats. It has been shown that the molecular genetic causes for polled phenotypes vary due to distinct genetic backgrounds across different species. In cows, the polled individuals are inherited as an autosomal monogenic dominant trait, where the P allele (hornless) exhibits dominance over the p allele (horned) [9]. Although the polled locus has been reported to be located on the proximal end of bovine chromosome 1 [10], the causal genetic variants linked to the polled phenotypes are distinct in cows from diverse breeds. For instance, a total of four distinct genetic variants have been identified to be linked to the polled phenotypes in different European breeds [11, 12]. In sheep, the researchers found that the insertion of a nucleotide sequence of about 1.78 kb in chromosome 10 was the factor leading to the emergence of hornless sheep [13, 14]. However, the insertion alone does not fully account for the polled phenotype observed in diverse sheep breeds [15]. In goats, the inheritance pattern of polled phenotypes resembles the monogenic autosomal dominant mode observed in cows. The absence of horn growth in goats is attributed to a complex structural genetic variant, characterized by the fusion of a large duplicated segment measuring 480 kb, located at position 150,334,286–150,818,098 bp on chromosome 1, into the deleted region spanning approximately 11.7 kb situated around 129 Mb on chromosome 1 [16, 17]. However, the specific elements of these intricate structural variations that contribute to the absence of horns remain indeterminate.

Polled intersexuality syndrome (PIS) goats are occasionally observed within polled flocks. The dominant inheritance of the polled phenotype is mainly linked to recessive intersexuality in goats [18]. Homozygous polled females typically exhibit infertility as intersex individuals, while homozygous polled males do not display any signs of intersexuality. Previous studies have demonstrated that a deletion spanning 11.7 Kb, located at approximately 129 Mb on chromosome 1, is postulated to underlie the etiology of PIS-affected Saanen goats [16]. The most significant observation is that a duplicated sequence of approximately 480 Kb, which was originally located at the genomic position of around 150.5 Mb on chromosome 1, underwent an inverse insertion event into the PIS locus [9, 16]. Although phenotypes are believed to arise from gene copy number amplification, PIS often cannot be fully explained by alterations in gene dosage alone [19].

Compared to low-density SNP microarrays, the utilization of whole genome sequencing (WGS) data enables precise localization of association signals in genome-wide association studies (GWAS), thereby significantly enhancing the likelihood of identifying candidate genes responsible for causal mutations affecting complex traits in animal and plant [20,21,22]. In this work, we conducted a GWAS to identify the genetic variants and key candidate genes associated with the polled phenotype in Xinong Saanen dairy goats. Our findings have unveiled that two genomic regions located on chromosome 1 play a pivotal role in determining the polled phenotypes observed in goats. The identification of key candidate genes and SNPs associated with PIS-affected goats enables the development of SNP chips for early prediction of infertility intersexuality, thereby facilitating advancements in the field.

Methods

Animals

In this study, a total of 106 Xinong Saanen dairy goats were randomly selected from the experimental farm situated at Northwest A&F University in Yangling, Shaanxi Province, China. Informed consent was obtained from the owner. The samples comprised 33 horned individuals, 70 polled individuals, and 3 PIS individuals (Table S1). The goats utilized in this study were born from 2019 to 2022. The animals were uniformly managed and provided with a mixed diet comprising corn, soybean meal, bran, rapeseed meal, and a mineral-vitamin premix. Additionally, the reproductive strategy employed was consistent. The types of horn phenotypes have been well documented. Three PIS goats were identified based on their genital phenotype (Figure S1). The experimental protocol, including blood sampling and PIS goat identification, was approved by the Animal Ethical and Welfare Committee of the College of Animal Science and Technology, Northwest A&F University, Yangling, China (protocol number DK2022008).

Samples collection, DNA extraction, and whole-genome sequencing

Blood samples were obtained via jugular venipuncture from selected goats. The DNA sample was extracted from blood and quantified using established experimental protocols to ensure accuracy and precision. Subsequently, the DNA samples were randomly fragmented into 350 bp fragments, followed by library construction through a series of sequential steps including end repair, polyA tailing, purification, PCR amplification, and so forth. Whole-genomic sequencing was conducted using the Illumina NovaSeq 6000 platform (Beijing, China).

Short-read alignment and variant calling

The Trimmomatic software (version 0.38) was employed for the removal of adaptor sequences and elimination of low-quality reads [23]. The high-quality reads were then aligned to the goat reference genome (Accession NO, GCA_001704415.1) using BWA software (version 0.7.17) with default settings [24]. The Picard software (v2.10.6) (http://broadinstitute.github.io/picard/) was utilized for the removal of duplicate reads and subsequently sort the bam files. The HaplotypeCaller tool in GATK (v3.8-1-0) was employed with default parameters to detect variants for each accession [25]. Subsequently, the CombineGVCFs and GenotypeGVCFs tools in GATK were utilized to call population-level variants.

Quality control and annotation for genotype data

To obtain high-quality variants, the raw variant calls (SNPs and Indels) were filtered based on quality and depth using VCFtools software with the following parameters: ‘--minQ 100 --minDP 3 --maxDP 100 --min-meanDP 3’ [26]. The biallelic SNPs used for subsequent analyses were obtained by applying a variant site filter using the PLINK software with the following parameters: ‘--geno 0.2 (excluding genotypes with more than 20% missing data), --maf 0.05 (removing variants with a minor allele frequency less than 0.05), and --biallelic-only’ [27]. In addition, the identified SNPs were annotated using the ANNOVAR software [28].

Population structure

In order to understand population structure and help control the false positive results of GWAS, we performed a population genetics analysis. The biallelic SNPs were filtered based on linkage disequilibrium (LD) and retained unlinked SNPS for population structure analysis and principal component analysis (PCA) using the PLINK software with the following parameters: ‘--indep-pairwise 50 10 0.2 (filtering out SNPs with a correlation coefficient greater than 0.2 in sliding windows of 50 Kb with a 10 SNP step)’. The genotypes of selected sites were extracted from the biallelic SNP vcf file and subsequently transformed into the ‘admixture’ format using PLINK. The genotypes of 106 animals were utilized for three replicate runs conducted with K = 2–5, where K represents the assumed number of subpopulations. For principal component analysis (PCA), we initially utilized the PLINK software to convert the genotype files in vcf format into *.bed, *.bim, and *.fam files (--make-bed). Subsequently, the genetic relationship matrix and PCA were calculated using gcta64 software with the parameters: ‘--make-grm’ and ‘--pca’, respectively [29].

Linkage disequilibrium analysis

To assess population characteristics, we performed linkage disequilibrium decay analysis on all samples, horned population, and polled population, respectively. Initially, PLINK software was employed to convert vcf files into ped and map formats (--recode ped). Subsequently, the following parameters were provided for LD analysis in PLINK: ‘--map, --ped, --noweb, --allow-no-sex, --allow-extra-chr, --missing-genotype 0, --maf 0.05, --r2, --ld-window 999999, --ld-window-kb 1000, and --ld-window-r2 0’. The results are summarized using the script decay_chrom.pl. After identifying a significant correlation interval through GWAS, we further utilized Local LD haplotype block maps to determine the presence of an LD relationship between the significantly associated single nucleotide polymorphisms (SNPs) and the target genes. The analysis of LD blocks was conducted using LDBlockShow [30].

Genome-wide association analysis

GWAS was performed using a mixed linear model (MLM) implemented in TASSEL software [31]. Genotypes (SNPs) and population structure (Q) were considered as fixed effects while accounting for kinship (K) as random effects in the model. The Q matrix was derived from the results of population structure analysis, and the K matrix was calculated using gcta64 software [29]. The genotype data was first sorted using the run_pipeline.pl script. Subsequently, the following command was executed: “perl run_pipeline.pl -fork1 -vcf sort.vcf -fork2 -r trait.txt -fork3 -q Q.txt -excludeLastTrait -fork4 -k K.txt -combine5 -input1 -input2 -input3 -intersect -combine6 -input5 -input4 -mlm -mlmVarCompEst P3D -mlmCompressionLevel None -runfork1 -runfork2 -runfork3”. Manhattan plots were generated using the CMplot package (https://github.com/YinLiLin/R-CMplot) within R v4.1.0.

Results

Sequencing variants and population stratification

The whole-genome sequencing (WGS) data were aligned to the goat reference genome (Accession NO, GCA_001704415.1), which was assembled using genomic DNA from a horned San Clemente goat. This alignment resulted in a mean genome coverage of 99.8% and an average sequencing depth of 13.21× (Table S1). After filtering out loci with deletion rates exceeding 20% and minor allele frequencies below 0.05, we identified a total of 15,944,055 biallelic single nucleotide polymorphisms (SNPs) and 1,391,598 insertions/deletions (indels) in the cohort comprising 106 individuals.

After filtering out SNPs with a correlation coefficient greater than 0.2 (r2 > 0.2) in sliding windows of 50 Kb with a 10 SNP step, we performed a population structure analysis. PCA revealed a slight population stratification among the 106 individuals, potentially suggesting familial relatedness (Fig. 1). However, the cumulative proportion of genetic variance explained by the first two principal components was relatively modest, accounting for only 3.96% and 3.06%, respectively. Similar findings were obtained from population structure analysis (Figure S2).

Fig. 1
figure 1

PCA plot showing the population structure of 106 sampled goats based on the identified biallelic SNPs. The red and green circles represented polled and horned goats, respectively

Full size image

Genome-wide association study

After undergoing rigorous quality control and passing the Hardy-Weinberg test, we successfully identified a set of 15,059,449 biallelic SNPs that are suitable for subsequent GWAS. To account for potential confounding factors related to population structure, we performed a robust mixed linear model (MLM) association analysis using the 106 goats. The threshold value was calculated using the Bonferroni multiple test correction (0.05 divided by the total number of SNAs, at a 5% level of significance) [32]. A total of 320 significant SNPs (Bonferroni-corrected P < 0.05, −log10P = 8.17) associated with the horned and polled phenotypes were identified, all of which were located on chromosome 1 (Table S2). The Manhattan plot revealed two association signals (Fig. 2). The first peak of association spanned from 128,817,052 to 133,005,441 bp on chromosome 1 and consisted of 189 genome-wide significant SNPs. Among these SNPs, the top ten most associated ones (-lg10P > 12.71) were located at positions ranging from 129,410,163 to 129,791,577 bp. The second association signal was identified at the genomic coordinates of 150,336,143 – 150,808,639 bp on chromosome 1. This region encompassed a total of 131 genome-wide significant SNPs, with the top 15 SNPs (−lg10P > 12.01) being most strongly associated within the range of 150,376,727 – 150,744,739 bp (Table S2).

Fig. 2
figure 2

Manhattan plots depicting the results of a genome-wide association study for single nucleotide polymorphisms located on autosomes 1–29 associated with the polled traits. The chromosomes are depicted in distinct colors and presented separately. The horizontal red dashed line denotes the genome-wide Bonferroni-corrected significance threshold

Full size image

Characterization of the genome-wide significant SNPs

The variant annotation revealed that out of a total of 320 genome-wide significant SNPs, 232 were located in intergenic regions, 48 in intronic regions, 39 in ncRNA intronic regions, and one in gene downstream regions (Table S3). The first associated region comprises nine genes, namely NMNAT3, RBP2, MRPS22, PIK3CB, ESYT3, MRAS, SOX14, PPP2R3A, and FOXL2 (ENSCHIG00000026523), and six ncRNAs (ENSCHIG00000006793, ENSCHIG00000006799, ENSCHIG00000006802, ENSCHIG0000001813, U1-201, and U6). The second associated region encompasses two genes (ERG and KCNJ15), one snRNA (U6), and one small nucleolar RNA (SNORA70).

Identification and characterization of variants associated with the horned and polled phenotypes

The significantly associated SNPs of the first region predominantly localize within the range of 129,000,000 to 130,005,000 bp on chromosome 1. According to variant annotation, the top 10 genome-wide most significantly associated SNPs in the first significant signal were all located within intergenic regions. Among these SNPs, we observed predominant heterozygosity in almost all polled goats, while all 33 horned goats exhibited homozygosity and carried reference alleles at 7 loci along with variant alleles at 3 loci. All 3 PIS goats displayed homozygosity with 3 loci displaying reference alleles and 7 loci exhibiting variant alleles - an opposite pattern compared to that observed in horned goats. Interestingly, the identical haplotype to that of the intersex goat was also found in 5 polled male goats (Fig. 3). Previous studies have demonstrated a dominant inheritance pattern in polled goats, wherein homozygous polled females (XX) are associated with infertile intersex individuals, while homozygous polled males do not exhibit intersexuality [18]. The genotypes of 7 SNPs (Chr1:129410163, Chr1:129564298, Chr1:129567274, Chr1:129710577, Chr1:129789816, Chr1:129791507, and Chr1:129791577) in our results were consistent with the aforementioned research conclusions, suggesting their potential association with the polled phenotype and intersexuality, thus indicating their suitability as marker SNPs. To screen tag SNPs and determine the potential candidate genes in this specific genomic region, we conducted LD block analysis using LDBlockShow. We identified a total of 23 LD blocks (Table S4). 3 out of 7 marker SNPs (including the leading associated SNP) exhibited a strong linkage disequilibrium with the FOXL2 gene (Fig. 4), which is widely recognized as a pivotal determinant of female sex development in mammals [33]. Therefore, 3 SNPs, Chr1:129789816, Chr1:129791507, and Chr1:129791577, can be regarded as potential indicators of the PIS-affected goats.

Fig. 3
figure 3

Summary of the genotype distributions for the top 10 genome-wide significantly associated SNPs in the first identified genomic region among horned, polled, and PIS goats. The three PIS goats exhibited homozygosity at three loci with reference alleles, while seven loci showed variant alleles - a contrasting pattern to that observed in horned goats. Notably, the same haplotype as the intersex goat was also identified in five polled male goats

Full size image
Fig. 4
figure 4

Linkage disequilibrium (LD) heatmap depicts the regions surrounding the strong peaks of the first significant association region. The genes MRPS22, FOXL2, and MRAS exhibit strong linkage disequilibrium relationships with the significant genome-wide association SNPs identified in the first association signal

Full size image

In the largest LD blocks, we observed a strong LD between 3 exonic SNPS in gene MRPS22 (Chr1:129294870, Chr1:129300519, and Chr1:129300647) and 39 significant genome-wide association SNPs (Fig. 4 and Table S4). Among the 39 SNPs examined, we observed that all three PIS goats exhibited homozygosity and carried variant alleles, while not all 33 horned goats displayed homozygosity and carried reference alleles. According to the variant annotation, we observed two synonymous mutations (Chr1:129294870, c.G > A, p.T150T and Chr1:129300647, c.G > C, p.T3T) and one nonsynonymous mutation (Chr1:129300519, c.G > A, p.A46V) in the exonic regions of MRPS22, but they did not show significant associations. To further determine whether there is a causative relationship between MRPS22 and polled intersex syndrome, we performed the LD analysis on the horned population and the polled population, respectively. The results revealed a significant LD between the 39 genome-wide significant association SNPs and MRPS22, forming an independent LD block in horned goats, while no such LD block was observed in polled goats. These findings showed that MRPS22 may play a role in the development of PIS-affected goats. Additionally, 29 genome-wide significant association SNPs located in intronic regions of MRAS were found to be enriched in the second-largest LD blocks (Fig. 4 and Figure S3). However, MRAS exhibited 4 synonymous mutations (Chr1:130287936, c.C > T, p.T198T, Chr1:130287942, c.C > A, p.R196R, Chr1:130296909, c.G > A, p.H102H, and Chr1:130296915, c.G > A, p.F100F) in the exonic regions.

The second association peak included two genes, namely ERG and KCNJ15. As shown in Fig. 5, we observed a strong linkage disequilibrium between the leading associated SNP (Chr1:150739498) and other genome-wide significant SNPs, with an r2 value exceeding 0.6. The LD analysis revealed two genome-wide significant association SNP (Chr1:150612293 and Chr1:150612295) in ERG located in a 12 SNPs LD block (Table S5), while the two other genome-wide significant SNPs did not form an LD block (Figure S4). We identified that a stop-gained mutation (Chr1:150634326, c.C90G: p.Y30X) and a nonsynonymous mutation (Chr1:150561787, c.G1070C: p.R357P) in ERG. Additionally, we found that KCNJ15 did not exhibit the significance associated with the genome-wide significant SNPs. These findings suggested that ERG is a potential candidate gene influencing the horned and polled phenotypes.

Fig. 5
figure 5

Linkage disequilibrium (LD) heatmap depicts the regions surrounding the strong peaks of the second significant association region. The genes ERG and KCNJ15 exhibit strong linkage disequilibrium relationships with the significant genome-wide association SNPs identified in the second association signal

Full size image

Discussion

In this study, we performed a genome-wide association study (GWAS) using a high-density single nucleotide polymorphism (SNP) map obtained from whole-genome sequencing (WGS) to identify the genetic loci underlying the polled phenotype in Xinong Saanen dairy goats. Our analysis revealed 320 genome-wide significant SNPs located on chromosome 1, which formed two distinct association peaks consistent with previous GWAS findings for the polled phenotype. While our investigation showed slight variation in the range of these associated peaks compared to prior studies, the primary regions remained consistent, with both peaks centered around 129 Mb and 150.5 Mb [17]. These discrepancies may be attributed to variances in the detection platforms or algorithms employed during the corresponding analysis, alterations in the genetic background of the examined goats, disparities in the size and composition of the study population, or stochastic or technical errors encountered during certain analyses.

The most important problem in GWAS is the control of false positives for significant SNP [34]. In our study, there was a slight population stratification among the 106 individuals, which is a risk for the reliability of the results of GWAS. Currently, there exist four commonly used methods for addressing the issue of population stratification, namely genomic control [35], structured association analysis [36], PCA [36], and mixed linear models MLM [35]. This study chose MLM to control population stratification. The fixed effects in the model encompass genotypes and population structure while accounting for kinship as random effects [37]. This made the GWAS results more reliable.

Previous studies have reported that the absence of horn growth in goats is attributed to a complex structural genetic variant, characterized by the fusion of a large duplicated segment measuring 480 kb [38], located at position 150,334,286–150,818,098 bp on chromosome 1, into the deleted region spanning approximately 11.7 kb situated around 129 Mb on chromosome 1 [16, 17]. Despite the phenotypes, which are believed to result from an increase in gene copy number, often eluding explanation through changes in gene dosage [19]. The specific elements of these intricate structural variants that contribute to the absence of horn growth in goats remain unclear.

To further ascertain the causative variations underlying the two association signals, we performed the LD analysis. In the first associated genomic region, several SNPs and genes were identified as marker variants and potential candidate genes, respectively. For the first time, we identified 3 SNPs (Chr1:129789816, Chr1:129791507, and Chr1:129791577) that can be considered as potential indicators of the PIS-affected goats. These genome-wide significant SNPs have a strong LD relationship with gene FOXL2. FOXL2 is widely recognized as a pivotal determinant of female sex development in mammals [33]. FOXL2 was also identified as the headgear-specific expressed gene in bovids and cervids based on transcriptome analyses [39], indicating its potential involvement in horn growth. In addition, the genotype distributions at these loci are consistent with previous studies that there is the absence of intersexuality in homozygous polled goats with a male genotype (XY) and the inheritance of polledness in goats is predominantly associated with recessive intersexuality.

The LD block showed that there was a strong LD between gene MRPS22 and 39 genome-wide significant association SNPs. The protein MRPS22, an essential component of the small subunit of mitochondrial ribosomes, plays a pivotal role in reproduction and ovarian development [40, 41]. According to the variant annotation, we identified two synonymous mutations and one nonsynonymous mutation (Chr1:129300519, G > A, p.A46V) within the exonic regions of MRPS22. Although these variants did not show significant associations in the genome-wide analysis, they formed an independent LD block with the 39 genome-wide markers specifically in horned goats. In contrast, no such association was observed in polled goats (including PIS goats). The present findings suggested that MRPS22 may play a potential role in the development of infertile intersex individuals. Additionally, our findings have demonstrated that all 3 PIS goats exhibited homozygosity and carried variant alleles for the 39 significant genome-wide association SNPs in the first association region, whereas not all 33 horned goats displayed homozygosity and carried reference alleles. The genotype distributions at these loci exhibited an incomplete correlation with the horned, polled, and PIS phenotypes in the 106 Xinong Saanen dairy goats, suggesting that these 39 SNPs can not be considered as potential markers of the PIS-affected goats.

Another gene MRAS has 29 genome-wide significant association SNPs located in intronic regions. However, MRAS exhibited 4 synonymous mutations (Chr1:130287936, c.C > T, p.T198T, Chr1:130287942, c.C > A, p.R196R, Chr1:130296909, c.G > A, p.H102H, and Chr1:130296915, c.G > A, p.F100F) in the exonic regions. Therefore, these single nucleotide polymorphisms (SNPs) and the gene MRAS may not be considered as causative variations or candidate genes for the polled phenotype in goats. The remaining genes in the first associated region, namely NMNAT3, RBP2, PIK3CB, ESYT3, SOX14, and PPP2R3A, were not assigned to any LD block and their SNPs did not exhibit significant associations yet. Despite a previous study demonstrating the association of both SOX14 and ESYT3 with horn bud development in sheep [39, 42], our study findings do not support their role in the formation of thorned/polled phenotype in goats.

The second associated genomic region includes two genes, namely ERG and KCNJ15. The gene ERG is implicated in chromosomal translocations, resulting in the formation of distinct fusion gene products, such as TMPSSR2-ERG [43], which regulate the expressions of three osteoblastic markers (i.e., ET-1, COL1A1, and ALPL) at the cellular level [44]. The SNPs located within or near the ERG gene demonstrated significant associations with body conformation traits in cows [45]. Our LD analysis identified two genome-wide significant association SNPs (Chr1:150612293 and Chr1:150612295) in ERG that exhibit a strong linkage disequilibrium with 10 nearby SNPs, providing support for the potential regulatory roles of ERG in horn growth and metabolism in Xinong Saanen dairy goats. The gene KCNJ15 is a protein-coding gene. The protein encoded by this gene functions as an integral membrane protein and exhibits characteristics of an inward-rectifier type potassium channel [46]. By analyzing the expression profile in otosclerosis patients, several studies have demonstrated the involvement of KCNJ15 in aberrant bone remodeling processes in humans [47, 48]. However, our study revealed that the associations of SNPs in KCNJ15 did not reach genome-wide significance.

The results presented in this study are based on theoretical causal variations derived from bioinformatics analysis. While these findings offer valuable insights into potential genetic markers for horn status and PIS, further experimental validation is necessary to confirm their functional significance. Future studies could focus on exploring the functional roles of these candidate genes through molecular biology techniques such as gene expression profiling or knockout experiments that could shed light on their precise contributions to horn development and PIS in goats. This work represents an important step towards unraveling the complex genetics behind horn status and PIS in goats. Three marker SNPs, which have been identified as closely linked to PIS, may contribute to the development of SNP chips for early prediction of PIS in goats. This advancement will facilitate breeding programs aimed at producing fertile herds with polled traits.

Conclusion

In conclusion, our findings suggested the horned/polled phenotype of Xinong Saanen dairy goats is strongly associated with two regions located on chromosome 1. Furthermore, using GWAS and LD analyses, we successfully identified three marker SNPs that are closely linked to PIS and several candidate genes associated with polled traits in Xinong Saanen dairy goats. The findings could contribute to the breeding of fertile herds with polled traits.

Data availability

The requisite data for assessing the findings can be located within the manuscript and/or Supplementary Materials. The raw sequencing data of 106 Xinong Saanen dairy goats, comprising 33 horned individuals, 70 polled individuals, and 3 individuals with polled intersexuality syndrome, were deposited in the NCBI SRA database (accession: PRJNA1077730, https://www.ncbi.nlm.nih.gov/bioproject/PRJNA1077730). Any further inquiries should be directed to the corresponding author.

Abbreviations

GWAS:

Genome-wide association study

LD:

Linkage disequilibrium

MLM:

Mixed linear model

PCA:

Principal component analysis

PIS:

Polled intersexuality syndrome

SNP:

Single nucleotide polymorphism

WGS:

Whole genome sequencing

References

  1. Saleh AA, Rashad AMA, Hassanine N, Sharaby MA, Zhao Y. Evaluation of morphological traits and physiological variables of several Chinese goat breeds and their crosses. Trop Anim Health Prod. 2021;53(1):74.

    Article  PubMed  Google Scholar 

  2. Lv FH, Agha S, Kantanen J, Colli L, Stucki S, Kijas JW, Joost S, Li MH. Ajmone Marsan P: adaptations to climate-mediated selective pressures in sheep. Mol Biol Evol. 2014;31(12):3324–43.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  3. Chen N, Cai Y, Chen Q, Li R, Wang K, Huang Y, Hu S, Huang S, Zhang H, Zheng Z, et al. Whole-genome resequencing reveals world-wide ancestry and adaptive introgression events of domesticated cattle in East Asia. Nat Commun. 2018;9(1):2337.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Zeder MA. Core questions in domestication research. Proc Natl Acad Sci U S A. 2015;112(11):3191–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  5. Allais-Bonnet A, Grohs C, Medugorac I, Krebs S, Djari A, Graf A, Fritz S, Seichter D, Baur A, Russ I, et al. Novel insights into the bovine polled phenotype and horn ontogenesis in Bovidae. PLoS ONE. 2013;8(5):e63512.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Reedman CN, Duffield TF, DeVries TJ, Lissemore KD, Winder CB. Graduate Student Literature Review: role of pain mitigation on the welfare of dairy calves undergoing disbudding. J Dairy Sci. 2022;105(8):6809–19.

    Article  CAS  PubMed  Google Scholar 

  7. Schafberg R. The history of breeding for polled cattle. Livest Sci. 2015;179:54–70.

    Article  Google Scholar 

  8. Grobler R, Marle-Kster EV, Visser CJLS. Challenges in selection and breeding of polled and scur phenotypes in beef cattle. 2021:104479.

  9. Simon R, Drögemüller C, Lühken G. The Complex and Diverse Genetic Architecture of the absence of horns (polledness) in domestic ruminants, including goats and Sheep. Genes 2022;13(5).

  10. Wiedemar N, Tetens J, Jagannathan V, Menoud A, Neuenschwander S, Bruggmann R, Thaller G, Drögemüller C. Independent polled mutations leading to complex gene expression differences in cattle. PLoS ONE. 2014;9(3):e93435.

    Article  PubMed  PubMed Central  Google Scholar 

  11. Utsunomiya YT, Torrecilha RBP, Milanesi M, Paulan SC, Utsunomiya ATH, Garcia JF. Hornless Nellore cattle (Bos indicus) carrying a novel 110 kbp duplication variant of the polled locus. Anim Genet. 2019;50(2):187–8.

    Article  PubMed  Google Scholar 

  12. Randhawa IAS, Burns BM, McGowan MR, Porto-Neto LR, Hayes BJ, Ferretti R, Schutt KM, Lyons RE. Optimized genetic testing for polledness in multiple breeds of cattle. G3 (Bethesda Md). 2020;10(2):539–44.

    Article  CAS  PubMed  Google Scholar 

  13. Duijvesteijn N, Bolormaa S, Daetwyler HD, van der Werf JHJ. Genomic prediction of the polled and horned phenotypes in Merino sheep. Genet Selection Evolution: GSE. 2018;50(1):28.

    Article  PubMed Central  Google Scholar 

  14. He X, Zhou Z, Pu Y, Chen X, Ma Y, Jiang L. Mapping the four-horned locus and testing the polled locus in three Chinese sheep breeds. Anim Genet. 2016;47(5):623–7.

    Article  CAS  PubMed  Google Scholar 

  15. Lühken G, Krebs S, Rothammer S, Küpper J, Mioč B, Russ I, Medugorac I. The 1.78-kb insertion in the 3’-untranslated region of RXFP2 does not segregate with horn status in sheep breeds with variable horn status. Genet Selection Evolution: GSE. 2016;48(1):78.

    Article  PubMed Central  Google Scholar 

  16. Pailhoux E, Vigier B, Chaffaux S, Servel N, Taourit S, Furet JP, Fellous M, Grosclaude F, Cribiu EP, Cotinot C, et al. A 11.7-kb deletion triggers intersexuality and polledness in goats. Nat Genet. 2001;29(4):453–8.

    Article  CAS  PubMed  Google Scholar 

  17. Guo J, Jiang R, Mao A, Liu GE, Zhan S, Li L, Zhong T, Wang L, Cao J, Chen Y, et al. Genome-wide association study reveals 14 new SNPs and confirms two structural variants highly associated with the horned/polled phenotype in goats. BMC Genomics. 2021;22(1):769.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Pannetier M, Elzaiat M, Thépot D, Pailhoux E. Telling the story of XX sex reversal in the goat: highlighting the sex-crossroad in domestic mammals. Sex Development: Genet Mol Biology Evol Endocrinol Embryol Pathol sex Determ Differ. 2012;6(1–3):33–45.

    Article  CAS  Google Scholar 

  19. Zlotorynski E. Genome organization: add a TAD of duplication. Nat Rev Mol Cell Biol. 2016;17(12):737.

    Article  CAS  PubMed  Google Scholar 

  20. Liang D, Zhao P, Si J, Fang L, Pairo-Castineira E, Hu X, Xu Q, Hou Y, Gong Y, Liang Z, et al. Genomic analysis revealed a convergent evolution of LINE-1 in Coat Color: a Case Study in Water buffaloes (Bubalus bubalis). Mol Biol Evol. 2021;38(3):1122–36.

    Article  CAS  PubMed  Google Scholar 

  21. Liu D, Zhang H, Yang Y, Liu T, Guo Z, Fan W, Wang Z, Yang X, Zhang B, Liu H, et al. Metabolome-based Genome-Wide Association Study of Duck Meat Leads to Novel genetic and biochemical insights. Adv Sci (Weinheim Baden-Wurttemberg Germany). 2023;10(18):e2300148.

    Google Scholar 

  22. Wang J, Yang W, Zhang S, Hu H, Yuan Y, Dong J, Chen L, Ma Y, Yang T, Zhou L, et al. A pangenome analysis pipeline provides insights into functional gene identification in rice. Genome Biol. 2023;24(1):19.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, et al. The genome analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20(9):1297–303.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJ, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81(3):559–75.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Wang K, Li M, Hakonarson H. ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res. 2010;38(16):e164.

    Article  PubMed  PubMed Central  Google Scholar 

  29. Yang J, Lee SH, Goddard ME, Visscher PM. GCTA: a tool for genome-wide complex trait analysis. Am J Hum Genet. 2011;88(1):76–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Dong SS, He WM, Ji JJ, Zhang C, Guo Y, Yang TL. LDBlockShow: a fast and convenient tool for visualizing linkage disequilibrium and haplotype blocks based on variant call format files. Brief Bioinform. 2021;22(4).

  31. Bradbury PJ, Zhang Z, Kroon DE, Casstevens TM, Ramdoss Y, Buckler ES. TASSEL: software for association mapping of complex traits in diverse samples. Bioinformatics. 2007;23(19):2633–5.

    Article  CAS  PubMed  Google Scholar 

  32. McCarthy MI, Abecasis GR, Cardon LR, Goldstein DB, Little J, Ioannidis JP, Hirschhorn JN. Genome-wide association studies for complex traits: consensus, uncertainty and challenges. Nat Rev Genet. 2008;9(5):356–69.

    Article  CAS  PubMed  Google Scholar 

  33. Uhlenhaut NH, Jakob S, Anlag K, Eisenberger T, Sekido R, Kress J, Treier AC, Klugmann C, Klasen C, Holter NI, et al. Somatic sex reprogramming of adult ovaries to testes by FOXL2 ablation. Cell. 2009;139(6):1130–42.

    Article  CAS  PubMed  Google Scholar 

  34. Jiang Y, Tang S, Wang C, Wang Y, Qin Y, Wang Y, Zhang J, Song H, Mi S, Yu F, et al. A genome-wide association study of growth and fatness traits in two pig populations with different genetic backgrounds. J Anim Sci. 2018;96(3):806–16.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Wang K, Liu D, Hernandez-Sanchez J, Chen J, Liu C, Wu Z, Fang M, Li N. Genome Wide Association Analysis Reveals New Production Trait Genes in a male Duroc Population. PLoS ONE. 2015;10(9):e0139207.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Fontanesi L, Schiavo G, Galimberti G, Calò DG, Scotti E, Martelli PL, Buttazzoni L, Casadio R, Russo V. A genome wide association study for backfat thickness in Italian large White pigs highlights new regions affecting fat deposition including neuronal genes. BMC Genomics. 2012;13:583.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. Yu J, Pressoir G, Briggs WH, Vroh Bi I, Yamasaki M, Doebley JF, McMullen MD, Gaut BS, Nielsen DM, Holland JB, et al. A unified mixed-model method for association mapping that accounts for multiple levels of relatedness. Nat Genet. 2006;38(2):203–8.

    Article  CAS  PubMed  Google Scholar 

  38. Simon R, Lischer HEL, Pieńkowska-Schelling A, Keller I, Häfliger IM, Letko A, Schelling C, Lühken G, Drögemüller C. New genomic features of the polled intersex syndrome variant in goats unraveled by long-read whole-genome sequencing. Anim Genet. 2020;51(3):439–48.

    Article  CAS  PubMed  Google Scholar 

  39. Wang Y, Zhang C, Wang N, Li Z, Heller R, Liu R, Zhao Y, Han J, Pan X, Zheng Z et al. Genetic basis of ruminant headgear and rapid antler regeneration. Science 2019, 364(6446).

  40. França MM, Mendonca BB. Genetics of ovarian insufficiency and defects of folliculogenesis. Best Pract Res Clin Endocrinol Metab. 2022;36(1):101594.

    Article  PubMed  Google Scholar 

  41. Tiosano D, Mears JA, Buchner DA. Mitochondrial dysfunction in primary ovarian insufficiency. Endocrinology. 2019;160(10):2353–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  42. Wilmore HP, Smith MJ, Wilcox SA, Bell KM, Sinclair AH. SOX14 is a candidate gene for limb defects associated with BPES and Möbius syndrome. Hum Genet. 2000;106(3):269–76.

    Article  CAS  PubMed  Google Scholar 

  43. Deplus R, Delliaux C, Marchand N, Flourens A, Vanpouille N, Leroy X, de Launoit Y, Duterque-Coquillaud M. TMPRSS2-ERG fusion promotes prostate cancer metastases in bone. Oncotarget. 2017;8(7):11827–40.

    Article  PubMed  Google Scholar 

  44. Delliaux C, Tian TV, Bouchet M, Fradet A, Vanpouille N, Flourens A, Deplus R, Villers A, Leroy X, Clézardin P, et al. TMPRSS2:ERG gene fusion expression regulates bone markers and enhances the osteoblastic phenotype of prostate cancer bone metastases. Cancer Lett. 2018;438:32–43.

    Article  CAS  PubMed  Google Scholar 

  45. Cole JB, Wiggans GR, Ma L, Sonstegard TS, Lawlor TJ Jr., Crooker BA, Van Tassell CP, Yang J, Wang S, Matukumalli LK, et al. Genome-wide association analysis of thirty one production, health, reproduction and body conformation traits in contemporary U.S. Holstein cows. BMC Genomics. 2011;12:408.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Zhang J, Han J, Li L, Zhang Q, Feng Y, Jiang Y, Deng F, Zhang Y, Wu Q, Chen B, et al. Inwardly rectifying potassium channel 5.1: structure, function, and possible roles in diseases. Genes Dis. 2021;8(3):272–8.

    Article  CAS  PubMed  Google Scholar 

  47. Karosi T. I Sziklai 2010 Etiopathogenesis of otosclerosis. Eur Archives oto-rhino-laryngology: Official J Eur Federation Oto-Rhino-Laryngological Soc (EUFOS): Affiliated German Soc Oto-Rhino-Laryngology - Head Neck Surg 267 9 1337–49.

    Google Scholar 

  48. Ealy M, Chen W, Ryu GY, Yoon JG, Welling DB, Hansen M, Madan A, Smith RJ. Gene expression analysis of human otosclerotic stapedial footplates. Hear Res. 2008;240(1–2):80–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

Heartfelt gratitude is extended to the diligent and dedicated research workers in the Dairy Goat Lab of Northwest A&F University. Their invaluable contributions have been instrumental in advancing our understanding of this vital field, and their insightful discussions have proven to be an indispensable resource.

Funding

This work was supported by Shaanxi Livestock and Poultry Breeding Double-chain Fusion Key Project of China (No. 2022GD-TSLD-46-0201) and Xinjiang Tianchi Outstanding Scholars Program (2023-068).

Author information

Authors and Affiliations

  1. Key Laboratory of Animal Genetics, Breeding and Reproduction of Shaanxi Province, College of Animal Science and Technology, Northwest A&F University, Yangling, 712100, P. R. China

    Fuhong Zhang, Qingqing Liu, Yaling Wang, Chenbo Shi, Lu Zhu, Jianqing Zhao, Weiwei Yao & Jun Luo

  2. Institute of Animal Husbandry Quality Standards, Xinjiang Academy of Animal Sciences, Urumchi, 830000, P. R. China

    Ping Gong

Authors
  1. Fuhong Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Qingqing Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Ping Gong
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Yaling Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Chenbo Shi
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Lu Zhu
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Jianqing Zhao
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Weiwei Yao
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Jun Luo
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

JL was in charge of designing the experiment. CBS, LZ, QQL, YLW, and WWY were tasked with collecting tissue samples and conducting verification tests. FHZ and JQZ curated the data and bioinformatics analysis. FHZ wrote the manuscript while FHZ, PG, and JL collaborated on interpreting the results and revising the manuscript. All authors read and approved the final version.

Corresponding author

Correspondence to Jun Luo.

Ethics declarations

Ethics approval and consent to participate

All animals received humane care as outlined in the Guide for the Care and Use of Experimental Animals of the National Institutes of Health. The animals used in study were under the ownership of the College of Animal Science and Technology, Northwest A&F University. Informed consent has been obtained from the owner. The experimental protocol, including blood sampling and PIS goat identification, was approved by the Animal Ethical and Welfare Committee of the College of Animal Science and Technology, Northwest A&F University, Yangling, China (protocol number DK2022008).

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Supplementary Material 2

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Zhang, F., Liu, Q., Gong, P. et al. Genome-wide association study provided insights into the polled phenotype and polled intersex syndrome (PIS) in goats. BMC Genomics 25, 661 (2024). https://doi.org/10.1186/s12864-024-10568-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-024-10568-9

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us