1- Department of Animal Science and Fisheries, Sari Agricultural Sciences and Natural Resources University (SANRU), Sari, Iran.
2- Center for immunity and immunotherapy, Seattle children's Research Institute, Seattle, USA.
2- Center for immunity and immunotherapy, Seattle children's Research Institute, Seattle, USA.
Abstract: (54 Views)
Introduction and Objective: Growth traits are among the most important economic indicators in livestock production, directly linked to meat yield, resource efficiency, and profitability. In yak, which is well adapted to high-altitude regions, growth performance is often constrained by harsh environmental conditions. Understanding the genetic basis of growth traits can provide insights that contribute to molecular breeding and genetic improvement. Genome-wide association studies (GWAS) are valuable tools for identifying genomic regions influencing economically important traits in livestock and for advancing genomic breeding programs. However, GWAS has several limitations: markers are tested individually, and stringent statistical thresholds often result in many causal variants being undetected, particularly for polygenic traits. On high-density SNP chips, the effect of a gene may be distributed across multiple markers, preventing any from reaching genome-wide significance. Moreover, extensive linkage disequilibrium in livestock populations complicates the precise identification of causal mutations. Additionally, GWAS does not capture the real interactions among genes within biological pathways or functional networks, and due to the biallelic nature of SNPs, it cannot fully represent the effects of multi-allelic QTLs. To overcome these challenges, it has been suggested that analyses should extend beyond SNP-level to gene-level or gene-set approaches. In this framework, genes harboring significant SNPs that participate in common pathways or biological processes are tested for statistical enrichment. Such approaches enhance the power to detect causal variants and provide deeper insights into the genetic architecture of complex traits. Based on this rationale, the present study utilized publicly available genomic data to functionally analyze candidate genes associated with growth traits in yak and to biologically interpret GWAS results in order to provide more practical implications for breeding programs. The main objective was to identify genomic loci associated with post-weaning growth in Ashidan yaks through GWAS and to investigate functional pathways by means of gene enrichment analysis.
Materials and Methods:
This study used previously published data from 354 female Ashidan yaks (Jia et al., 2020). Genotyping had been performed using the Illumina BovineHD BeadChip, which included 777,962 SNP markers. Data quality control was conducted in PLINK. SNPs with more than 0.05 missing data, minor allele frequency (MAF) < 0.05, or significant deviation from Hardy–Weinberg equilibrium (p < 10-6) were excluded, as were individuals with >0.05 missing data. Missing genotypes were imputed using the LD-kNNi algorithm in TASSEL. Population structure was evaluated by principal component analysis (PCA), and the first five PCs were included as covariates in the mixed linear model. GWAS was carried out in TASSEL v5.2 using a mixed linear model (MLM) incorporating both a kinship matrix (K) and population structure (Q). To identify biological pathways associated with growth traits, Gene Ontology (GO) analysis was performed. All genes located within ±15 kb of SNPs with p-value < 0.05 were extracted based on the Bos taurus UMD3.1 genome assembly. These candidate genes were then submitted to the g:Profiler tool for enrichment analysis in three GO domains: Biological Process (BP), Molecular Function (MF), and Cellular Component (CC), as well as in KEGG pathways. Significance thresholds were set at corrected p-value < 0.05.
Results:
Following quality control, one individual was excluded due to poor genotyping quality, leaving 353 yaks for further analysis. From the initial SNP dataset, 2,540 SNPs were removed due to high missingness, 65,122 SNPs were removed due to low MAF, and four SNPs deviating from Hardy–Weinberg equilibrium were excluded. Ultimately, 32,573 high-quality SNPs remained. The average genotyping rate of the retained individuals was 99.6%. For the first trait (average daily gain from 6 to 12 months), no genomic regions reached genome-wide significance (p < 1.5 × 10⁻⁶). In contrast, for the second trait (average daily gain from 6 to 30 months), 22 significant SNPs (p < 1.5 × 10⁻⁶) were detected on chromosome 1, with the strongest signals appearing between 1.8 and 2.1 Mb. Within this region, the gene C1H21orf62 emerged as the main candidate, harboring several significant SNPs. Downstream of this region, additional candidate genes such as GCFC and SYNJ1 were also identified. For the third trait (average daily gain from 12 to 30 months), a highly significant region was again detected on chromosome 1, overlapping with the same region as the second trait, with all significant SNPs closely associated with the C1H21orf62 gene. Functional enrichment analysis further revealed that the identified growth-associated genes were significantly enriched in multiple categories. At the molecular function level, significant enrichment was observed for ion binding (GO:0043167), small molecule binding (GO:0036094), transmembrane transporter activity (GO:0022857), cytoskeletal motor activity (GO:0003774), and hydrolase activity (GO:0016787). At the biological process level, candidate genes were enriched in cell adhesion (GO:0007155), monoatomic ion transport (GO:0006811), and nervous system development (GO:0007399). At the cellular component level, genes were predominantly enriched in cell junctions (GO:0030054), cell periphery (GO:0071944), and plasma membrane (GO:0005886). KEGG pathway analysis revealed that the identified genes were mainly involved in motor protein-related pathways (KEGG:04814), emphasizing the regulatory complexity of growth-related genes.
Conclusion:
This study successfully identified significant SNPs and genomic regions associated with post-weaning growth traits in Ashidan yaks. While no significant associations were detected for early growth (6–12 months), extended growth intervals (6–30 and 12–30 months) revealed several key loci and biologically relevant pathways. The integration of high-density genotyping, genotype imputation, mixed linear modeling, and gene enrichment analysis provided a comprehensive overview of the genetic architecture underlying growth in yaks. These findings lay a valuable foundation for future molecular studies and provide practical markers for use in yak breeding and selection programs.
Materials and Methods:
This study used previously published data from 354 female Ashidan yaks (Jia et al., 2020). Genotyping had been performed using the Illumina BovineHD BeadChip, which included 777,962 SNP markers. Data quality control was conducted in PLINK. SNPs with more than 0.05 missing data, minor allele frequency (MAF) < 0.05, or significant deviation from Hardy–Weinberg equilibrium (p < 10-6) were excluded, as were individuals with >0.05 missing data. Missing genotypes were imputed using the LD-kNNi algorithm in TASSEL. Population structure was evaluated by principal component analysis (PCA), and the first five PCs were included as covariates in the mixed linear model. GWAS was carried out in TASSEL v5.2 using a mixed linear model (MLM) incorporating both a kinship matrix (K) and population structure (Q). To identify biological pathways associated with growth traits, Gene Ontology (GO) analysis was performed. All genes located within ±15 kb of SNPs with p-value < 0.05 were extracted based on the Bos taurus UMD3.1 genome assembly. These candidate genes were then submitted to the g:Profiler tool for enrichment analysis in three GO domains: Biological Process (BP), Molecular Function (MF), and Cellular Component (CC), as well as in KEGG pathways. Significance thresholds were set at corrected p-value < 0.05.
Results:
Following quality control, one individual was excluded due to poor genotyping quality, leaving 353 yaks for further analysis. From the initial SNP dataset, 2,540 SNPs were removed due to high missingness, 65,122 SNPs were removed due to low MAF, and four SNPs deviating from Hardy–Weinberg equilibrium were excluded. Ultimately, 32,573 high-quality SNPs remained. The average genotyping rate of the retained individuals was 99.6%. For the first trait (average daily gain from 6 to 12 months), no genomic regions reached genome-wide significance (p < 1.5 × 10⁻⁶). In contrast, for the second trait (average daily gain from 6 to 30 months), 22 significant SNPs (p < 1.5 × 10⁻⁶) were detected on chromosome 1, with the strongest signals appearing between 1.8 and 2.1 Mb. Within this region, the gene C1H21orf62 emerged as the main candidate, harboring several significant SNPs. Downstream of this region, additional candidate genes such as GCFC and SYNJ1 were also identified. For the third trait (average daily gain from 12 to 30 months), a highly significant region was again detected on chromosome 1, overlapping with the same region as the second trait, with all significant SNPs closely associated with the C1H21orf62 gene. Functional enrichment analysis further revealed that the identified growth-associated genes were significantly enriched in multiple categories. At the molecular function level, significant enrichment was observed for ion binding (GO:0043167), small molecule binding (GO:0036094), transmembrane transporter activity (GO:0022857), cytoskeletal motor activity (GO:0003774), and hydrolase activity (GO:0016787). At the biological process level, candidate genes were enriched in cell adhesion (GO:0007155), monoatomic ion transport (GO:0006811), and nervous system development (GO:0007399). At the cellular component level, genes were predominantly enriched in cell junctions (GO:0030054), cell periphery (GO:0071944), and plasma membrane (GO:0005886). KEGG pathway analysis revealed that the identified genes were mainly involved in motor protein-related pathways (KEGG:04814), emphasizing the regulatory complexity of growth-related genes.
Conclusion:
This study successfully identified significant SNPs and genomic regions associated with post-weaning growth traits in Ashidan yaks. While no significant associations were detected for early growth (6–12 months), extended growth intervals (6–30 and 12–30 months) revealed several key loci and biologically relevant pathways. The integration of high-density genotyping, genotype imputation, mixed linear modeling, and gene enrichment analysis provided a comprehensive overview of the genetic architecture underlying growth in yaks. These findings lay a valuable foundation for future molecular studies and provide practical markers for use in yak breeding and selection programs.
Type of Study: Research |
Subject:
ژنتیک و اصلاح نژاد دام
Received: 2025/09/2 | Accepted: 2025/12/21
Received: 2025/09/2 | Accepted: 2025/12/21
Send email to the article author
| Rights and permissions | |
 |
This work is licensed under a Creative Commons Attribution-NonCommercial 4.0 International License. |
