Unveiling the common loci for six body measurement traits in Chinese Wenshan cattle

Introduction: Body measurement traits are integral in cattle production, serving as pivotal criteria for breeding selection. Wenshan cattle, a local breed in China’s Yunnan province, exhibit remarkable genetic diversity. However, the molecular mechanisms regulating body measurement traits in Wenshan cattle remain unexplored. Methods: In this study, we performed a genome-wide association method to identify genetic architecture for body height body length hip height back height (BAH), waist height and ischial tuberosity height using the Bovine 50 K single nucleotide polymorphism Array in 1060 Wenshan cattles. Results: This analysis reveals 8 significant SNPs identified through the mixed linear model (MLM), with 6 SNPs are associated with multiple traits and 4 SNPs are associated with all 6 traits. Furthermore, we pinpoint 21 candidate genes located in proximity to or within these significant SNPs. Among them, Scarb1, acetoacetyl-CoA synthetase and HIVEP3 were implicated in bone formation and rarely encountered in livestock body measurement traits, emerge as potential candidate genes regulating body measurement traits in Wenshan cattle. Discussion: This investigation provides valuable insights into the genetic mechanisms underpinning body measurement traits in this unique cattle breed, paving the way for further research in this domain.


Intruduction
The nutritional value of beef, characterized by its high-quality animal protein content, underscores its significance in the meat market (Wu G, 2020).With an escalating demand for beef due to the increased consumption, understanding and enhancing the economic traits associated with beef cattle have become paramount (MacDonald, 2003).Body measurement traits, such as body length (BL), body height (BH) and back height (BAH), are widely used as important indicators for selection criterion to improve the production of the beef cattle (Gritsenko et al., 2023).According to previous study, these body measurement traits have demonstrated correlations with various other economically significant characteristics, such as reproductive traits (Bessa et al., 2021;Shin et al., 2021), longevity (Forabosco et al., 2004;Setati et al., 2004), feed efficiency (Benfica et al., 2020), carcass traits (Munim, et al., 2013) and growth traits (Naserkheil et al., 2020).Furthermore, body measurement traits can also be used to predict cattle's health status (Schmidtmann et al., 2023).In the past few decades, the improvement of beef cattle body measurements has relied on conventional breeding techniques.Body measurement traits, as a kind of quantitative traits which are controlled by a few major genes and numerous minor genes, are also influenced by various environmental variables, such as diet, prenatal conditions and living environment (Larson et al., 2023;Dawson et al., 2022).Thus, the efficacy of traditional methods is constrained (Ring et al., 2019).Fortunately, due to the reduction in sequencing costs, the completion of the cattle genome sequence, and the continuous improvement of the genetic evaluation methods, molecular marker-assisted breeding has become an effective and reliable way to improve beef cattle's body measurement traits (Naserkheil et al., 2020).
Wenshan cattle, colloquially referred to as Wenshan-humped yellow cattle, represent a distinctive beef cattle breed primarily found in the Wenshan Zhuang and Miao Autonomous Prefecture of Yunnan Province in China.Wenshan cattle has excellent features like adaptation to rough feed, stress resistance capacity, early mature (China National Commission of Animal Genetic Resources, 2011), tolerance to hot and humid environments (Liu et al., 2022;Yan et al., 2022), superior meat quality, and a wealth of genetic diversity (Nie et al., 1999;Chen et al., 2020), and has important value of research.Chen, Zhi et al. (2019)showed that the contents of fat, myristic acid, palmitic acid and palmitoleic acid in the longissimus muscle of 1-year old Wenshan bull were higher than those of Simmental cattle.However, akin to other beef cattle breeds, Wenshan cattle have their inherent limitations.As to 18-month old cattle, the average body weight of Wenshan cattle was significantly lower than Simmental cattle (p = 0.0004) (Li et al., 2019).The growth performance and carcass traits of Wenshan cattle were also lower than those of Yunling cattle and Simmental cattle, respectively (Fan et al., 2020).At 12 months, Wenshan cattle display notably smaller body measurements, including withers height, body slanting length, chest circumference, and hip and rump length, compared to Yunling cattle and Simmental cattle (p < 0.05) (Meng et al., 2020).These deficiencies, such as smaller size, long growth cycle and lower growth rate, along with diminished meat production performance and higher feeding costs in comparison to widely used breeds, which restrict its further development.Therefore, improvement in the growth traits is the key to the development of Wenshan cattle industry.
With the rapid development of single nucleotide polymorphism (SNP) arrays and the reduction of the genotyping costs, genome-wide association study (GWAS) has been widely used to identify genomic regions that control economic traits in domestic animals (Zhang et al., 2019;Song et al., 2023).GWAS employs thousands of SNP markers to efficiently associate economic traits in beef cattle, such as body measurement traits, and many candidate genes has been identified in recent years (van den Berg et al., 2022;Niu et al., 2021;Bouwman et al., 2018).In a recent study, 463 Wagyu beef cattles were genotyped by 770K SNP array to identify candidate genes for body measurement traits, and 18, 5 and 1 SNPs associated with hip height, body height and body length were detected respectively through GWAS (An et al., 2019).Vanvanhossou et al. (2020) detected that PIK3R6 and PIK3R1 had direct functional associations with height and body size for 4 Beninese indigenous cattle breeds.Terakado, et al. (2018) identified that the gene TMEM68 had associated with yearling height in Nelore cattle using GWAS based on the Illumina BovineHD BeadChip.Recently, An et al. (2020) conducted GWAS on 1,217 Chinese Simmental beef cattle to study the changes of body size traits such as heart size (HS), abdominal size (AS), body height (BH), body length (BL), hip height (HH) and cannon bone size (CS) at three different growth stages (6, 12 and 18 months).In the above research, Illumina Bovine HD 770 K BeadChip was used through GWAS, and revealed 58 significant SNPs compatible with 21 genes correlated to body size in Simmental beef cattle.Due to the difficulty in collecting phenotype data and the high cost of large population sequencing, studies on the molecular mechanisms of body size traits in beef cattle remain limited.However, GWAS analyses face a significant challenge in controlling false positives and false negatives that may arise due to population structure and kinship.To solve this problem, mixed linear models (MLMs) were frequently used, which incorporate covariates for population structure and principal component analysis (PCA).Linkage disequilibrium (LD) serves as the foundation of GWAS, which provided insight in identifying genomic regions to improve economically important traits (McKay et al., 2007;Porto-Neto et al., 2014).LD is mainly determined by the physical distance between markers, and may also be influenced by several other demographic and evolutionary factors including genetic bottleneck, population stratification, selection, inbreeding, genetic drift, effective population size, migration, mutation, and recombination rate (Reich et al., 2011;Karimi et al., 2015;Singh et al., 2021).Consequently, GWAS is suitable for studying body measurement traits.
In this study, the GWAS was conducted using the Bovine 50K SNP Array to identify significant SNPs associated with six body measurement traits of 1060 Wenshan cattles, including body height (BH), body length (BL), hip height (HH), back height (BAH), waist height (WH) and ischial tuberosity height (ITH).To enhance the accuracy of our analysis and control for population structure, we employed the Mixed Linear Model (MLM) and used population structure as covariates.The primary objective of this study was to map significant SNPs, with a particular focus on common SNPs, and to identify candidate genes that play a role in body measurement-related traits.Besides, Mixed Linear Model (MLM) was used for GWAS analysis, and population structure was used as covariates.The objective of this study was to map significantly associated SNPs, especially common SNPs, and to identify candidate genes involved in body measurement related traits.This study represents the first comprehensive genome-wide investigation of body measurement-related traits in Wenshan cattle, with the aim of identifying candidate genes and potential markers.The findings of this research provide valuable insights and pave the way for further exploration of the genetic mechanisms underlying body measurement traits in Wenshan cattle.This study offers valuable insights for the further investigation of potential genetic mechanism of body measurement traits in Wenshan cattle.

Ethics statement
All of the animals employed in this research were handled following the guidelines established by the Ministry of Agriculture and Rural Affairs of China for use of experimental animals.The ethics committee of Yunnan Agricultural University (YNAU, Kunming, China) approved of the entire research.

Animals, phenotypic collection and statistical analysis
According to the date of birth, 1,060 Wenshan cattles (900 females and 160 males) aged around 12-monthold (±15 days) were collected from Yunnan Guduo Agriculture and Animal Husbandry Co., Ltd. in Yunnan Province, China.These cattles were reared under the homologous management conditions and similar environments.In this study, 6 phenotypic traits including BH, BL, WH, BAH, ITH and HH were measured simultaneously for each individual by the method as in Ref. (An et al., 2020).Ear tissue samples were collected from the up 1,060 Wenshan cattles, and were stored at −80 °C.
SPSS25 software was used to make descriptive statistics, including the number, minimum, maximum, mean, standard deviation, and coefficient of variation for our 6 body measurement traits in Wenshan cattle.Phenotypic correlations among the body measurement traits were calculated within the R statistical environment and used to determine whether they reflected the relationships between the GWAS results.

Genotyping and quality control
Ear samples of 1,060 Wenshan cattle were used to extract the DNA using the Tissues Genomic DNA (Omega Bio-Tek, Norcross, GA, United States) kit according to the manufacturer's instructions.After that, the DNA was quantified and genotyped using the Bovine 50K SNP Array and the PLINK v1.90 Software was used for quality control (Purcell et al., 2007).Briefly, The raw genotypic data were controlled according to call rates >0.95 and markers with call rates >0.95, minor allele frequency (MAF) > 0.01, and Hardy-Weinberg (HWE) p > 10 -6 were retained.All markers located on sex chromosomes or in unmapped regions were excluded.Missing genotypes were imputed using the Beagle software (Browning and Browning, 2009).After that, 46,284 high-quality SNPs were used for subsequent analyses.

Population structure
To investigate population structure of Wenshan cattle, population stratification was detected using eigenvalues and eigenvectors calculated by GCTA software (Li et al., 2016) based on principal component analysis (PCA).The SNPs density distributions plot was drawn using the R package "CMplot" (Yin, L. et al., 2021).

Association analysis
To test for normality, the Shapiro test was performed on our phenotypic data of 6 body measurement traits.Since the original data are not normally distributed, a log transformation was performed to normalize it.The normalized data on the 6 body measurement traits were then subjected to a one-way analysis of variance (ANOVA), followed by post hoc Fisher's least significant difference (LSD) test using the Agricolae package (Steel and Torrie, 1997) in R 4.2 (https://www.R-project.org/).The log-transformed phenotypic data were eventually used in the GWAS analysis.The single variable mixed linear model (MLM) of GEMMA software (Zhou and Stephens, 2012) was employed for association analysis between body measurement traits and filtered SNP markers (n = 46,284).The first three principal components of principal component analysis were used as covariables to conduct association analysis.The univariate linear mixing model is as follows: Where y is the phenotypic value vector; γ is the regression coefficient; PCA is a covariate vector; β is the fixed effect vector of birth year, birth month and sex; α is SNP effect value vector; μ is the residual polygenic effect value vector, which follows the [a~MVN(0, A σ 2 a )] distribution, and A represents the molecular affinity matrix.e is the residual effect value vector, following the [e ~MVN(0, I σ 2 e )] distribution, I represents the unit vector; X, Z and W are the correlation matrices of β, α and μ, respectively.
Using the Bonferroni correction, the genome-wide significant thresholds and suggestion threshold were set as p < 0.05/N, p < 1/N, respectively, where N was the number of SNPs tested in the analyses.Genome-wide significant and suggestive levels were set as 0.05⁄ 46,284 = 1.08 × 10 −6 and 1⁄ 46,284 = 2.16 × 10 −5 , respectively.Finally, the manhattan and Quantile-Quantile (Q-Q) plots of GWAS were drawn using the R package "qqman".

Candidate gene selection and functional annotations
Candidate genes were selected within 0.25 Mb upstream or downstream of the significant SNPs based on LD value for Wenshan cattle population.The major genome browser, Ensembl (https://uswest.ensembl.org/index.html), was used to annotate the significant SNPs identified for the 6 body measurement traits in Wenshan cattle.Candidate genes were then selected according to their biological function.Frontiers in Genetics frontiersin.org3 Results

Phenotype data statistics
The statistical information on the six body measurement traits is shown in Table 1.The mean values for BL, BH, WH, ITH, BAH and HH were 101.58 cm, 104.33 cm, 104.30cm, 92.23 cm, 100.63 cm and 88.44 cm, respectively.Coefficients of variation for the 6 body measurement traits were 9.6, 8.0, 7.7, 7.9, 7.8 and 7.8, respectively.
The phenotypic correlation coefficients for the 6 body measurement traits were shown in Table 2.The results revealed that significant positive correlations were found between all of these 6 traits (p < 0.01), ranging from 0.603 to 0.961.The highest correlation was found between BH and BAH (r = 0.961, p < 0.001).

Population structure
After elimination of the monomorphic loci and loci with minor allele frequency (MAF) below 1%, from the total of 63,791 SNP markers, 46,284 SNPs remained for associations with body measurement of Wenshan cattle.The density distributions of the filtered SNPs are shown in Figure 1.The highest number of markers was on chromosome 1 (3,168 SNPs), and the lowest number of markers was on chromosome 28 (824 SNPs).SNPs were found in almost all of the non-overlapping 1 Mb regions of the genome, indicating the reliability of the data.
In our study, in order to investigate the genetic structure of 1,060 Wenshan cattle, quality-controlled SNP markets were utilized for principal component analysis (PCA).The analysis revealed that there were no distinct genetic clusters among the samples, as shown in Figure 2A.Furthermore, the heatmap and dendrogram of the kinship matrix confirmed the absence of clear clusters in the population, indicating that the Wenshan cattle population in present study is irrelevant to family (Figure 2B).

FIGURE 1
The filtered SNPs density distributions on chromosomes.The horizontal axis (x-axis) shows the chromosome length (Mb).Color index indicates the number of labels.
For BH trait, a total of 6 SNPs were detected, with two located on chromosome 17 and the others on chromosomes 3, 13, 16 and 22, respectively.Among them, the most significant SNP was Affx-257,095,832 (p = 1.53 × 10 −64 ), which located on BTA13: 5,348,767 bp.For BL trait, a total of 6 SNPs were detected, with two located on chromosome 17 and the others on chromosomes 2, 13, 22 and 24, respectively.Among them, the most significant SNP was Affx-41,315,554 (p = 1.58 × 10 −45 ), which located on BTA17: 52,029,736 bp.For HH, WH and BAH trait, 6 identical SNPs were detected for every one of them, with two located on chromosome 17 and the others on chromosomes 2, 3, 13 and 22, respectively.For HH, WH and BAH, Affx-41,315,554 was the most significant one among their selected SNPs, which located on BTA17: 52,029,736 bp, and the p-values were 2.52 × 10 −64 , 6.26E×10 −68 and 2.3 × 10 −68 , respectively.For ITH trait, a total of 4 SNPs were detected, with two located on chromosome 17 and the others on chromosomes 13 and  22, respectively.Among them, the most significant SNP was Affx-41,315,554 (p = 8.96 × 10 −64 ), which located on BTA17: 52,029,736 bp.The GWAS results for the 6 body measurement traits were illustrated by Manhattan plots and Quantile-Quantile (Q-Q) plots in Figure 3. QQ plots of BH, BL, WH, BAH, ITH and HH traits showed no evidence of population stratification in our study.

Identification of candidate genes
To identify candidate genes associated with the 6 body measurement traits of Wenshan cattle, we searched for genes within a 0.5 Mb of the identified SNPs using the UCSC genome browser.In total, 21 genes were identified within or overlapping the candidate regions of 8 SNPs (Table 4).And, 1, 6, 13 and 1 genes were mapped for 095,832,873,673,315,062,550,521,055,744,044,208,517,222 did not map any gene.Among these genes, 20 genes were associated with all the 6 traits, and mainly distributed on BTA13, 17 and 22. HIVEP3, which located on BTA3, was associated with WH, BAH, BH and HH, respectively.Through reviewing relevant research reports, we found that Scarb1, AACS and HIVEP3 as candidate genes for body measurement traits according to their biological functions.

Discussion
Population stratification is a crucial factor to consider in GWAS, as it can lead to an increased rate of false positives (Devlin and Roeder, 1999;Price et al., 2006;Kang et al., 2010;Zhou and Stephens, 2012;An et al., 2019).Many studies have indicated that adding population stratification to GWAS models can significantly improve the accuracy of analysis results (Devlin and Roeder, 1999;Price et al., 2006;Kang et al., 2010;An et al., 2019).MLM with principal component analysis (PCA) model was widely used to simulate population structure, kinship, and family structure, which is currently the most effective method to reduce population stratification (Price et al., 2006;Zhang et al., 2010).In the study of a commercial broiler chicken population, Mebratie, Wossenie et al. (2019) conducted GWAS using MLM approach, and found that using models that account for the population structure may reduce bias and increase accuracy of the estimated SNP effects in the association analysis.From the QQ plots (Figure 3), it can be seen that using the MLM model in our study can effectively reduce the impact of population stratification.Meanwhile, one of the major challenges in GWAS is multiple hypothesis testing (Joo et al., 2016).GWAS involves the assessment of millions of statistical tests, and as such, the p-value threshold for statistical significance must be adjusted to control the overall false positive rate.In order to control false positives or ensure the accuracy of analysis results, the Bonferroni correction approach which may be a little strict, and is often applied in GWAS (Al-Mamun et al., 2015;Puglisi et al., 2016;Abdalla et al., 2021).In the study for body measurement traits in Chinese Holstein cattle based on GWAS, 11 significant SNPs were identified using the Bonferroni correction approach, and identified 6 related genes (Abdalla et al., 2021).Thus, the Bonferroni was used as multiple hypothesis testing in this study and mapped 8 significant SNPs.
Body measurement traits in cattle are not only indicators of their physical development and size but are also important for assessing breed quality and production performance, such as weight.We observed a strong positive correlation among various body measurement traits, which is consistent with findings from other related studies.GUO et al. (2023) recently found Jinnan cattle exhibit a high positive correlation in body measurement traits (0.47-0.96), with a correlation of about 0.93 between WH and BH, which was aligns with our study (0.92).A study conducted by Afolayan et al. (2007) found that correlations between body measurements including HH and BL, at 400 days of age in Australian, and reported the values for phenotypic correlations were also highly positive (0.56-0.78).Another study (Kamprasert et al., 2019) found that the correlation between HH and BL of Brahman cattle at 400 days (0.80) was higher than our correlation of 0.629, but fell to 0.64 at 600 days, which indicated that the correlation between body measurement traits was influenced by both genetics and environment.Considering the strong correlation between multi body measurement traits, their genetic regulatory mechanisms should also have some similarities.In this study, we identified 8 significant SNPs, of which 6 appeared in more than one of the 6 body measurement traits, which was uniform with their phenotype correlations.
More importantly, sample size is the most critical factor that affects statistical power and limited sample size is hence a major hurdle in GWAS for traits that are difficult or expensive to measure (Gebreyesus et al., 2019).In this study, 1,060 Wenshan cattle were used for the GWAS of body measurement traits, with a much larger sample size than other similar studies in beef cattle (Buzanskas et al., 2014;An J. et al., 2019;Zepeda-Batista et al., 2021), thus ensuring the accuracy of our results.In this study, we identified three important genes including Scarb1, AACS and HIVEP3, which have been indicated to be related to bone metabolism, but there have been no previous studies on their association with livestock body measurement traits.These results indicate that those three genes may be specific to body measurement traits in Chinese Wenshan cattle and have high research value.
Scavenger receptor type I (Scarb1) on BTA17 was related to all the 6 body measurement traits.Scarb1 produced by the SCARB1 gene, is the major receptor for high-density lipoprotein (HDL).The role of Scarb1 in bone formation has been the subject of conflicting findings in recent research.Martineau et al.(2014a;2014b;2014c) reported that Scarb1 knockout (KO) mice have increased trabecular bone and with no changes in osteoclasts parameters.In contrast, Tourkova et al. (2019) showed that Scarb1 KO mice was osteopenic relative to the wild type.Morever, osteoblast and osteoclast-related mRNAs of KO mice greatly decreased compared to WT mice, suggesting that Scarb1 is required for normal bone differentiation.More strangely, the latest research from Palmieri M et al. (2023a;2023b) reported that bone mass was not affected in Scarb1 KO mice and thought Scarb1 did not contribute to bone homeostasis.However, it is established that HDL, via its main protein component apolipoprotein A1, is essential for normal bone formation by affecting osteoblast (Blair et al., 2016;Papachristou et al., 2017).Then, as the major receptor for HDL, we speculate that Scarb1 may play an important role in bone formation, and relevant molecular mechanism needs further research.
HIVEP3, also known as Schnurri-3 (SHN3), is a big zinc finger protein and an essential regulator of bone formation.Jones C et al. (2006) found that HIVEP3 −/− mice reveled osteosclerosis due to increased osteoblast activity and increased bone mass.Further research found that, HIVEP3 regulated the bone mass by controlling the protein levels of Runx2 which was the principal transcriptional regulator of osteoblast differentiation (Otto et al., 1997;Jones et al., 2006;Jones et al., 2007).Further investigation into the role of HIVEP3 in osteoblasts has revealed that HIVEP3 functions as a dampener of ERK (extracellular signal-regulated kinase) activity, particularly downstream of WNT signaling in osteoblasts.Mutations in the HIVEP3 gene can disrupt this interaction, leading to abnormal activation of ERK and hyperactivity of osteoblasts in vivo (Shim et al., 2013).A genome-wide association analysis in Chinese populations suggested that the HIVEP3 gene might be potentially correlated with Femoral Neck Bone Mineral Content (BMC) and Hip Geometry (Hu et al., 2020).By searching related researches, this study was the first to find that the HIVEP3 gene had a potential relationship with BH, HH, WH, and BAH in the livestock population, especially in beef cattle.Thus, HIVEP3 can be regarded as a candidate gene for Wenshan cattle body measurement traits.
The acetoacetyl-CoA synthetase (AACS) gene is known for the synthesis of biologically important lipogenic substances (Endemann et al., 1982).Previous reports have suggested that high level of serum lipids can trigger bone metabolic disorders (Bredella et al., 2013;Alsahli et al., 2016).In 1999, Mundy et al. (1949) found that statins are powerful cholesterol lowering medications that can improve new bone formation in rodents.Another study (Hasegawa et al., 2012) demonstrated that knockdown of AACS in vivo decreases peroxisome proliferator-activated receptor γ (PPARγ) and CCAAT/enhancer binding protein α (C/EBP α), which can enhance bone resorption and have a critical role in relationship between obesity and bone loss (Motyl et al., 2011;Akune et al.,  Frontiers in Genetics frontiersin.org08 2004).In obesity, AACS was highly expressed in the differentiated osteoclasts, but did not in osteoblast differentiation (Yamasaki et al., 2016).In situ hybridization, AACS was observed in several regions of the embryo, including the backbone region (especially the somite) and adult femur epiphysis.Collectively, these findings indicated that AACS may be involved in bone homeostasis through its impact on adipogenic transcription factors such as the C/EBP family and PPARγ.Then, we speculate that AACS plays an important role associated with the bone homeostasis, which is worthy of further study and can be selected as a candidate gene for body measurement traits.

Conclusion
To sum up, 8 SNPs were detected to be associated with 6 body measurement traits.Interestingly, 6 of these SNPs were detected in more than one of our 6 body measurement traits, which was consistent with their phenotype correlations.Meanwhile, we found 21 candidate genes located nearby or within the associated SNPs.Among them, Scarb1, AACS and HIVEP3 have been indicated to be related to bone formation and are also discovered rarely in livestock body measurement traits, which can be used as candidate genes for Wenshan cattle.This study offers valuable insights for the further investigation of potential genetic mechanism of body measurement traits in Wenshan cattle.

FIGURE 2
FIGURE 2 Population structure analysis.(A) The principal component analysis (PCA) for the third principal component (PC3) against the first principal component (PC1) and the second principal component (PC2).(B) The hierarchical clustering and heat map of the pairwise kinship matrix values.

FIGURE 3
FIGURE 3 Manhattan plots and QQ plots for six body measurement traits using MLM.(A) BL, (B) BH, (C) ITH, (D) WH, (E) BAH, (F) HH.Negative log10 p-values of the filtered high-quality SNPs were plotted against their genomic positions.The solid lines of orange and blue correspond to the Bonferroni-corrected thresholds of p = 2.16 × 10 −5 and p = 1.08 × 10 −6 , respectively.

FIGURE 4
FIGURE 4Venn diagram representation of the common significant loci in 6 body measurement traits for Chinese Wenshan cattle.

TABLE 1
Descriptive statistics of 6 body measurement traits in Chinese Wenshan cattle.

TABLE 2
Correlation coefficients of 6 body measurement traits in Chinese Wenshan cattle.

TABLE 3 Significant
SNPs identified for 6 body measurement traits in Chinese Wenshan cattle.

TABLE 4
List of candidate genes associated with 6 body measurement traits in Chinese Wenshan cattle.