All cows considered for this study belonged to the Holstein breed, the worldwide dominating dairy breed, and were from commercial herds in northeastern Germany. Phenotypes were collected as scores based on visual inspection and thus completely non-invasive. Cows were not specifically restrained while scoring but rather the possible occurrence of BDD lesions was scored while cows were milked, according to the usual daily routine of a commercial dairy farm. This daily routine was completely unchanged for the scoring and not altered in any way, for example, timewise or by the order of animals for scoring. SNP array genotypes were available for all cows included in the study and had been based on blood samples taken for projects of the respective herdbook association. Thus, the results of SNP array genotyping were available from a central database. Likewise, DNA samples that had been stored as retention samples after SNP array genotyping were available from the two laboratories that had carried out the genotyping with the consent of the two herdbook associations involved, namely, RinderAllianz GmbH and RBB Rinderproduktion Berlin-Brandenburg GmbH. The collection of samples was approved by the Lower Saxony State Office for Consumer Protection and Food Safety (33.19-42502-05-17A196) according to §8a Abs. 1 Nr. 2 of the German Animal Protection Law.

Scoring of phenotypes was conducted in two batches. The first batch consisted of seven large dairy herds visited between October 2015 and April 2016. Herds were selected based on the type of the milking parlor. External rotary parlors were preferred to enable BDD assessment from the outside of the milking carousel. Six of the herds had external rotary parlors, while one farm was equipped with a side-by-side milking parlor. Herds were selected for phenotyping from a pool of 56 herds participating in a scientific project covering a broad range of phenotypic traits for examinations of feasibility of genomic selection based on individual cow phenotypes. In these herds, 20,000 cows, focusing on young cows in their first lactation, had been genotyped based on a protocol consisting of random sampling. Cows were drawn to be included in the study before phenotyping started.

The second batch of data consisted of cows in a further set of six large dairy farms visited between April 2018 and July 2018. Participating herds in this case had enrolled in a genotyping program organized by their herdbook association with the goal of obtaining genomic breeding values for their livestock. The genotyping program was organized in such a way that young stock and cows during first lactation were genotyped by birth year cohorts without any exemptions within birth year. Thus, no selection was practiced. Depending on the time when the herd had enrolled, genotypes were available at least for first lactation and partly also for cows in second lactation.

Every farm in both batches of data was visited three times by a team of three trained persons at intervals of 3 weeks. The team leader and third author is a veterinarian specifically trained for assessing BDD status and was present during all farm visits. All cows present at the time of milking were examined in each herd. As cows in a commercial dairy herd may enter or leave the herd based on herd management, not all cows could be assessed during all three visits. The study group included only cows that were inspected at least twice, with a clear majority of 81.7% being scored three times. After phenotype collection, the resulting data were merged with the genotypes available thus forming the data set useable for analysis. By mere coincidence, after merging phenotypes with genotypes, the first and the second batch of data resulted in 2,520 cows each.

Every cow received a score for the M-stage (M0, M2, or M4) (Döpfer et al., 1997) and a score for signs of chronicity, be it hyperkeratotic or proliferative. A detailed illustration of the scoring system with corresponding symptoms is shown in Supplementary Figure S1. Phenotypes used in our study were defined as follows: for trait TBIN, individuals were coded as 0 if all repetitive scores were M0 and coded as 1 otherwise. For trait TBINA, code 0 was defined as mentioned earlier while a code of 1 denoted an M2 case in any of the repetitive scores. Trait TCHRONA was defined as code = 0 if no chronic proliferative signs were detected in any of the repeated scores and coded = 1 otherwise.

DNA extraction and genotyping was carried out by the IFN Schönow e.V (Institut für Fortpflanzung landwirtschaftlicher Nutztiere Schönow e.V., Bernau, Germany) and the Zentrum für molekulare Diagnostik ZMD at the Institute of Veterinary Medicine (Georg-August-University Göttingen, Germany). The majority of all cows in batches 1 and 2 were genotyped using the Illumina EuroG10K BeadChip versions 4 and 5 (Illumina Inc., San Diego, CA, United States). A subset of 516 animals within batch 1 was analyzed using the Illumina BovineSNP50 BeadChip versions 1 and 2. Quality control was set as call rate >90% and minor allele frequencies >1%. Imputation from the 10-K SNP chip information to the 50-K variant via the bioinformatics tool FImpute version 2.2 (Sargolzaei et al., 2014) with an accuracy of 99.5% was performed by Vereinigte Informationssysteme Tierhaltung (VIT, Verden/Aller, Germany). Finally, a total of 45,613 SNP markers could be employed for a genome-wide association approach.

For GWAS, the software based on the BLUPF90 software family (Misztal et al., 2002) adapted for genomic analyses (Aguilar et al., 2011) was employed in the form of GBLUP analysis. Modeling was performed according to a threshold model as applicable for binary traits. THRGIBBSF90 was used to estimate the variance components and then to predict the GEBV. The SNP effects were calculated using postGSf90 software. The single-trait threshold animal model for GBLUP/GWAS is as follows: P R ( Y i j k l = 1 ) = φ ( μ + h e r d i + p a r i t y j + D I M k + a n i m a l l ) , (1) where PR is the probability of occurrence, y _ijkl is the vector of phenotypic observations coded as binary traits across 2 or 3 observations, φ is the probit link function, μ is the overall mean, herd _i is the fixed effect of the farm (i = 1, …, 7 for batch 1 and i = 1, …, 13 for batch 1 + 2 combined), parity _j is the fixed effect for parity (primipar or multipar), DIM _k is the fixed effect for lactation stage grouped into classes (k = 1, …, 8 for DIM <50, 50–99, 100–149, 150–199, 200–249, 250–299, 300–349, and >350), and animal _l is the random additive genetic effect of animals (l =1, …,5040). For random additive genetic effects, the covariance structure is given by the matrix G. G was set up comprising all cows with phenotypes and genotypes and their genotyped sires as available from the central database. To control for multiple testing, the false discovery rate (FDR) was used (Benjamini and Hochberg 1995). An initial GWAS analysis was conducted after sampling of batch 1 had been completed, and a final analysis included all individuals from batch 1 and 2. As a main result from GWAS, two SNPs with significant effects were obtained for chromosomes 11 and 19 after analysis of batch 1. For identification of animals to be sequenced for candidate genes, the two positions of these SNPs subsequently were subject to an analysis of phased haplotypes. The procedure of the haplotype analysis along with detailed results is presented in detail in Supplementary Tables S1, S2. In brief, phased haplotype data inferred from genotype information of the study animals along with their sires were provided by the Vereinigte Informationssysteme Tierhaltung (VIT, Verden/Aller, Germany), having employed the software FImpute version 2.2 (Sargolzaei et al., 2014). For all 2,520 genotyped cows in batch 1 and both regions, a window of 10 SNPs surrounding the target SNP, (HapMap60356-rs29024194, BTA11, and BTA-45551-no-rs, BTA19) was selected, and haplotypes of all animals were determined.

In a first run of haplotype analysis, all individuals homozygous for a given haplotype formed one effect class while the vast majority of animals were assigned to a class “heterozygous” as they showed some form of heterozygosity for the haplotype. After analysis of traits TBIN, TBINA, and TCHRONA, the results provided an indication on the status of haplotypes being favorable, intermediate, or unfavorable. In a second run, for each allele, the status of each individual was defined based on whether the haplotype allele occurred twice (homozygous), once (heterozygous), or not at all (all other combinations). Only haplotype alleles appearing in at least 10 animals were included in further analysis. The impact of each haplotype allele was estimated applying the SAS procedure GLIMMIX under the following threshold model: P R ( Y i j k l = 1 ) = φ ( μ + h e r d i + p a r i t y j + D I M k + h a p l o t y p e l ) , (2) where PR is the probability of occurence, Y _ijkl is the binary trait (TBIN, TBINA, and TCHRONA; 1 = affected and 0 = unaffected), φ is the probit link function, μ is the overall mean, herd _i is the fixed effect of the farm (i = 1, …, 7), parity _j is the fixed effect for parity (primipar or multipar), DIM _k is the fixed effect for the lactation stage grouped into classes (k = 1, …, 8 for DIM <50, 50–99, 100–149, 150–199, 200–249, 250–299, 300–349, and >350), and haplotype _l is the fixed effect for haplotype (l = 0, 1, and 2 for none, heterozygous, and homozygous, respectively).

Haplotype alleles that would warrant further analysis were chosen on the basis of the p-values of significance as well as the occurrence of homozygous and heterozygous individuals, the number of homozygous individuals (N ≥ 5), and a clear differentiation between LSM estimates of homozygous and heterozygous animals.

For position 90,100,118 (HapMap60356-rs29024194, BTA11), three haplotypes were identified with one being unfavorable in the homozygous status and two being favorable. For each of the three haplotypes, sires being heterozygous were identified with one sire exhibiting haplotype 669/701 (both haplotypes favorable), another sire showing 701/164 haplotypes (favorable/unfavorable), and a further sire with 669/164 haplotypes (favorable/unfavorable). For all three sires, two homozygous daughters for each haplotype were selected for sequencing of the target gene. For position 44,597,888 (rs41603040—BTA-45551-no-rs, BTA19), two haplotypes were identified with one being favorable and one being unfavorable in the homozygous status. For five heterozygous sires, a total of 11 homozygous daughters were sequenced for the target gene.

High-throughput examination of rs208894039 (CMPK2) and rs109521151 (ASB16) was carried out using fluorescence resonance energy transfer (FRET) for 2,485 (CMPK2) and 2,471 (ASB16) animals from batch 1 with phenotype information regarding the established binary BDD traits.

Since in the genome database Ensembl (Genome assembly ARS-UCD1.2) (Zerbino et al., 2018), A and G are, respectively, specified as reference and alternative alleles for rs208894039, genotypes were denoted as homozygous reference (A/A), heterozygous variant (A/G), or homozygous variant (G/G). For SNP rs109521151, A is denoted as the reference allele and G as the variant allele (genome assembly ARS-UCD1.2). Genotype information was statistically tested against the traits using the SAS procedure GLIMMIX under a threshold model considering herd, parity, and days in milk (DIM) as fixed effects.

Phenotype data were collected from cows in a total of 13 large dairy herds. Table 1 displays the results of descriptive statistics for the binary phenotypes as defined across multiple visits. TBIN = 0 denotes a healthy cow throughout visits, while TBIN is coded as 1 for any lesion detected, be it an infectious stage (M2) or a chronic case of M4. TBINA refers to infectious cases scored as 1, while all other stages are coded as 0. Finally, TCHRONA refers to scores given for signs of chronicity with proliferations, coded as 1, while all other cases are scored as 0.

Initially, seven herds were included in the study and phenotyped in 2015 and 2016, and these data are denoted as batch 1. In 2018, cows in further six herds were phenotyped forming batch 2 of the data. Herds included in batch 1 are numbered from one to seven, while second batch herds are numbered from 8 to 13. Differences between herds appear to be of high magnitude when considering infectious cases only as indicated by means for TBINA. Herds 3, 4, 10, and 12 exhibit rather low numbers of M2 lesions and in herds 8 and 9 infectious cases are even scarcer. In contrast, herds 2 and 13 show rather high values for infectious cases. Phenotypes, TBIN and TCHRONA, however, show that in all herds BDD is abundant. It should be noted that phenotypes TBIN and TCHRONA are closely related, since many M4 cases exhibit chronic signs, that is, proliferation and/or hyperkeratotic lesions.

Results for GWAS are presented in Figure 1 (batch 1) and Figure 2 (batch 1 + 2). Details on the significance of individual SNP are given in Supplementary Tables S3A, S3B. Figure 1A indicates one clearly interesting chromosomal region for phenotype TBIN on BTA19, based on chromosome-wise FDR values for three SNPs. Another region of interest is located on BTA11 at position 90,100,118 with one SNP chromosome-wise significant FDR and trait TCHRONA focusing on chronic proliferative cases of BDD (Figure 1B). No significant p-values were found for any chromosomal region and phenotype TBINA (data not shown). Preliminary analyses (not shown) for this phenotype applying simpler approaches for GWAS neglecting co-variances between animals arising from genetic relationships, yielded p-values that differed among each other but provided no indications for any chromosomal region with a suspected effect on the trait. It has to be concluded that given the data shown in Table 1 with most herds exhibiting only low numbers of active M2 lesions, the definition of phenotype TBINA is not useful. For phenotype TCHRONA (Figure 1B), two SNPs with the chromosome-wise FDR value were also found on BTA19. For TBIN and TCHRONA, additional chromosomal regions of interest were found on BTA13 and BTA29. However, as these were not detectable in a consistent manner for both phenotypes, these regions were not considered for further investigation.

Figure 2 shows GWAS results for the full data of batches 1 and 2 combined for phenotypes TBIN and TCHRONA. Due to the low incidence of M2 cases in herds 3, 4, 8, 9, 10, and 12, again no significant estimates of SNP effects could be obtained for phenotype TBINA. Hence, the results for TBINA are not shown. For TBIN (Figure 2A), nine chromosomal regions exhibit a chromosome-wise significance using FDR values. On BTA11, a non-significant -log10 p-value of 4.23 is obtained for HAPMAP60356-rs29024194 at position 90,100,118. Other SNPs around this position show -log10 p-values between 2.0 and 3.0. For BTA19, a total of 11 SNPs show chromosome-wise significant FDR values; among them, two SNPs show a Bonferroni-corrected genome-wide significance for trait TCHRONA (Figure 2B). Considering phenotype TCHRONA (Figure 2B), the general picture for chromosomal regions of interest appears to be similar to TBIN.

Figure 3 displays the GWAS results (Ensembl genome assembly UMD3.1) for traits TBIN and BTA11, focusing on the region between 80 and 100 Mbp. The signal for HAPMAP60356-rs29024194 at position 90,100,118 is flanked by less pronounced signals. The closest gene to HAPMAP60356-rs29024194 is CMPK2 (cytidine/uridine monophosphate kinase 2; Ensembl ID ENSBTAG00000019979) (Zerbino et al., 2018), which lies approximately 41.1 kb upstream of this polymorphism. Since studies in different species indicate an important role of CMPK2 regarding immunomodulatory signaling pathways (Lee and O’Brien, 1995; Zhang et al., 2020), especially in association with bacterial infections (Zhong et al., 2018; Feng et al., 2021), this gene was selected as a candidate gene for a potential influence on BDD.

In Figure 4, the GWAS results (Ensembl genome assembly UMD3.1) for TBIN and BTA19 are shown for the region between 35 and 55 Mbp. Of particular note is BTA-45551-no-rs at position 44,597,888, chromosome-wise significant for both traits and both datasets. From this SNP, the gene ASB16 (ankyrin repeat and SOCS box containing 16; Ensembl ID ENSBTAG000019658) (Zerbino et al., 2018) is located 16.2 kb downstream. For ASB proteins, interaction with numerous effector molecules within inflammatory signaling pathways has already been shown (Andresen et al., 2014). In addition, ASB proteins play an important role in the regulation of protein activities through their function as the substrate recognition unit in a subset of E3 ubiquitin ligases, which is significant in maintaining a balanced immunological response (Anasa et al., 2018; Hou et al., 2021). Considering this, ASB16 was selected as a potential BDD candidate gene.

The two chromosomal regions identified in the GWAS and the respective candidate genes were subsequently used as a basis for a haplotype analysis and subsequent identification of animals for sequencing of the candidate genes. This analysis revealed a remarkable difference between samples at genomic position 90,083,819 (BTA11), which lies 307 bp downstream of CMPK2. This distinct position with rs208894039 is not included on the Illumina Bovine BeadChip used in this study. Sequencing of ASB16 was performed analogously and uncovered a clear variation between samples at genomic position 43,995,946 (BTA19) within exon 5. This position has dbSNP number rs109521151. The alternative G allele of this SNP causes a missense variation leading to incorporation of arginine instead of histidine at amino acid position 447 of ASB16 protein.

Genotypes for rs208894039 and rs109521151 were the results from a FRET analysis of batch 1, with a total of N = 2,485 of evaluable genotypes among the 2,520 total number of samples. Table 2 summarizes the results for the threshold model analysis of rs208894039 (A/G, BTA11) downstream of CMPK2 at position 90,083,819 that was identified after sequencing of this genetic region. Apart from the genotype, as fixed effects in the model herd, parity, and classes for the lactation stage were included. For TBIN, substantial differences exist between LSMeans (LSM) of the three genotypes. WT/WT animals show lowest values for the status = 1 of being diseased while MT/MT animals show highest values. Heterozygous animals appear to react similar to MT/MT animals, thus pointing to an incomplete dominance of the G allele for the disease status. Differences comparing MT/MT with WT/WT and WT/WT with Het are highly significant (p < 0.0001). Comparisons of MT/MT against Het were not significant. For TBINA, including only M2 lesions for status = 1 (diseased), none of the comparisons for the three genotypes were significant. In the bottom part of Table 2, the results for TCHRONA are displayed, thus focusing on chronic proliferative cases vs. other animals taken as healthy. The results shown are highly analogous to the comparisons for TBIN although less pronounced. In summary, for the polymorphism rs208894039, substantial differences between genotypes can be identified despite the lower frequency of the MT allele (q = 0.24). More than half of the animals possess the advantageous genotype WT/WT. However, only one copy of the MT allele drastically changes susceptibility to BDD suggesting a dominant effect of the alternative allele over the reference allele.

Table 3 displays the results for rs109521151 (A/G, BTA19) at position 43,995,946. For this polymorphism within ASB16, N = 2,471 animals had evaluable genotypes. For phenotype TBIN, drastic differences can be observed between MT/MT and WT/WT genotypes, with WT/WT genotypes being advantageous. MT/MT animals are less frequent, with an allele frequency of p (MT) = 0.128. Heterozygous animals exhibit a diseases status of almost intermediate incidence when compared to MT/MT and WT/WT animals. Contrast of MT/MT vs. WT/WT animals was significant with p = 0.0011. Comparing WT/WT animals to heterozygous animals yielded significant differences with p < 0.0001. Contrasts between MT/MT and heterozygous animals were not significant. The results for phenotype TBINA follow the same pattern as shown for TBIN; however, all comparisons between genotypes were not significant. For TCHRONA, that is, when contrasting chronic proliferative cases vs. all other animals, elevated values for the status of phenotype equal to 1 can be found for MT/MT animals, while WT/WT animals show substantially lower values and heterozygous animals have intermediate values. Contrasts between MT/MT and WT/WT animals as well as WT/WT and heterozygous animals were highly significant (p < 0.001), while contrasts between MT/MT and heterozygous animals were significant at a level of p = 0.03. Again, for phenotype TCHRONA, a similar pattern was observed in comparison to phenotype TBIN. In summary of the results in Tables 2, 3, in both cases, the alleles associated with increased susceptibility for BDD are substantially less frequent.

In Table 4, the combination of analysis of the SNP rs208894039 and rs109521151 is shown. A total of N = 2,461 animals had evaluable genotypes for a combination of both polymorphisms. Genotypes are denoted as such that the BTA11 polymorphism is given first, followed by the BTA19 genotype. For interpretation of Table 4, it has to be taken into consideration that the desirable and more common allele for both the CMPK2 variant and the ASB16 variant is WT (see Tables 2, 3). This is confirmed as shown in Table 4 as the WT/WT—WT/WT genotype shows lowest values for TBIN and TCHRONA for the status of being diseased. The combined WT/WT—WT/WT genotype is the most common of all combinations, thus underpinning possible effects of a positive indirect selection. In contrast, the MT/MT—MT/MT genotype exhibits most undesirable values for all phenotypes. For all combinations of combined genotypes, contrasts were estimated and tested for their significance. In the case of phenotype TBINA, not a single contrast approached significance, thus emphasizing again that scoring of the M2 stage animals vs. all other animals does not appear to reveal genomic effects. Hence, the results for TBINA may be dropped from further interpretation. For phenotypes, TBIN and TCHRONA, many LSMeans as well as many contrasts showed significance (Supplementary Table S4) and thus indeed warrant a detailed presentation for the disease status. In Figure 5, the LSMeans for both phenotypes, TBIN (Figure 5A) and TCHRONA (Figure 5B), were ordered from the highest to the lowest value, and the resulting bars were grouped according to the number of copies of the two variants. In both cases, rs109521151 (ASB16) is shown to have a dominant effect and rs208894039 (CMPK2) a modulatory effect.