Genome-Wide Selective Analysis of Boer Goat to Investigate the Dynamic Heredity Evolution under Different Stages

Simple Summary Boer goat is a world-famous meat-type goat breed that underwent long-term artificial selection from African indigenous animals. The current study displayed the genome-wide selection signature analyses of South African indigenous goat (AF), African Boer (BH), and Australian Boer (AS), to investigate the hereditary basis of artificial selection in different stages. Moreover, the θπ, FST, and XP-CLR methods were used to search for the candidate signatures of positive selection in Boer goats. Ten genes (e.g., BMPR1B, DNER, ITGAL, and KIT) under selection in both groups were identified and are potentially responsible for reproduction, metabolism, growth, and development. This study provided a comprehensive overview of genomic variations in Boer goat, which may provide a basis for further resource protection and improvement of this breed. Abstract Boer goats, as kemp in meat-type goats, are selected and bred from African indigenous goats under a long period of artificial selection. Their advantages in multiple economic traits, particularly their plump growth, have attracted worldwide attention. The current study displayed the genome-wide selection signature analyses of South African indigenous goat (AF), African Boer (BH), and Australian Boer (AS) to investigate the hereditary basis of artificial selection in different stages. Four methods (principal component analysis, nucleotide diversity, linkage disequilibrium decay, and neighbor-joining tree) implied the genomic diversity changes with different artificial selection intensities in Boer goats. In addition, the θπ, FST, and XP-CLR methods were used to search for the candidate signatures of positive selection in Boer goats. Consequently, 339 (BH vs. AF) and 295 (AS vs. BH) candidate genes were obtained from SNP data. Especially, 10 genes (e.g., BMPR1B, DNER, ITGAL, and KIT) under selection in both groups were identified. Functional annotation analysis revealed that these genes are potentially responsible for reproduction, metabolism, growth, and development. This study used genome-wide sequencing data to identify inheritance by artificial selection. The results of the current study are valuable for future molecular-assisted breeding and genetic improvement of goats.


Introduction
Following the advances of high-throughput sequencing technology, the scientific community has widely explored relative heredity basis with adaptability and economic phenotype of animals using genome wide selective analysis strategy (GWSA) [1][2][3]. In particular, a series of candidate genes related to domestication and artificial selection of domestic animal was determined. For example, the mutation genotype pattern located in the promoter region of PDGFD may be the cause of fat deposition in the tail of sheep [4]; the AHR gene related to female reproduction is under strong positive selection during the domestication of pigs [5]. An 8-kb sequence spanning the AMY2B locus showed signals of selection of key roles in starch digestion and fat metabolism [6].
Furthermore, many studies displayed the GWSA to investigate the key genes of environmental adaptability [7] and multiple economic traits [8,9] in goats, as well as their population migration [10], phylogeny, and domestic history [11]. For instance, a study found that MUC6 gene mutation in goats can improve the antiparasitic ability of their gastrointestinal tract to ensure that they have more adaptability to the human environment [12].
According to relevant reports, the indigenous goat accounted for more than 63% of total stock in South Africa, and it displayed a remarkable contribution to local livestock husbandry. Moreover, high-performance breeds under artificial selection appeared with the transformation of small-scale production systems to large-scale commercial agriculture [13]. Boer goat is a world-famous meat-type goat breed that underwent long-term artificial selection from African indigenous animals, and Australian Boer goats are introduced from Africa and further artificially selected [14]. To date, many countries have introduced Boer goats to improve local breeds and achieved excellent benefit. However, there is very little research literature on genomic changes in Boer goats. Therefore, the genetic basis of artificial selection in the improvement of economic traits can be elucidated by investigating the genomic genetic divergence between Boer goats and African indigenous goats with different stages.
In this study, we displayed the GWSA among African indigenous, South African Boer, and Australian Boer goats to identify the inheritance shake by artificial selection. Results were valuable for future molecular-assisted breeding and genetic improvement of goats.

Animals, Genome Sequencing, and Data Acquisition
The experiment was carried out in strict accordance with the guidelines of the International Cooperation Committee of Animal Welfare on the care and use of experimental animals.
Genomic DNA of 10 African Boer goats (BH) were supplied by Sokoine University of Agriculture (SUA) of Tanzania. Sequencing libraries of all samples were constructed by using NEBNext ® MLtra DNA library preparation kit (Illumina ® , 15,026,486 Rev. C, San Diego, CA, US) and whole genome sequencing was performed by Illumina NovaSeq 6000 × Ten platform (BGI, Shenzhen, China), the sequencing depth, for BH animals, was expected around 10×. The genome datasets of total 30 Australian Boer goats (AS) from our previous study [15] and 30 African indigenous goats (AF) were obtained from NCBI SRA database (PRJNA671542), respectively (Table S1).

Genome-Wide Selective Sweep Analysis and Gene Annotation
GWSA of SNPs were performed with two groups, as follows: (1) 10 African Boer goats (case) versus 30 African indigenous individuals (control); and (2) 30 Australian Boer goat (case) versus 10 African Boer goats (control). For the SNPs dataset, the pairwise fixation index (F ST ) and π ratio (π case /π control ) were calculated with 40 kb sliding windows and 20 kb step size using vcftools [16]. The cross-population composite likelihood ratio (XP-CLR) is a method for the determination of the selection signal, can avoid ascertainment biases and successfully detect older signals and the selections on standing variation, which means the allele frequency changed very rapidly due to random drift between two populations. Non-overlapping sliding windows of 50 kb were used, and the maximum number of SNPs within each window was 600 [17]. The top 5% was chosen as the significance threshold for XP-CLR. The variants were annotated with ANNOVAR (ANNOate VARiation), which can utilize annotation databases from annotation data set conforming to Generic Feature Format version 3 (GFF3). Candidate genes were annotated by the intersection of the common parameters with top 5% (F ST, πratio, XP-CLR). We can estimate the fraction of intronic, exonic, missense, stop-gained, etc. variants after annotated by Variant Effect Predictor (VEP). Additionally, candidate genes were subjected to gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) with KOBAS 3.0 (http://kobas.cbi. pku.edu.cn/ (accessed on 28 July 2021)) rely on human.

Results
Obtained from 70 animals were 18,321,865 SNPs. The frequency of most SNPs is located in the intergenic (66.79%) and intron (26.66%) regions. The Ti/TvRatio is 2.38. On the contrary, the least frequency of SNPs is distributed in the 3 UTR (0.32%) and 5 UTR (0.07%; Table S2).
The phylogenetic network revealed that all populations were divided into single linkage ( Figure 1A), which corresponded to that of the PCA pattern ( Figure 1B). Analysis result of LD (r 2 ) showed the highest mean LD and the slowest decay in the AS population. In contrast, the lowest mean LD and the fastest decay in the AF population ( Figure 1C) are shown. Coincidentally, the highest genome-wide nucleic acid diversity (ND) is AF in comparison with BH and AS ( Figure 1D). The GWSA result of group 1 (BH vs. AF) revealed that 339 candidate genes (Table  S3) were extracted from intersected windows by the top 5% of common parameters of FST, θπ rate, and XP-CLR ( Figure 2A). The GO analysis results (Table S4) showed that a high growth trait-related GO terms were identified, such as tissue development (GO:0009888), muscle fiber development (GO:0048747), and fibroblast growth factor binding (GO:0017134). In addition, several reproduction-related terms were also identified, such as positive regulation of intracellular estrogen receptor signaling pathway (GO:0033148), gonad development (GO:0008406), and sperm-egg recognition (GO:0035036). In addition, the KEGG analysis results (Table S4) showed that a significant number of genes are enriched in the top 10 (corrected p value < 0.01) pathways (i.e., transcriptional misregulation in cancer, pathways in cancer, glutamatergic synapse, oocyte meiosis, human T-cell leukemia virus 1 infection, Rap1 signaling pathway, Ras signaling pathway, and progesterone-mediated oocyte maturation). The GWSA result of group 1 (BH vs. AF) revealed that 339 candidate genes (Table S3) were extracted from intersected windows by the top 5% of common parameters of F ST , θπ rate, and XP-CLR ( Figure 2A). The GO analysis results (Table S4) showed that a high growth trait-related GO terms were identified, such as tissue development (GO:0009888), muscle fiber development (GO:0048747), and fibroblast growth factor binding (GO:0017134). In addition, several reproduction-related terms were also identified, such as positive regulation of intracellular estrogen receptor signaling pathway (GO:0033148), gonad development (GO:0008406), and sperm-egg recognition (GO:0035036). In addition, the KEGG analysis results (Table S4) showed that a significant number of genes are enriched in the top 10 (corrected p value < 0.01) pathways (i.e., transcriptional misregulation in cancer, pathways in cancer, glutamatergic synapse, oocyte meiosis, human T-cell leukemia virus 1 infection, Rap1 signaling pathway, Ras signaling pathway, and progesterone-mediated oocyte maturation).
Subsequently, the result of GWAS in group 2 (AS vs. BH) using the SNP dataset revealed that 295 candidate genes (Table S5) were obtained by the top 5% of three parameters ( Figure 2B) at the gene overlapping region. It was similar to the GO and KEGG (Table S6) enrichment in the rich multiple categories of group 1, such as those with well-known pathways (i.e., PI3K-Akt signaling pathway, Rap1 signaling pathway, MAPK signaling pathway, and Ras signaling pathway). Numerous GO terms were related to growth development (e.g., striated muscle tissue development, cardiac muscle contraction, and positive regulation of smooth muscle cell proliferation), reproduction (e.g., intracellular estrogen receptor signaling pathway, spermatogenesis, and male gonad development), immunity (e.g., T cell costimulation, T cell migration, and B cell apoptotic process), and metabolism (e.g., insulin secretion, glycogen catabolic process, and phosphate-containing compound metabolic process). Subsequently, the result of GWAS in group 2 (AS vs. BH) using the SNP dataset revealed that 295 candidate genes (Table S5) were obtained by the top 5% of three parameters ( Figure 2B) at the gene overlapping region. It was similar to the GO and KEGG (Table  S6) enrichment in the rich multiple categories of group 1, such as those with well-known pathways (i.e., PI3K-Akt signaling pathway, Rap1 signaling pathway, MAPK signaling pathway, and Ras signaling pathway). Numerous GO terms were related to growth development (e.g., striated muscle tissue development, cardiac muscle contraction, and positive regulation of smooth muscle cell proliferation), reproduction (e.g., intracellular estrogen receptor signaling pathway, spermatogenesis, and male gonad development), immunity (e.g., T cell costimulation, T cell migration, and B cell apoptotic process), and metabolism (e.g., insulin secretion, glycogen catabolic process, and phosphate-containing compound metabolic process).
Interestingly, outstanding visible divergence in digestion and metabolism-related KEGG pathways was found according to the results of the gene functional enrichment of selected genes at two different stages (Table 1). Notably, 10 genes (BMPR1B, DNER, EPHA6, ERCC6L2, ITGAL, KIT, PROM1, RBFOX2, RFX3, and U6) under selection in both groups were identified.  Interestingly, outstanding visible divergence in digestion and metabolism-related KEGG pathways was found according to the results of the gene functional enrichment of selected genes at two different stages (Table 1). Notably, 10 genes (BMPR1B, DNER, EPHA6, ERCC6L2, ITGAL, KIT, PROM1, RBFOX2, RFX3, and U6) under selection in both groups were identified.

Discussion
In this study, the ND of both Boer goat populations was lower than African indigenous goats. The higher LD linkage within the AS population indicated that the frequency of genotype may be gradually changed at different stages as humans expected during the artificial selection process [18][19][20]. The phylogenetic and PCA pattern relationship of the three populations was inconsistent with their geographic location, which implied that artificial selection may contribute to the changes in animal genome diversity greatly and rapidly.
Boer goat was initially selected and bred from AF goats and has a superior economic phenotype and fecundity [14,21]. Interestingly, 10 genes were identified to undergo selection in both breeding stages of the Boer goat, which was widely confirmed to be involved in various biological functions (e.g., reproduction, growth, and metabolism).
Despite a series of metabolism-related (MRD) genes found in the two stages of BH and AS, the genes and their MRD pathways were different. For example, in terms of amino acid metabolism, six pathways related to amino acid metabolism were enriched and involved in the digestion and absorption of 17 essential and nonessential amino acids in the AF to BH stages. In the BH to AS stage, only three MRD pathways were enriched including eight kinds of amino acids. The AF feeding conditions are generally single and backward [22]. In addition, BH has undergone artificial selection and breeding in the commercial system with improved living conditions and medical security. Corresponding changes have been made in the environmental adaptation, metabolism, and other hereditary traits of Boer goats. Furthermore, the differences in the ecological environment of the habitat, breeding standards, and improved raising level in AS led to the rapid coevolution of the metabolism and environment of Boer goats. This was also confirmed by a large number of candidate genes related to the metabolism of fatty acids, minerals, and vitamins, as identified in AS vs. BH. Thus, the coevolution between the metabolic inheritance and the integrated living environment in animals has been universally confirmed [23][24][25].
Particularly, some candidate genes related to metabolism have been attracting attention. Ectonucleotide pyrophosphatase/phosphodiesterase3 (ENPP3), as a regulator of N-acetylglucosaminyltransferase GnT-IX (GnT-Vb), would have widespread and significant impacts on glycosyltransferase activities [26]. Short-chain enoyl-CoA hydratase (ECHS1) is involved in amino and fatty acid catabolism in mitochondria and its deficiency results in Leigh syndrome or exercise-induced dystonia [27]. Niemann-Pick C1-Like 1 (NPC1L1) was widely verified to be a cholesterol import protein and mediated intestinal cholesterol absorption [28]. It also participated in the intestinal absorption of fat-soluble vitamins [29]. In addition, genes that have been continuously selected in terms of immune resistance exist due to changes in the living environment. The integrin LFA-1 (also known as ITGAL), also known as CD11a, is upregulated in metastatic melanoma, highly expressed in most immune cell populations [30]. ITGAL is closely linked to the pathogenesis of diverse immune-related diseases [31]. Therefore, the capture of these genes implied that the metabolic and physiological functions of Boer goats are continuously affected by the selective breeding process.
A series of reproduction-related genes were identified in different stages. For example, many studies have shown that BMPR1B is associated with litter size in animals [32][33][34]. It is the main gene that affects the ovulation rate and litter size in sheep [35][36][37]. In addition, a series of studies have confirmed that the BMPR1B can affect follicular development, ovulation [38], and cause ovarian insufficiency [39]. These genes may directly affect the fertility of Boer goats under strong artificial selective pressure at different stages of Boer goat. In addition, KIT is verified to be involved in mammalian oocyte growth and follicular development [40]. A lack of KIT can lead to the loss of oocytes. In the ovaries of rodents, KIT displayed an important role in the migration, proliferation, and survival of primordial germ cells [41].
Among melanin-related genes, KIT is pivotal in the melanogenesis signaling pathway, and mutations or deletion of KIT can cause different hair and skin colors in mammals [42,43].
For example, a single A-G base missense mutation in exon 13 has been reported to cause differential expression of KIT in Liaoning cashmere goats, resulting in different coat colors [44]. This may be the genomic evidence for the changes in coat color in Boer goats at different stages.
Subsequently, excellent growth performance is one of the important goals of Boer goat breeding. Therefore, investigating the selective genes between different breeding stages of Boer goats may help elucidate the genetic basis of muscle development [15]. Notably, some genes identified in this study attracted the attention of the authors. For example, as an essential for maintaining skeletal muscle mass and protein homeostasis, the RBFOX2 was confirmed to control the fusion of myoblasts in myogenesis by coordinating the alternative splicing of Mef2d and Rock2 [45]. Notch activators (EDNR) are putatively involved in the interconnected signaling networks that control satellite cell function [46]. Studies have shown that EDNR gene expression is closely related to marbling levels in larger beef cattle populations [47]. In addition, it has also been shown that the delivery of a negative estrogen receptor gene (EDNR) abrogates estrogen-and progesterone-regulated gene expression [48].
Finally, the evolution process of animals started from a very distant era, some of the inferior traits of these animals were eliminated, and the superior traits were improved by humans to provide a stable breeding state according to their requirements [49]. Artificial selection resulted in outstanding physiologically and phenotypically changes in these animals [49,50]. This is powerful proof of the coevolution between animal hereditary and feeding environment.
Of course, this study still has a limitation, e.g., the sample size of African Boer goats (10) is observably less than other sample groups, which is also a common flaw in a large number of previous related studies [51][52][53] due to the high sequencing cost. However, high sequencing depth with tenfold genome coverage of each individual and strict quality control conditions were performed in this study can ensure the detection rate of genetic variation enrichment and genotype accuracy of the whole genome.

Conclusions
This study provided a comprehensive overview of genomic variations in Boer goat genomes under different stages. The characterization of population structure and genomic diversity will point out the direction for genetic assessment and the development for reasonable breeding strategies of Boer goats. Moreover, a series of candidate genes that may be important for the meat quality traits as well growth and fertility of this breed were identified. In addition, some genes also implied that the metabolic and physiological functions of Boer goats are continuously affected by the selective breeding process. These results provide a basis for further resource protection and breeding improvement of this breed.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ani12111356/s1, Table S1: Basic information and data readsdepth of 70 individuals; Table S2: SNP mutation types and numbers from 70 goats individuals; Table  S3: 339 candidate genes were extracted from top 5% of common parameters of F ST , θπ rate and XP-CLR from BH vs. AF; Table S4: Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of 339 candidate genes from top 5% of common parameters of F ST , θπ rate and XP-CLR from BH vs. AF based on SNP datas; Table S5: 295 candidate genes were extracted from top 5% of common parameters of F ST , θπ rate and XP-CLR from AS vs. BH; Table S6: Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of 295 candidate genes from top 5% of common parameters of F ST , θπ rate and XP-CLR from AS vs. BH based on SNP datas.

Data Availability Statement:
The sequence data has been deposited at NCBI with the accession number PRJNA671542, PRJNA770516.

Conflicts of Interest:
The authors declare that they have no known competing financial interest or personal relationship that could have appeared to influence the work reported in this paper.