This study collected two populations of Xinong Saanen Dairy Goat (XNG, n = 40) and Guanzhong Dairy Goat (GZG, n = 10) from Shaanxi Province. Each goat had 5 mL of blood collected from the jugular vein, and DNA was extracted using the standard phenol-chloroform method. Samples with determined DNA concentrations were subjected to whole-genome sequencing by Huada Company, with paired-end libraries constructed using the Huada T7 platform having an average insert size of 500 bp per individual and an average read length of 150 bp. In addition, the resequencing data of 84 individuals were downloaded from public databases, including 8 breeds as follows: Australian Alpine Dairy Goat (AAG, n = 2), Australian Saanen Dairy Goat (ASG, n = 6), Guanzhong Dairy Goat (GZG, n = 26), New Zealand Saanen Dairy Goat (NSG, n = 6), New Zealand Alpine Dairy Goat (NAG, n = 4), Laoshan Dairy Goat (LSG, n = 9), Nubian Dairy Goat (NBY, n = 15), Tugenburg Dairy Goat (TGB, n = 16). Resulting in a total of 134 samples in the study.

The study utilized Trimmomatic v0.38 (21) for filtering the paired-end sequences. Next, BWA-MEM (v0.7.15-r1140) (22) aligned the clean data with the goat reference genome ARS1.2 (GCF_001704415.2, https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_001704415.2/). Subsequently, samtools (23) was utilized to construct a BAM file index for mapping. The BAM files were then sorted and potential duplicate reads were removed utilizing the Picard¹ tool. After mapping, SNP calling was performed using the “Haplotype Caller,” “Genotype GVCFs,” and “Select Variants” modules in the GATK genomic analysis tool package (GATK, version 3.8-1-0-gf15c1c3ef) (24). The initial SNPs were filtered with the “Variant Filtration” module based on the parameters: “QD < 2.0, FS > 60.0, MQ < 40.0, MQRankSum < −12.5, ReadPosRankSum < −8.0, and SOR > 3.0,” and an average sequencing depth of variants within the range “<1/3× and >3×” for all individuals. Lastly, the ANNOVAR (25) software was utilized for functional annotation of the SNPs.

After applying linkage disequilibrium (LD) to prune SNPs, the PLINK (26) program is utilized to compute runs of homozygosity (ROH) with specific parameters as follows: (a) 100 consecutive homozygous SNPs, (b) minimum of 50 SNPs per window, (c) 500 kb of contiguous homozygous length, (d) minimum density of 1 SNP per 50 kb, (e) window overlap rate of 0.05, and (f) each window containing 1 heterozygous and 2 missing calls. The subsequent analysis results are categorized into 0.5–1 Mb, 1–2 Mb, and 2–4 Mb for visualization. The heterozygosity of the SNPs is calculated to estimate the inbreeding coefficient (Fhom). VCFtools (27) are used to analyze nucleotide diversity for each breed in non-overlapping 50 kb windows. Next, the PopLDdecay (28) software is used to calculate the decay of LD to assess the physical distance between SNPs within haplotype blocks of different breeds. Finally, the frequencies of runs of homozygosity (ROH) islands for Xinong Saanen dairy goats and Guanzhong dairy goats are calculated, and regions with frequencies above 20% are determined as ROH hotspot islands.

Three methods are used for Variant Call Format (VCF) filtering and population structure estimation: (a) constructing a neighbor-joining tree with MEGA v10.2.6 (29) and visualizing it using iTOL (30); (b) principal component analysis (PCA) with EIGENSOFT v5.0 (31) software; and (c) conducting population structure analysis using ADMIXTURE v1.3.0 (32). Cross-validation is employed to calculate the cross-validation error and determine the optimal K value (assuming the ancestral population K value ranges from 2 to 8).

Based on the characteristics and origins of the dairy goat breeds, we categorized 134 dairy goats into native dairy goat (NDG) and imported dairy goat (IDG) breeds. We then compared the genomes of these two dairy goat breeds and estimated the signal scanning regions using a combination of nucleotide diversity (θπ NDG/IDG) and fixation index (FST) with a VCFtools, employing 50 kb sliding windows and 25 kb sliding steps. Additionally, Tajima’s D (33) statistic and cross-population composite likelihood ratio test (XP-CLR) (34) were employed to identify potential regional differences between different breeds. XP-CLR is a likelihood method for detecting selective sweeps by jointly modeling multi-locus allele frequency differences between two groups. The scanned regions detected from the intersection of the two parameters with the highest 5% threshold were annotated to identify candidate genes. Lastly, Bedtools (35) was used to annotate the selected regions for subsequent analysis.

The enrichment module in KOBAS (36) was utilized to identify pathways and Gene Ontology (GO) terms that showed statistically significant associations (at a significance level of p ≤ 0.05). Subsequently, the results were visualized through the OmicShare tool.²

To identify the genomic difference between local (XNG, GZG, LSG) and imported breeds of dairy goats (AAG, ASG, NSG, NAG, TGB, NBY), we compared 134 sets of whole-genome resequencing data, including 40 genomes of XNG and 10 genomes of GZG, combined with 84 published genome data that includes both local and imported breeds (GZG, LSG, AAG, ASG, NSG, NAG, TGB, NBY) (Figure 1). The data was aligned to the goat reference genome ARS1.2 (GCF_001704415), achieving an average alignment rate of 97.74% and a mean sequencing depth of 11.3×. After the removal of low-quality sequences, the average clean reads for each sample were 373,568,004 bp, with GC content ranging from 41.82 to 42.33%. The average alignment rate for the studied populations (XNG, GZG) is 97.5%, accompanied by an average sequencing depth of 17.01×. The imported breed populations (AAG, ASG, NSG, NAG, TGB) exhibited an average alignment rate of 97.04% and an average sequencing depth of 22.51× (Supplementary Table S1).

After SNP variation detection and statistical analysis of the obtained 134 dairy goat genome data sets, a total of 33,011,806 SNPs were detected. Among them, the proportion of intronic SNPs was the highest (32,705,301), accounting for 55.04% of the total; the proportion of exonic SNPs was approximately 1.142%. There were 283,685 synonymous SNPs and 227,235 missense SNPs. The transition/transversion ratio (Ti/Tv) was computed and analyzed, revealing a ratio of 2.344, consistent with the general value for goat populations. These findings indicate a high quality of sequencing, thus the obtained data meet the requirements for subsequent analysis (Supplementary Table S2).

The LD analysis shows that GZG and XNG exhibit the fastest decay rate (Figure 2A), suggesting lower domestication levels and higher genetic diversity. Nucleotide diversity analysis (Figure 2B) reveals higher diversity in local breeds compared to imported ones. ROH is an indicator of kinship proximity and inbreeding history, and the analysis reveals the genomic patterns of recent demographic history. Short ROH indicates that ancient inbreeding is much more prevalent in the XNG and GZG compared to other imported breeds. In addition, XNG exhibits a higher value of ROH, indicating a longer breeding history (Figure 2C). Meanwhile, the coefficient of inbreeding due to recent inbreeding is highest for each individual genome in TGB, with a range from 0.1 to 0.51 (Figure 2D).

This study conducted ROH island analysis on XNG and GZG. In the XNG samples, 16 ROH islands with a frequency of 20% were identified on chromosomes 12, 13, and 19 (Figure 3A), annotating to 31 genes. Similarly, on chromosomes 1, 7, 12, 13, and 18 of the GZG, 15 ROH islands with a frequency of 20% were identified (Figure 3B), annotating 25 genes. The Venn diagram (Supplementary Figure S1) illustrates 4 common genes between the two breeds (MPHOSPH8, LOC108637296, PSPC1, NBEA). In addition, the GO and KEGG enrichment analyses (Supplementary Figure S2) indicate enrichment of those identified genes in 23 biological processes (BP), 9 molecular functions (MF), and 2 cellular components (CC). The enriched pathways primarily include biological regulation, metabolic processes, immune system processes, and transcriptional regulation. Furthermore, the KEGG analysis reveals enrichment in valine, leucine, and isoleucine degradation, apoptosis, PI3K-Akt signaling pathway, and Rap1 signaling pathway.

To comprehend the population structure of local dairy goats in China, we looked into the SNP results derived from the resequencing data. Employing imported breeds as outgroups, we constructed SNP differentiation clustering based on various populations, the results of which can be used to mutually verify with other clustering methods. The percentages of PC1, PC2, and PC3 are 7.91, 5.61, and 3.60%, respectively, (Figure 4A; Supplementary Figure S3). There are three major clusters that AAG, ASG, LSG, NAG, NSG, and TGB are closely clustered, XNG and GZG form another cluster, while NBY is away from both of them. The result indicates that local breeds can be differentiated from imported breeds, and there are differences in clustering patterns among the 9 populations. The genetic data of the populations were utilized to calculate the degree of kinship between individuals, construct a genetic distance matrix, and build a phylogenetic tree utilizing the distance matrix. Although GZG and XNG are partially overlapped in the phylogenetic tree, the XNG group is clustered on its own, and imported breeds can be differentiated from local breeds, which is consistent with the PCA results (Figure 4B). In addition, for the 9 population, the number of subgroups (value of K) was preset to 2–9 for clustering, and the clustering results were cross-validated. The optimal number of subgroups was identified based on the minimum cross-validation error rate. K = 4 was identified as the optimal number of subgroups, and there were significant genetic differences among the XNG, NBY, NAG, and TGB populations. Regardless of the value of K, it consistently displayed genetic differentiation between local and imported breeds, indicating different genetic compositions among the populations, and relatively consistent genetic backgrounds among the different 9 populations (Figure 4C).

To identify the regions of selection among the native dairy goat breeds and imported dairy goat breeds, we further conducted a comparative analysis using the resequencing data. Next, using the top 5% of both FST (Figure 5A) and θπ ratio (Figure 5B) cutoffs (Figure 5C), we observed that the native dairy goat breeds exhibited a total of 41.25 Mb across 607 selectively scanned regions. These regions encompassed 549 genes (Figure 5D) and showed significant enrichment in pathways such as the calcium signaling pathway, MAPK signaling pathway, fatty acid metabolism, and the PI3K-Akt signaling pathway (Figure 5E). On the other hand, imported dairy goat breeds displayed 25.05 Mb in 385 selective scanning regions, which encompassed 387 genes (Figure 5D). Pathway enrichment analysis revealed that those genes are involved in oocyte meiosis, cell senescence, calcium signaling pathway, GnRH signaling pathway, and the synthesis and secretion of growth hormone pathway (Figure 5F). These findings emphasize the genetic functional differences between local and imported breeds (Supplementary Figure S4).

To further explore the traits corresponding to selection signals, we performed FST, PI, and XP-CLR assays to discover positive selection in dairy goats. Meanwhile, for a more comprehensive annotation of genes We integrated and annotated the QTLdb database and a previously published genome-wide association analysis on important traits of dairy goats into our selection signaling results (Figure 6A). We found a large number of genes associated with milk production traits on chromosomes 5, 13, and 19 (Figure 6A), which were mainly involved in the monoacylglycerol biosynthetic process (GO:0006640), the long-chain fatty-acyl-CoA metabolic process (GO:00353336), monoacylglycerol metabolic process (GO:0046462), and fatty acid homeostasis (GO:0055089) (Figure 6B). In addition, these genes were also significantly enriched in the PI3K-Akt signaling pathway, glycerolipid metabolism, Fat digestion and absorption, and Regulation of lipolysis in adipocyte signaling pathways (Supplementary Figure S5). Integrative analysis (Fst, PI, XP-CLR) revealed two overlapping genes (Supplementary Figure S6), ASIP and STK3, respectively. Notably, integrative analysis (Fst, PI, XP-CLR, ROH_island) of the genes obtained from the above analytical methods revealed seven overlapping genes, including PRELID3B, ATP5E, CTSZ, NELFCD, LOC108637417, TRNAF-GAA, and TRNAN-GUU (Figure 6C). We then analyzed the allele frequency distribution of selected regions for STK3, GHR, and PRELID3B since many methods have detected these three genes, which showed that the gene frequencies of the selected genes differed between the local and imported breeds (Figure 6D; Supplementary Figure S7). Through the selection feature, we found that some regions were strongly selected. Some candidate genes were associated with milk production traits, growth and development, and immune traits.