Genome-wide association and functional genomic analyses for various hoof health traits in North American Holstein cattle

Hoof diseases are a major welfare and economic issue in the global dairy cattle production industry, which can be minimized through improved management and breeding practices. Optimal genetic improvement of hoof health could benefit from a deep understanding of the genetic background and biological underpinning of indicators of hoof health. Therefore, the primary objectives of this study were to perform genome-wide association studies, using imputed high-density genetic markers data from North American Holstein cattle, for 8 hoof-related traits: digital dermatitis, sole ulcer, sole hemorrhage, white line lesion, heel horn erosion, interdigital dermatitis, interdigital hyperplasia, and toe ulcer, and a hoof health index. De-regressed estimated breeding values from 25,580 Holstein animals were used as pseudo-phenotypes for the association analyses. The genomic quality control, genotype phasing, and genotype imputation were performed using the PLINK (version 1.9), Eagle (version 2.4.1), and Minimac4 software, respectively. The functional genomic analyses were performed using the GALLO R package and the DAVID platform. We identified 22, 34, 14, 22, 28, 33, 24, 43, and 15 significant markers for digital dermatitis, heel horn erosion, interdigital dermatitis, interdigital hyperplasia, sole hemorrhage, sole ulcer, toe ulcer, white line lesion disease, and the hoof health index, respectively. The significant markers were located across all autosomes, except BTA10, BTA12, BTA20, BTA26, BTA27, and BTA28. Moreover, the genomic regions identified overlap with various previously re-ported quantitative trait loci for exterior, health, meat and carcass, milk, production, and reproduction traits. The enrichment analyses identified 44 significant gene ontology terms. These enriched genomic regions harbor various candidate genes previously associated with bone development, metabolism, and infectious and immunological diseases. These findings indicate that hoof health traits are highly polygenic and influenced by a wide range of biological processes.

Appropriate management practices can reduce the prevalence of hoof diseases in dairy herds (Meskini et al., 2023).Previous reports have indicated that hoof health-related traits are heritable, with heritability estimates ranging from 0.01 to 0.14 (Malchiodi et al., 2017,

Genotype Imputation and Quality Control
As a first step, genotype imputation was performed from a MD SNP panel containing 44,315 SNPs to a HD SNP panel containing 311,725 SNPs.A total of 25,580 animals (11,182 females and 14,398 males) had MD data, whereas the reference HD population contained 2,507 North American Holstein animals (562 females and 1945 males).Before genotype imputation, the SNPs from the MD SNP panel that were not present in the HD panel were excluded.Moreover, quality control (QC) was performed using the PLINK 1.9 software (Purcell et al., 2007) to exclude SNPs: (1) with a call rate <0.95; (2) with an extreme departure from Hardy-Weinberg equilibrium (P < 10 −8 ), as an indication of genotyping errors; (3) located in nonautosomal chromosomes; and (4) with unknown position based on the ARS-UCD1.2cattle genome assembly.After QC,41,612 and 296,456 SNPs remained in the MD and HD panels, respectively.The genotype phasing step was performed using Eagle 2.4.1 software (Loh et al., 2016), and the genotype imputation was done by Minimac4 software (Das et al., 2016).After genotype imputation, an additional QC was performed to exclude individuals or genotypes with call rate <0.90, SNPs with minor allele frequency < 0.01, as well as SNPs with extreme departure from Hardy-Weinberg equilibrium (P < 10 −8 ).Finally, the number of SNPs that remained for the GWAS analyzes was 282,434 for DD, 280,898 for HHE, 282,827 for HH, 282,829 for ID, 282,492 for IH, 280,959 for SH, 280,844 for SU, 280,854 for TU, and 280,947 for WL.

Genome-Wide Association Analyses
A mixed linear model was used to estimate the SNP effects using the GCTA package (Yang et al., 2011).The univariate model can be described as follows: y = μ + Xβ + Za + e, where y is the vector of dEBVs for each trait; μ is the overall mean; β is the fixed effect of the SNP being tested for statistical association with each trait; a is a vector of random polygenic effects where G is the genomic-based relationship matrix (GRM; VanRaden, 2008) and σ a 2 is the additive genetic variance; X and Z are incidence matrices for the effects in β and a, respectively; and e is a vector of residuals e N e ∼ 0 2 , , Iσ ( )       where I is an identity matrix and σ e 2 is the residual variance.When using dE-BVs (with different reliabilities) as phenotype in GWAS, it is recommended to account for the differential residual variances (Jiang et al., 2019).However, in this specific study, EBVs with reliability lower than 0.50 were filtered out and the dEBVs used for the GWAS analyses had high average reliabilities with a reasonably small average coefficient of variation across traits of 10.55% (Table 1).Under these circumstances, differential weighting of the residuals was assumed to have a small effect on the results.
Mixed linear model-based association analyses were performed using the GCTA software, excluding the chromosome harboring the candidate SNP being tested before calculating the GRM (Yang et al., 2014).This method is computationally less efficient but more powerful as compared with the mixed linear model analysis including the chromosome of the candidate SNP (Yang et al., 2014).Therefore, a total of 29 GRM were alternatively constructed by randomly sampling 50,000 SNP from all chromosomes (after QC) to avoid overcorrection for population structure.After the GWAS analyses, all SNPs were ranked based on their P-values and clumped according to their linkage disequilibrium (LD) values (r 2 threshold of 0.9), which has been suggested as a preferable method compared with the traditional LD pruning strategy (Privé et al., 2018).The genomic inflation factor (λ) was calculated as λ = median(χ 2 )/0.456 (Bacanu et al., 2000), for which a 95% confidence interval of λ value was further derived.

Multiple Testing Correction
As a large number of variants were tested, the traditional Bonferroni correction would be highly conservative because not all tests are independent due to LD among markers (Johnson et al., 2010).Thus, to avoid excessive false negative results, we applied a modified Bonferroni correction by using the number of independent chromosome segments (Me) at the genome-wide level (Li et al., 2015) instead of the total number of tests.The Me is a function of both effective population size (Ne) and genome length (L) in Morgans, which can be calculated as Me = (2 × Ne × L)/log(Ne × L) (Goddard et al., 2011).Based on current literature reports across mammal species and the lack of more accurate conversion values for cattle, 1 cM was considered to be equivalent to 1 Mb (Wang et al., 2016), and the Ne value used was 66, as it is the most conservative Ne value reported for the same Holstein cattle population (Makanjuola et al., 2020).The average level of LD in the studied population based on |D'| and r 2 for markers 40 to 60 kb apart was 0.72 and 0.20, respectively (Bohmanova et al., 2010).A SNP effect was considered to be statistically significant if its −log 10 (P-value) was higher than the genome-wide threshold, which was calculated by dividing 0.05 by Me.The value of Me and the corresponding significance threshold used were 2,044.73 and 4.61, respectively.

Functional Genomic Analyses
The SNP coordinates were based on the ARS-UCD1.2assembly of the cattle reference available genome in the GenBank accession (GCA_002263795.2).The GALLO R package (Fonseca et al., 2020) was used to detect positional genes and quantitative trait loci (QTL) located 100 Kb upstream or downstream from the significant SNPs.The QTL database used was the Animal QTLdb release 49 (Hu et al., 2019).Subsequently, functional enrichment analyses were performed for the group of genes found for each trait using the DAVID platform (Huang et al., 2009).Multiple testing was accounted for based on the false discovery rate of 0.05.
The list of significant SNPs and candidate genes found for each trait are shown in Tables 2, 3, 5, and 6.Detailed information about the GWAS, including Pvalues, SNP effects, minor allele frequency, QTL, and candidate genes found within an interval of 100 Kb upstream and downstream from the significant SNPs are provided in Supplemental File S2 (https: / / figshare .com/s/ 8fb6aba738fdc041cd05).The number of significant (P-value < 2.45 × 10 −5 ) SNPs was 15 for HH, 22 for DD, 34 for HHE, 14 for ID, 22 for IH, 28 for SH, 33 for SU, 24 for TU, and 43 for WLD.We did not find any significant SNPs located on BTA10, BTA12, BTA20, BTA26, BTA27, and BTA28 for any of the traits evaluated, but BTA14 had significant SNPs for all traits.The rs110059642 SNP (BTA5) was the only one associated with more than one trait (DD and IH) and 7 genes are located close to this SNP.In addition, 28 genes were associated with 2 other traits, 5 with DD and HHE, 4 with DD and ID, 14 with IH and SH, and 5 with SU and WLD (Supplemental File S2).Nine of these genes are located on BTA3, 7 on BTA5, 5 on BTA13 and 14 on BTA19.
A total of 2,473 QTL were previously reported to overlap with the genomic regions identified with significant markers for HH index, infectious and noninfectious diseases on many chromosomes, with the exception of BTA10, BTA12, BTA20, BTA26, BTA27, and BTA28 and they can be classified into different QTL groups (exterior, health, meat and carcass, milk, production,    and reproduction).Additional details about these QTL are presented in the Supplemental File S2.

Hoof Health Index
The GWAS for HH found 15 significant SNPs distributed across the chromosomes BTA3, BTA5, BTA8, BTA11, BTA14, BTA15, and BTA17 (Supplemental File S2).The most significant SNP (rs3423383926) for HH was found on BTA14: 76 ,259 ,291 bp.Eleven candidate genes with known biological functions (Table 2) were found near the significant SNPs, but no genes were found close to the significant SNPs located on BTA11, BTA15, and BTA17.The functional enrichment analyses did not identify any significant gene ontology (GO) terms for HH.

Hoof Health Index
The HH index is a weighted linear combination of hoof health traits used as a selection index for genetically improving hoof health in Canadian Holstein cattle.Eleven candidate genes were found harboring or near the 7 chromosomes associated with this index.The NTNG1 (Netrin G1) gene located on BTA3 encodes a member of the Netrin protein family.The netrins play key roles in organogenesis and tissue morphogenesis, including directing cell migration, cell-cell interactions, and cell-extracellular matrix adhesion, processes that are important in wound repair (Iorio et al., 2015).A SNP (rs135185485; BTA3: 36,106,918 bp) located within this gene lies close to a previously reported QTL (BTA3: 36 ,159 ,159 ,482) associated with digital cushion thickness (Stambuk et al., 2020), indicating a high LD between this SNP and the QTL.The gene SLC2A13 (Solute Carrier Family 2 Member 13) located on BTA5 is expressed in the basal layer of the epidermis and encodes transporters of myoinositol, which influences osmoregulation of primary epidermal keratinocytes (Foster et al., 2020).The genes GLIPR1 (GLI Pathogenesis Related 1) in BTA5 and CPNE3 (Copine 3) on BTA14 are associated with important  health-related pathways such as immune system, innate immune system, and neutrophil degranulation.
The gene KLF9 (KLF Transcription Factor 9), located on BTA8, positively regulates TRIM33 to inhibit abnormal fibroblast proliferation, migration, and inflammation (Huang et al., 2022).Fibroblasts play a key role in synthesizing new connective tissue and contracting it to promote wound healing (Shephard et al., 2004).The WWP1 (WW Domain Containing E3 Ubiquitin Protein Ligase 1) gene may also be involved in innate immunity (Zhi and Chen, 2012), and the ubiquitin proteasome system is a key regulatory factor in the proliferation, differentiation, and collagen secretion of skin fibroblasts (Shen et al., 2021).
Various QTL (Supplemental File S2) were located in the genomic regions around the SNPs significantly associated with HH trait.Heel depth is an exterior QTL type that was associated with a much higher than average incidence of lesions when heel depth was extremely shallow in Canadian Holstein cows (Chapinal et al., 2013).Other QTL types were related to health (bovine tuberculosis susceptibility, clinical mastitis, gastrointestinal nematode burden, immunoglobulin G level, somatic cell count, and tick resistance), meat and carcass (connective tissue amount), and production (average daily gain, body height, body weight, dry matter intake, feed conversion ratio, length of productive life, and residual feed intake).The fact that the genomic regions associated with hoof health diseases have also been related to other relevant traits indicate the biological complexity of the mechanisms involved in resistance to hoof diseases.
Health traits may help explain the relationship between these QTL and hoof diseases.Susceptibility to bovine tuberculosis may be associated with the efficiency of the innate and adaptive immune systems in the host defense against tuberculosis and mycobacteria (Koul et al., 2004;Raman et al., 2010).Furthermore, cows with lameness have an increased incidence of mastitis compared with healthy cows (Souza et al., 2006).It was previously reported that the prevalence of gastrointestinal parasites in donkeys was associated with lameness (Ayele et al., 2006) and the IgG level in the plasma of dairy cows with high levels of lameness significantly increased, indicating that humoral immunity is activated by the stressor lameness (Sun et al., 2015).Krpálková et al. (2019) reported that the prevalence of hoof disease detected in the first month of first and second lactations was associated with increases of 58,000 and 45,000 cell/mL of somatic cell count in dairy cows, respectively.Various genomic regions identified in this study were previously associated with tick resistance in cattle (Supplemental File S2).Ticks are among the most relevant vectors of diseases affecting livestock, and selection for resistance may be a complementary tool to assist parasite control in Nellore cattle (Passafaro et al., 2015).
Two SNPs (rs134623665-BTA8: 46 ,514 ,143 bp and rs3423383926-BTA14: 76 ,259 ,291 bp) associated with the hoof health index are located near to the QTL regions previously reported to be associated with the amount of connective tissue.This previously reported region is located on BTA8 (46,557,557,877) and on BTA14 (76,171,71,344), respectively (Supplemental File S2).The connective tissue is essential for normal growth and development, but many diseases have long been associated with the breakdown of the collagenous matrix of bone, cartilage, and related tissues in humans (Reynolds et al., 1994).With respect to production trait QTL, Schöpke et al. (2013) reported that heifers between 50 and 99 DIM had the highest frequencies of laminitis and increased biomechanical stress caused by different factors (weight, social rank, standing time), which was a presumed cause of the increased susceptibility for hoof diseases in Holstein cows.This study found QTL regions (BTA3: 31 ,970 ,767 ,592;BTA11: 15 ,887 ,668 ,573 and BTA15: 17 ,567 ,415 ,621) previously associated with weaning weight in beef cattle (Kim et al., 2003;McClure et al., 2010).Heavy animals tend to have greater probability of hoof diseases (Pérez-Cabal and Charfeddine, 2016).In a survey to identify the reasons for hoof disorders in crossbred dairy cattle, Kumar et al. (2019) reported that body weight greater than 410 kg, injured hocks, and animals with abnormal claws were important risk groups for hoof disorders with higher incidence (43.0%, 51.7%, and 32.6%, respectively).Thus, body weight can influence the incidence of hoof diseases, as cows have a limited ability to shift weight from front to back, and is an important follow-up measure for farmers (Neveux et al., 2006).

Infectious Hoof Diseases
Digital dermatitis, ID, and HHE are 3 infectious hoof-related diseases (Knappe-Poindecker et al., 2013).Digital dermatitis is transmitted via the environment by an infectious agent that might infect keratinized cells in the epidermal layers, resulting in the production of a keratolytic toxin, which stimulates the epidermal proliferation and hyperplasia (Biemans, 2018;Bay et al., 2023).Digital dermatitis has a high prevalence (13.3%) in Canadian Holstein cattle (Lactanet, 2022).Heel horn erosion and ID are usually more frequent in freestall system herds with wet and unhygienic environments (Knappe-Poindecker et al., 2013).The present study identified various SNPs associated with these traits (Table 3), including clear Manhattan plot peaks (Figure 2) located on BTA3 (103.7-113.1 Mb) for all the 3 traits, on BTA14 (46.1-54.5 Mb) for HHE, and on BTA23 (24.04-27.4Mb) for DD.
Various SNPs located on BTA3 (105.4-108.4Mb), which were associated with DD, overlap with previously reported QTL for clinical mastitis in dairy cattle (Lund et al., 2008;Sahana et al., 2008;Sodeland et al., 2011).This is an interesting result as DD and mastitis are moderately genetically correlated (0.49; Lai et al., 2021a).This genomic region also contains another SNP (rs136841499; 108.4 Mb) associated with HHE, and 6 SNPs associated with ID are located between 103.7 and 113.1 Mb on BTA3 (Supplemental File S2).Therefore, further study is needed to investigate this genomic region through a finer screening using whole-genome sequence data.Some genes on this BTA3 region were previously associated with diseases caused by bacteria.For instance, ZC3H12A (zinc finger CCCH-type containing 12A) has been linked with mastitis in dairy cattle (Sharifi et al., 2018;Umesh et al., 2022) and 2 close-by SNPs (i.e., rs134482817 and rs137184511) overlap with a previously reported QTL for HHE (Supplemental File S2).The genes CITED4 (Cbp/p300 interacting transactivator with Glu/Asp rich carboxy-terminal domain 4) and NFYC (nuclear transcription factor Y subunit gamma) were found near variants associated with Mycobacterium avium subspecies paratuberculosis infection status in Holstein cattle (Settles et al., 2009).Variants in the CITED4 gene were also associated with monocyte count in humans (Sakaue et al., 2021).RIMS3 (Regulating synaptic membrane exocytosis 3) is also associated with host-pathogen adaptation to Mycobacterium tuberculosis (Phelan et al., 2023).The genes MFSD2A, MYCL, TRIT1, and PPIE were found near a SNP (rs110237486) associated with mastitis prevalence in dual-purpose German Black Pied cattle (Meier et al., 2020), and PPIE (Peptidylprolyl Isomerase E) encodes a member of PPIases family, which are associated with regulation of innate and adaptive immunity functions (Nath and Isakov, 2015).
Another important genomic region (BTA3: 24.05-27.41Mb) was associated with DD (Table 3).This region harbors genes from the major histocompatibility complex, which is known as the Bovine Leukocyte Antigen (Takeshima and Aida, 2006) and is related to immune regulation against infectious diseases.SNPs located in this region of BTA23 were previously associated with Holstein cow livability and infectious diseases in cattle (Freebern et al., 2020).IL17A (Interleukin 17A) and IL17F (Interleukin 17F) are members of the IL-17 family and influence host defense, cell trafficking, immune modulation, and tissue repair, with important functions in the induction of innate immune defenses (Ge et al., 2020).The activation of the proinflammatory IL-17 signaling pathway may influence DD progression in Holstein-Friesian cattle (Vermeersch et al., 2022).Moreover, IL17A has been linked with cattle tuberculosis (Waters et al., 2016) and somatic cell count in Chinese Holstein and Inner-Mongolia Sanhe cattle (Usman et al., 2017).SNPs located in the IL17A and IL17F genes were also associated with milk production in Chinese Holstein (Ghulam-Mohyuddin et al., 2022), indicating potential pleiotropic QTL in this region.The C2 and C4A genes were previously associated with inflammation response (Alper et al., 2003;Korwin-Kossakowska et al., 2022) and C4A was linked to gene expression of inflammatory mechanisms in the parenchyma of the mammary gland of dairy cows (Korwin-Kossakowska et al., 2022).CFB is associated with mastitis in cattle (Sharifi et al., 2018) and bovine respiratory disease in commercial crossbred cattle (Scott et al., 2021).
Two pathways (complement and coagulation cascades, and Staphylococcus aureus infection) were significantly enriched for DD.Complement and coagulation cascades is a plasma protein system that is activated by pathogen infection (Markiewski et al., 2008).The 3 main consequences of this complement activation are: 1) opsonization of pathogens; 2) recruitment of inflammatory and immune competent cells; and 3) direct killing of pathogens (Bajic et al., 2015).Staphylococcus aureus is a pathogen that often causes skin and soft tissue infections.Its virulence is a consequence of suppression of innate and adaptive immune responses (DeDent et al., 2012).Staphylokinase, secreted by S. aureus, promotes the establishment of skin infections in humans and increased bacterial penetration through skin barriers by activating plasminogen (Kwiecinski et al., 2013).Therefore, the significant biological processes for DD (innate immunity and complement pathway) may influence the cow's susceptibility to hoof health problems.
Twenty-two SNPs located between 46.1 and 54.5 Mb on BTA14 were associated with HHE.This important genomic region overlaps with previously reported SNPs located in a region between 48.1 and 49.9 Mb that were associated with Mycobacterium avium ssp.paratuberculosis infection in Holstein cattle (Mallikarjunappa et al., 2018).Moreover, other candidate genes were found in this genomic region.For instance, TRPS1 (Transcriptional Repressor GATA Binding 1) is a key regulatory gene that influences wound regeneration (Mascharak et al., 2022).EXT1 (Exostosin Glycosyltransferase 1) is related to fibroblast growth factor, which can mediate various processes, including wound healing (Liu et al., 2021).EXT1 was found to be a candidate gene for feet and legs malformation in Nellore cattle (Silva et al., 2022).RAD21 (RAD21 Cohesin Complex Component) participates in cellular senescence in dermal fibroblast and its expression was associated with the proliferative ability of fibroblasts in the dermis (Okuno et al., 2022).SAMD12 (sterile α motif domain containing 12) is involved in transmembrane receptor protein tyrosine kinase signaling pathway, which has immune response functions, including mastitis resistance in Holstein cattle (Wagner et al., 2021).CSMD3 (CUB and Sushi Multiple Domains 3) is a candidate gene for body size in cattle (An et al., 2019;Ghoreishifar et al., 2020).This is relevant as a previous study reported association between body weight and claw disorders in Spanish Holstein cows (Pérez-Cabal and Charfeddine, 2016).
The PLG (Plasminogen) gene located on BTA9 contains significant markers associated with ID.PLG encodes an active serine protease plasmin, the principal enzymatic mediator of fibrinolysis that plays a fundamental role in the dissolution of key proteins involved in immunity and tissue repair (Seillier et al., 2022).Various QTL were found exclusively for infectious hoof diseases (Supplemental File S2).The QTL for DD (BTA21: 60 ,781 ,781 ,779) and HHE (BTA3: 108 ,428 ,428 ,619) were also found by Croué et al. (2019) supporting the association of the QTL regions found in the current study with the infectious hoof health diseases.Other health traits were previously associated with SNPs located in the genomic regions identified in this study, including antibody-mediated immune response, bovine leukemia virus susceptibility, bovine spongiform encephalopathy, final packed red blood cell volume, general disease susceptibility, initial packed red blood cell volume, IL-4 level, mean corpuscular hemoglobin concentration, and white blood cell number (Supplemental File S2).The possibility of transmission of bovine leukemia virus through equipment from the hoof trimmer has been reported (Kuczewski et al., 2022).Finally, the QTL found in the highlighted genomic regions in the current study are also associated with diseases that affect the immune response (Waldvogel et al., 2000) and hematological changes resulting from infections (Bengaly et al., 2002).

Noninfectious Hoof Diseases
Five noninfectious hoof diseases (SH, SU, TU, DH, and WLD) were included in this study.The hoof capsule is formed by cornified epidermal tissue known as horn, which is described in descending order of hardness, as wall, sole, heel, and white line.Horn is produced by the keratinocytes of the epidermis by the specialized process known as keratinization (Hoblet and Weiss, 2001).For SH, 13 keratin family genes (LOC100295015, LOC618495, LOC100294937, LOC788284, LOC787225, KRTAP1-1, KRTAP3-3, LOC540285, KRTAP3-1, LOC618938, KRT40, KRT23, and KRT39) were found near the SNP (rs29021711; BTA19: 41 .17Mb).This gene family was enriched for 2 significant cellular components ("Keratin filament" and "Keratin") and one molecular function ("Structural molecule activity"; Table 7).Keratins are the characteristic structural proteins of the highly cornified epidermis of the skin, feathers, and hoof (Tomlinson et al., 2004).Keratin is responsible for forming fibrous proteins rich in cysteine, which are the main constituent of the horn, nails, hair, Sousa et al.: GENOMICS OF HOOF HEALTH TRAITS epidermis, and feathers (Patrucco et al., 2019).The keratin family genes regulate various human clinical disorders, including epidermolytic ichthyosis, superficial epidermolytic ichthyosis, epidermolysis bullosa simplex, palmoplantar keratoderma, and white-sponge nevus (Ho et al., 2022).In cattle, the keratin family genes were suggested as candidate genes for susceptibility and resilience to paratuberculosis (Alonso-Hearn et al., 2022), and a mutation in the bovine keratin-5 gene (KRT5) was found as the causal variation for epidermolysis bullosa simplex (Ford et al., 2005), a condition presented as loss of skin and mucosa from contact areas and inflammation in Friesian-Jersey crossbred bulls.In an earlier study, Lai et al. (2021a) reported genes related to keratinization as candidate genes for various claw lesions in Holstein cattle.
Another group of genes (STAT5A, STAT5B, and STAT3) located near the SNP on BTA19 (rs3423466061; 42.4 Mb) are candidate genes for SH and IH.The STAT (Signal Transducers and Activators of Transcription) family of proteins contains transcription factors that are specifically activated to regulate gene transcription when cells encounter cytokines and growth factors (Banerjee and Resat, 2016).STAT proteins regulate host innate and acquired immune responses (Kisseleva et al., 2002).STAT3 is a key transcriptional regulator of genes associated with innate immunity and wound healing (Hughes et al., 2012), and a previous study reported an interaction between STAT3 and keratin-17 to increase keratinocyte proliferation in humans (Yang et al., 2018).A variant in STAT5A was associated with the serum biochemical assays of IL-6 in Holstein cattle (Usman et al., 2014), a cytokine that acts in innate and adaptive immune responses and is an important inflammatory biomarker.Another variant (rs134993207; BTA19: 43 ,038 ,655) in the STAT5 gene was associated with mastitis in Nordic Holstein cattle (Cai et al., 2018).The biological process JAK-STAT cascade involved in the growth hormone signaling pathway was associated with the STAT5A, STAT5B, and STAT3 genes (Table 7).In research carried out by Brown et al., (2004), a clinical improvement in hoof pathologies of Holstein cows was detected when regulation of growth hormone secretion was used, through the GHRH gene (Growth hormone-releasing hormone).It was also found the association of growth hormone with some traits in Raini Cashmere Goats (Gooki et al., 2018) such as toe height of hand hoof, toe height of foot hoof and heel height of foot hoof.
Four genes (BPIFA2A, BPIFB3, BPIFB4, BPIFB6) of the BPI-fold (BPIF) superfamily located near the same SNP (rs110484302, BTA13: 62 .4Mb) were associated with SU.The BPIF superfamily has been associated with key health-related pathways, such as antimicrobial peptides, immune system, and innate immune system in humans (Theprungsirikul et al., 2021).Three other candidate genes were found to be associated with SU.POFUT1 (Protein O-Fucosyltransferase 1), located on BTA13 (61.69-61.74Mb), influences mammalian notch signaling (Stahl et al., 2008), in which the activation or inhibition of notch signaling may alter the functions of vascular endothelial cells, keratinocytes, and fibroblasts in a wound healing model (Chigurupati et al., 2007).The rs109265528 SNP (BTA13: 58 ,741 ,419) was also found to be associated with a previously reported QTL for SU in Holstein cattle (Supplemental File S2).CD36 (CD36 Molecule), located on , is a multifunctional glycoprotein that acts as receptor for a broad range of ligands with important functions in TLR2/CD36-p38 MAPK signaling and can improve the defense against bacterial infection (Li et al., 2013).ALOXE3 (Arachidonate Lipoxygenase 3) is a member of the lipoxygenases family, which are key enzymes in the biosynthesis of a variety of highly active oxylipins.Oxylipins regulate the modulation of epithelial proliferation and differentiation and are involved in inflammation, wound healing, inflammatory skin diseases, and in the maintenance of the skin permeability barrier (Krieg and Fürstenberger, 2014).
Enrichment analyses found no significant GO terms for TU, but functional annotation analyses identified candidate genes for this trait.Three genomic regions located on BTA8 (Table 6) were associated with TU.The first of them (22.0-22.1 Mb) overlaps with the genes CDKN2A (Cyclin Dependent Kinase Inhibitor 2A) and CDKN2B (Cyclin Dependent Kinase Inhibitor 2B), which are checkpoint proteins for cell cycle arrest (Wiesner et al., 2010).An accumulation of CDKN2A protein in keratinocytes is observed in normal-aged skin as well as neoplastic and inflammatory skin disorders (Hashimoto-Hachiya et al., 2019).The second region of BTA8 (31.4 to 31.6 Mb) harbors the TYRP1 (Tyrosinase Related Protein 1) gene, which is related to melanin biosynthesis (Zhang et al., 2022) and peri-wound melanin levels may influence wound healing (Rowledge et al., 2016).The third region (BTA8: 110.6-112.6Mb) harbors the PXDN (Peroxidasin) gene, which participates of collagen-related pathways, including assembly of collagen fibrils and other multimeric structures, collagen formation, crosslinking of collagen fibril, and extracellular matrix organization.Collagen plays critical role in the regulation of the phases of wound healing in its native fibrillar conformation or as soluble components in the wound milieu (Mathew-Steiner et al., 2021).
Five genes (HBZ, HBM, HBQ1, HBA1, HBA) had significant GO terms for WLD (Table 7).These genes are associated with 2 biological processes (hydrogen peroxide catabolic and cellular oxidant detoxification) and 2 cellular components, including haptoglobinhemoglobin complex and hemoglobin complex.Hemoglobin has the function of transporting oxygen in the blood plasma, but an inadequate intake of minerals can affect the synthesis of hemoglobin and the increase of toxic effects of oxygen metabolites (NRC, 2001).The inadequate intake of minerals is associated with malformation of keratins in cattle, contributing to hoof disorders (Tomlinson et al., 2004).Human patients with diabetes mellitus and foot disorders had elevated levels of hemoglobin (Li et al., 2021).The AXIN1 (Axin 1) gene, located on BTA25, encodes a component of the β-catenin destruction complex, which contributes to protection against various human pathogenesis.Moreover, the degradation and plasma sequestration of Axin-1 increase bacterial invasiveness and inflammatory responses (Jin et al., 2017).β-catenin is involved with protection of the epidermis from mechanical stresses (Ray et al., 2013).DNMT3B (DNA Methyltransferase 3 Beta) encodes a primary methyltransferase enzyme that catalyzes DNA methylation to establish cell fates during development and plays important roles in various adult tissues, including the epidermis (Rinaldi et al., 2016).Moreover, DNMT genes regulates DNA methylation in inflammatory cells, keratinocytes, and fibroblasts with consequent effects on wound healing (Pastar et al., 2021).PSMA6 (Proteasome 20S Subunit Alpha 6) is involved in adaptive immune system pathways (Chiao et al., 2021).It encodes a protein involved in the major histocompatibility complex class I antigen processing, and its disruption generates autoantigens or other proinflammatory antigens for presentation to pathogenic CD8+ T cells (Sukhov et al., 2016).This gene is located near a QTL previously associated with SH in Holstein cattle (Supplemental File S2).SLC38A2 (Solute Carrier Family 38 Member 2) encodes an amino acid transmembrane protein, and the amino acids influences the structural integrity of the hooves in cattle (Langova et al., 2020).BRINP1 (BMP/Retinoic Acid Inducible Neural Specific 1) inhibits cell proliferation, mediates cell death, and regulates expression of components of the plasminogen pathway.BRINP1 expression fell as the urea concentration increase in Holstein-Friesian cows (Cheng et al., 2015), and high concentrations of urea in the blood can damage sensitive lamellae and corium in the hoof (Langova et al., 2020).Various QTL were found exclusively for noninfectious hoof diseases (Supplemental File S2) for the exterior type (IH, SH, SU, body condition score, and bone quality).In addition, QTL were found for health (ConA-induced cell proliferation), meat and carcass (bone percentage and bone weight), and production (body size).Concana-valin A (ConA) is an antigen-independent mitogen and functions as a signal one inducer, leading T cells to polyclonal proliferation (Ando et al., 2014), being an important immunological response pathway.Fiore et al. (2019) identified that osteolysis affecting the pedal bone is associated with hoof lesions.Additionally, cows selected for larger body size were more often found to have leg and foot problems in Holstein cattle (Hansen et al., 1999).
Only a single marker was the same between infectious and noninfectious diseases (rs110059642; BTA5: 30 ,490 ,039), which was associated with DD and IH and located close to 7 genes (Tables 3 and 5).Among these genes, C1QL4 and TUBA1C (Table 4) showed significance for GO terms for infectious disease.These results demonstrate that the traits studied are highly polygenic and influenced by a wide range of biological mechanisms.

Implications and Next Steps
This study identified numerous genomic regions and candidate genes associated with hoof health traits in Holstein cattle.Some genomic regions overlap across traits, such as BTA3 (105.5-108.4Mb) for DD, HHE, and ID; BTA13 (62.1-64.2Mb) for TU, SU, and WL; BTA14 (45.2-47.8Mb) for HH, HHE,ID,and IH;and, for IH and SH.The GO terms and biological pathways indicate that the results found are biologically relevant.They have been associated with pathways of immune regulation, inflammatory activation, tissue repair, and collagen formation.Thus, they are directly related to the mechanisms of expression for symptoms of infectious and noninfectious hoof health diseases.The SNPs identified could be added to lower-density SNP panels to potentially increase genomic prediction accuracies.The next steps would be to perform GWAS based on whole-genome sequence genotypes, focusing on the regions identified in this study and using additional omics approaches, such as transcriptomics, for validating the role of the SNPs and candidate genes identified.In addition, performing validation in independent populations using larger data sets is recommended to evaluate the consistency of the results.

CONCLUSIONS
The genome-wide association analyses identified SNPs significantly associated with infectious and noninfectious hoof diseases and with the overall hoof health index trait on most chromosomes.The rs110059642 SNP (BTA5: 30,490,039) was associated with DD and Sousa et al.: GENOMICS OF HOOF HEALTH TRAITS IH, and 28 genes were associated with 2 trait pairs.The regions around the significant SNPs overlapped with various previously reported QTL for exterior, health, meat and carcass, milk, production and reproduction traits.A total of 235 genome-wide significant SNPs and 44 gene ontology terms and biological pathways were enriched.Many candidate genes involved in biological pathways associated with bone development, metabolism, and infectious and immunological diseases were identified.These findings indicate that hoof health traits are highly polygenic and influenced by a wide range of biological processes.Further studies with other independent populations and with whole-genome sequence data are warranted to validate the identified candidate genes, especially those in the regions of overlapping markers (located on BTA3, BTA5, BTA13, BTA14, and BTA19) and with more than one gene associated with the traits studied.

Figure 2 .
Figure 1.Manhattan plot of the genome-wide association analysis for hoof health index using imputed high-density SNP genotypes in Canadian Holstein cattle.The statistically significant SNPs after a genome-wide Bonferroni correction are colored in red.The significance threshold is indicated by the horizontal black line.

Figure 3 .
Figure 3. Manhattan plot of the genome-wide association analysis for interdigital hyperplasia (a), sole hemorrhage (b), and white line disease (c) using imputed high-density SNP genotypes in Canadian Holstein cattle.The statistically significant SNPs after a genome-wide Bonferroni correction are colored in red.The significance threshold is indicated by the horizontal black line.

Figure 4 .
Figure 4. Manhattan plot of the genome-wide association analysis for sole ulcer (a) and toe ulcer (b) using imputed high-density SNP genotypes in Canadian Holstein cattle.The statistically significant SNPs after a genome-wide Bonferroni correction are colored in red.The significance threshold is indicated by the horizontal black line.
chromosome; MAF = minor allele frequency.Table 5 (Continued).Single nucleotide polymorphisms associated with interdigital hyperplasia (IH), sole hemorrhage (SH), and white line disease (WLD), and the corresponding positional genes located 100 Kb upstream or downstream from the significant SNPs for each trait 1 Sousa et al.: GENOMICS OF HOOF HEALTH TRAITS

Table 2 .
The SNPs associated with hoof health index and the corresponding positional genes located 100 Kb upstream or downstream from the significant SNPs

Table 4 .
Sousa et al.:GENOMICS OF HOOF HEALTH TRAITS Significant terms in the functional annotation analysis for group of genes located 100 Kb upstream or downstream from the significant SNPs for digital dermatitis in 1 KEGG_PATHWAY = (Kyoto Encyclopedia of Genes and Genomes; http: / / www .genome.ad.jp/kegg/ ) is a collection of databases dealing with genomes, biological pathways, diseases, drugs, and chemical substances; UP_KW_BIOLOGICAL_PROCESS = (UniprotKB Keywords; https: / / www .uniprot.org/help/keywords)keywordsassigned to proteins because they are involved in a particular biological process; UP_KW_DOMAIN = (UniprotKB Keywords; https: / / www .uniprot.org/help/keywords)keywordsassigned to proteins domain; UP_KW_MOLECULAR_FUNCTION = (UniprotKB Keywords; https: / / www .uniprot.org/help/keywords)keywordsassigned to proteins due to their particular molecular function; FDR = false discovery rate.Sousa et al.: GENOMICS OF HOOF HEALTH TRAITSTable5.Single nucleotide polymorphisms associated with interdigital hyperplasia (IH), sole hemorrhage (SH), and white line disease (WLD), and the corresponding positional genes located 100 Kb upstream or downstream from the significant SNPs for each trait 1 Trait

Table 6 .
Single nucleotide polymorphisms associated with sole ulcer (SU) and toe ulcer (TU), and the corresponding positional genes located 100 Kb upstream or downstream from the significant SNPs for each trait 1 Trait