Genomic Characteristics and Selection Signatures in Indigenous Chongming White Goat (Capra hircus)

The Chongming white goat (CM) is an indigenous goat breed exhibits unique traits that are adapted to the local environment and artificial selection. By performing whole-genome re-sequencing, we generated 14–20× coverage sequences from 10 domestic goat breeds to explore the genomic characteristics and selection signatures of the CM breed. We identified a total of 23,508,551 single-nucleotide polymorphisms (SNPs) and 2,830,800 insertion–deletion mutations (indels) after read mapping and variant calling. We further specifically identified 1.2% SNPs (271,713) and 0.9% indels (24,843) unique to the CM breed in comparison with the other nine goat breeds. Missense (SIFT < 0.05), frameshift, splice-site, start-loss, stop-loss, and stop-gain variants were identified in 183 protein-coding genes of the CM breed. Of the 183, 36 genes, including AP4E1, FSHR, COL11A2, and DYSF, are involved in phenotype ontology terms related to the nervous system, short stature, and skeletal muscle morphology. Moreover, based on genome-wide FST and pooled heterozygosity (Hp) calculation, we further identified selection signature genes between the CM and the other nine goat breeds. These genes are significantly associated with the nervous system (C2CD3, DNAJB13, UCP2, ZMYND11, CEP126, SCAPER, and TSHR), growth (UCP2, UCP3, TSHR, FGFR1, ERLIN2, and ZNF703), and coat color (KITLG, ASIP, AHCY, RALY, and MC1R). Our results suggest that the CM breed may be differentiated from other goat breeds in terms of nervous system owing to natural or artificial selection. The whole-genome analysis provides an improved understanding of genetic diversity and trait exploration for this indigenous goat breed.


INTRODUCTION
Global livestock genetic diversity encompasses species, breeds, strains, and their variations. The diversity between and within breeds rather than species is always considered to be essential for the purposes of conservation biology (Yaro et al., 2017). Goats are important husbandry animals with an ancient domestication history and economic value. It has been reported that their domestication occurred approximately 10,000 years ago and spread worldwide following human migrations and trade routes (Benjelloun et al., 2015;Kim et al., 2019). China is the global leader in terms of goat production (Skapetas and Bampidis, 2016), including many commercial, indigenous, and composite breeds. Owing to the weaker production potential compared with that of some global commercial breeds, indigenous goat breeds are facing the problem of genetic invasion and resource degradation in recent years (Kim et al., 2019;Monau et al., 2020).
Indigenous goats are unique groups that have developed under the forces of natural selection, domestication, and local environmental conditions. Their genetic background, unique traits, and biodiversity, as well as their environmental adaptability, serve as a kind of biological heritage that can be used for the genetic improvement of animal husbandry (Kim et al., 2019). The Chongming white goat (CM) is an indigenous island-type goat breed distributed and bred on Shanghai Chongming Island. Although there have not been sufficient systematic studies on the selection signal genes for island climate adaptation, this may be a direction of research. The CM goats were also marked as special agricultural products with geographical indications registered by the State Administration for Industry and Commerce of China. Moreover, the relatively small size and ease of husbandry make it a suitable animal model for some specific human diseases such as acute spinal cord compression injuries . However, even with rich domestication history, studies on the genome-wide characteristics of this indigenous goat have not been reported yet to the best of our knowledge.
Whole-genome sequencing (WGS) information bypasses the study limitations of maternal mitochondrial DNA and paternal Y chromosome inheritance in species evolution and population history dynamics. Moreover, WGS can also ameliorate the current problems of breed bias and insufficient markers in genotyping chips, and the new genetic variation information also provides research materials for the further production of highdensity chips. The use of WGS and reference genome alignment can identify more comprehensive, genome-wide variation, such as the presence of single-nucleotide polymorphisms (SNPs), insertion-deletion mutations (indels), and structural variants (SVs). Although high-throughput sequencing technology enables the comparison of patterns of polymorphisms at wholegenome levels, characterizing genetic structure from individual sequencing data still remains expensive for non-model species. The DNA pooling strategy (Pool-seq) represents an attractive and cost-effective alternative for SNP discovery (Futschik and Schlötterer, 2010;Hivert et al., 2018) and have been successfully applied to many domestic animal studies such as chicken (Rubin et al., 2010), pig (Rubin et al., 2012), goat Guo et al., 2018), and sheep .
Molecular genetic markers, such as high-density SNPs, are often used to analyze genetic variation at the genome-wide level caused by natural and artificial selection of species. For example, the Chinese Tibetan cashmere goat adapted to a highaltitude area and hypoxic environment. In which, 339 genes were identified through exome sequencing, which are potentially under high-altitude selection (Song et al., 2016). Another study (Kim et al., 2019) analyzed the genome of the Korean indigenous goat (KNG) in comparison with that of nine other goat breeds and revealed that the KNG has selection signatures for Salmonella infection and cardiomyopathy pathways. Wang et al. (2016) generated 9-13× coverage sequences from eight domesticated goat breeds and identified 22 genomic regions that may be associated with specific phenotypic traits, including coat color, body size, cashmere traits, and high-altitude adaptation.
In this study, we aimed to explore the genomic characteristics in the CM goat. Utilizing WGS and genomic approaches, we conducted the first comparative genomic study to reveal the genome-wide variation and selection signatures in the CM breed.

Animals
We collected 100 goat samples from 10 breeds, with 10 individuals per breed (five bucks and five does). Samples of each breed were randomly collected from two to three local breeding farms in their province ( Figure 1A). The 10 breeds included an indigenous CM breed from the farms of Shanghai Academy of Agricultural Sciences (SAAS), which were located on the Chongming Island of Shanghai. The other nine breeds were from the following regions of China: commercial Boer goats (Bor) from Jiangsu Province; Haimen goats (HM) from Jiangsu Province; Xuhuai goats (XH) from Jiangsu Province; Matou goats (MT) from Hubei Province; Da'er goats (DE) from Sichuan Province; Nanjiang yellow goats (NJ) from Sichuan Province; Macheng black goats (MB) from Hubei Province; Yichang white goats (YC) from Yichang City, Hubei Province; and Qin goats (QS) from Jining City, Shandong Province. The phenotypic descriptions of these 10 breeds can be found in Supplementary Data S1.

Genomic DNA Extraction and Whole-Genome Sequencing
The genomic DNA was extracted from goat ear punches using the cetyl trimethylammonium bromide (CTAB) method (Thomas et al., 1997). DNA isolated from 10 animals per breed was pooled evenly (by µg of DNA) into a single sample. We ensured that the qualities of all the DNA samples were recorded as given OD260/280 = 1.8∼2.0 and OD260/230 2.0. One microgram of pooling DNA was used for library preparation using the TruSeq TM Nano DNA Library Prep Kit (Illumina, San Diego, CA, United States) according to the manufacturer's guidelines. Briefly, the qualified DNA was randomly fragmented using ultrasoundbased strategy (Larguinho et al., 2010). The fragments were then subjected to end-repair followed by adaptor ligation to the fragments. Fragments of ∼400 bp in length were selected, enriched by magnetic beads, and further amplified by polymerase chain reaction. High-throughput sequencing was performed using the Illumina HiSeq TM platform (2 × 150 bp read length). The library preparation and sequencing were completed by Shanghai Majorbio Company, China 1 .

Raw Data Quality Control and Filtering
The raw reads underwent the following filter to get the clean data: (1) remove adaptor sequences; (2) remove bases containing non-AGCT at the 5' end; (3) trim the end of the reads with the sequencing quality score less than Q20; (4) remove reads containing N > 10%; and (5) discard reads with length less than 25 bp after above filtering.
After the quality control and filtering, the clean reads were summarized with the following statistics: (1) number of the reads; (2) number of the bases; (3) percentages of Q30 bases; (4) average and median insertion size; and (5) average sequencing depth.

Genomic Alignment and Variant Calling
We used the San Clemente goat's WGS (ARS1) (Bickhart et al., 2017) as the mapping reference genome 2 , which is the latest de novo assembled sequence (Capra hircus). The clean reads were mapped against the ARS1 by using Burrows Wheeler Alignment (BWA v0.7.15) (Li and Durbin, 2009) with default parameters. We used the Picard v2.23.1 3 to remove duplicate reads and the GATK (v3.8.0) (McKenna et al., 2010) to call variants (SNPs and indels). The VCFtools v0.1.15 ) was used for variant filtering. More precisely, SNPs and indels were filtered by mapping quality >50 and sequencing depth 10. The BreakDancer v1.3.6 (Chen et al., 2009) was applied to detect the SVs. Variant annotation was conducted with the Ensembl Variant Effect Predictor (VEP 4 ) (McLaren et al., 2016). The Sorting Intolerant from Tolerant (SIFT) was used for predicting whether an amino acid substitution affects protein function on the basis of sequence homology and the physical properties of amino acids (Sim et al., 2012).

Gene Functional Enrichment Analysis
Gene set enrichment analysis was performed at g: Profiler 5 , a web server for functional enrichment analysis and conversions of gene lists (Raudvere et al., 2019). For gene ontology (GO; biological processes) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, the goat genomic annotation was used. For phenotype ontology, since this information is not available for goat, we used human phenotype ontology (HP) instead to perform the analysis. The P-value generated from the test was automatically adjusted to account for multiple comparisons using the Benjamini and Hochberg correction (Benjamini and Hochberg, 1995). The false discovery rate (FDR) < 0.05 was required to determine the genes that are significantly overrepresented in those categories.

Principal Component Analysis and Detection of Selection Signature Genes
The SNP variant file was transformed to PLINK PED format via PLINK v1.9 (Purcell et al., 2007) to examine the genetic relationships with principal component analysis (PCA) among the 10 goat breeds.
We used two approaches include F ST and pooled heterozygosity (Hp) to identify the selection signature differences between the genomes of the CM and the other goat breeds. The weighted population pairwise F ST values were calculated for each SNP according to the formula of F ST = S 2 p(1−p)+S 2 /r , in which S 2 represents the sampling variance of allele frequencies between two breeds andp represents the overall average allele frequency across breeds (Weir and Cockerham, 1984;Guo et al., 2018). F ST values were further averaged over SNPs using a 100-kb sliding window with a 50-kb step size and Z-transformed as follows: where µF ST and σF ST are the mean and standard deviation for the F ST , respectively. This analysis was done with popoolation2 v1.201 (Kofler et al., 2011) and R package v3.5.1 (R Core Team, 2013).
The identification of selective sweeps was performed by calculating the ZHp for each pool and SNP. The numbers of reads corresponding to the major (n Maj ) and minor (n Min ) allele frequencies and for each window in each breed pool, the Hp score is calculated according to the formula of Hp = 2 n Maj n Min / n Maj + n Min 2 , where n Maj , and n Min are the sums of n Maj and n Min for all SNPs in the window. Individual Hp values were Z transformed as follows: ZHp = (Hp -µHp)/σHp, where µHp and σHp are the mean and standard deviation for the Hp, respectively. Consistent with the F ST , we calculated the ZHp value using a 100-kb sliding window with a 50-kb step size.
Windows with extremely high F ST values (99% percentile) and low ZHp scores (1% percentile) were proposed to be under selection. The genome-wide scan results of ZF ST and ZHp can be found in Supplementary Data S2, S3.

Read Mapping and Variant Identification
In this study, genome sequencing yielded a total of 652.25 Gb of raw data from the 10 goat breeds. Approximately 195.6 (Bor) to 253 (MT) million clean reads were obtained after quality control and filtration. More than 97.5% clean reads were mapped against to the reference genome ARS1 with 5× coverage of ∼94%. The average sequencing depth ranges from 14.24 × (YC) to 20.3 × (QS). These results indicated that high-quality sequences were obtained (Table 1).
A total of 23,508,551 SNPs and 2,830,800 indels were detected, out of which 12,939,052 SNPs and 1,571,878 indels were identified in the CM breed. The number of SNPs in the other nine breeds ranged from 11,871,653 (Bor) to 13,597,084 (MT). Similarly, the lowest number (1,486,903) of indels was detected in the Bor breed and the largest number (1,704,999) in the MT breed. By comparing the other nine breeds, we further identified 1.2% SNPs (271,713) and 0.9% indels (24,843) unique to the CM breed. The other nine breeds' specific variants ranged from 259,098 (MB) to 540,055 (QS) for SNPs, and 26,031 (NJ) to 59,999 (QS) for indels ( Table 2).
Based on the PCA results ( Figure 1B), the CM breed showed the closest relationship with the HM breed. The Bor breed demonstrated the most divergent genetic relationships with the CM breed. This is consistent with the fact that the Bor goat is originated from South African regions and was introduced domestically as an excellent commercial meat goat (Malan, 2000). The DE and NJ breeds are found in the Sichuan Province of China and are clustered together on the PCA map. This was also in line with our expectations.

Candidate Coding Functional Genes and Related Phenotypes in the Chongming White Goat Breed
The CM breed is our research of interest. Thus, we aimed to further investigate the phenotypes related to the CM-specific genes with LoF variants. In total, 351 LoF variants in 183 coding genes were identified (Supplementary Data S5). We performed the phenotype enrichment analysis for those genes on g: Profiler and obtained 338 HP terms with FDR < 0.05. Next, we filtered the genes with those involved in more than 50 HP terms. This resulted in 36 out of 183 genes. The top 50 phenotype terms related to those 36 genes are listed in Figure 5. Of note, all the genes participated in regulation of the nervous system (HP:0000707, FDR = 0.0002). Furthermore, most of the genes were associated with several other phenotype terms. For instance, 20 genes (e.g., SYNE2, AP4E1, and KBTBD13) were enriched in the "skeletal muscle morphology" (HP:0011805, FDR = 0.0016); 18 genes (e.g., AP4E1, and FSHR, and Col11a2) involved in the "short stature" (HP:0004322, FDR = 0.0025) and "body height" (HP:000002, FDR = 0.0025); and 15 (e.g., IFT172, NALCN, and IARS2) genes are found in the terms of "ear morphology" (HP:0031703, FDR = 0.001106) and "hearing" (HP:0000364, FDR = 0.00098), simultaneously.

Selection Signature Genes
Co-selection Signature Genes The CM goats are domesticated under the natural conditions of Chongming Island in Shanghai. In order to detect the genomewide selection signatures, we calculated F ST values between the CM and each of the other nine breeds using a 100-kb sliding widow with a step size of 50 kb. The 99% percentile Z score transformed F ST values for each comparison was calculated and averaged (ZF ST = 3.27). This value was used as the selection signature cutoff threshold. On the basis of this, we identified hundreds of genes that underwent selection, ranging from 228 between CM and YC to 330 between CM and DE. Of those, 37 genes were observed having co-selection signatures in more than four comparisons. In addition to F ST , we also calculated the ZHp scores for each breed and identified 16 out of the 37 co-selection genes with low ZHp scores < −3 (1% percentile, Figure 6). Of those, seven genes (C2CD3, DNAJB13, UCP2, ZMYND11, CEP126, SCAPER, and TSHR) are related to six nervous system-related terms (HP:0012758, HP:0012759, HP:0000707, HP:0012639, HP:0012638, and HP:0002011). Four genes (UCP2, UCP3, SCAPER, and TSHR) are related to the body weight (HP:0004324 and HP:0004323), and four genes (C2CD3, UCP2, ZMYND11, and TSHR) are involved in modulating the phenotype of globe developmental delay (HP:0001263).   Notably, both F ST and ZHp indicate that the 29.46-29.65 Mb region on Chr 15 has a strong selection signature in both the CM and HM breeds (Figure 7). This region contains five genes, with C2CD3, DNAJB13, and UCP2 involved in six phenotypes of the nervous system including abnormality of nervous system physiology (HP:0012638), abnormality of nervous system (HP:0000707), neurodevelopmental abnormality (HP:0012759), neurodevelopmental delay (HP:0012758), morphological abnormality of the central nervous system (HP:0002011), and abnormality of nervous system morphology (HP:0012639).

Unique Selection Signature Genes in Chongming White Goat Breed
Through pairwise F ST comparisons between the CM and the other breeds, we also identified several unique selection genes such as SLC22A15, KITLG, GALNT10, FGFR1, ERLIN2, and ZNF703 in CM vs Bor (Figure 7A). Among the notable ones are FGFR1, ERLIN2, and ZNF703 located on Chr 27 at 12-12.7 Mb, which is a unique strong selection signal with high F ST value (>7) and low ZHp scores (<-7.4) only observed in the Bor breed. Moreover, SDCCAG8 and STAG1 with high ZF ST values and low ZHp scores were detected in the CM and MB breeds. Five coat color-related genes (KITLG, ASIP, AHCY, RALY, and MC1R) were identified in the Bor, NJ, and MB breeds (Figure 7). The

Genetic and Genomic Studies of the Goat
In this study, we used high-throughput sequencing technologies to perform WGS for 10 goat breeds, pooling libraries from 100 goat individuals. After genome-wide mapping and variant calling, we obtained ∼23.5 million SNPs and 2.8 million indels from the 10 goat breeds. By calculating the distribution of genome-wide SNP polymorphisms in the 100-kb nonoverlapping window, we identified five common chromosomal regions with high polymorphism on goat genome. The results are consistent with those reported in previous goat genome study (Guo et al., 2018). Moreover, several breed-specific, highly polymorphic regions were identified. Most of the annotated genes in these genomic regions are uncharacterized LOC symbols, and the biological importance of them remains to be elucidated.

Genomic Characteristics of the Chongming White Goat Breed
Many indigenous breeds, which might represent highly valuable biodiversity resources, are threatened with extinction owing to their substitution by cosmopolitan breeds. The CM goat also faces this dilemma. In our study, we identified 271,713 CM-specific SNPs and 24,843 indels, which is the first whole-genome study to the best of our knowledge to characterize genetic polymorphisms in this breed. The results clearly demonstrate that the CM breed has the closest phylogenetic relationship with the HM and the most divergent relationships with the Bor breed. It has also been argued that both the CM and HM breeds belong to the Yangtze River Delta white goat species in China. The phylogenetic results among these breeds are also consistent with the expectation from the geographical distribution.
LoF variants are typically thought to have deleterious effects on protein-coding genes. The effects of these variants upon protein structure and expression in turn lead to phenotypic differences. In this study, we identified 183 genes with the CM breed-specific LoF variants. Out of these 183, 36 genes are enriched in more than 50 phenotype terms, FIGURE 6 | Phenotype overview of the 16 co-selection signature genes. Blue square indicates genes appearing in more than four of 10 breeds [false discovery rate (FDR) < 0.05]. Red squares indicate the nine genes enriched in the top 30 phenotype terms. including "nervous system, " "short stature, " and "skeletal muscle morphology" (Figure 5). The genes significantly enriched in the "nervous system" phenotypes, possibly reflecting the fact that the most critical phenotypic changes during the initial steps of animal domestication affect brain and neuronal development. Variants in these genes may produce behavioral traits that ease the domestication process or be a product of domestication itself (Carneiro et al., 2014).
The body size of the CM breed is relatively small (mature weights for bucks and does are 23-32 and 18-25 kg, respectively) compared with the Bor breed (mature weights for bucks and does are 90-130 and 80-100 kg, respectively (Lu, 2001). It is worth noting that 18 genes (e.g., AP4E1, FSHR, and COL11A2) are related to the terms of "short stature" and "body height." Mutations in the AP4E1 gene are associated with neurodevelopmental disorders in humans, which are characterized by intellectual disability, global developmental delay, and short stature (Jamra et al., 2011;Moreno-De-Luca et al., 2011). Moreover, Ap4e1 mutant mice display a lean body mass phenotype (Eppig et al., 2005). Follicle-stimulating Frontiers in Genetics | www.frontiersin.org hormone receptor (FSHR) plays a prominent role in mammalian reproduction (Papadimitriou et al., 2016), and ovarian follicle FSHR expression levels in prolific Lezhi black goats were reported to be significantly higher than those in non-prolific Tibetan goats (Zi et al., 2013). In mice, Fshr mutant females have body weights 20% greater than those of wild-type control females (Danilovich et al., 2000). Similarly, homozygote Col11a2 mutant mice are 25% smaller at birth than their wild-type counterparts (Eppig et al., 2005).
In addition, skeletal muscle is a significant contributor to overall body mass and an important structural component of the mammalian body (Berridge et al., 2018). Dysferlin (encoded by DYSF) is a skeletal muscle protein important for organization and structure of the sarcolemma (Spuler et al., 2008). DYSF mutations are associated with muscular dystrophy in humans (Bashir et al., 1998) and mice (Bittner et al., 1999). DYSF has also been identified as an alternatively spliced gene in goat leg muscle, likely playing crucial roles in skeletal muscle development (Xu et al., 2018). A study has revealed several candidate genes associated with goat ear morphology (Brito et al., 2017). In our study, 15 genes were related to the terms of "ear morphology" and "hearing" simultaneously. Therefore, these 36 genes may play important roles in the CM breed-specific phenotypic traits, and study of the biological underpinnings of these genes in relation to the traits is necessary and important.

Co-selection Signature Genes
Through calculated F ST values, we filtered 37 genes that have co-selection signals, 16 of which had ZHp scores < -3, which indicate that these genes have significant selection signatures with a high degree of confidence. Moreover, seven genes (C2CD3, DNAJB13, UCP2, ZMYND11, CEP126, SCAPER, and TSHR) were significantly enriched in six nervous system-related phenotypes (Figure 6). Some of them were reported to have clear functions in human or mouse studies. For example, ZMYND11 variants were identified in mice with anxiety-like behavior (Parker et al., 2016). Pathogenic variants in ZMYND11 have been associated with intellectual disability, behavioral abnormalities, and seizures in humans (Yates et al., 2020).
The region on Chr 15 at 29.46-29.65 Mb contains five genes (C2CD3, DNAJB13, UCP2, UCP3, and PAAF1), which is a strong selective signal region only in the CM and HM breeds. Given their close phylogenetic relationship, and the fact that three genes (C2CD3, DNAJB13, and UCP2) are implicated in six nervous system-related traits, we have reason to suppose that the nervous system in the CM and HM breeds may differ from that in the other goat breeds. Owing to the lack of sufficient phenotypic information on the nervous system of goats, the contribution of these genes is still unclear and remains to be further elucidated.
As mentioned above, it is of interest to identify genes important for goat growth and development and how these genes differ between breeds. Five selection signature genes (UCP2, UCP3, ZMYND11, SCAPER, and TSHR) are contributing to the phenotype term of "body weight" or "global developmental delay." Our results suggested that thyroid-stimulating hormone receptor (TSHR) is one of the co-selective signature genes between the CM and the DE, XH, NJ, and YC breeds (Figure 6). TSHR has a pivotal role in metabolic regulation and photoperiod control of reproduction in vertebrates. Studies had identified TSHR under selective pressure in other domestic animals, such as chicken (Rubin et al., 2010) and sheep . Accumulating evidence suggests that the uncoupling proteins UCP2 and UCP3 have a role in thermogenesis and energy efficiency (Nagy et al., 2004;de Almeida Brondani et al., 2014). Thus, these five genes, which we identified in this study, may play vital roles in the growth, performance, and metabolism of the goat, but the molecular mechanisms have not been fully elucidated yet.

Unique Selection Signature Genes in Chongming White Goat Breed
Some selection signatures were only identified between one and two breeds. There may be a connection between these unique selective signature genes and the traits of a specific breed. For example, Bor goats are excellent for meat production, grow quickly, and have high fertility rates. We detected a unique selective region located on Chr 27 at 12-12.7 Mb, which includes three genes (FGFR1, ERLIN2, and ZNF703) (Figures 7A,D). Genetic studies have placed the FGFR1 gene at the top of major ontogenic pathways that enable gastrulation, tissue development, and organogenesis (Terranova et al., 2015). Moreover, the genomic region that includes ERLIN2 and ZNF703 was reported to be associated with average daily gain (ADG) in Nellore cattle (Olivieri et al., 2016). Although the function of these unique genes and the relationship with the characteristics of the goat breeds are still unclear, these unique strongly selective genes that include SDCCAG8 and STAG1 in the CM and MB breeds are worthy of further study.

Coat Color Genes
Animal coat color is an important trait that is easy to observe. Studies have shown that some coat color genes have been underwent selections. We identified five coat color-related genes that have selection signatures in this study. The KITLG gene shows a strong selection signature in the Bor breed with unique high F ST and low ZHp scores (Figures 7A,D). KITLG encodes for the ligand of c-Kit, which is associated with coat color in pigs (Hadjiconstantouras et al., 2008) and UV-protective eye area pigmentation in cattle (Pausch et al., 2012). KITLG genomic region also reported has a significant effect on skin and hair color of humans (Miller et al., 2007;Sulem et al., 2007). RALY was reported to be associated with saddle tan and black-and-tan phenotypes in the dogs (Dreger et al., 2013). The RALY-ELF2S2-ASIP locus has influence on skin and hair pigmentation in the Nanjiang yellow goat (Guo et al., 2018). In the present study, we also detected the selective signal of RALY in the NJ breed, which further emphasizes that the RALY gene may contribute to the tan color of goats. Agouti signaling protein gene (ASIP) result in various coat patterns in domestic mammals 6 including mouse (Voisey and Van Daal, 2002), dog (Kerns et al., 2004), cat (Eizirik et al., 2003), rabbit (Fontanesi et al., 2010), horse (Rieder et al., 2001), sheep (Norris and Whan, 2008), and donkeys (Abitbol et al., 2015). It was also studied as a strong candidate gene associated with the black coat color in a Chinese Taihang black goat . In this study, only the MB breed has an all-black coat, and the results suggest that the ASIP and AHCY genes may contribute to the black coat color of the MB breed. For the MC1R gene, we observed the low ZHp values in the MB (black) and QS breeds (blue gray). The melanocortin one receptor (MC1R) is a melanocytic Gs protein-coupled receptor that regulates skin pigmentation, UV responses, and melanoma risk (Wolf Horrell et al., 2016). Studies have shown that this gene is important for regulation of eumelanin (black/brown) and phaeomelanin (red/yellow) syntheses in mammalian melanocytes . The coat color of animals may be determined by the combined action of several genes. These five genes are present in varying degrees of selective signatures in our 10 goat breeds, which highlighted that they play important roles in determining the coat color phenotypes of goats.

CONCLUSION
In this study, by using the WGS method and genomic approaches, we conducted a comparative genomic study for the first time to reveal the genomic characteristics and selection signatures between the CM and the other nine goat breeds. Our results identified 36 genes with the CM breed-specific LoF variants that were significantly enriched in the phenotype terms of "nervous system, " "short stature, " and "skeletal muscle morphology." From another perspective, hundreds of selection signature genes were identified in this study that are possibly associated with important traits, such as growth (UCP2, UCP3, ZMYND11, SCAPER, TSHR, FGFR1, ERLIN2, and ZNF703), nervous system (C2CD3, DNAJB13, UCP2, ZMYND11, CEP126, SCAPER, and TSHR), and coat color (KITLG, ASIP, AHCY, RALY, and M1CR). The genes with functional mutations and the selection signatures explored in this study both reflect the fact that the nervous system of the CM breed may differ from the other breeds. Although the trait-related genes require further validation through biological experiments, our findings provide an improved understanding for biodiversity conservation and trait exploration for the domestic goats.

DATA AVAILABILITY STATEMENT
Raw sequence data of 10 sheep breeds in the bam format from this study were deposited in the NCBI Short Read Archive (SRA; http://www.ncbi.nlm.nih.gov/sra) under accession number PRJNA608539.

ETHICS STATEMENT
The animal study was reviewed and approved by Experimental Animal Management Committee of SAAS.

AUTHOR CONTRIBUTIONS
JG and JD managed the project. FX and JG carried out the bioinformatics analyses. YL, DZ, FS, JY, CL, HL, and HY collected the samples. JG and FX wrote the manuscript. KR reviewed and edited the English writing. All authors read and approved the final manuscript.