Genome-Wide Association Study Identifies Candidate Genes Associated with Feet and Leg Conformation Traits in Chinese Holstein Cattle

Simple Summary Feet and leg problems are among the major reasons for dairy cows leaving the herd, as well as having direct association with production and reproduction efficiency, health (e.g., claw disorders and lameness) and welfare. Hence, understanding the genetic architecture underlying feet and conformation traits in dairy cattle offers new opportunities toward the genetic improvement and long-term selection. Through a genome-wide association study on Chinese Holstein cattle, we identified several candidate genes associated with feet and leg conformation traits. These results could provide useful information about the molecular breeding basis of feet and leg traits, thus improving the longevity and productivity of dairy cattle. Abstract Feet and leg conformation traits are considered one of the most important economical traits in dairy cattle and have a great impact on the profitability of milk production. Therefore, identifying the single nucleotide polymorphisms (SNPs), genes and pathways analysis associated with these traits might contribute to the genomic selection and long-term plan selection for dairy cattle. We conducted genome-wide association studies (GWASs) using the fixed and random model circulating probability unification (FarmCPU) method to identify SNPs associated with bone quality, heel depth, rear leg side view and rear leg rear view of Chinese Holstein cows. Phenotypic measurements were collected from 1000 individuals of Chinese Holstein cattle and the GeneSeek Genomic Profiler Bovine 100 K SNP chip was utilized for individual genotyping. After quality control, 984 individual cows and 84,906 SNPs remained for GWAS work; as a result, we identified 20 significant SNPs after Bonferroni correction. Several candidate genes were identified within distances of 200 kb upstream or downstream to the significant SNPs, including ADIPOR2, INPP4A, DNMT3A, ALDH1A2, PCDH7, XKR4 and CADPS. Further bioinformatics analyses showed 34 gene ontology terms and two signaling pathways were significantly enriched (p ≤ 0.05). Many terms and pathways are related to biological quality, metabolism and development processes; these identified SNPs and genes could provide useful information about the genetic architecture of feet and leg traits, thus improving the longevity and productivity of Chinese Holstein dairy cattle.


Introduction
Body conformation traits are considered economically important traits in dairy cattle [1]. Improving the accuracy of selection for body conformation traits would enhance

Ethics Statement
All procedures for collecting hair follicle samples and measuring phenotypic traits were carried out strictly in accordance with the guidelines proposed by the China Council on Animal Care and Ministry of Agriculture of the People's Republic of China. The study was also approved by the Institutional Animal Care and Use Committee of School of the Yangzhou University Animal Experiments Ethics Committee (License Number: SYXK (Su) IACUC 2012-0029), Yangzhou University.

Body Measurements and DNA Samples Collection
The experimental population consisted of 1000 Chinese Holstein cows raised on four farms, (199 individuals from Sihong Farm, 214 from Xuyi Farm, 224 from Xuzhou Farm and 363 from Huaxia Farm). All these farms are located in the northern part of Jiangsu province, China. Phenotypic descriptive traits (BQ, bone quality; HD, heel depth; RLSV, rear leg side view; RLRV, rear leg rear view) were measured individually and scored on a 1-9 scale score according to the National Standards of People's Republic of China Code of practice of type classification in Chinese Holstein (GBT35568-2017). The measurement of traits for each cow were performed by three trained technicians and the average of the measurements taken by the different professionals was used as the phenotype for each trait to ensure the accuracy of the data. For the genotyping analysis, hair follicle samples were collected individually and stored in special paper envelopes for each cow to prevent it from environmental and DNA contamination; each paper envelope must contain enough hair samples, not less than 50 hairs.

Phenotypic and Genetic Parameters
Descriptive statistics and pairwise Pearson correlation coefficients of phenotypic traits were determined using the computer-based software IBM-SPSS, version 25.
Genetic analysis was carried out using the Derivative-free approach to MUltivariate analysis (DMU) software [26] to estimate heritability and genetic correlation between pairs of traits with the animal model as follows: y ijklm = u + Herd i + Year j + Season k + Parity 1 + a m + e ijklm where y ijklm is the phenotype in the jth year, kth season and lth parity of the mth individual from the ith herd; u is overall mean of the population, Herd i is the herd effect according to a cow's origin from one of the four herds; Year j is the jth year effect, Season k is the kth season effect and parity is the effect of lth parity; a is the additive effect of the mth individual, which was evaluated by the pedigree information, and e is the residual in the jth year, kth season and lth parity of the mth individual from the ith herd. All effects were treated as random except the overall mean. The pedigree of the cows could be traced back at least three generations (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020), the parities of cows were between 1 and 4 and four seasons of were defined, December-February, March-May, June-August and September-November.

Genotyping and Quality Control
Each individual from the experimental population was genotyped using the GeneSeek Genomic Profiler Bovine 100 k SNP chip (Neogen Corporation, http://www.neogenchina. com.cn/ (accessed on 28 June 2020)) based on ARC-UCD1.2/bosTau9 as the genome reference. Genomic DNA was extracted from hair follicle samples. The GeneSeek Genomic Profiler Bovine 100 K SNP chip containing 100,000 SNPs was utilized for individual genotyping.
Marker intervals and linkage disequilibrium (LD) were calculated to estimate R square for all markers and plotted the marker distribution as showed in Figure 1d. LD decay fell off quickly within 100 kb physical distance, then decreased slowly afterward.

Population Structure Analysis
We used the Plink 1.90 software [27] to implement principal component analysis (PCA) on 1000 cows genotyping with 84,406 SNPs in Chinese Holstein cattle herds raised at four breeding herds to investigate the population structure and PCA was plotted using the ggplot package In R 4.0.2. The admixture software [28] was also used to study the population structure and correct population stratification.

Genome-Wide Association Studies
In this study, the fixed and random model circulating probability unification (Farm-CPU) method was used to carry out the genome-wide association analysis [29]. The FarmCPU method uses a fixed effect model and a random effect model iteratively. In GWASs, population stratification is the main cause of false positive correlations. Therefore, the fixed effect model tests SNPs one at a time. The significant SNPs are evaluated in the random effect model and the validated SNPs are fitted as covariates in the fixed effect model to control population structure. The model can be written as follows: (1) where y i is the observation of the ith individual; M i1 , M i2 . . . M it are the genotypes of t pseudo QTNs, initiated as an empty set; b 1 , b 2 , . . . , b t are the corresponding effects of the pseudo QTNs; S ij is the genotype of the ith individual and jth genetic marker; d j is the corresponding effect of the jth genetic marker; e i is the residual having a distribution with zero mean and variance of σ 2 e .

Population Structure Analysis
We used the Plink 1.90 software [27] to implement principal component analysis (PCA) on 1000 cows genotyping with 84,406 SNPs in Chinese Holstein cattle herds raised at four breeding herds to investigate the population structure and PCA was plotted using the ggplot package In R 4.0.2. The admixture software [28] was also used to study the population structure and correct population stratification.

Genome-Wide Association Studies
In this study, the fixed and random model circulating probability unification (Farm-CPU) method was used to carry out the genome-wide association analysis [29]. The Farm-CPU method uses a fixed effect model and a random effect model iteratively. In GWASs, population stratification is the main cause of false positive correlations. Therefore, the fixed effect model tests SNPs one at a time. The significant SNPs are evaluated in the random effect model and the validated SNPs are fitted as covariates in the fixed effect model After substitution, every marker has its own p value. The p values and the associated marker map are used to update the selection of pseudo QTNs using the SUPER algorithm (Settlement of MLM Under Progressively Exclusive Relationship) [30] in a REM as follows: where y i and e i are the same as in Equation (1) and u i is the total genetic effect of the ith individual. The total type 1 error (false positive) rate was controlled at 5% and the significance threshold of the GWAS was determined according to this formula (0.05/Nsnp), where Nsnp is the number of SNPs remaining after quality control [31]. Subsequently, the significant threshold for the GWAS was 5.9 × 10 −7 (0.05/84,406) after Bonferroni correction.
Quantile-quantile (Q-Q) and Manhattan plots were drawn using the CMplot package in the R 3.1.1 software [32].

Function and Pathway Enrichment and Network Analysis
In our study, we submitted the candidate genes obtained by GWAS into the Database for Annotation, Visualization and Integrated Discovery (DAVID) [33] for the Gene Ontology (GO) terms [34] and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis [35]. The statistically significant p value for functional analysis and pathway analysis was defined at p ≤ 0.05. Protein-protein interactions among genes were performed using the online Search Tool for the Retrieval of Interacting Genes (STRING) database v11.0 [36] with the Cytoscape software.v3.8.2 to visualize the resultant PPI network.

Descriptive Statistics and Heritability Estimation of Feet and Legs Traits
The feet and leg traits for the 1000 Chinese Holstein cows from four farms measured in this study included heel depth (HD), bone quality (BQ), rear leg side view (RLSV) and the rear leg rear view (RLRV). The descriptive statistics (mean, maximum and minimum values and standard deviation) of phenotypic measurements for these traits are showed in Table 1, where the mean of the HD trait had the greater score (7.04) and the RLSV trait revealed the lowest score, with 3.92. The pairwise genetic and Pearson phenotypic correlation between the traits are provided in Table 2. BQ phenotypically positively correlated with RLSV, and the other traits were low or moderate negatively correlated with RLSV. However, BQ was genetically positively correlated with HD, RLSV and RLRV. HD was negatively correlated with all traits, both genetically and phenotypically, except BQ had a highly positive genetic correlation. RLSV was highly positively genetically correlated with BQ, while it was negative with HD and RLRV, whereas it had a negative phenotypic correlation with BQ, HD and RLRV. Estimates of the heritability results were 0.15, 0.05, 0.17 and 0.15 for heel depth, bone quality, rear leg side view and rear leg rear view, respectively ( Table 2).

Information of SNPs
After quality control, 984 individual cows and 84,406 SNPs were used to conduct the GWAS. The filtered SNPs were distributed on all 29 chromosomes. Chromosome 1 has shown the greatest number of SNPs, whereas chromosome 25 contained the fewest. The minor allele frequency (MAF) for all SNPs was re-calculated after quality control; only MAF above 5% remained. The LD decay line dramatically decreased at the beginning then tended to be slow after 100 kb distance. (Figure 1).

Population Stratification
Principal component analysis (PCA) and the ADMIXTURE program were used to visualize family structure in this study. The results reveal all individuals were grouped into two unequal sized clusters, as shown in Figure 2a,b; as for the admixture population, using ADMIXURE with k value ranging from 1 to 7 (Figure 2c), based on the cross-validation errors, K = 4 was identified to be the optimal number of genetic clusters defining the population structure among the four herds of Chinese Holstein cattle. Also, the results show that the first two principal components were about 21% of the variation (Figure 2a), and they were fitted as covariate variables in the association analysis using the FarmCPU model. Population stratification based on the PCA analysis results was considered and incorporated into the mixed linear model.

GWAS Results
The FarmCPU model was used to conduct the genome-wide association analysis in the present study. In Figure 3, The quantile-quantile (Q-Q) plots illustrated the model used in this study for GWAS analysis was reasonable. The lambda (inflation factor (λ)) was close to 1 (0.9 < 1.01) and the point at the upper right corner of (Q-Q) plots are the significant markers associated with the traits under study ( Figure 3); therefore, the population stratification was adequately controlled. Meanwhile, Manhattan plots are used to visualize GWAS significance level (−log10 of p value of each SNP) by chromosome loca-

GWAS Results
The FarmCPU model was used to conduct the genome-wide association analysis in the present study. In Figure 3, The quantile-quantile (Q-Q) plots illustrated the model used in this study for GWAS analysis was reasonable. The lambda (inflation factor (λ)) was close to 1 (0.9 < 1.01) and the point at the upper right corner of (Q-Q) plots are the significant markers associated with the traits under study ( Figure 3); therefore, the population stratification was adequately controlled. Meanwhile, Manhattan plots are used to visualize GWAS significance level (−log10 of p value of each SNP) by chromosome location (Figure 4).
In the present study, 20 SNPs (Table 3) were past the threshold and significantly associated with four traits of feet and legs (heel depth, bone quality, rear leg side view and rear leg rear view).

Identified Candidate Genes
The results of the linkage disequilibrium (LD) analysis indicate the LD decay (R 2 ), as showed in Figure 1d, tends to be stable when the distance is 200 kb. Therefore, genes located within this region (200 kb) of the significant SNP are defined as candidate genes.

Identified Candidate Genes
The results of the linkage disequilibrium (LD) analysis indicate the LD decay (R 2 ), as showed in Figure 1d, tends to be stable when the distance is 200 kb. Therefore, genes located within this region (200 kb) of the significant SNP are defined as candidate genes.

Functional Analysis
Using the UCSC Genome browser and NCBI database through the Asian server for cow assembly April 2018 (ARC-UCD1.2/bosTau9), 105 genes were obtained within the region of 200 kb up/downstream of the significant SNPs for the feet and legs traits; these genes were used for further enrichment and pathway analysis (gene ontology and KEGG). Gene ontology enrichment analysis revealed 34 GO terms were significantly enriched (p < 0.05), which consist of 18 biological process terms, 1 cellular component and 15 molecular function terms (Table S1 and Figure 5).
The STRING database was used to conduct protein-protein interaction network (PPIN) analysis using all the genes previously used in the functional analysis, to identify the interactions between these genes. Figure 6 shows that there is several interacting between genes (containing 78 nodes connected via 139 edges); the proportional interaction strength between these genes was shown by the intensity of staining between the lines that linked one gene to another ( Figure 6).

Discussion
The esitmated genetic correlations in our study ranged between −0.35 and 0.84 ( Table  2); these results are in agreement with the study by Olasege et al. [37], who reported that between moderate negative and strong postive genetic correlation existed among feet and leg traits, ranging from −0.35 (between rear leg view and set of rear legs) to 0.74 (between foot angle and bone quality) in Chinese Holstein cattle. Among the feet and leg traits, the

Discussion
The esitmated genetic correlations in our study ranged between −0.35 and 0.84 (Table 2); these results are in agreement with the study by Olasege et al. [37], who reported that between moderate negative and strong postive genetic correlation existed among feet and leg traits, ranging from −0.35 (between rear leg view and set of rear legs) to 0.74 (between foot angle and bone quality) in Chinese Holstein cattle. Among the feet and leg traits, the strongest genetic correlation was found between RLSV and BQ (0.84) ( Table 2); if two traits have a high and postive genetic correlation, it implies that most QTL affect them both in the same direction [38]. For intance, most SNPs with significant influences affect weight and height in the same direction, thus helping to explain the known high genetic correlation between these two traits [39].
Population stratification or cryptic relatedness are considered major challenges that face gnome-wide association analysis. Thus, the presence of population stratification can cause spurious association due to systematic ancestry differences [40] that can lead to false positives in GWAS results [41]; therefore, stratification must be well controlled across the experimental population to maximize power to detect true associations [42,43].
There are several methods to correct population stratification; carefully choosing the statistical model can be a useful method to correct and minimize the chances of type 1 error (false positive associations) [44].
To eliminate false positives, the most effective strategy is either (1) fitting population structure as covariates in a General Linear Model (GLM) [45] or (2) fitting both population structure and each individual's total genetic effect as covariates in a Mixed Linear Model (MLM) [46] to make adjustments for testing markers. As MLM can miss some potentially important discoveries and lead to false negatives, due to the confounding among population structures, quantitative trait nucleotides (QTNs) and kinship. Therefore, in the present study, we conducted Fixed and random model Circulating Probability Unification (FarmCPU) model analysis because of its advantages to completely control false positives, eliminate confounding and improve computational efficiency by using the fixed effect model and random effect model iteratively [29].
The inflation factor (λ) should be close to 1, after correcting the population stratification [45]. In the present study, the Q-Q plot in Figure 3 shows that the deviation of the observed value from the expected value is near to 1 and the inflation factor (λ) is 0.9 < 1.01, both indicating that the population stratification was properly corrected by using an appropriate model.
The principal component analysis (PCA) method [47] and ADMIXTURE program analysis [28] are able to detect the population structure by classifying individuals into grouped ancestry based on their genetic makeup. Figure 2a shows all the experimental population was clustered into two groups (one large and one small). This indicates that each group clustered closely together has a genetic relationship. The division in two groups might be attributed to using Holstein semen from different countries or the cows in these farms may contain blood from other breeds, as it is known one of the requirements for registering Holstein cattle in China is that the cattle have at least 87.5 percent blood of Holstein (Chinese Holstein, GB/T 3157 2008).
In the present study, genome-wide association analysis identified 20 significant SNPs associated with four traits belonging to feet and legs trait in Chinese Holstein cattle using the FarmCPU model; among them, the most three significant SNPs were BovineHD2000006450-rs137022628, located near the ACTBL2 gene, ARS-BFGL-NGS-91167-rs41565304, located within the ADIPOR2 gene and no genes were found within 200 kb upstream or downstream of the BovineHD0300022142-rs43350216 SNP.
Until now, only one report has discussed the functional role of the ACTBL2 gene. Hödebeck et al. [48] reported that the silencing of ACTBL2 leads to diminished motility of human arterial smooth muscle cells. These authors also demonstrated that the expression of ACTBL2 in smooth muscle cells under stretch conditions depends on the nuclear factor 5 of activated T-cells (NFAT5).
Adiponectin receptor 2 (ADIPOR2) is a member of the adipocytokines gene family. Generally, adiponectin receptor genes play an important role in bone and whole-body energy homeostasis that could be induced by activation of the AMPK signaling pathway; ADIPOR1 and ADIPOR2 also play a crucial role in metabolic pathways, which lead to the regulation of lipid and glucose metabolism, oxidative stress and inflammation [49]. Ouyang et al. [50] reported ADIPOR2 has a role in intramuscular fat content in the longissimus dorsi (LD) muscle of different pig breeds. Lewis et al. [49] discussed in detail the role of adiponectin signaling in bone homeostasis and its effect on the lifestyle of human healthy aging and disease. Moreover, ADIPOR2 is considered to be one of the genes associated with conformation and reproductive performance traits (hip height, rump length, calving ease, height, ovulation, type and rump angle) in one of the South African cattle breeds [51].
In addition to the ADIPOR2 gene, our study identified other genes that were significantly associated with rear-leg side view traits (INPP4A, DNMT3A, ALDH1A2 and PCDH7).
Some studies concerning the inositol polyphosphate-4-phosphatase type I A gene (INPP4A) reported it to have an association with the thoracic vertebrae number in sheep [52] and a role in regulating the normal functioning of the spinal cord neurons [53]. Additionally, the metabolic process [54] also presented different isoforms in transcriptomic investigation of meat tenderness in two Italian cattle breeds [55].
The DNA methyltransferase 3 alpha gene (DNMT3A) was associated to obese adipose tissue in transgenic mice [56] and significantly associated with beef quality traits [57]. Meanwhile, the aldehyde dehydrogenase 1 family member A2 gene (ALDH1A2) was associated with the growth trait in the Blanco Orejinegro cattle breed from Colombia [58], as well as intramuscular fatty acid composition in rabbits [59]. ALDH1A2 was also associated with lactation persistency in Canadian Holstein cattle [60], carcass traits in Chinese Simmental beef cattle [61] and litter traits in Landrace and Large White pigs [62]. Moreover, ALDH1A2 has been reported to play critical roles in the synthesis of retinoic acid, the active derivative of vitamin A, which is important for limb and organ development [63].
The protocadherin 7 gene (PCDH7) might be related to the body size trait of Chinese Simmental beef cattle [64], residual feed intake in Nelore cattle steers [65], feed efficiency trait in pig [66], internal organ traits in chickens [67] and sperm traits in Assaf rams [68].
The BovineHD4100014792-rs134139959 SNP located within F-box and leucine rich repeat protein 7 gene (FBXL7) was associated with the rear leg rear view trait in the present study, but, in previous studies, was associated with body length in water buffaloes (Bubalus bubalis) [69], clinical mastitis in first lactating US Holstein dairy cattle [70], subclinical ketosis in second and later lactations in Canadian Holstein dairy cattle [71] and also litter traits in Landrace and Large White pigs [62].
Furthermore, the results indicate the TMEM229A, POLE, XKR4 and CADPS genes were harboring significant SNPs correlated to the bone quality trait. In relation to the BovineHD0400024774-rs133088614 SNP, located close to the transmembrane protein 229A gene (TMEM229A), Juan Carlos [72] described the polymorphism of this gene as having a great effect on dairy yield per lactation in Colombian Holstein cattle. There are no further publications concerning this gene in cattle.
Functional analyses of list candidate genes showed a DNA polymerase epsilon (POLE) gene potentially associated with reproductive traits [73]; they might be implicated in milk yield regulation by performing certain biological functions [74], while our findings suggest that POLE is associated with the bone quality trait.
The BovineHD1400007035-rs136017102 SNP on Chr14, located close to the XKR4 gene, in previous studies, was reported to be associated with feed intake and growth phenotypes in cattle [75], residual feed intake in Australian Angus cattle [76], carcass weight (CWT) and eye muscle area (EMA) in Korean Hanwoo cattle [77], carcass trait [78], fat deposition in genome-wide association study for carcass trait in Brazilian Nellore cattle [79] and subcutaneous rump fat thickness in indicine and composite cattle [80].
The significant SNP BovineHD2200011035-rs110949452 on Chr22 is positioned within the calcium dependent secretion activator gene (CADPS); this gene is a member of the Ca +2 -dependent activator for the secretion protein family [81]. Interestingly, our study detected this gene is associated with the bone quality trait, while Vargas et al. [25] annotated the same gene within the top windows related to the overall quality of feet and leg in Nellore cattle. Kim et al. [82] reported this gene as selection signatures in Korean Brindle Hanwoo cattle obtained from genome-wide SNP analysis. Another study in two pig breeds identified the CADPS gene as associated with body weight [83].
In the present study, we found the ARS-BFGL-BAC-36389-rs109652453 SNP, which is associated with heel depth trait, to be located within synaptonemal complex protein 2, such as SYCP2L. Previously, most SYCP2L gene studies related to fertility. He et al. [84] reported a novel gene responsible for human premature ovarian insufficiency. For instance, SYCP2L was demonstrated to play a significant role in the survival of oocytes [85]. This gene is a paralogue of the synaptonemal complex protein gene (SYCP2), regulating females' reproductive aging. Therefore, it was suggested to have an association with the age of natural menopause in humans [85]. In addition, this gene is involved in the goat fecundity trait of Chinese Guangfeng goat [86].
Gene ontology enrichment analysis and KEGG pathways, in our study, revealed a number of GO terms and KEGG pathways. These involved significant genes related to the traits under examination, for example, the regulation of biological quality term (GO:0065008) in the biological process category, which has 22 genes ( Figure 5 and Table S1); among them, the CADPS, BARHL2, ADIPOR2 and ALDH1A2 genes were considered as nearest genes to our significant SNPs. According to the Mouse Genome Database (MGD) http://www.informatics.jax.org/ (accessed on 25 May 2021), this term and its sub-term (homeostatic process; GO:0042592) are related to any process that modulates a qualitative or quantitative trait of a biological quality, such as size, mass, shape, color, etc. As discussed earlier, some of the candidate genes involved in this term might play an important role in the biological regulation of body conformation traits and, subsequently, contribute to feet and legs traits.
Moreover, there were two terms (embryonic forelimb morphogenesis, GO:0035115, and forelimb morphogenesis, GO:0035136) that were related to the limb development term.
On another hand, the P2RX3 and P2RX2 genes were noticed in many biological process terms, such as urinary bladder smooth muscle contraction, GO:0014832, urinary tract smooth muscle contraction, GO:0014848, response to organophosphorus, GO:0046683, response to purine-containing compound, GO:0014074, peristalsis, GO:0022803, and in the molecular function category, such as channel activity (GO:0015267) and terms related with it (cation channel activity, gated channel activity, ATP-gated ion channel activity, extracellular ATP-gated cation channel activity and substrate-specific channel activity), and passive transmembrane transporter activity, GO:0022803. The P2RX3 and P2RX2 genes are subunits of the P2X gene family. P2X receptors are non-selective cation channels gated by extracellular ATP and exhibit relatively high Ca 2+ permeability [87]; the P2rx2 receptor mainly functions in sensory neurons, neuromuscular junction formation, whereas P2RX3, with multiple copies, were found in bony fishes [88].
The top two KEGG analysis pathways were bta00531: Glycosaminoglycan degradation, which has three genes, and bta01100: Metabolic pathways, which comprise 12 genes.
The functional analysis yielded many terms and pathways related to biological quality, metabolism and development. Therefore, it is reasonable to presume that all significant SNPs and candidate genes might be associated with feet and legs conformation.

Conclusions
In conclusion, this study identified 20 significant SNPs associated with feet and legs traits (bone quality, heel depth, rear leg side view and rear leg rear view) in Chinese Holstein cattle. Serval genes harbor SNPs (ADIPOR2, INPP4A, DNMT3A, ALDH1A2, PCDH7, XKR4 and CADPS), identified mainly to participate in biological quality, metabolism and development processes. Our findings provide useful information to understand genetic architecture and give some fundamentals for molecular-based breeding, leading to genetic improvement programs on feet and leg traits. Further investigations to validate the biological functions of these genes are recommended.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.