We performed single-step genome-wide association (ssGWA) analyses, which simultaneously combine all genotypic, phenotypic and pedigree information available13. Our genomic heritability estimates for BINP (0.21–0.24) were slightly higher than Schöpke et al.10 (0.19), who investigated a similar binary trait but for only one observation period using data on 729 Holstein heifers in the USA. For PROP, our genomic heritability estimates (PROP: 0.21–0.25) were lower than Biemans et al.11 (0.37), who also examined the fraction of observations a cow was DD free by scoring 1513 Holstein–Friesian cows in the Netherlands. Other heritability estimates in previous studies5,11,14 fell within the wide range of 0.01–0.40, reflective of differences in data, populations, trait definition, and methods of statistical analysis. Nevertheless, most studies agreed that development of DD can be attributed to both genetic contributions and environmental influences.
Our single locus analyses highlighted six SNPs on autosomes 7, 15, 16, 17, 19 and 21 associated with the studied traits. Two of these markers on chromosomes 7 and 15 affected both DD phenotypes, namely BINP at the LATE stage and PROP at the PEAK and LATE stages of the lactation cycle. Despite reaching a suggestive statistical significance level after multiple testing adjustments, the proportion of the additive genetic variance of the trait accounted for by each SNP was below 0.18%, suggesting a likely polygenic control of the DD phenotypes of study.
According to Aguilar et al.15, a p-value determines if an individual marker has an effect that is seemingly different from zero with statistical assessment, while the actual SNP effect size mainly associates with the attributable part of the genetic variance of the trait, with no statistical assessment. Therefore, in the present study, we examined further the genomic regions harbouring the significant SNPs, together with the top 10 genomic windows by amount of variance explained.
For the window-based ssGWA, we utilised the sliding windows approach which has been proposed as a robust method to detect combined significant patterns of clustered SNPs associated with complex polygenic traits16. Two key points emerged from our results. Firstly, the top 10 genomic windows accounting for the highest proportion of genetic variance of the trait were located in multiple autosomes across the entire genome. Secondly, regions harbouring the suggestive significant individual SNPs from single locus ssGWA did not overlap with the top 10 genomic windows by proportion of genetic variance explained. These observations further attest to the polygenic nature of the studied phenotypes, where multiple genes controlling DD are widely distributed across the genome, each with a relatively modest effect. This profile fits the infinitesimal model assumptions for complex trait analyses.
The identified genomic windows of the present study harboured positional candidate genes, whose functional relevance to DD, was explored with functional enrichment analysis. Based on the results of this analysis, the ALOXE3 (arachidonate lipoxygenase 3), ALOX12B (arachidonate 12-lipoxygenase, 12R type), ALOX12E (arachidonate lipoxygenase, epidermal), ALOX15 (arachidonate 15-lipoxygenase), ALOX12 (arachidonate 12-lipoxygenase, 12S type), CMKLR1 (chemerin chemokine-like receptor 1), C1QBP (complement C1q binding protein) and ARRB2 (arrestin beta 2) genes could be associated with DD development. Among them, the ALOXE3, ALOX12B, ALOX12E, ALOX15 and ALOX12 genes belong to the family of lipoxygenases, which play a known role in the modulation of epithelial proliferation and/or differentiation as well as in inflammation, wound healing, inflammatory skin diseases17. Indeed, the ALOXE3 gene has been reported in a recent study of sole ulcer, a non-contagious hoof condition in cattle18. The latter study reported the same genomic region in chromosome 19 that was identified here, encompassing the ALOXE3 (arachidonate lipoxygenase 3), HES7 (hes family bHLH transcription factor 7), U6 (U6 spliceosomal RNA), PER1 (period circadian regulator 1), VAMP2 (vesicle associated membrane protein 2), U8 (U8 small nucleolar RNA), TMEM107 (transmembrane protein 107), BORCS6 (BLOC-1 related complex subunit 6), AURKB (aurora kinase B), CTC1 (CST telomere replication complex component 1), PFAS (phosphoribosylformylglycinamidine synthase), RANGRF (RAN guanine nucleotide release factor), SLC25A35 (solute carrier family 25 member 35) and ARHGEF15 (Rho guanine nucleotide exchange factor 15) genes.
Amongst the other genes identified in the present study, CMKLR1 reportedly influences adipose tissue development, inflammation, and glucose homeostasis in mice19. In cattle, chemerin is a novel regulator of lactogenesis via its own receptor in mammary epithelial cells20. Moreover, C1QBP is involved in regulation of T cell proliferation21 and has been reportedly associated with milk protein percentage in cattle22. For ARRB2, a mouse study23 found that β-arrestin 2 inhibits the production of proinflammatory chemokines in cultured epidermal keratinocytes in response to T cell–derived cytokines in vitro and in T cell–driven allergic skin inflammation in vivo.
In addition to the aforementioned genes, the DDX10 (DEAD-box helicase 10) gene participating in the anterior head development has been previously associated with the foot-and-mouth disease in cattle24. A recent study25 reported functional variants in this gene potentially involved in innate immune response. Furthermore, six genes (KHDRBS3, RBL2, PDPN, RTTN, BSND and LCORL) reported here were also identified in a previous study on non-infectious claw horn disruptive lesions in the same cattle population, where KHDRBS3 (KH RNA binding domain containing, signal transduction associated 3) and RBL2 (RB transcriptional corepressor like 2) were associated with sole ulcer, PDPN (podoplanin), RTTN (rotatin) and BSND (barttin CLCNK type accessory subunit beta) with sole haemorrhage, and LCORL (ligand dependent nuclear receptor corepressor like) with white line disease26.
In the present study, DD phenotypes were defined at four distinct timepoints in a cow’s lactation. Each timepoint is associated with different metabolic and physiological state, conditions, and challenges for the animal. Certain overlap in suggestive significant SNPs (Supplementary Tables S1 and S2) and genomic windows (Supplementary Table S3) were observed across these timepoints. Moreover, genomic correlations between the DD phenotypes at different stages of lactation were sufficient high to imply a similar genomic profile controlling the disease development throughout the cow’s production cycle. Given the polygenic nature of the traits, the latter is the strongest result in this respect, suggesting that genetic selection to improve DD resistance at any timepoint would have an intertemporal benefit to foot health.
Instead of farmer or hoof trimming records used widely in previous studies5,27, we deployed here a large number of phenotypic records based on thorough inspections of lifted cow feet by a trained veterinarian. Although, admittedly, we involved a limited number of four farms, the focus of our study was on the accuracy and detail of individual clinical phenotyping for the in-depth analysis of DD.
In conclusion, results from the present study contribute to the understanding of the genetic mechanism underlying the DD lesion development in dairy cattle across distinct stages of the lactation cycle. Sizeable heritability estimates, combined with the largely polygenic mode of genetic control revealed in ssGWA analyses suggest that genomic selection based estimated breeding values would be the best approach to enhancing the inherent resistance of animals to DD infection. Genomic variants identified in the present study may be incorporated in and inform the genetic evaluation processes leading to increased accuracy of genomic selection. Our findings could contribute to future meta-analyses or multivariate analyses with larger samples size in order to detect genome-wide significant markers. Results may also serve as a basis for functional gene expression studies, aiming to identify causal genetic entities impacting on DD.
