In the experiment, the breeding environment complied with the relevant standards for general animal experimental facilities in the Chinese national standard “Laboratory Animal Environment and Facilities” (GB 14925–2023). The feeding and experimental operations of the animals met the requirements of animal welfare. In this study, no anesthesia or euthanasia was performed on experimental animals.

The Inner Mongolia cashmere goats (n = 461) used in this study were all obtained from Inner Mongolia Yiwei White Cashmere Goat Co., Ltd. All the goats were reared under consistent feeding environments and nutritional conditions, being provided with the same commercial diet and having unrestricted access to water. All Inner Mongolia cashmere goats used in this study were born from 2010 to 2018, and all of the data were recorded in detail. Before weaning, some nutritional supplements were provided in addition to the ewe’s milk. The lambs were weaned when they were 4 months old. Birth weight (BW), weaning weight (WW) and weight at 12 months (YW) were measured using electronic scales in the same environment. Additionally, the weaning date was recorded. Average daily gain before weaning (ADG) was calculated as (weaning weight - birth weight) /days to weaning. BW was measured 0.5 h after birth, and WW and YW were measured 12 h after fasting (2). Ear tissue samples were carefully collected using ear-notch forceps. Immediately after collection, they were swiftly transferred into pre-prepared cryotubes filled with 75% alcohol. Subsequently, these cryotubes were stored at −80°C until DNA extraction was carried out.

Total Genomic DNA was extracted from ear tissue using the traditional phenol-chloroform method, yielding a concentration of 50 ng/μL. The quality of DNA in all samples (461 DNA samples) was evaluated based on light absorption ratios (A260 / 280 and A260/230) and gel electrophoresis (18). Following DNA extraction, qualified samples were fragmented using a Covaris crusher. The DNA fragments were then end-repaired, polyA-tailed, ligated with sequencing adapters, purified, and PCR-amplified to construct the sequencing library. The library was preliminarily quantified using Qubit 3.0, with insert size distribution verified by the Agilent 5,300 system. After passing quality control, the library was subjected to paired-end sequencing (PE150) on the DNBSEQ-T7 platform.

The raw data were disconnected, filtered and subjected to various quality control steps to obtain clean data for subsequent bioinformatics analysis. In order to ensure the accuracy of the analysis, FastQC v0.11.5 was used to perform strict quality inspection on the original reading length (19). The original readings were then compared with the Capra hircus reference genome (Genome assembly ASM4082201v1)¹ using the Burrows-Wheeler Aligner (BWA aln) v0.7.8 with default parameters (20). On average, 99% of the reads were successfully mapped to the reference genome, achieving a final average sequencing coverage of 9.46 × (ranging from 7 × to 12×) per individual. The initial BAM files, containing sequence alignment data for each individual, were generated using BWA and subsequently sorted and indexed using SAMtools v1.0 (21). Then the genome features such as GC content, repeat and gap content, read count and absolute copy number were calculated using the sliding window method (1 kb window, 800 bp step) (22). For the initial CNVnator output, quality control was performed by filtering copy number variations (CNVs) based on length (50 bp to 5 Mb). Further CNVs filtering was performed using Plink v1.90 beta to remove CNVs with minor allele frequency (MAF < 0.01). CNV regions (CNVRs) were defined by merging quality-controlled CNVs with ≥1 bp overlap at identical genomic positions across all samples using HandyCNV (R v4.4.1). HandyCNV defines CNVR as three types: amplification, deletion and mixing (when deletion and amplification occur in the same region) (23).

The general linear model was used to perform the GWAS analysis based on CNV for early growth traits in IMCGs with Plink v1.90 software. The formula for this model is as follows:

Where y is the phenotype vector, α is the fixed effect vector, the fixed effects included herds, sex, maternal age, and year of measurement. X is the structure matrix of fixed effect. β is CNV effect, Z is the structure matrix of CNV effect; e is the residual effect, and the distribution is e ~ N ( 0 , σ 2 ) .

In the CNV-based GWAS, the Bonferroni method was used to determine the genome-wide significant (0.05/N) threshold, where N represents the number of CNVs. Given that is a stringent criterion, a more lenient threshold was also used for detecting the suggestive (1/N) CNVs (24, 25). The qqman package in R software was used to plot the Manhattan and Q-Q plots (26).

Eight candidate CNVs associated with early growth traits in IMCGs were selected for validation, including two for BW traits (CNV_DEL_17406 and CNV_DEL_18821), three for ADG and WW traits (CNV_DEL_11189, CNV_DEL_17895, and CNV_DUP_18956), and three for YW traits (CNV_DEL_4359, CNV_DEL_4552, and CNV_DUP_20170). These CNVs were subsequently verified using real-time quantitative polymerase chain reaction (qPCR). The primers for these eight CNVs were designed using Primer Premier 5 software (Supplementary Table S1). Following the methodology established by Sonika et al. (27), the ACTB were selected as the reference gene, because the gene is highly conserved in goats and exists in the reference genome in the form of a single copy. A total of 64 DNA samples were selected for verification, including 32 samples containing the target CNV, and 32 normal samples without copy number variation in the test area were used as controls. The qPCR experiment was conducted using the 2 × SG Green qPCR Mix (with ROX; SINOGENE, Beijing, China). The PCR reaction was performed in a total volume of 15 μL, containing 1 μL DNA template (50 ng/μL), 0.25 μL each of forward and reverse primers (10 pM/μL), 5 μL of 2 × Blue-SYBR Green mixture, and 6 μL of nuclease-free water. The thermal cycling protocol consisted of an initial denaturation at 95°C for 10 min, followed by 40 cycles of amplification (95°C for 20 s, 60°C for 30 s), and a final dissociation curve analysis (95°C for 15 s, 60°C for 30 s, and 95°C for 15 s). All reactions were carried out in triplicate on a 96-well transparent reaction plate, and the average cycle threshold (Ct) values were calculated for subsequent copy number analysis. The relative copy number variation in the target region was determined using the 2(^-ΔΔCt) method, where ΔΔCt was calculated as follows: [(mean Ct of target gene in test sample)–(mean Ct of reference gene (GCG) in test sample)]–[(mean Ct of target gene in reference sample)–(mean Ct of GCG in reference sample)]. Based on this calculation, a value approximating 2 was considered normal, while values ≥ 3 and ≤ 1 indicated copy number gain and loss, respectively.

After using qPCR to verify CNVs, the phenotypic data of individuals with significant CNVs in GWAS analysis were tested by t-test to verify the authenticity of GWAS results.

The physical position information was retrieved from the Capra hircus reference genome (Genome assembly ASM4082201v1, see Footnote 1). Candidate genes located within 300 kb upstream and downstream of the significant CNVs were identified using the bedtools software (28). Subsequently, the overlapping genes were subjected to enrichment analysis using the Database for Annotation, Visualization, and Integrated Discovery (DAVID) tool, which included Gene Ontology (GO) functional annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. Enrichment terms with statistical significance (p < 0.05), determined by Fisher’s exact test, were selected to further investigate genes involved in relevant biological pathways and processes (17, 29). Gene functions were queried using Ensembl Biomart.²

In this study, the descriptive statistics of our early growth traits, including birth weight (BW), weaning weight (WW), average daily gain (ADG), and yearling weight (YW), were summarized in Table 1. The mean values (±SD) were determined as 2.66 ± 0.45 kg for BW, 21.90 ± 3.28 kg for WW, 0.16 ± 0.03 kg for ADG, and 31.21 ± 3.60 kg for YW. The coefficients of variation (CV) for the four traits were calculated at 17.10, 14.98, 18.75, and 11.53%, respectively. Prior to statistical analysis, data quality control were implemented, including the exclusion of phenotypic records with missing values and the removal of outliers exceeding the threshold of mean ± 3 standard deviations. Consequently, the sample sizes presented in Table 1 for each trait are slightly reduced from the initial cohort of 461 individuals. Furthermore, normality assessment through statistical tests confirmed that all traits followed a normal distribution pattern (Figure 1), validating the suitability of parametric statistical methods for subsequent analyses.

Following quality control procedures, 11.87 TB of high-quality reads were retained from the initial 13.41 TB raw sequencing data obtained from 461 Inner Mongolian Cashmere Goats (IMCGs). Genome-wide CNV detection was performed across all 29 autosomes using CNVnator (v0.4.1). The initial analysis identified 315,891 CNV events, categorized as either losses (copy number = 0 or 1) or gains (copy number ≥ 2). To ensure data reliability, stringent filtering criteria were implemented. CNVs with population frequencies below 0.01 (observed in fewer than 5 individuals) were excluded as potential false positives. Furthermore, overlapping CNVs detected at identical genomic positions across multiple individuals were merged and systematically reannotated. This rigorous filtering process yielded a final set of 26,003 non-redundant CNVs, comprising 10,216 gain-type and 15,787 loss-type events (Figures 2a,b), each present in at least five individuals. The cumulative length of identified CNVs spanned 3,469.5 Mb, exceeding the total length of the goat autosomes (2,605.7 Mb). This apparent discrepancy arises from the overlapping nature of CNVs across individuals during the detection process, preventing direct calculation of chromosomal coverage. To address this, we generated two complementary visual representations: illustrates the density distribution of CNVs across chromosomes was shown in Figure 2c, while Figure 2d presents CNV coverage using a 1 Mb sliding window approach along the chromosomes.

The genomic distribution and characteristics of CNVs across the 29 autosomes are summarized in Table 2. On average, each chromosome contained 897 CNVs, consisting of 352 gain-type and 545 loss-type events. While CNV distribution generally correlated with chromosome length, notable exceptions to this trend were observed. Chromosome 17 exhibited the highest CNV count (1,638). In contrast, chromosome 26 showed the lowest CNV frequency (359).

The identified CNVs exhibited a size range from 0.8 kb to 4.99 Mb, with a mean length of 133 kb. Specifically, the largest CNV identified was CNV_DUP_13234 located on chromosome 16, spanning 4.99 Mb, while the smallest was CNV_DEL_7892 on chromosome 10, measuring 0.8 kb.

Of the total 26,003 CNVs detected, gain-type CNVs accounted for 39.3% (n = 10,216), ranging in size from 4.8 kb to 4.99 Mb. These gain events collectively spanned 2,285.9 Mb, with an average length of 224 kb. In contrast, loss-type CNVs represented 60.7% (n = 15,787), ranging from 0.8 kb to 4.85 Mb in length. The cumulative length of deletion CNVs was 1,183.5 Mb, with an average size of 75 kb.

The CNV distribution patterns across chromosomes in 461 IMCG individuals are illustrated in Figure 3, where copy numbers of 0 or 1 represent losses and 3 or 4 represent gains. The individual CNV count distribution, revealing that the top three individuals contained 2,097, 2,004, and 1,955 CNVs respectively, while the bottom three individuals exhibited only 43, 34, and 31 CNVs. The average population was 685 CNVs per individual (Figure 3a). The size distribution of CNVs across different copy number categories was shown in Figure 3b. The majority of CNVs (less than 0.05 Mb) were significantly more abundant than larger CNVs (2–5 Mb). Loss-type CNVs (copy number = 0 or 1) substantially outnumbered gain-type CNVs (copy number = 3 or 4), with most CNVs concentrated below 0.3 Mb in length.

Chromosomal distribution patterns are shown in Figure 3c. Loss-type CNVs (copy number = 0) were most prevalent on chromosome 18 (9,155 events) and least frequent on chromosome 26 (2,925 events). Similarly, CNVs with copy number = 1 showed maximal distribution on chromosome 13 (22,394 events) and minimal on chromosome 26 (2,641 events). For gain-type CNVs, those with copy number = 3 were most abundant on chromosome 10 (3,448 events) and least on chromosome 26 (183 events), while copy number = 4 CNVs peaked on chromosome 1 (3,520 events) and were rarest on chromosome 14 (32 events). This chromosomal distribution analysis confirms the relatively low CNV frequency on chromosome 26, consistent with the merged and reannotated CNV data presented in Table 2. The observed patterns suggest potential chromosome-specific mechanisms influencing CNV formation and maintenance in IMCGs.

Through integration of the 26,003 quality-controlled CNVs, we identified 5,014 CNV regions (CNVRs; Supplementary Table 2) by merging adjacent CNVs with overlapping regions exceeding 1 bp. These CNVRs spanned all autosomes, ranging from 75 to 344 regions per chromosome, with a cumulative length of 1,015.43 Mb, representing 38.97% of the goat autosomal genome.

The CNVRs were classified into three categories: 1,085 gain-type CNVRs with a total length of 175 Mb (average length: 161.3 kb), 3,406 loss-type CNVRs spanning 328.76 Mb (average length: 96.5 kb), and 523 mixed-type CNVRs covering 511.67 Mb (average length: 978.3 kb). Chromosome 19 exhibited the highest CNVR coverage (65.58%), while chromosome 9 showed the lowest coverage (27.20%).

Analysis of the 20 longest CNVRs revealed that 50% were mixed-type, with 6 located in chromosomal telomeric regions. This distribution pattern suggests that CNVs frequently occur in highly repetitive telomeric regions, which are known hotspots for large-scale genomic recombination and replication events.

To investigate the functional significance of CNVs in goat growth development, Genome-wide association studies (GWAS) were conducted for four early growth traits to investigate the functional significance of CNVs in goat growth development. The results are presented in Table 3. Using the Bonferroni correction method, A genome-wide significance threshold of 3.85 × 10⁻⁵ (corresponding to 1/N, where N = 26,003) was established for identifying significant CNV-trait associations. The GWAS analysis revealed 11 CNVs that surpassed the significance threshold, demonstrating associations with early growth traits. Specifically, two genome-wide significant CNVs associated with BW on chromosomes 18 and 20, three CNVs associated with both WW and ADG on chromosomes 13, 19, and 20, and six CNVs associated with YW distributed across chromosomes 3, 5, 21, 22 and 23 were identified. Consistent with the previously reported high genetic correlation between WW and ADG, our analysis identified two overlapping CNVs (CNV_DEL_11189 and CNV_DEL_17895) that showed significant associations with both traits. These CNV regions, along with all other identified variants, were precisely mapped and validated using bedtools software. Functional annotation of the 11 significant CNVs revealed their association with 91 protein-coding genes, as detailed in Table 4.

The GWAS analysis for BW identified two significant loss CNVs: CNV_DEL_17406 on chromosome 18 and CNV_DEL_18821 on chromosome 20 (Figure 4a). These two CNVs were functionally annotated to four candidate genes, including Zinc finger protein 845 (ZN845), Endogenous retrovirus group K member 25 Env polyprotein (ENK25), Zinc finger protein 160 (ZN160), and Zinc finger protein 208 (ZN208) among others (Table 3).

The GWAS analysis for ADG and WW identified three significant CNVs: two loss-type CNVs (CNV_DEL_11189 on chromosome 13 and CNV_DEL_17895 on chromosome 19) showing associations with both traits, and one gain-type CNV (CNV_DUP_18956 on chromosome 20) specifically associated with ADG (Figure 5). The overlapping associations of CNV_DEL_11189 and CNV_DEL_17895 with both ADG and WW support the previously reported high genetic correlation between these traits. Functional annotation revealed that these three CNVs encompassed 66 candidate genes, including Transcription factor SOX-15 (SOX15) and Thialysine N-epsilon-acetyltransferase (SAT2; Table 3).

The GWAS analysis for YW identified six significant CNVs distributed across chromosomes 3,5,21,22 and 23: CNV_DUP_3067, CNV_DEL_4359, CNV_DEL_4552, CNV_DUP_19581, CNV_DUP_20170, and CNV_DUP_21558 (Figures 6a,b). Functional annotation of these six CNVs revealed their association with 32 candidate genes, including Inositol 1,4,5-trisphosphate-gated calcium channel (ITPR1) and ADP-ribosylation factor-like protein 8B (ARL8B), among others (Table 3).

The qPCRs were performed to verify 8 CNVs in 32 samples. As shown in Figure 4, more than 78.0% of the results are consistent with the type of CNV predicted using CNVnator (Figure 7).

Comparative analysis revealed that individuals carrying either CNV_DEL_17406 or CNV_DEL_18821 exhibited significantly higher average BW (p < 0.05) compared to the overall IMCG population mean of 2.69 kg (Figure 8a). This finding suggests a potential positive association between these CNVs and increased birth weight in the IMCG population.

Comparative analysis demonstrated significant phenotypic effects: individuals carrying CNV_DEL_11189 or CNV_DUP_18956 exhibited higher average ADG (p < 0.05) than the population mean (0.16 kg), while those with CNV_DEL_17895 showed lower ADG (p < 0.05; Figure 8b). Similarly, for WW traits, CNV_DEL_11189 carriers displayed significantly higher WW (p < 0.05) than the population mean (21.90 kg), whereas CNV_DEL_17895 carriers showed reduced WW (p < 0.05; Figure 8c).

Comparative phenotypic analysis demonstrated distinct effects of these CNVs on YW: individuals carrying CNV_DEL_4359, CNV_DEL_4552, or CNV_DUP_20170 exhibited significantly higher YW (p < 0.05) than the population mean (31.21 kg). While CNV_DUP_3067 and CNV_DUP_19581 carriers showed elevated YW compared to the population mean, these differences did not reach statistical significance. In contrast, individuals with CNV_DUP_21558 displayed significantly lower YW (p < 0.05) than the population average (Figure 8d). Notably, CNV_DEL_4359 and CNV_DEL_4552 showed particularly strong positive associations with YW, with carriers exhibiting significantly higher weights than the population mean (Figure 8e).

Comparative analysis using independent samples revealed significant phenotypic divergence (p < 0.05) between CNV carriers and the Inner Mongolia Cashmere Goats (IMCGs) population baseline in 9 out of 11 candidate copy number variations (CNVs) associated with early growth traits. While two CNV loci (CNV_DEL_17406 and CNV_DUP_20170) showed non-significant p-values (0.19 and 0.23, respectively), their effect sizes (Cohen’s d > 0.4) demonstrated meaningful biological divergence from population means. This validation framework achieved 81.8% concordance between GWAS-identified loci and measurable phenotypic effects, with the remaining variants showing directional consistency in trait modulation, thereby providing orthogonal validation for the GWAS findings.

Genomic annotation was performed on genes overlapping with the 11 significant CNVs, including their upstream and downstream 300 kb flanking regions. Based on the Ensembl annotation of the Capra hircus reference genome (Genome assembly ASM4082201v1), a total of 91 genes were identified. To elucidate the functional relevance of these genes to early growth traits in IMCGs, comprehensive enrichment analyses were conducted using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO) terms. The gene set enrichment analysis revealed several functionally relevant terms associated with early growth traits in IMCGs. Specifically, the analysis identified 20 significant GO terms, comprising 6 cellular component terms, 7 biological process terms, and 7 molecular function terms. KEGG demonstrated that these genes are predominantly involved in crucial biological processes, including cellular development, hormonal biosynthesis, protein phosphorylation, and related metabolic pathways (Figure 9).

Through integrative analysis combining data from the GeneCards database and published literature, we identified several genes implicated in key biological pathways and processes. Seven candidate genes were selected from this analysis based on their overlap with significant CNVs and enrichment in gene set analysis (Table 4). The candidate gene analysis revealed ZN845 was as the primary candidate for BW, while SOX15, FGF11, GPS2, and DVL2 emerged as candidate genes for WW and ADG. For YW, SPRY4 and STAT2 were identified as potential candidate genes. These candidate genes associated traits are comprehensively summarized in Table 4.

Functional enrichment analysis revealed distinct biological roles for the identified candidate genes. For BW traits, ZN845 did not exhibit significant enrichment in any specific pathways within our analysis framework. For ADG and WW traits, the SOX15 demonstrated significant associations with critical cellular processes, including cell differentiation, DNA-binding transcription factor activity, and RNA polymerase II-specific regulation. FGF11 was prominently enriched in pathways governing protein phosphorylation modulation and cellular differentiation processes. GPS2 showed specific involvement in transcriptional regulation mediated by RNA polymerase II, particularly in the context of human T-cell leukemia virus 1 infection. DVL2 was significantly enriched in multiple pathways, including the regulation of cell population proliferation, signaling pathways controlling stem cell pluripotency, and positive regulation of protein phosphorylation. For YW traits, STAT2 was associated with fundamental biological processes such as cell proliferation, same-protein binding interactions, and the regulation of osteoclast differentiation. SPRY4 is significantly enriched in cell solutes (Table 4).