Two significant SNPs were found associated with the litter size: an SNP located in the genomic region on chromosome 6 (casein kappa [CSN3]) at the genome-wide level, and another SNP located in the genomic region on chromosome 24 (transcription factor 4 [TCF4]) at the chromosome-wide level (Figure 1A). The CSN3 gene, which encodes a major milk protein known as kappa-casein, plays a vital role in protein synthesis and stability of casein micelle. Another study suggests this gene is an essential putative candidate gene for milk-producing traits in different dairy cow and goat breeds [18]. Lan et al [19] demonstrated that CSN3-TaqI locus was associated with litter size in Saanen dairy goats, which is consistent with our result. In addition, several candidate genes involved in this trait have been identified in other goat populations. Previously, Lan et al [20] reported two SNPs to be significantly associated with the litter size in Yunnan black goats, located downstream of the solute carrier family 4 member 10 (SLC4A10) gene and upstream of the T-box brain transcription factor 1 (TBR1) gene, respectively. Islam et al [21] identified several candidate genes like meiosis arrest female 1 (MARF1), synaptonemal complex protein 2 (SYCP2), transmembrane protein 200C (TMEM200C), steroidogenic factor-1 (SF1), adenylate cyclase 1 (ADCY1), and bone morphogenetic protein 5 (BMP5) involved in goat fecundity in Chinese goat populations. In another study, the KiSS-1 metastasis suppressor (KISS1) gene was related to the litter size in the Chinese Arbas Cashmere goat [3], coincident with the findings from the Black Bengal goat and Jining grey goat. Genome-wide screening in the Yunshang black goat identified some loci like gap junction protein beta 2 (GJB2), gap junction protein beta 6 (GJB6), stabilin 2 (STAB2), and heparan sulfate-glucosamine 3-sulfotransferase 2 (HS3ST2) to significantly affect the litter size [22]. Most recently, a genome-wide selective sweep analysis has revealed that the insertion/deletion variants within the insulin-like growth factor 2 mRNA-binding protein 2 (IGF2BP2) gene were correlated with the litter size in Shaanbei white Cashmere goat [23].

The above research results demonstrated various efforts for identifying the causal genes responsible for litter size in goats, but the results vary among breeds. These inconsistent findings might suggest that the complex genetic mechanisms underpin litter size variation. Our results provided a helpful reference for explaining the genetic mechanism of litter size in goats.

Nipples (or teats) are epidermal appendages on the udders or breasts of mammals. Nipple number is a valuable trait crucial to female goat reproduction performance and offspring growth [6]. Sheep and goats normally have only two functional nipples, whereas a small number of female goats under investigation possess one or two supernumerary nipples. Studies suggested that supernumerary nipples might be susceptible to bacterial ingress through the teat canal and consequently impede their production efficiency considerably [24]. Therefore, the supernumerary nipples in goats and sheep are undesirable.

In this study, no SNPs were identified for nipple number trait (Figure 1B). The BBX high mobility group box domain containing (BBX) and CD47 molecule (CD47) gene has been reported to be significantly linked to the supernumerary nipple phenotype in Wadi sheep [5]. Additionally, another recent genome-wide comparative analysis has shown the γ-aminobutyric acid type A receptor subunit γ3 (GABRG3), LDL receptor related protein 1B (LRP1B), and mono-ADP ribosylhydrolase 2 (MACROD2) genes to be associated with the nipple numbers in Luzhong mutton sheep [25]. Furthermore, there was noteworthy relevance of nipple number and reproduction traits, e.g., litter size in mammals such as pigs and goats [6], suggesting that further research could combine both the nipple number and litter size trait for seeking candidate genes. Nevertheless, a GWAS on Saanen dairy goat illustrated 17 regions on ten chromosomes to be significantly associated with the supernumerary nipple trait at the chromosome-wide level. In contrast, there were no SNPs at the genome-wide level. It suggested that supernumerary nipple trait might be inherited in a polygenic fashion, consistent with the recent findings in cattle [7]. Therefore, the failure in identifying the SNPs for the nipple number trait implied that other causal loci and feasible mechanisms might exist in the three Chinese native goat populations.

Congenital appendages (wattles) are commonly found on the throat of goats and are reported to be under the genetic control by the W gene with an inheritance. Previous studies reported that female goats with wattle have 7% higher fecundity than those without wattle, and highlighted the importance of this trait in marker-assisted breeding.

At the chromosome-wide level, two significant SNPs were identified for the wattle trait located in the genomic regions on chromosomes 9 (katanin catalytic subunit A1 [KATNA1]) and 1 (ubiquitin associated and SH3 domain containing A [UBASH3A]) (Figure 1C). Katanin is a heterodimer microtubule-severing enzyme comprising a catalytic p60 subunit A1 and an 80 kDa regulatory p80 subunit B1, in which KATNA1 encodes the p60 subunit A. It is also involved in regulating neuronal progenitor proliferation in vivo during embryonic development and adult neurogenesis. On the other hand, the USBSH3A gene encodes one of the two family members belonging to the T-cell ubiquitin ligand which can negatively regulate the T-cell signaling. Previous studies have indicated that genetic variations in the UBASH3A gene were associated with several autoimmune diseases such as diabetes, rheumatoid arthritis, and systemic lupus erythematosus.

Earlier, a GWAS performed in nine different Swiss goat breeds identified two significant SNPs on chromosome 9 which contained the formin 1 (FMN1) and gremlin 1 (GREM1) genes [13], responsible for limb development and outgrowth. These findings supported earlier assumptions that wattles are rudimentary developed extremities. Also, this is the first study demonstrating the associations between the KATNA1 and UBASH3A genes with wattle in goats based on our current knowledge.

The genetic determinants of goat coat color are not only important for improving the artificial selection efficiency of the desired coat color but also crucial for understanding the different evolutionary dynamics in shaping the phenotypic variation in farm and wild animals. Two significant SNPs were identified for coat color trait (white vs brown) at the genome-wide level and were located in the genomic regions on chromosomes 26 (sortilin-related VPS10 domain-containing receptor 3 [SORCS3]) and 24 (Golgi-associated protein Dymeclin [DYM]). On the other hand, at the chromosome-wide level, one significant SNP was identified for the coat color trait located in the genomic regions on chromosome 20 (phosphodiesterase 4D [PDE4D]) (Figure 1D).

The SORCS3 gene encodes a type-I receptor transmembrane protein having a vital role as a sorting agency within the cells and transports various intracellular proteins. The DYM gene encodes a Golgi-associated protein (Dymeclin), whose deficiency typically leads to Dyggve-Melchior-Clausen syndrome (DMC). The melanosomes are discontinuous membranous organelles formed by the Golgi apparatus-endoplasmic reticulum-lysosome, and they are the only organelles synthesizing and depositing melanin in the eukaryotic cells. In addition, the studies in mice and broilers revealed the SORCS3 gene to affect feed efficiency, energy balance, and orexigenic peptide production and suggested a relationship between pigmentation and energy metabolism [26,27]. Thus, it is most probable that SORCS3 and DYM have essential roles in the variation of coat color in goats in the present study. Moreover, the PDE4D gene, encodes the phosphodiesterase 4D having a 3′,5′-cyclic-AMP phosphodiesterase activity, with an essential role in signaling pathways of melanin synthesis. Although well-known genes like agouti signaling protein (ASIP), tyrosinase related protein 1 (TYRP1), melanocyte inducing transcription factor (MITF), melanocortin 1 receptor (MC1R), and KIT proto-oncogene (KIT) have been associated with coat color in various goat populations [28], our results provided a new variation source of coat color in goats.

The color pattern is important for both wild and domesticated animals. The black dorsal line is a unique coat color pattern in goats and many animals such as Chinese hamsters, Sugar gliders, and slow loris. The genetic basis of this trait mechanism remains unclear.

At the genome-wide level and chromosome-wide level, two significant SNPs were identified for the black dorsal line trait located in the genomic regions on chromosome 18 (transcription factor 25 [TCF25]) and 15 (CXADR like membrane protein [CLMP]), respectively (Figure 1E). The T cell transcription factors (TCF) are essential modulators for Wnt-induced transcriptional activation, with essential roles in cell fate decisions during embryonic development and stem cell homeostasis in somatic niches of normal tissues. The TCF transcription factors participate in melanoma progression by regulating the versican expression. GWAS by Wang et al [29] demonstrated that TCF25 might be a putative candidate gene for coat color in Tan sheep. In mouse embryos, TCF25 is intensely expressed in the dorsal root ganglia. For this reason, we assume that the regional differential expression of TCF25 might modulate the subsequent melanocyte migration leading to the dark middorsal strip in goats. The previous investigation on the CLMP gene has reported this gene to be responsible for encoding the CXADR-like membrane protein localized to the junctional complexes between the endothelial and epithelial cells and might have a role in cell-cell adhesion. The expression of the CLMP gene in the human white adipose tissue is implicated in adipocyte maturation and obesity development. Therefore, the association between the black dorsal line and the CLMP gene remains yet to be verified.

Beard is a qualitative trait that significantly benefits the survivability, climate adaptability, and productivity of goats in tropical regions. The beard traits have been evaluated in some goat breeds [2]. However, the underlying molecular mechanism needs more clarification.

At the genome-wide level, one significant SNP was identi fied for beard trait on chromosome 8, but no candidate gene was found. While at the chromosome-wide level, three significant SNPs were identified for beard trait in the genomic regions on chromosome 2 (PAX3), 5 (PIK3C2G), and 28 (PLA2G12B and OIT3) (Figure 1F). The PAX3 gene, which encodes a member of the paired box family of transcription factor 3 (PAX3), is critical for developing melanoblasts. The PIK3C2G gene encodes phosphatidylinositol-4-phosphate 3-kinase catalytic subunit type 2 gamma, associated with cell migration. According to other reports, the protein encoded by the PLA2G12B gene belongs to the phospholipase A2 (PLA2) group of enzymes and might function as a negative regulator of tumor progression. While the OIT3 gene encoding the oncoprotein-induced transcript 3 is a novel biomarker of alternatively activated macrophages and facilitates hepatocellular carcinoma metastasis.

Previously, a GWAS performed on the Beijing fatty chicken showed that the copy number variations in the homeobox B7 (HOXB7) and homeobox B8 (HOXB8) genes were involved in forming beard traits [30]. While in goats, evaluating the beard traits has been carried out in some goat breeds in Southwest China recently [2], but no previous studies have so far investigated the goat beard traits by genome selection approach. To our knowledge, this study attempted to study the candidate genes or SNPs associated with goat beard traits using GWAS for the first time, and our results might provide a potential phenotypic selection strategy for goat breeding in Southwest China.

Some of the goats under investigation were found to have long beards while others have short beards. A study in anestrus dairy goats evaluated some morphometric traits including beard length, for exploring the possible role of these traits in goat reproductive success [31]. Such morphometric advantages shown by the high social rank goats emphasized the need for comprehending the biological basis of relevant animal traits and yet, the genetic mechanism of beard length remains unclear. Therefore, exploring the genetic basis of beard length traits could facilitate the sustainable utilization, conservation, and improvement of the genetic resources in goats.

At the genome-wide level, one significant SNP was identi fied for the beard length trait and was located in the genomic regions on chromosome 18 (KCNG4) (Figure 1G). The KCNG4 gene encodes the potassium voltage-gated channel modifier subfamily G member 4, functioning as a modulatory subunit. A previous study also identified the KCNG4 gene as a candidate gene for multiple sclerosis in humans. This study could be the first one to report that the KCNG4 gene is associated with the beard length trait in goats, and the results could also help in revealing the genetic mechanism of beard length traits in goats.

Genetic improvement of hair/wool-related features is a major goal in the goat and sheep industry, while hind leg hair is associated with the hair/wool quality and quantity in goats. Although the genetic variation of hair/wool traits has been explored in sheep and goats [8], the genetic architecture controlling the wool traits has not been fully elucidated. Therefore, a comprehensive understanding of the molecular genetic mechanisms of wool traits is needed for improving the breeding efficiency of goats.

At the chromosome-wide level, three significant SNPs were identified for the hind leg hair traits. One of the identified SNPs is located on chromosome 17 (between GLRB and GRIA2); the other SNP is located on chromosome 28 (PGBD5); the third SNP is located on chromosome 4 but not annotated to any known gene (Figure 1H). The GLRB gene encodes the glycine receptor β subunit with an important role in binding agonists and conducting ion channels. The missense mutations in the glycine receptor β subunit might underlie disease pathology in the human startle disease. The GRIA2 gene encodes glutamate ionotropic receptor AMPA type subunit 2, which is a predominant excitatory neurotransmitter receptor in the mammalian brain. A recent study showed that the de-novo variants in the GRIA2 gene could cause neurodevelopmental disorders in humans. The PGBD5 gene encodes the piggyBac transposable element derived protein 5, which belongs to the group of transposases. According to other reports, the PGBD5 gene could also promote site-specific oncogenic mutations in human solid tumors. Furthermore, a previous GWAS performed on yearling Merino sheep identified that the fin bud initiation factor homolog (FIBIN), hydroxysteroid 17-beta dehydrogenase 11 (HSD17B11), and protein inhibitor of activated STAT 1 (PIAS1) genes are related to the development and growth of hair/wool, which is inconsiderate with the genes identified in our study. To date, our study could be the first to report that the GLRB, GRIA2, and PGBD5 genes are associated with the hind leg hair traits in goats, and the results could also help to reveal the genetic mechanism of hair/wool traits in goats.

In some sheep breeds, black skin spots (skin pigmentation) can also produce pigmented fibers which significantly decrease the fiber quality. Therefore, elucidating the genetic mechanism of skin color in goats and sheep is of great economic value.

In this study, no SNPs were identified for the skin color traits (Figure 1I). Earlier studies have reported that the TYRP1, FBJ murine osteosarcoma viral oncogene homolog (C-FOS), Krüppel-like factor 4 (KLF4), ubiquitin-fold modifier conjugating enzyme 1 (UFC1), and ASIP genes are associated with skin pigmentation in sheep and goats [32,33]. However, the genetic mechanism underlying skin color traits remains unclear. For livestock, a majority of complex traits are likely controlled by a larger number of genes that individually confer a small effect. Furthermore, a GWAS performed on the Hereford cattle provided evidence that eyelid pigmentation is a heritable trait influenced by many loci [12]. Moreover, given that most tests for association proceed in a simple SNP-by-trait fashion using stringent significance levels to offset multiple testing, many of the genes contributing to trait variation remain undetected [33]. However, in a study performed on the pigmentation phenotype (piebald) in Merino sheep, insulin-like growth factor binding protein 7 (IGFBP7), platelet-derived growth factor alpha (PDGFRA), and CD9 molecule (CD9) were identified as gene drivers for the piebald phenotype by combing GWAS with gene expression data [10]. This study was implicated in exploring approaches that sought to incorporate GWAS into systems biology, such as considering merging SNP variation with the analyses of differential gene expression for investigating the basis of a skin pigmentation phenotypic trait.